source: trunk/Clp/src/ClpSimplexOther.cpp @ 1862

Last change on this file since 1862 was 1862, checked in by forrest, 8 years ago

keep parametrics flexible for advanced users

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 355.2 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1862 2012-06-08 15:09:20Z forrest $ */
2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9
10#include "CoinHelperFunctions.hpp"
11#include "ClpSimplexOther.hpp"
12#include "ClpSimplexDual.hpp"
13#include "ClpSimplexPrimal.hpp"
14#include "ClpEventHandler.hpp"
15#include "ClpHelperFunctions.hpp"
16#include "ClpFactorization.hpp"
17#include "ClpDualRowDantzig.hpp"
18#include "ClpNonLinearCost.hpp"
19#include "ClpDynamicMatrix.hpp"
20#include "CoinPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
22#include "CoinBuild.hpp"
23#include "CoinMpsIO.hpp"
24#include "CoinFloatEqual.hpp"
25#include "ClpMessage.hpp"
26#include <cfloat>
27#include <cassert>
28#include <string>
29#include <stdio.h>
30#include <iostream>
31#ifdef HAS_CILK
32#include <cilk/cilk.h>
33#else
34#define cilk_for for
35#define cilk_spawn
36#define cilk_sync
37#endif
38#ifdef INT_IS_8
39#define COIN_ANY_BITS_PER_INT 64
40#define COIN_ANY_SHIFT_PER_INT 6
41#define COIN_ANY_MASK_PER_INT 0x3f
42#else
43#define COIN_ANY_BITS_PER_INT 32
44#define COIN_ANY_SHIFT_PER_INT 5
45#define COIN_ANY_MASK_PER_INT 0x1f
46#endif
47/* Dual ranging.
48   This computes increase/decrease in cost for each given variable and corresponding
49   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
50   and numberColumns.. for artificials/slacks.
51   For non-basic variables the sequence number will be that of the non-basic variables.
52
53   Up to user to provide correct length arrays.
54
55*/
56void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
57                                  double * costIncreased, int * sequenceIncreased,
58                                  double * costDecreased, int * sequenceDecreased,
59                                  double * valueIncrease, double * valueDecrease)
60{
61     rowArray_[1]->clear();
62     columnArray_[1]->clear();
63     // long enough for rows+columns
64     assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_);
65     rowArray_[3]->clear();
66     int * backPivot = rowArray_[3]->getIndices();
67     int i;
68     for ( i = 0; i < numberRows_ + numberColumns_; i++) {
69          backPivot[i] = -1;
70     }
71     for (i = 0; i < numberRows_; i++) {
72          int iSequence = pivotVariable_[i];
73          backPivot[iSequence] = i;
74     }
75     // dualTolerance may be zero if from CBC.  In fact use that fact
76     bool inCBC = !dualTolerance_;
77     if (inCBC)
78          assert (integerType_);
79     dualTolerance_ = dblParam_[ClpDualTolerance];
80     double * arrayX = rowArray_[0]->denseVector();
81     for ( i = 0; i < numberCheck; i++) {
82          rowArray_[0]->clear();
83          //rowArray_[0]->checkClear();
84          //rowArray_[1]->checkClear();
85          //columnArray_[1]->checkClear();
86          columnArray_[0]->clear();
87          //columnArray_[0]->checkClear();
88          int iSequence = which[i];
89          if (iSequence < 0) {
90               costIncreased[i] = 0.0;
91               sequenceIncreased[i] = -1;
92               costDecreased[i] = 0.0;
93               sequenceDecreased[i] = -1;
94               continue;
95          }
96          double costIncrease = COIN_DBL_MAX;
97          double costDecrease = COIN_DBL_MAX;
98          int sequenceIncrease = -1;
99          int sequenceDecrease = -1;
100          if (valueIncrease) {
101               assert (valueDecrease);
102               valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
103               valueDecrease[i] = valueIncrease[i];
104          }
105
106          switch(getStatus(iSequence)) {
107
108          case basic: {
109               // non-trvial
110               // Get pivot row
111               int iRow = backPivot[iSequence];
112               assert (iRow >= 0);
113#ifndef COIN_FAC_NEW
114               double plusOne = 1.0;
115               rowArray_[0]->createPacked(1, &iRow, &plusOne);
116#else
117               rowArray_[0]->createOneUnpackedElement( iRow, 1.0);
118#endif
119               factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
120               // put row of tableau in rowArray[0] and columnArray[0]
121               matrix_->transposeTimes(this, -1.0,
122                                       rowArray_[0], columnArray_[1], columnArray_[0]);
123#if COIN_FAC_NEW
124               assert (!rowArray_[0]->packedMode());
125#endif
126               double alphaIncrease;
127               double alphaDecrease;
128               // do ratio test up and down
129               checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease,
130                               costDecrease, sequenceDecrease, alphaDecrease);
131               if (!inCBC) {
132                    if (valueIncrease) {
133                         if (sequenceIncrease >= 0)
134                              valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence);
135                         if (sequenceDecrease >= 0)
136                              valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence);
137                    }
138               } else {
139                    int number = rowArray_[0]->getNumElements();
140#if COIN_FAC_NEW
141                    const int * index = rowArray_[0]->getIndices();
142#endif
143                    double scale2 = 0.0;
144                    int j;
145                    for (j = 0; j < number; j++) {
146#ifndef COIN_FAC_NEW
147                         scale2 += arrayX[j] * arrayX[j];
148#else
149                         int iRow=index[j];
150                         scale2 += arrayX[iRow] * arrayX[iRow];
151#endif
152                    }
153                    scale2 = 1.0 / sqrt(scale2);
154                    //valueIncrease[i] = scale2;
155                    if (sequenceIncrease >= 0) {
156                         double djValue = dj_[sequenceIncrease];
157                         if (fabs(djValue) > 10.0 * dualTolerance_) {
158                              // we are going to use for cutoff so be exact
159                              costIncrease = fabs(djValue / alphaIncrease);
160                              /* Not sure this is good idea as I don't think correct e.g.
161                                 suppose a continuous variable has dj slightly greater. */
162                              if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) {
163                                   // can improve
164                                   double movement = (columnScale_ == NULL) ? 1.0 :
165                                                     rhsScale_ * inverseColumnScale_[sequenceIncrease];
166                                   costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
167                              }
168                         } else {
169                              costIncrease = 0.0;
170                         }
171                    }
172                    if (sequenceDecrease >= 0) {
173                         double djValue = dj_[sequenceDecrease];
174                         if (fabs(djValue) > 10.0 * dualTolerance_) {
175                              // we are going to use for cutoff so be exact
176                              costDecrease = fabs(djValue / alphaDecrease);
177                              if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
178                                   // can improve
179                                   double movement = (columnScale_ == NULL) ? 1.0 :
180                                                     rhsScale_ * inverseColumnScale_[sequenceDecrease];
181                                   costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
182                              }
183                         } else {
184                              costDecrease = 0.0;
185                         }
186                    }
187                    costIncrease *= scale2;
188                    costDecrease *= scale2;
189               }
190          }
191          break;
192          case isFixed:
193               break;
194          case isFree:
195          case superBasic:
196               costIncrease = 0.0;
197               costDecrease = 0.0;
198               sequenceIncrease = iSequence;
199               sequenceDecrease = iSequence;
200               break;
201          case atUpperBound:
202               costIncrease = CoinMax(0.0, -dj_[iSequence]);
203               sequenceIncrease = iSequence;
204               if (valueIncrease)
205                    valueIncrease[i] = primalRanging1(iSequence, iSequence);
206               break;
207          case atLowerBound:
208               costDecrease = CoinMax(0.0, dj_[iSequence]);
209               sequenceDecrease = iSequence;
210               if (valueIncrease)
211                    valueDecrease[i] = primalRanging1(iSequence, iSequence);
212               break;
213          }
214          double scaleFactor;
215          if (rowScale_) {
216               if (iSequence < numberColumns_)
217                    scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
218               else
219                    scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
220          } else {
221               scaleFactor = 1.0 / objectiveScale_;
222          }
223          if (costIncrease < 1.0e30)
224               costIncrease *= scaleFactor;
225          if (costDecrease < 1.0e30)
226               costDecrease *= scaleFactor;
227          if (optimizationDirection_ == 1.0) {
228               costIncreased[i] = costIncrease;
229               sequenceIncreased[i] = sequenceIncrease;
230               costDecreased[i] = costDecrease;
231               sequenceDecreased[i] = sequenceDecrease;
232          } else if (optimizationDirection_ == -1.0) {
233               costIncreased[i] = costDecrease;
234               sequenceIncreased[i] = sequenceDecrease;
235               costDecreased[i] = costIncrease;
236               sequenceDecreased[i] = sequenceIncrease;
237               if (valueIncrease) {
238                    double temp = valueIncrease[i];
239                    valueIncrease[i] = valueDecrease[i];
240                    valueDecrease[i] = temp;
241               }
242          } else if (optimizationDirection_ == 0.0) {
243               // !!!!!! ???
244               costIncreased[i] = COIN_DBL_MAX;
245               sequenceIncreased[i] = -1;
246               costDecreased[i] = COIN_DBL_MAX;
247               sequenceDecreased[i] = -1;
248          } else {
249               abort();
250          }
251     }
252     rowArray_[0]->clear();
253     //rowArray_[1]->clear();
254     //columnArray_[1]->clear();
255     columnArray_[0]->clear();
256     //rowArray_[3]->clear();
257     if (!optimizationDirection_)
258          printf("*** ????? Ranging with zero optimization costs\n");
259}
260/*
261   Row array has row part of pivot row
262   Column array has column part.
263   This is used in dual ranging
264*/
265void
266ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
267                                 CoinIndexedVector * columnArray,
268                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
269                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
270{
271     double acceptablePivot = 1.0e-9;
272     double * work;
273     int number;
274     int * which;
275     int iSection;
276
277     double thetaDown = 1.0e31;
278     double thetaUp = 1.0e31;
279     int sequenceDown = -1;
280     int sequenceUp = -1;
281     double alphaDown = 0.0;
282     double alphaUp = 0.0;
283
284     int addSequence;
285
286     for (iSection = 0; iSection < 2; iSection++) {
287
288          int i;
289          if (!iSection) {
290               work = rowArray->denseVector();
291               number = rowArray->getNumElements();
292               which = rowArray->getIndices();
293               addSequence = numberColumns_;
294          } else {
295               work = columnArray->denseVector();
296               number = columnArray->getNumElements();
297               which = columnArray->getIndices();
298               addSequence = 0;
299          }
300
301          for (i = 0; i < number; i++) {
302               int iSequence = which[i];
303               int iSequence2 = iSequence + addSequence;
304#ifndef COIN_FAC_NEW
305               double alpha = work[i];
306#else
307               double alpha = !addSequence ? work[i] : work[iSequence];
308#endif
309               if (fabs(alpha) < acceptablePivot)
310                    continue;
311               double oldValue = dj_[iSequence2];
312
313               switch(getStatus(iSequence2)) {
314
315               case basic:
316                    break;
317               case ClpSimplex::isFixed:
318                    break;
319               case isFree:
320               case superBasic:
321                    // treat dj as if zero
322                    thetaDown = 0.0;
323                    thetaUp = 0.0;
324                    sequenceDown = iSequence2;
325                    sequenceUp = iSequence2;
326                    break;
327               case atUpperBound:
328                    if (alpha > 0.0) {
329                         // test up
330                         if (oldValue + thetaUp * alpha > dualTolerance_) {
331                              thetaUp = (dualTolerance_ - oldValue) / alpha;
332                              sequenceUp = iSequence2;
333                              alphaUp = alpha;
334                         }
335                    } else {
336                         // test down
337                         if (oldValue - thetaDown * alpha > dualTolerance_) {
338                              thetaDown = -(dualTolerance_ - oldValue) / alpha;
339                              sequenceDown = iSequence2;
340                              alphaDown = alpha;
341                         }
342                    }
343                    break;
344               case atLowerBound:
345                    if (alpha < 0.0) {
346                         // test up
347                         if (oldValue + thetaUp * alpha < - dualTolerance_) {
348                              thetaUp = -(dualTolerance_ + oldValue) / alpha;
349                              sequenceUp = iSequence2;
350                              alphaUp = alpha;
351                         }
352                    } else {
353                         // test down
354                         if (oldValue - thetaDown * alpha < -dualTolerance_) {
355                              thetaDown = (dualTolerance_ + oldValue) / alpha;
356                              sequenceDown = iSequence2;
357                              alphaDown = alpha;
358                         }
359                    }
360                    break;
361               }
362          }
363     }
364     if (sequenceUp >= 0) {
365          costIncrease = thetaUp;
366          sequenceIncrease = sequenceUp;
367          alphaIncrease = alphaUp;
368     }
369     if (sequenceDown >= 0) {
370          costDecrease = thetaDown;
371          sequenceDecrease = sequenceDown;
372          alphaDecrease = alphaDown;
373     }
374}
375/** Primal ranging.
376    This computes increase/decrease in value for each given variable and corresponding
377    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
378    and numberColumns.. for artificials/slacks.
379    For basic variables the sequence number will be that of the basic variables.
380
381    Up to user to provide correct length arrays.
382
383    When here - guaranteed optimal
384*/
385void
386ClpSimplexOther::primalRanging(int numberCheck, const int * which,
387                               double * valueIncreased, int * sequenceIncreased,
388                               double * valueDecreased, int * sequenceDecreased)
389{
390     rowArray_[0]->clear();
391     rowArray_[1]->clear();
392     lowerIn_ = -COIN_DBL_MAX;
393     upperIn_ = COIN_DBL_MAX;
394     valueIn_ = 0.0;
395     for ( int i = 0; i < numberCheck; i++) {
396          int iSequence = which[i];
397          double valueIncrease = COIN_DBL_MAX;
398          double valueDecrease = COIN_DBL_MAX;
399          int sequenceIncrease = -1;
400          int sequenceDecrease = -1;
401
402          switch(getStatus(iSequence)) {
403
404          case basic:
405          case isFree:
406          case superBasic:
407               // Easy
408               valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
409               valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
410               sequenceDecrease = iSequence;
411               sequenceIncrease = iSequence;
412               break;
413          case isFixed:
414          case atUpperBound:
415          case atLowerBound: {
416               // Non trivial
417               // Other bound is ignored
418#ifndef COIN_FAC_NEW
419               unpackPacked(rowArray_[1], iSequence);
420#else
421               unpack(rowArray_[1], iSequence);
422#endif
423               factorization_->updateColumn(rowArray_[2], rowArray_[1]);
424               // Get extra rows
425               matrix_->extendUpdated(this, rowArray_[1], 0);
426               // do ratio test
427               checkPrimalRatios(rowArray_[1], 1);
428               if (pivotRow_ >= 0) {
429                    valueIncrease = theta_;
430                    sequenceIncrease = pivotVariable_[pivotRow_];
431               }
432               checkPrimalRatios(rowArray_[1], -1);
433               if (pivotRow_ >= 0) {
434                    valueDecrease = theta_;
435                    sequenceDecrease = pivotVariable_[pivotRow_];
436               }
437               rowArray_[1]->clear();
438          }
439          break;
440          }
441          double scaleFactor;
442          if (rowScale_) {
443               if (iSequence < numberColumns_)
444                    scaleFactor = columnScale_[iSequence] / rhsScale_;
445               else
446                    scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
447          } else {
448               scaleFactor = 1.0 / rhsScale_;
449          }
450          if (valueIncrease < 1.0e30)
451               valueIncrease *= scaleFactor;
452          else
453               valueIncrease = COIN_DBL_MAX;
454          if (valueDecrease < 1.0e30)
455               valueDecrease *= scaleFactor;
456          else
457               valueDecrease = COIN_DBL_MAX;
458          valueIncreased[i] = valueIncrease;
459          sequenceIncreased[i] = sequenceIncrease;
460          valueDecreased[i] = valueDecrease;
461          sequenceDecreased[i] = sequenceDecrease;
462     }
463}
464// Returns new value of whichOther when whichIn enters basis
465double
466ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
467{
468     rowArray_[0]->clear();
469     rowArray_[1]->clear();
470     int iSequence = whichIn;
471     double newValue = solution_[whichOther];
472     double alphaOther = 0.0;
473     Status status = getStatus(iSequence);
474     assert (status == atLowerBound || status == atUpperBound);
475     int wayIn = (status == atLowerBound) ? 1 : -1;
476
477     switch(getStatus(iSequence)) {
478
479     case basic:
480     case isFree:
481     case superBasic:
482          assert (whichIn == whichOther);
483          // Easy
484          newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
485          break;
486     case isFixed:
487     case atUpperBound:
488     case atLowerBound:
489          // Non trivial
490     {
491          // Other bound is ignored
492#ifndef COIN_FAC_NEW
493          unpackPacked(rowArray_[1], iSequence);
494#else
495          unpack(rowArray_[1], iSequence);
496#endif
497          factorization_->updateColumn(rowArray_[2], rowArray_[1]);
498          // Get extra rows
499          matrix_->extendUpdated(this, rowArray_[1], 0);
500          // do ratio test
501          double acceptablePivot = 1.0e-7;
502          double * work = rowArray_[1]->denseVector();
503          int number = rowArray_[1]->getNumElements();
504          int * which = rowArray_[1]->getIndices();
505
506          // we may need to swap sign
507          double way = wayIn;
508          double theta = 1.0e30;
509          for (int iIndex = 0; iIndex < number; iIndex++) {
510
511               int iRow = which[iIndex];
512#ifndef COIN_FAC_NEW
513               double alpha = work[iIndex] * way;
514#else
515               double alpha = work[iRow] * way;
516#endif
517               int iPivot = pivotVariable_[iRow];
518               if (iPivot == whichOther) {
519                    alphaOther = alpha;
520                    continue;
521               }
522               double oldValue = solution_[iPivot];
523               if (fabs(alpha) > acceptablePivot) {
524                    if (alpha > 0.0) {
525                         // basic variable going towards lower bound
526                         double bound = lower_[iPivot];
527                         oldValue -= bound;
528                         if (oldValue - theta * alpha < 0.0) {
529                              theta = CoinMax(0.0, oldValue / alpha);
530                         }
531                    } else {
532                         // basic variable going towards upper bound
533                         double bound = upper_[iPivot];
534                         oldValue = oldValue - bound;
535                         if (oldValue - theta * alpha > 0.0) {
536                              theta = CoinMax(0.0, oldValue / alpha);
537                         }
538                    }
539               }
540          }
541          if (whichIn != whichOther) {
542               if (theta < 1.0e30)
543                    newValue -= theta * alphaOther;
544               else
545                    newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
546          } else {
547               newValue += theta * wayIn;
548          }
549     }
550     rowArray_[1]->clear();
551     break;
552     }
553     double scaleFactor;
554     if (rowScale_) {
555          if (whichOther < numberColumns_)
556               scaleFactor = columnScale_[whichOther] / rhsScale_;
557          else
558               scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
559     } else {
560          scaleFactor = 1.0 / rhsScale_;
561     }
562     if (newValue < 1.0e29)
563          if (newValue > -1.0e29)
564               newValue *= scaleFactor;
565          else
566               newValue = -COIN_DBL_MAX;
567     else
568          newValue = COIN_DBL_MAX;
569     return newValue;
570}
571/*
572   Row array has pivot column
573   This is used in primal ranging
574*/
575void
576ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
577                                   int direction)
578{
579     // sequence stays as row number until end
580     pivotRow_ = -1;
581     double acceptablePivot = 1.0e-7;
582     double * work = rowArray->denseVector();
583     int number = rowArray->getNumElements();
584     int * which = rowArray->getIndices();
585
586     // we need to swap sign if going down
587     double way = direction;
588     theta_ = 1.0e30;
589     for (int iIndex = 0; iIndex < number; iIndex++) {
590
591          int iRow = which[iIndex];
592#ifndef COIN_FAC_NEW
593          double alpha = work[iIndex] * way;
594#else
595          double alpha = work[iRow] * way;
596#endif
597          int iPivot = pivotVariable_[iRow];
598          double oldValue = solution_[iPivot];
599          if (fabs(alpha) > acceptablePivot) {
600               if (alpha > 0.0) {
601                    // basic variable going towards lower bound
602                    double bound = lower_[iPivot];
603                    oldValue -= bound;
604                    if (oldValue - theta_ * alpha < 0.0) {
605                         pivotRow_ = iRow;
606                         theta_ = CoinMax(0.0, oldValue / alpha);
607                    }
608               } else {
609                    // basic variable going towards upper bound
610                    double bound = upper_[iPivot];
611                    oldValue = oldValue - bound;
612                    if (oldValue - theta_ * alpha > 0.0) {
613                         pivotRow_ = iRow;
614                         theta_ = CoinMax(0.0, oldValue / alpha);
615                    }
616               }
617          }
618     }
619}
620/* Write the basis in MPS format to the specified file.
621   If writeValues true writes values of structurals
622   (and adds VALUES to end of NAME card)
623
624   Row and column names may be null.
625   formatType is
626   <ul>
627   <li> 0 - normal
628   <li> 1 - extra accuracy
629   <li> 2 - IEEE hex (later)
630   </ul>
631
632   Returns non-zero on I/O error
633
634   This is based on code contributed by Thorsten Koch
635*/
636int
637ClpSimplexOther::writeBasis(const char *filename,
638                            bool writeValues,
639                            int formatType) const
640{
641     formatType = CoinMax(0, formatType);
642     formatType = CoinMin(2, formatType);
643     if (!writeValues)
644          formatType = 0;
645     // See if INTEL if IEEE
646     if (formatType == 2) {
647          // test intel here and add 1 if not intel
648          double value = 1.0;
649          char x[8];
650          memcpy(x, &value, 8);
651          if (x[0] == 63) {
652               formatType ++; // not intel
653          } else {
654               assert (x[0] == 0);
655          }
656     }
657
658     char number[20];
659     FILE * fp = fopen(filename, "w");
660     if (!fp)
661          return -1;
662
663     // NAME card
664
665     if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
666          fprintf(fp, "NAME          BLANK      ");
667     } else {
668          fprintf(fp, "NAME          %s       ", strParam_[ClpProbName].c_str());
669     }
670     if (formatType >= 2)
671          fprintf(fp, "FREEIEEE");
672     else if (writeValues)
673          fprintf(fp, "VALUES");
674     // finish off name
675     fprintf(fp, "\n");
676     int iRow = 0;
677     for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
678          bool printit = false;
679          if( getColumnStatus(iColumn) == ClpSimplex::basic) {
680               printit = true;
681               // Find non basic row
682               for(; iRow < numberRows_; iRow++) {
683                    if (getRowStatus(iRow) != ClpSimplex::basic)
684                         break;
685               }
686               if (lengthNames_) {
687                    if (iRow != numberRows_) {
688                         fprintf(fp, " %s %-8s       %s",
689                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
690                                 columnNames_[iColumn].c_str(),
691                                 rowNames_[iRow].c_str());
692                         iRow++;
693                    } else {
694                         // Allow for too many basics!
695                         fprintf(fp, " BS %-8s       ",
696                                 columnNames_[iColumn].c_str());
697                         // Dummy row name if values
698                         if (writeValues)
699                              fprintf(fp, "      _dummy_");
700                    }
701               } else {
702                    // no names
703                    if (iRow != numberRows_) {
704                         fprintf(fp, " %s C%7.7d     R%7.7d",
705                                 getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
706                                 iColumn, iRow);
707                         iRow++;
708                    } else {
709                         // Allow for too many basics!
710                         fprintf(fp, " BS C%7.7d", iColumn);
711                         // Dummy row name if values
712                         if (writeValues)
713                              fprintf(fp, "      _dummy_");
714                    }
715               }
716          } else  {
717               if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
718                    printit = true;
719                    if (lengthNames_)
720                         fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
721                    else
722                         fprintf(fp, " UL C%7.7d", iColumn);
723                    // Dummy row name if values
724                    if (writeValues)
725                         fprintf(fp, "      _dummy_");
726               }
727          }
728          if (printit && writeValues) {
729               // add value
730               CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
731               fprintf(fp, "     %s", number);
732          }
733          if (printit)
734               fprintf(fp, "\n");
735     }
736     fprintf(fp, "ENDATA\n");
737     fclose(fp);
738     return 0;
739}
740// Read a basis from the given filename
741int
742ClpSimplexOther::readBasis(const char *fileName)
743{
744     int status = 0;
745     bool canOpen = false;
746     if (!strcmp(fileName, "-") || !strcmp(fileName, "stdin")) {
747          // stdin
748          canOpen = true;
749     } else {
750          FILE *fp = fopen(fileName, "r");
751          if (fp) {
752               // can open - lets go for it
753               fclose(fp);
754               canOpen = true;
755          } else {
756               handler_->message(CLP_UNABLE_OPEN, messages_)
757                         << fileName << CoinMessageEol;
758               return -1;
759          }
760     }
761     CoinMpsIO m;
762     m.passInMessageHandler(handler_);
763     *m.messagesPointer() = coinMessages();
764     bool savePrefix = m.messageHandler()->prefix();
765     m.messageHandler()->setPrefix(handler_->prefix());
766     status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
767                          status_,
768                          columnNames_, numberColumns_,
769                          rowNames_, numberRows_);
770     m.messageHandler()->setPrefix(savePrefix);
771     if (status >= 0) {
772          if (!status) {
773               // set values
774               int iColumn, iRow;
775               for (iRow = 0; iRow < numberRows_; iRow++) {
776                    if (getRowStatus(iRow) == atLowerBound)
777                         rowActivity_[iRow] = rowLower_[iRow];
778                    else if (getRowStatus(iRow) == atUpperBound)
779                         rowActivity_[iRow] = rowUpper_[iRow];
780               }
781               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
782                    if (getColumnStatus(iColumn) == atLowerBound)
783                         columnActivity_[iColumn] = columnLower_[iColumn];
784                    else if (getColumnStatus(iColumn) == atUpperBound)
785                         columnActivity_[iColumn] = columnUpper_[iColumn];
786               }
787          } else {
788               memset(rowActivity_, 0, numberRows_ * sizeof(double));
789               matrix_->times(-1.0, columnActivity_, rowActivity_);
790          }
791     } else {
792          // errors
793          handler_->message(CLP_IMPORT_ERRORS, messages_)
794                    << status << fileName << CoinMessageEol;
795     }
796     return status;
797}
798/* Creates dual of a problem if looks plausible
799   (defaults will always create model)
800   fractionRowRanges is fraction of rows allowed to have ranges
801   fractionColumnRanges is fraction of columns allowed to have ranges
802*/
803ClpSimplex *
804ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
805{
806     const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
807     bool changed = false;
808     int numberChanged = 0;
809     int iColumn;
810     // check if we need to change bounds to rows
811     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
812          if (columnUpper_[iColumn] < 1.0e20 &&
813                    columnLower_[iColumn] > -1.0e20) {
814               changed = true;
815               numberChanged++;
816          }
817     }
818     int iRow;
819     int numberExtraRows = 0;
820     if (numberChanged <= fractionColumnRanges * numberColumns_) {
821          for (iRow = 0; iRow < numberRows_; iRow++) {
822               if (rowLower_[iRow] > -1.0e20 &&
823                         rowUpper_[iRow] < 1.0e20) {
824                    if (rowUpper_[iRow] != rowLower_[iRow])
825                         numberExtraRows++;
826               }
827          }
828          if (numberExtraRows > fractionRowRanges * numberRows_)
829               return NULL;
830     } else {
831          return NULL;
832     }
833     if (changed) {
834          ClpSimplex * model3 = new ClpSimplex(*model2);
835          CoinBuild build;
836          double one = 1.0;
837          int numberColumns = model3->numberColumns();
838          const double * columnLower = model3->columnLower();
839          const double * columnUpper = model3->columnUpper();
840          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
841               if (columnUpper[iColumn] < 1.0e20 &&
842                         columnLower[iColumn] > -1.0e20) {
843                    if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
844                         double value = columnUpper[iColumn];
845                         model3->setColumnUpper(iColumn, COIN_DBL_MAX);
846                         build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
847                    } else {
848                         double value = columnLower[iColumn];
849                         model3->setColumnLower(iColumn, -COIN_DBL_MAX);
850                         build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
851                    }
852               }
853          }
854          model3->addRows(build);
855          model2 = model3;
856     }
857     int numberColumns = model2->numberColumns();
858     const double * columnLower = model2->columnLower();
859     const double * columnUpper = model2->columnUpper();
860     int numberRows = model2->numberRows();
861     double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
862     double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
863
864     const double * objective = model2->objective();
865     CoinPackedMatrix * matrix = model2->matrix();
866     // get transpose
867     CoinPackedMatrix rowCopy = *matrix;
868     const int * row = matrix->getIndices();
869     const int * columnLength = matrix->getVectorLengths();
870     const CoinBigIndex * columnStart = matrix->getVectorStarts();
871     const double * elementByColumn = matrix->getElements();
872     double objOffset = 0.0;
873     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
874          double offset = 0.0;
875          double objValue = optimizationDirection_ * objective[iColumn];
876          if (columnUpper[iColumn] > 1.0e20) {
877               if (columnLower[iColumn] > -1.0e20)
878                    offset = columnLower[iColumn];
879          } else if (columnLower[iColumn] < -1.0e20) {
880               offset = columnUpper[iColumn];
881          } else {
882               // taken care of before
883               abort();
884          }
885          if (offset) {
886               objOffset += offset * objValue;
887               for (CoinBigIndex j = columnStart[iColumn];
888                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
889                    int iRow = row[j];
890                    if (rowLower[iRow] > -1.0e20)
891                         rowLower[iRow] -= offset * elementByColumn[j];
892                    if (rowUpper[iRow] < 1.0e20)
893                         rowUpper[iRow] -= offset * elementByColumn[j];
894               }
895          }
896     }
897     int * which = new int[numberRows+numberExtraRows];
898     rowCopy.reverseOrdering();
899     rowCopy.transpose();
900     double * fromRowsLower = new double[numberRows+numberExtraRows];
901     double * fromRowsUpper = new double[numberRows+numberExtraRows];
902     double * newObjective = new double[numberRows+numberExtraRows];
903     double * fromColumnsLower = new double[numberColumns];
904     double * fromColumnsUpper = new double[numberColumns];
905     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
906          double objValue = optimizationDirection_ * objective[iColumn];
907          // Offset is already in
908          if (columnUpper[iColumn] > 1.0e20) {
909               if (columnLower[iColumn] > -1.0e20) {
910                    fromColumnsLower[iColumn] = -COIN_DBL_MAX;
911                    fromColumnsUpper[iColumn] = objValue;
912               } else {
913                    // free
914                    fromColumnsLower[iColumn] = objValue;
915                    fromColumnsUpper[iColumn] = objValue;
916               }
917          } else if (columnLower[iColumn] < -1.0e20) {
918               fromColumnsLower[iColumn] = objValue;
919               fromColumnsUpper[iColumn] = COIN_DBL_MAX;
920          } else {
921               abort();
922          }
923     }
924     int kRow = 0;
925     int kExtraRow = numberRows;
926     for (iRow = 0; iRow < numberRows; iRow++) {
927          if (rowLower[iRow] < -1.0e20) {
928               assert (rowUpper[iRow] < 1.0e20);
929               newObjective[kRow] = -rowUpper[iRow];
930               fromRowsLower[kRow] = -COIN_DBL_MAX;
931               fromRowsUpper[kRow] = 0.0;
932               which[kRow] = iRow;
933               kRow++;
934          } else if (rowUpper[iRow] > 1.0e20) {
935               newObjective[kRow] = -rowLower[iRow];
936               fromRowsLower[kRow] = 0.0;
937               fromRowsUpper[kRow] = COIN_DBL_MAX;
938               which[kRow] = iRow;
939               kRow++;
940          } else {
941               if (rowUpper[iRow] == rowLower[iRow]) {
942                    newObjective[kRow] = -rowLower[iRow];
943                    fromRowsLower[kRow] = -COIN_DBL_MAX;;
944                    fromRowsUpper[kRow] = COIN_DBL_MAX;
945                    which[kRow] = iRow;
946                    kRow++;
947               } else {
948                    // range
949                    newObjective[kRow] = -rowUpper[iRow];
950                    fromRowsLower[kRow] = -COIN_DBL_MAX;
951                    fromRowsUpper[kRow] = 0.0;
952                    which[kRow] = iRow;
953                    kRow++;
954                    newObjective[kExtraRow] = -rowLower[iRow];
955                    fromRowsLower[kExtraRow] = 0.0;
956                    fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
957                    which[kExtraRow] = iRow;
958                    kExtraRow++;
959               }
960          }
961     }
962     if (numberExtraRows) {
963          CoinPackedMatrix newCopy;
964          newCopy.setExtraGap(0.0);
965          newCopy.setExtraMajor(0.0);
966          newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
967          rowCopy = newCopy;
968     }
969     ClpSimplex * modelDual = new ClpSimplex();
970     modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
971                            fromColumnsLower, fromColumnsUpper);
972     modelDual->setObjectiveOffset(objOffset);
973     modelDual->setDualBound(model2->dualBound());
974     modelDual->setInfeasibilityCost(model2->infeasibilityCost());
975     modelDual->setDualTolerance(model2->dualTolerance());
976     modelDual->setPrimalTolerance(model2->primalTolerance());
977     modelDual->setPerturbation(model2->perturbation());
978     modelDual->setSpecialOptions(model2->specialOptions());
979     modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
980     delete [] fromRowsLower;
981     delete [] fromRowsUpper;
982     delete [] fromColumnsLower;
983     delete [] fromColumnsUpper;
984     delete [] newObjective;
985     delete [] which;
986     delete [] rowLower;
987     delete [] rowUpper;
988     if (changed)
989          delete model2;
990     modelDual->createStatus();
991     return modelDual;
992}
993// Restores solution from dualized problem
994int
995ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
996                                 bool checkAccuracy)
997{
998     int returnCode = 0;;
999     createStatus();
1000     // Number of rows in dual problem was original number of columns
1001     assert (numberColumns_ == dualProblem->numberRows());
1002     // If slack on d-row basic then column at bound otherwise column basic
1003     // If d-column basic then rhs tight
1004     int numberBasic = 0;
1005     int iRow, iColumn = 0;
1006     // Get number of extra rows from ranges
1007     int numberExtraRows = 0;
1008     for (iRow = 0; iRow < numberRows_; iRow++) {
1009          if (rowLower_[iRow] > -1.0e20 &&
1010                    rowUpper_[iRow] < 1.0e20) {
1011               if (rowUpper_[iRow] != rowLower_[iRow])
1012                    numberExtraRows++;
1013          }
1014     }
1015     const double * objective = this->objective();
1016     const double * dualDual = dualProblem->dualRowSolution();
1017     const double * dualDj = dualProblem->dualColumnSolution();
1018     const double * dualSol = dualProblem->primalColumnSolution();
1019     const double * dualActs = dualProblem->primalRowSolution();
1020#if 0
1021     ClpSimplex thisCopy = *this;
1022     thisCopy.dual(); // for testing
1023     const double * primalDual = thisCopy.dualRowSolution();
1024     const double * primalDj = thisCopy.dualColumnSolution();
1025     const double * primalSol = thisCopy.primalColumnSolution();
1026     const double * primalActs = thisCopy.primalRowSolution();
1027     char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
1028     printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
1029     for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
1030          printf("%d at %c primal %g dual %g\n",
1031                 iRow, ss[dualProblem->getRowStatus(iRow)],
1032                 dualActs[iRow], dualDual[iRow]);
1033     printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
1034     for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
1035          printf("%d at %c primal %g dual %g\n",
1036                 iColumn, ss[dualProblem->getColumnStatus(iColumn)],
1037                 dualSol[iColumn], dualDj[iColumn]);
1038     printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
1039     for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
1040          printf("%d at %c primal %g dual %g\n",
1041                 iRow, ss[thisCopy.getRowStatus(iRow)],
1042                 primalActs[iRow], primalDual[iRow]);
1043     printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
1044     for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
1045          printf("%d at %c primal %g dual %g\n",
1046                 iColumn, ss[thisCopy.getColumnStatus(iColumn)],
1047                 primalSol[iColumn], primalDj[iColumn]);
1048#endif
1049     // position at bound information
1050     int jColumn = numberRows_;
1051     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1052          double objValue = optimizationDirection_ * objective[iColumn];
1053          Status status = dualProblem->getRowStatus(iColumn);
1054          double otherValue = COIN_DBL_MAX;
1055          if (columnUpper_[iColumn] < 1.0e20 &&
1056                    columnLower_[iColumn] > -1.0e20) {
1057               if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1058                    otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1059               } else {
1060                    otherValue = columnLower_[iColumn] + dualDj[jColumn];
1061               }
1062               jColumn++;
1063          }
1064          if (status == basic) {
1065               // column is at bound
1066               if (otherValue == COIN_DBL_MAX) {
1067                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1068                    if (columnUpper_[iColumn] > 1.0e20) {
1069                         if (columnLower_[iColumn] > -1.0e20) {
1070                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1071                                   setColumnStatus(iColumn, atLowerBound);
1072                              else
1073                                   setColumnStatus(iColumn, isFixed);
1074                              columnActivity_[iColumn] = columnLower_[iColumn];
1075                         } else {
1076                              // free
1077                              setColumnStatus(iColumn, isFree);
1078                              columnActivity_[iColumn] = 0.0;
1079                         }
1080                    } else {
1081                         setColumnStatus(iColumn, atUpperBound);
1082                         columnActivity_[iColumn] = columnUpper_[iColumn];
1083                    }
1084               } else {
1085                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1086                    //printf("other dual sol %g\n",otherValue);
1087                    if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1088                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1089                              setColumnStatus(iColumn, atLowerBound);
1090                         else
1091                              setColumnStatus(iColumn, isFixed);
1092                         columnActivity_[iColumn] = columnLower_[iColumn];
1093                    } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1094                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1095                              setColumnStatus(iColumn, atUpperBound);
1096                         else
1097                              setColumnStatus(iColumn, isFixed);
1098                         columnActivity_[iColumn] = columnUpper_[iColumn];
1099                    } else {
1100                         abort();
1101                    }
1102               }
1103          } else {
1104               if (otherValue == COIN_DBL_MAX) {
1105                    // column basic
1106                    setColumnStatus(iColumn, basic);
1107                    numberBasic++;
1108                    if (columnLower_[iColumn] > -1.0e20) {
1109                         columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1110                    } else if (columnUpper_[iColumn] < 1.0e20) {
1111                         columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1112                    } else {
1113                         columnActivity_[iColumn] = -dualDual[iColumn];
1114                    }
1115                    reducedCost_[iColumn] = 0.0;
1116               } else {
1117                    // may be at other bound
1118                    //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1119                    if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1120                         // column basic
1121                         setColumnStatus(iColumn, basic);
1122                         numberBasic++;
1123                         //printf("Col %d otherV %g dualDual %g\n",iColumn,
1124                         // otherValue,dualDual[iColumn]);
1125                         columnActivity_[iColumn] = -dualDual[iColumn];
1126                         columnActivity_[iColumn] = otherValue;
1127                         reducedCost_[iColumn] = 0.0;
1128                    } else {
1129                         reducedCost_[iColumn] = objValue - dualActs[iColumn];
1130                         if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1131                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1132                                   setColumnStatus(iColumn, atLowerBound);
1133                              else
1134                                   setColumnStatus(iColumn, isFixed);
1135                              columnActivity_[iColumn] = columnLower_[iColumn];
1136                         } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1137                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1138                                   setColumnStatus(iColumn, atUpperBound);
1139                              else
1140                                   setColumnStatus(iColumn, isFixed);
1141                              columnActivity_[iColumn] = columnUpper_[iColumn];
1142                         } else {
1143                              abort();
1144                         }
1145                    }
1146               }
1147          }
1148     }
1149     // now rows
1150     int kExtraRow = jColumn;
1151     int numberRanges = 0;
1152     for (iRow = 0; iRow < numberRows_; iRow++) {
1153          Status status = dualProblem->getColumnStatus(iRow);
1154          if (status == basic) {
1155               // row is at bound
1156               dual_[iRow] = dualSol[iRow];;
1157          } else {
1158               // row basic
1159               setRowStatus(iRow, basic);
1160               numberBasic++;
1161               dual_[iRow] = 0.0;
1162          }
1163          if (rowLower_[iRow] < -1.0e20) {
1164               if (status == basic) {
1165                    rowActivity_[iRow] = rowUpper_[iRow];
1166                    setRowStatus(iRow, atUpperBound);
1167               } else {
1168                    assert (dualDj[iRow] < 1.0e-5);
1169                    rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1170               }
1171          } else if (rowUpper_[iRow] > 1.0e20) {
1172               if (status == basic) {
1173                    rowActivity_[iRow] = rowLower_[iRow];
1174                    setRowStatus(iRow, atLowerBound);
1175               } else {
1176                    rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1177                    assert (dualDj[iRow] > -1.0e-5);
1178               }
1179          } else {
1180               if (rowUpper_[iRow] == rowLower_[iRow]) {
1181                    rowActivity_[iRow] = rowLower_[iRow];
1182                    if (status == basic) {
1183                         setRowStatus(iRow, isFixed);
1184                    }
1185               } else {
1186                    // range
1187                    numberRanges++;
1188                    Status statusL = dualProblem->getColumnStatus(kExtraRow);
1189                    //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1190                    //     iRow,status,kExtraRow,statusL, dualSol[iRow],
1191                    //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1192                    if (status == basic) {
1193                         assert (statusL != basic);
1194                         rowActivity_[iRow] = rowUpper_[iRow];
1195                         setRowStatus(iRow, atUpperBound);
1196                    } else if (statusL == basic) {
1197                         numberBasic--; // already counted
1198                         rowActivity_[iRow] = rowLower_[iRow];
1199                         setRowStatus(iRow, atLowerBound);
1200                         dual_[iRow] = dualSol[kExtraRow];;
1201                    } else {
1202                         rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1203                         assert (dualDj[iRow] < 1.0e-5);
1204                         // row basic
1205                         //setRowStatus(iRow,basic);
1206                         //numberBasic++;
1207                         dual_[iRow] = 0.0;
1208                    }
1209                    kExtraRow++;
1210               }
1211          }
1212     }
1213     if (numberBasic != numberRows_) {
1214          printf("Bad basis - ranges - coding needed\n");
1215          assert (numberRanges);
1216          abort();
1217     }
1218     if (optimizationDirection_ < 0.0) {
1219          for (iRow = 0; iRow < numberRows_; iRow++) {
1220               dual_[iRow] = -dual_[iRow];
1221          }
1222     }
1223     // redo row activities
1224     memset(rowActivity_, 0, numberRows_ * sizeof(double));
1225     matrix_->times(1.0, columnActivity_, rowActivity_);
1226     // redo reduced costs
1227     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1228     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1229     checkSolutionInternal();
1230     if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
1231          returnCode = 1;
1232#ifdef CLP_INVESTIGATE
1233          printf("There are %d dual infeasibilities summing to %g ",
1234                 numberDualInfeasibilities_, sumDualInfeasibilities_);
1235          printf("and %d primal infeasibilities summing to %g\n",
1236                 numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
1237#endif
1238     }
1239     // Below will go to ..DEBUG later
1240#if 1 //ndef NDEBUG
1241     if (checkAccuracy) {
1242       // Check if correct
1243       double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1244       double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1245       double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1246       double * dual = CoinCopyOfArray(dual_, numberRows_);
1247       this->dual(); //primal();
1248       CoinRelFltEq eq(1.0e-5);
1249       for (iRow = 0; iRow < numberRows_; iRow++) {
1250         assert(eq(dual[iRow], dual_[iRow]));
1251       }
1252       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1253         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1254       }
1255       for (iRow = 0; iRow < numberRows_; iRow++) {
1256         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1257       }
1258       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1259         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1260       }
1261       delete [] columnActivity;
1262       delete [] rowActivity;
1263       delete [] reducedCost;
1264       delete [] dual;
1265     }
1266#endif
1267     return returnCode;
1268}
1269/* Does very cursory presolve.
1270   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1271*/
1272ClpSimplex *
1273ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1274                        int & nBound, bool moreBounds, bool tightenBounds)
1275{
1276     //#define CHECK_STATUS
1277#ifdef CHECK_STATUS
1278     {
1279          int n = 0;
1280          int i;
1281          for (i = 0; i < numberColumns_; i++)
1282               if (getColumnStatus(i) == ClpSimplex::basic)
1283                    n++;
1284          for (i = 0; i < numberRows_; i++)
1285               if (getRowStatus(i) == ClpSimplex::basic)
1286                    n++;
1287          assert (n == numberRows_);
1288     }
1289#endif
1290
1291     const double * element = matrix_->getElements();
1292     const int * row = matrix_->getIndices();
1293     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1294     const int * columnLength = matrix_->getVectorLengths();
1295
1296     CoinZeroN(rhs, numberRows_);
1297     int iColumn;
1298     int iRow;
1299     CoinZeroN(whichRow, numberRows_);
1300     int * backColumn = whichColumn + numberColumns_;
1301     int numberRows2 = 0;
1302     int numberColumns2 = 0;
1303     double offset = 0.0;
1304     const double * objective = this->objective();
1305     double * solution = columnActivity_;
1306     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1307          double lower = columnLower_[iColumn];
1308          double upper = columnUpper_[iColumn];
1309          if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
1310               backColumn[iColumn] = numberColumns2;
1311               whichColumn[numberColumns2++] = iColumn;
1312               for (CoinBigIndex j = columnStart[iColumn];
1313                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1314                    int iRow = row[j];
1315                    int n = whichRow[iRow];
1316                    if (n == 0 && element[j])
1317                         whichRow[iRow] = -iColumn - 1;
1318                    else if (n < 0)
1319                         whichRow[iRow] = 2;
1320               }
1321          } else {
1322               // fixed
1323               backColumn[iColumn] = -1;
1324               solution[iColumn] = upper;
1325               if (upper) {
1326                    offset += objective[iColumn] * upper;
1327                    for (CoinBigIndex j = columnStart[iColumn];
1328                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1329                         int iRow = row[j];
1330                         double value = element[j];
1331                         rhs[iRow] += upper * value;
1332                    }
1333               }
1334          }
1335     }
1336     int returnCode = 0;
1337     double tolerance = primalTolerance();
1338     nBound = 2 * numberRows_;
1339     for (iRow = 0; iRow < numberRows_; iRow++) {
1340          int n = whichRow[iRow];
1341          if (n > 0) {
1342               whichRow[numberRows2++] = iRow;
1343          } else if (n < 0) {
1344               //whichRow[numberRows2++]=iRow;
1345               //continue;
1346               // Can only do in certain circumstances as we don't know current value
1347               if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
1348                    // save row and column for bound
1349                    whichRow[--nBound] = iRow;
1350                    whichRow[nBound+numberRows_] = -n - 1;
1351               } else if (moreBounds) {
1352                    // save row and column for bound
1353                    whichRow[--nBound] = iRow;
1354                    whichRow[nBound+numberRows_] = -n - 1;
1355               } else {
1356                    whichRow[numberRows2++] = iRow;
1357               }
1358          } else {
1359               // empty
1360               double rhsValue = rhs[iRow];
1361               if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
1362                    returnCode = 1; // infeasible
1363               }
1364          }
1365     }
1366     ClpSimplex * small = NULL;
1367     if (!returnCode) {
1368       //printf("CRUNCH from (%d,%d) to (%d,%d)\n",
1369       //     numberRows_,numberColumns_,numberRows2,numberColumns2);
1370          small = new ClpSimplex(this, numberRows2, whichRow,
1371                                 numberColumns2, whichColumn, true, false);
1372#if 0
1373          ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1374          if (rowCopy) {
1375               assert(!small->rowCopy());
1376               small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1377                                    numberColumns2, whichColumn));
1378          }
1379#endif
1380          // Set some stuff
1381          small->setDualBound(dualBound_);
1382          small->setInfeasibilityCost(infeasibilityCost_);
1383          small->setSpecialOptions(specialOptions_);
1384          small->setPerturbation(perturbation_);
1385          small->defaultFactorizationFrequency();
1386          small->setAlphaAccuracy(alphaAccuracy_);
1387          // If no rows left then no tightening!
1388          if (!numberRows2 || !numberColumns2)
1389               tightenBounds = false;
1390
1391          int numberElements = getNumElements();
1392          int numberElements2 = small->getNumElements();
1393          small->setObjectiveOffset(objectiveOffset() - offset);
1394          handler_->message(CLP_CRUNCH_STATS, messages_)
1395                    << numberRows2 << -(numberRows_ - numberRows2)
1396                    << numberColumns2 << -(numberColumns_ - numberColumns2)
1397                    << numberElements2 << -(numberElements - numberElements2)
1398                    << CoinMessageEol;
1399          // And set objective value to match
1400          small->setObjectiveValue(this->objectiveValue());
1401          double * rowLower2 = small->rowLower();
1402          double * rowUpper2 = small->rowUpper();
1403          int jRow;
1404          for (jRow = 0; jRow < numberRows2; jRow++) {
1405               iRow = whichRow[jRow];
1406               if (rowLower2[jRow] > -1.0e20)
1407                    rowLower2[jRow] -= rhs[iRow];
1408               if (rowUpper2[jRow] < 1.0e20)
1409                    rowUpper2[jRow] -= rhs[iRow];
1410          }
1411          // and bounds
1412          double * columnLower2 = small->columnLower();
1413          double * columnUpper2 = small->columnUpper();
1414          const char * integerInformation = integerType_;
1415          for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1416               iRow = whichRow[jRow];
1417               iColumn = whichRow[jRow+numberRows_];
1418               double lowerRow = rowLower_[iRow];
1419               if (lowerRow > -1.0e20)
1420                    lowerRow -= rhs[iRow];
1421               double upperRow = rowUpper_[iRow];
1422               if (upperRow < 1.0e20)
1423                    upperRow -= rhs[iRow];
1424               int jColumn = backColumn[iColumn];
1425               double lower = columnLower2[jColumn];
1426               double upper = columnUpper2[jColumn];
1427               double value = 0.0;
1428               for (CoinBigIndex j = columnStart[iColumn];
1429                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1430                    if (iRow == row[j]) {
1431                         value = element[j];
1432                         break;
1433                    }
1434               }
1435               assert (value);
1436               // convert rowLower and Upper to implied bounds on column
1437               double newLower = -COIN_DBL_MAX;
1438               double newUpper = COIN_DBL_MAX;
1439               if (value > 0.0) {
1440                    if (lowerRow > -1.0e20)
1441                         newLower = lowerRow / value;
1442                    if (upperRow < 1.0e20)
1443                         newUpper = upperRow / value;
1444               } else {
1445                    if (upperRow < 1.0e20)
1446                         newLower = upperRow / value;
1447                    if (lowerRow > -1.0e20)
1448                         newUpper = lowerRow / value;
1449               }
1450               if (integerInformation && integerInformation[iColumn]) {
1451                    if (newLower - floor(newLower) < 10.0 * tolerance)
1452                         newLower = floor(newLower);
1453                    else
1454                         newLower = ceil(newLower);
1455                    if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1456                         newUpper = ceil(newUpper);
1457                    else
1458                         newUpper = floor(newUpper);
1459               }
1460               newLower = CoinMax(lower, newLower);
1461               newUpper = CoinMin(upper, newUpper);
1462               if (newLower > newUpper + tolerance) {
1463                    //printf("XXYY inf on bound\n");
1464                    returnCode = 1;
1465               }
1466               columnLower2[jColumn] = newLower;
1467               columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1468               if (getRowStatus(iRow) != ClpSimplex::basic) {
1469                    if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1470                         if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1471                              // can only get here if will be fixed
1472                              small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1473                         } else {
1474                              // solution is valid
1475                              if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1476                                        fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1477                                   small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1478                              else
1479                                   small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1480                         }
1481                    } else {
1482                         //printf("what now neither basic\n");
1483                    }
1484               }
1485          }
1486          if (returnCode) {
1487               delete small;
1488               small = NULL;
1489          } else if (tightenBounds && integerInformation) {
1490               // See if we can tighten any bounds
1491               // use rhs for upper and small duals for lower
1492               double * up = rhs;
1493               double * lo = small->dualRowSolution();
1494               const double * element = small->clpMatrix()->getElements();
1495               const int * row = small->clpMatrix()->getIndices();
1496               const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1497               //const int * columnLength = small->clpMatrix()->getVectorLengths();
1498               CoinZeroN(lo, numberRows2);
1499               CoinZeroN(up, numberRows2);
1500               for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1501                    double upper = columnUpper2[iColumn];
1502                    double lower = columnLower2[iColumn];
1503                    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1504                    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1505                         int iRow = row[j];
1506                         double value = element[j];
1507                         if (value > 0.0) {
1508                              if (upper < 1.0e20)
1509                                   up[iRow] += upper * value;
1510                              else
1511                                   up[iRow] = COIN_DBL_MAX;
1512                              if (lower > -1.0e20)
1513                                   lo[iRow] += lower * value;
1514                              else
1515                                   lo[iRow] = -COIN_DBL_MAX;
1516                         } else {
1517                              if (upper < 1.0e20)
1518                                   lo[iRow] += upper * value;
1519                              else
1520                                   lo[iRow] = -COIN_DBL_MAX;
1521                              if (lower > -1.0e20)
1522                                   up[iRow] += lower * value;
1523                              else
1524                                   up[iRow] = COIN_DBL_MAX;
1525                         }
1526                    }
1527               }
1528               double * rowLower2 = small->rowLower();
1529               double * rowUpper2 = small->rowUpper();
1530               bool feasible = true;
1531               // make safer
1532               for (int iRow = 0; iRow < numberRows2; iRow++) {
1533                    double lower = lo[iRow];
1534                    if (lower > rowUpper2[iRow] + tolerance) {
1535                         feasible = false;
1536                         break;
1537                    } else {
1538                         lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1539                    }
1540                    double upper = up[iRow];
1541                    if (upper < rowLower2[iRow] - tolerance) {
1542                         feasible = false;
1543                         break;
1544                    } else {
1545                         up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1546                    }
1547               }
1548               if (!feasible) {
1549                    delete small;
1550                    small = NULL;
1551               } else {
1552                    // and tighten
1553                    for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1554                         if (integerInformation[whichColumn[iColumn]]) {
1555                              double upper = columnUpper2[iColumn];
1556                              double lower = columnLower2[iColumn];
1557                              double newUpper = upper;
1558                              double newLower = lower;
1559                              double difference = upper - lower;
1560                              if (lower > -1000.0 && upper < 1000.0) {
1561                                   for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1562                                        int iRow = row[j];
1563                                        double value = element[j];
1564                                        if (value > 0.0) {
1565                                             double upWithOut = up[iRow] - value * difference;
1566                                             if (upWithOut < 0.0) {
1567                                                  newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1568                                             }
1569                                             double lowWithOut = lo[iRow] + value * difference;
1570                                             if (lowWithOut > 0.0) {
1571                                                  newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1572                                             }
1573                                        } else {
1574                                             double upWithOut = up[iRow] + value * difference;
1575                                             if (upWithOut < 0.0) {
1576                                                  newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1577                                             }
1578                                             double lowWithOut = lo[iRow] - value * difference;
1579                                             if (lowWithOut > 0.0) {
1580                                                  newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1581                                             }
1582                                        }
1583                                   }
1584                                   if (newLower > lower || newUpper < upper) {
1585                                        if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1586                                             newUpper = floor(newUpper);
1587                                        else
1588                                             newUpper = floor(newUpper + 0.5);
1589                                        if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1590                                             newLower = ceil(newLower);
1591                                        else
1592                                             newLower = ceil(newLower - 0.5);
1593                                        // change may be too small - check
1594                                        if (newLower > lower || newUpper < upper) {
1595                                             if (newUpper >= newLower) {
1596                                                  // Could also tighten in this
1597                                                  //printf("%d bounds %g %g tightened to %g %g\n",
1598                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1599                                                  //     newLower,newUpper);
1600#if 1
1601                                                  columnUpper2[iColumn] = newUpper;
1602                                                  columnLower2[iColumn] = newLower;
1603                                                  columnUpper_[whichColumn[iColumn]] = newUpper;
1604                                                  columnLower_[whichColumn[iColumn]] = newLower;
1605#endif
1606                                                  // and adjust bounds on rows
1607                                                  newUpper -= upper;
1608                                                  newLower -= lower;
1609                                                  for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1610                                                       int iRow = row[j];
1611                                                       double value = element[j];
1612                                                       if (value > 0.0) {
1613                                                            up[iRow] += newUpper * value;
1614                                                            lo[iRow] += newLower * value;
1615                                                       } else {
1616                                                            lo[iRow] += newUpper * value;
1617                                                            up[iRow] += newLower * value;
1618                                                       }
1619                                                  }
1620                                             } else {
1621                                                  // infeasible
1622                                                  //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1623                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1624                                                  //     newLower,newUpper);
1625#if 1
1626                                                  delete small;
1627                                                  small = NULL;
1628                                                  break;
1629#endif
1630                                             }
1631                                        }
1632                                   }
1633                              }
1634                         }
1635                    }
1636               }
1637          }
1638     }
1639#if 0
1640     if (small) {
1641          static int which = 0;
1642          which++;
1643          char xxxx[20];
1644          sprintf(xxxx, "bad%d.mps", which);
1645          small->writeMps(xxxx, 0, 1);
1646          sprintf(xxxx, "largebad%d.mps", which);
1647          writeMps(xxxx, 0, 1);
1648          printf("bad%d %x old size %d %d new %d %d\n", which, small,
1649                 numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
1650#if 0
1651          for (int i = 0; i < numberColumns_; i++)
1652               printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1653          for (int i = 0; i < numberRows_; i++)
1654               printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
1655#endif
1656     }
1657#endif
1658#ifdef CHECK_STATUS
1659     {
1660          int n = 0;
1661          int i;
1662          for (i = 0; i < small->numberColumns(); i++)
1663               if (small->getColumnStatus(i) == ClpSimplex::basic)
1664                    n++;
1665          for (i = 0; i < small->numberRows(); i++)
1666               if (small->getRowStatus(i) == ClpSimplex::basic)
1667                    n++;
1668          assert (n == small->numberRows());
1669     }
1670#endif
1671     return small;
1672}
1673/* After very cursory presolve.
1674   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1675*/
1676void
1677ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1678                             const int * whichRow,
1679                             const int * whichColumn, int nBound)
1680{
1681#ifndef NDEBUG
1682     for (int i = 0; i < small.numberRows(); i++)
1683          assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1684     for (int i = 0; i < small.numberColumns(); i++)
1685          assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
1686#endif
1687     getbackSolution(small, whichRow, whichColumn);
1688     // and deal with status for bounds
1689     const double * element = matrix_->getElements();
1690     const int * row = matrix_->getIndices();
1691     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1692     const int * columnLength = matrix_->getVectorLengths();
1693     double tolerance = primalTolerance();
1694     double djTolerance = dualTolerance();
1695     for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1696          int iRow = whichRow[jRow];
1697          int iColumn = whichRow[jRow+numberRows_];
1698          if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1699               double lower = columnLower_[iColumn];
1700               double upper = columnUpper_[iColumn];
1701               double value = columnActivity_[iColumn];
1702               double djValue = reducedCost_[iColumn];
1703               dual_[iRow] = 0.0;
1704               if (upper > lower) {
1705                    if (value < lower + tolerance && djValue > -djTolerance) {
1706                         setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1707                         setRowStatus(iRow, ClpSimplex::basic);
1708                    } else if (value > upper - tolerance && djValue < djTolerance) {
1709                         setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1710                         setRowStatus(iRow, ClpSimplex::basic);
1711                    } else {
1712                         // has to be basic
1713                         setColumnStatus(iColumn, ClpSimplex::basic);
1714                         reducedCost_[iColumn] = 0.0;
1715                         double value = 0.0;
1716                         for (CoinBigIndex j = columnStart[iColumn];
1717                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1718                              if (iRow == row[j]) {
1719                                   value = element[j];
1720                                   break;
1721                              }
1722                         }
1723                         dual_[iRow] = djValue / value;
1724                         if (rowUpper_[iRow] > rowLower_[iRow]) {
1725                              if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1726                                        fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1727                                   setRowStatus(iRow, ClpSimplex::atLowerBound);
1728                              else
1729                                   setRowStatus(iRow, ClpSimplex::atUpperBound);
1730                         } else {
1731                              setRowStatus(iRow, ClpSimplex::isFixed);
1732                         }
1733                    }
1734               } else {
1735                    // row can always be basic
1736                    setRowStatus(iRow, ClpSimplex::basic);
1737               }
1738          } else {
1739               // row can always be basic
1740               setRowStatus(iRow, ClpSimplex::basic);
1741          }
1742     }
1743     //#ifndef NDEBUG
1744#if 0
1745     if  (small.status() == 0) {
1746          int n = 0;
1747          int i;
1748          for (i = 0; i < numberColumns; i++)
1749               if (getColumnStatus(i) == ClpSimplex::basic)
1750                    n++;
1751          for (i = 0; i < numberRows; i++)
1752               if (getRowStatus(i) == ClpSimplex::basic)
1753                    n++;
1754          assert (n == numberRows);
1755     }
1756#endif
1757}
1758/* Tightens integer bounds - returns number tightened or -1 if infeasible
1759 */
1760int
1761ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1762{
1763     // See if we can tighten any bounds
1764     // use rhs for upper and small duals for lower
1765     double * up = rhsSpace;
1766     double * lo = dual_;
1767     const double * element = matrix_->getElements();
1768     const int * row = matrix_->getIndices();
1769     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1770     const int * columnLength = matrix_->getVectorLengths();
1771     CoinZeroN(lo, numberRows_);
1772     CoinZeroN(up, numberRows_);
1773     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1774          double upper = columnUpper_[iColumn];
1775          double lower = columnLower_[iColumn];
1776          //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1777          for (CoinBigIndex j = columnStart[iColumn];
1778                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1779               int iRow = row[j];
1780               double value = element[j];
1781               if (value > 0.0) {
1782                    if (upper < 1.0e20)
1783                         up[iRow] += upper * value;
1784                    else
1785                         up[iRow] = COIN_DBL_MAX;
1786                    if (lower > -1.0e20)
1787                         lo[iRow] += lower * value;
1788                    else
1789                         lo[iRow] = -COIN_DBL_MAX;
1790               } else {
1791                    if (upper < 1.0e20)
1792                         lo[iRow] += upper * value;
1793                    else
1794                         lo[iRow] = -COIN_DBL_MAX;
1795                    if (lower > -1.0e20)
1796                         up[iRow] += lower * value;
1797                    else
1798                         up[iRow] = COIN_DBL_MAX;
1799               }
1800          }
1801     }
1802     bool feasible = true;
1803     // make safer
1804     double tolerance = primalTolerance();
1805     for (int iRow = 0; iRow < numberRows_; iRow++) {
1806          double lower = lo[iRow];
1807          if (lower > rowUpper_[iRow] + tolerance) {
1808               feasible = false;
1809               break;
1810          } else {
1811               lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
1812          }
1813          double upper = up[iRow];
1814          if (upper < rowLower_[iRow] - tolerance) {
1815               feasible = false;
1816               break;
1817          } else {
1818               up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1819          }
1820     }
1821     int numberTightened = 0;
1822     if (!feasible) {
1823          return -1;
1824     } else if (integerType_) {
1825          // and tighten
1826          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1827               if (integerType_[iColumn]) {
1828                    double upper = columnUpper_[iColumn];
1829                    double lower = columnLower_[iColumn];
1830                    double newUpper = upper;
1831                    double newLower = lower;
1832                    double difference = upper - lower;
1833                    if (lower > -1000.0 && upper < 1000.0) {
1834                         for (CoinBigIndex j = columnStart[iColumn];
1835                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1836                              int iRow = row[j];
1837                              double value = element[j];
1838                              if (value > 0.0) {
1839                                   double upWithOut = up[iRow] - value * difference;
1840                                   if (upWithOut < 0.0) {
1841                                        newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1842                                   }
1843                                   double lowWithOut = lo[iRow] + value * difference;
1844                                   if (lowWithOut > 0.0) {
1845                                        newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1846                                   }
1847                              } else {
1848                                   double upWithOut = up[iRow] + value * difference;
1849                                   if (upWithOut < 0.0) {
1850                                        newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1851                                   }
1852                                   double lowWithOut = lo[iRow] - value * difference;
1853                                   if (lowWithOut > 0.0) {
1854                                        newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1855                                   }
1856                              }
1857                         }
1858                         if (newLower > lower || newUpper < upper) {
1859                              if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1860                                   newUpper = floor(newUpper);
1861                              else
1862                                   newUpper = floor(newUpper + 0.5);
1863                              if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1864                                   newLower = ceil(newLower);
1865                              else
1866                                   newLower = ceil(newLower - 0.5);
1867                              // change may be too small - check
1868                              if (newLower > lower || newUpper < upper) {
1869                                   if (newUpper >= newLower) {
1870                                        numberTightened++;
1871                                        //printf("%d bounds %g %g tightened to %g %g\n",
1872                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1873                                        //     newLower,newUpper);
1874                                        columnUpper_[iColumn] = newUpper;
1875                                        columnLower_[iColumn] = newLower;
1876                                        // and adjust bounds on rows
1877                                        newUpper -= upper;
1878                                        newLower -= lower;
1879                                        for (CoinBigIndex j = columnStart[iColumn];
1880                                                  j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1881                                             int iRow = row[j];
1882                                             double value = element[j];
1883                                             if (value > 0.0) {
1884                                                  up[iRow] += newUpper * value;
1885                                                  lo[iRow] += newLower * value;
1886                                             } else {
1887                                                  lo[iRow] += newUpper * value;
1888                                                  up[iRow] += newLower * value;
1889                                             }
1890                                        }
1891                                   } else {
1892                                        // infeasible
1893                                        //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1894                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1895                                        //     newLower,newUpper);
1896                                        return -1;
1897                                   }
1898                              }
1899                         }
1900                    }
1901               }
1902          }
1903     }
1904     return numberTightened;
1905}
1906/* Parametrics
1907   This is an initial slow version.
1908   The code uses current bounds + theta * change (if change array not NULL)
1909   and similarly for objective.
1910   It starts at startingTheta and returns ending theta in endingTheta.
1911   If reportIncrement 0.0 it will report on any movement
1912   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1913   If it can not reach input endingTheta return code will be 1 for infeasible,
1914   2 for unbounded, if error on ranges -1,  otherwise 0.
1915   Normal report is just theta and objective but
1916   if event handler exists it may do more
1917   On exit endingTheta is maximum reached (can be used for next startingTheta)
1918*/
1919int
1920ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
1921                             const double * lowerChangeBound, const double * upperChangeBound,
1922                             const double * lowerChangeRhs, const double * upperChangeRhs,
1923                             const double * changeObjective)
1924{
1925     bool needToDoSomething = true;
1926     bool canTryQuick = (reportIncrement) ? true : false;
1927     // Save copy of model
1928     ClpSimplex copyModel = *this;
1929     int savePerturbation = perturbation_;
1930     perturbation_ = 102; // switch off
1931     while (needToDoSomething) {
1932          needToDoSomething = false;
1933          algorithm_ = -1;
1934
1935          // save data
1936          ClpDataSave data = saveData();
1937          // Dantzig
1938          ClpDualRowPivot * savePivot = dualRowPivot_;
1939          dualRowPivot_ = new ClpDualRowDantzig();
1940          dualRowPivot_->setModel(this);
1941          int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1942          int iRow, iColumn;
1943          double * chgUpper = NULL;
1944          double * chgLower = NULL;
1945          double * chgObjective = NULL;
1946
1947
1948          if (!returnCode) {
1949               // Find theta when bounds will cross over and create arrays
1950               int numberTotal = numberRows_ + numberColumns_;
1951               chgLower = new double[numberTotal];
1952               memset(chgLower, 0, numberTotal * sizeof(double));
1953               chgUpper = new double[numberTotal];
1954               memset(chgUpper, 0, numberTotal * sizeof(double));
1955               chgObjective = new double[numberTotal];
1956               memset(chgObjective, 0, numberTotal * sizeof(double));
1957               assert (!rowScale_);
1958               double maxTheta = 1.0e50;
1959               if (lowerChangeRhs || upperChangeRhs) {
1960                    for (iRow = 0; iRow < numberRows_; iRow++) {
1961                         double lower = rowLower_[iRow];
1962                         double upper = rowUpper_[iRow];
1963                         if (lower > upper) {
1964                              maxTheta = -1.0;
1965                              break;
1966                         }
1967                         double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0;
1968                         double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0;
1969                         if (lower > -1.0e20 && upper < 1.0e20) {
1970                              if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1971                                   maxTheta = (upper - lower) / (lowerChange - upperChange);
1972                              }
1973                         }
1974                         if (lower > -1.0e20) {
1975                              lower_[numberColumns_+iRow] += startingTheta * lowerChange;
1976                              chgLower[numberColumns_+iRow] = lowerChange;
1977                         }
1978                         if (upper < 1.0e20) {
1979                              upper_[numberColumns_+iRow] += startingTheta * upperChange;
1980                              chgUpper[numberColumns_+iRow] = upperChange;
1981                         }
1982                    }
1983               }
1984               if (maxTheta > 0.0) {
1985                    if (lowerChangeBound || upperChangeBound) {
1986                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1987                              double lower = columnLower_[iColumn];
1988                              double upper = columnUpper_[iColumn];
1989                              if (lower > upper) {
1990                                   maxTheta = -1.0;
1991                                   break;
1992                              }
1993                              double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0;
1994                              double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0;
1995                              if (lower > -1.0e20 && upper < 1.0e20) {
1996                                   if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1997                                        maxTheta = (upper - lower) / (lowerChange - upperChange);
1998                                   }
1999                              }
2000                              if (lower > -1.0e20) {
2001                                   lower_[iColumn] += startingTheta * lowerChange;
2002                                   chgLower[iColumn] = lowerChange;
2003                              }
2004                              if (upper < 1.0e20) {
2005                                   upper_[iColumn] += startingTheta * upperChange;
2006                                   chgUpper[iColumn] = upperChange;
2007                              }
2008                         }
2009                    }
2010                    if (maxTheta == 1.0e50)
2011                         maxTheta = COIN_DBL_MAX;
2012               }
2013               if (maxTheta < 0.0) {
2014                    // bad ranges or initial
2015                    returnCode = -1;
2016               }
2017               if (maxTheta < endingTheta) {
2018                 char line[100];
2019                 sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
2020                         endingTheta,maxTheta);
2021                 handler_->message(CLP_GENERAL,messages_)
2022                   << line << CoinMessageEol;
2023                 endingTheta = maxTheta;
2024               }
2025               if (endingTheta < startingTheta) {
2026                    // bad initial
2027                    returnCode = -2;
2028               }
2029          }
2030          double saveEndingTheta = endingTheta;
2031          if (!returnCode) {
2032               if (changeObjective) {
2033                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2034                         chgObjective[iColumn] = changeObjective[iColumn];
2035                         cost_[iColumn] += startingTheta * changeObjective[iColumn];
2036                    }
2037               }
2038               double * saveDuals = NULL;
2039               reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
2040               assert (!problemStatus_);
2041               for (int i=0;i<numberRows_+numberColumns_;i++)
2042                 setFakeBound(i, noFake);
2043               // Now do parametrics
2044               handler_->message(CLP_PARAMETRICS_STATS, messages_)
2045                 << startingTheta << objectiveValue() << CoinMessageEol;
2046               while (!returnCode) {
2047                    //assert (reportIncrement);
2048                 parametricsData paramData;
2049                 paramData.startingTheta=startingTheta;
2050                 paramData.endingTheta=endingTheta;
2051                 paramData.maxTheta=COIN_DBL_MAX;
2052                 paramData.lowerChange = chgLower;
2053                 paramData.upperChange = chgUpper;
2054                    returnCode = parametricsLoop(paramData, reportIncrement,
2055                                                 chgLower, chgUpper, chgObjective, data,
2056                                                 canTryQuick);
2057                 startingTheta=paramData.startingTheta;
2058                 endingTheta=paramData.endingTheta;
2059                    if (!returnCode) {
2060                         //double change = endingTheta-startingTheta;
2061                         startingTheta = endingTheta;
2062                         endingTheta = saveEndingTheta;
2063                         //for (int i=0;i<numberTotal;i++) {
2064                         //lower_[i] += change*chgLower[i];
2065                         //upper_[i] += change*chgUpper[i];
2066                         //cost_[i] += change*chgObjective[i];
2067                         //}
2068                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
2069                           << startingTheta << objectiveValue() << CoinMessageEol;
2070                         if (startingTheta >= endingTheta)
2071                              break;
2072                    } else if (returnCode == -1) {
2073                         // trouble - do external solve
2074                         needToDoSomething = true;
2075                    } else if (problemStatus_==1) {
2076                      // can't move any further
2077                      if (!canTryQuick) {
2078                        handler_->message(CLP_PARAMETRICS_STATS, messages_)
2079                          << endingTheta << objectiveValue() << CoinMessageEol;
2080                        problemStatus_=0;
2081                      }
2082                    } else {
2083                         abort();
2084                    }
2085               }
2086          }
2087          reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2088
2089          delete dualRowPivot_;
2090          dualRowPivot_ = savePivot;
2091          // Restore any saved stuff
2092          restoreData(data);
2093          if (needToDoSomething) {
2094               double saveStartingTheta = startingTheta; // known to be feasible
2095               int cleanedUp = 1;
2096               while (cleanedUp) {
2097                    // tweak
2098                    if (cleanedUp == 1) {
2099                         if (!reportIncrement)
2100                              startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2101                         else
2102                              startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2103                    } else {
2104                         // restoring to go slowly
2105                         startingTheta = saveStartingTheta;
2106                    }
2107                    // only works if not scaled
2108                    int i;
2109                    const double * obj1 = objective();
2110                    double * obj2 = copyModel.objective();
2111                    const double * lower1 = columnLower_;
2112                    double * lower2 = copyModel.columnLower();
2113                    const double * upper1 = columnUpper_;
2114                    double * upper2 = copyModel.columnUpper();
2115                    for (i = 0; i < numberColumns_; i++) {
2116                         obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2117                         lower2[i] = lower1[i] + startingTheta * chgLower[i];
2118                         upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2119                    }
2120                    lower1 = rowLower_;
2121                    lower2 = copyModel.rowLower();
2122                    upper1 = rowUpper_;
2123                    upper2 = copyModel.rowUpper();
2124                    for (i = 0; i < numberRows_; i++) {
2125                         lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2126                         upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2127                    }
2128                    copyModel.dual();
2129                    if (copyModel.problemStatus()) {
2130                      char line[100];
2131                      sprintf(line,"Can not get to theta of %g\n", startingTheta);
2132                      handler_->message(CLP_GENERAL,messages_)
2133                        << line << CoinMessageEol;
2134                         canTryQuick = false; // do slowly to get exact amount
2135                         // back to last known good
2136                         if (cleanedUp == 1)
2137                              cleanedUp = 2;
2138                         else
2139                              abort();
2140                    } else {
2141                         // and move stuff back
2142                         int numberTotal = numberRows_ + numberColumns_;
2143                         CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2144                         CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2145                         CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2146                         cleanedUp = 0;
2147                    }
2148               }
2149          }
2150          delete [] chgLower;
2151          delete [] chgUpper;
2152          delete [] chgObjective;
2153     }
2154     perturbation_ = savePerturbation;
2155     char line[100];
2156     sprintf(line,"Ending theta %g\n", endingTheta);
2157     handler_->message(CLP_GENERAL,messages_)
2158       << line << CoinMessageEol;
2159     return problemStatus_;
2160}
2161/* Version of parametrics which reads from file
2162   See CbcClpParam.cpp for details of format
2163   Returns -2 if unable to open file */
2164int 
2165ClpSimplexOther::parametrics(const char * dataFile)
2166{
2167  int returnCode=-2;
2168  FILE *fp = fopen(dataFile, "r");
2169  char line[200];
2170  if (!fp) {
2171    handler_->message(CLP_UNABLE_OPEN, messages_)
2172      << dataFile << CoinMessageEol;
2173    return -2;
2174  }
2175
2176  if (!fgets(line, 200, fp)) {
2177    sprintf(line,"Empty parametrics file %s?",dataFile);
2178    handler_->message(CLP_GENERAL,messages_)
2179      << line << CoinMessageEol;
2180    fclose(fp);
2181    return -2;
2182  }
2183  char * pos = line;
2184  char * put = line;
2185  while (*pos >= ' ' && *pos != '\n') {
2186    if (*pos != ' ' && *pos != '\t') {
2187      *put = static_cast<char>(tolower(*pos));
2188      put++;
2189    }
2190    pos++;
2191  }
2192  *put = '\0';
2193  pos = line;
2194  double startTheta=0.0;
2195  double endTheta=0.0;
2196  double intervalTheta=COIN_DBL_MAX;
2197  int detail=0;
2198  bool good = true;
2199  while (good) {
2200    good=false;
2201    // check ROWS
2202    char * comma = strchr(pos, ',');
2203    if (!comma)
2204      break;
2205    *comma = '\0';
2206    if (strcmp(pos,"rows"))
2207      break;
2208    *comma = ',';
2209    pos = comma+1;
2210    // check lower theta
2211    comma = strchr(pos, ',');
2212    if (!comma)
2213      break;
2214    *comma = '\0';
2215    startTheta = atof(pos);
2216    *comma = ',';
2217    pos = comma+1;
2218    // check upper theta
2219    comma = strchr(pos, ',');
2220    good=true;
2221    if (comma)
2222      *comma = '\0';
2223    endTheta = atof(pos);
2224    if (comma) {
2225      *comma = ',';
2226      pos = comma+1;
2227      comma = strchr(pos, ',');
2228      if (comma)
2229        *comma = '\0';
2230      intervalTheta = atof(pos);
2231      if (comma) {
2232        *comma = ',';
2233        pos = comma+1;
2234        comma = strchr(pos, ',');
2235        if (comma)
2236          *comma = '\0';
2237        detail = atoi(pos);
2238        if (comma) 
2239        *comma = ',';
2240      }
2241    }
2242    break;
2243  }
2244  if (good) {
2245    if (startTheta<0.0||
2246        startTheta>endTheta||
2247        intervalTheta<0.0)
2248      good=false;
2249    if (detail<0||detail>1)
2250      good=false;
2251  }
2252  if (intervalTheta>=endTheta)
2253    intervalTheta=0.0;
2254  if (!good) {
2255    sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2256    handler_->message(CLP_GENERAL,messages_)
2257      << line << CoinMessageEol;
2258    fclose(fp);
2259    return -2;
2260  }
2261  if (!fgets(line, 200, fp)) {
2262    sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2263    handler_->message(CLP_GENERAL,messages_)
2264      << line << CoinMessageEol;
2265    fclose(fp);
2266    return -2;
2267  }
2268  double * lowerRowMove = NULL;
2269  double * upperRowMove = NULL;
2270  double * lowerColumnMove = NULL;
2271  double * upperColumnMove = NULL;
2272  double * objectiveMove = NULL;
2273  char saveLine[200];
2274  saveLine[0]='\0';
2275  std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2276  int gotRow[] = { -1, -1, -1, -1, -1};
2277  int orderRow[5];
2278  assert(sizeof(gotRow) == sizeof(orderRow));
2279  int nAcross = 0;
2280  pos = line;
2281  put = line;
2282  while (*pos >= ' ' && *pos != '\n') {
2283    if (*pos != ' ' && *pos != '\t') {
2284      *put = static_cast<char>(tolower(*pos));
2285      put++;
2286    }
2287    pos++;
2288  }
2289  *put = '\0';
2290  pos = line;
2291  int i;
2292  good = true;
2293  if (strncmp(line,"column",6)) {
2294    while (pos) {
2295      char * comma = strchr(pos, ',');
2296      if (comma)
2297        *comma = '\0';
2298      for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2299        if (headingsRow[i] == pos) {
2300          if (gotRow[i] < 0) {
2301            orderRow[nAcross] = i;
2302            gotRow[i] = nAcross++;
2303          } else {
2304            // duplicate
2305            good = false;
2306          }
2307          break;
2308        }
2309      }
2310      if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2311        good = false;
2312      if (comma) {
2313        *comma = ',';
2314        pos = comma + 1;
2315      } else {
2316        break;
2317      }
2318    }
2319    if (gotRow[0] < 0 && gotRow[1] < 0)
2320      good = false;
2321    if (gotRow[0] >= 0 && gotRow[1] >= 0)
2322      good = false;
2323    if (gotRow[0] >= 0 && !lengthNames())
2324      good = false;
2325    if (gotRow[4]<0) {
2326      if (gotRow[2] < 0 && gotRow[3] >= 0)
2327        good = false;
2328      else if (gotRow[3] < 0 && gotRow[2] >= 0)
2329        good = false;
2330    } else if (gotRow[2]>=0||gotRow[3]>=0) {
2331      good = false;
2332    }
2333    if (good) {
2334      char ** rowNames = new char * [numberRows_];
2335      int iRow;
2336      for (iRow = 0; iRow < numberRows_; iRow++) {
2337        rowNames[iRow] =
2338          CoinStrdup(rowName(iRow).c_str());
2339      }
2340      lowerRowMove = new double [numberRows_];
2341      memset(lowerRowMove,0,numberRows_*sizeof(double));
2342      upperRowMove = new double [numberRows_];
2343      memset(upperRowMove,0,numberRows_*sizeof(double));
2344      int nLine = 0;
2345      int nBadLine = 0;
2346      int nBadName = 0;
2347      bool goodLine=false;
2348      while (fgets(line, 200, fp)) {
2349        goodLine=true;
2350        if (!strncmp(line, "ENDATA", 6)||
2351            !strncmp(line, "COLUMN",6))
2352          break;
2353        goodLine=false;
2354        nLine++;
2355        iRow = -1;
2356        double upper = 0.0;
2357        double lower = 0.0;
2358        char * pos = line;
2359        char * put = line;
2360        while (*pos >= ' ' && *pos != '\n') {
2361          if (*pos != ' ' && *pos != '\t') {
2362            *put = *pos;
2363            put++;
2364          }
2365          pos++;
2366        }
2367        *put = '\0';
2368        pos = line;
2369        for (int i = 0; i < nAcross; i++) {
2370          char * comma = strchr(pos, ',');
2371          if (comma) {
2372            *comma = '\0';
2373          } else if (i < nAcross - 1) {
2374            nBadLine++;
2375            break;
2376          }
2377          switch (orderRow[i]) {
2378            // name
2379          case 0:
2380            // For large problems this could be slow
2381            for (iRow = 0; iRow < numberRows_; iRow++) {
2382              if (!strcmp(rowNames[iRow], pos))
2383                break;
2384            }
2385            if (iRow == numberRows_)
2386              iRow = -1;
2387            break;
2388            // number
2389          case 1:
2390            iRow = atoi(pos);
2391            if (iRow < 0 || iRow >= numberRows_)
2392              iRow = -1;
2393            break;
2394            // lower
2395          case 2:
2396            upper = atof(pos);
2397            break;
2398            // upper
2399          case 3:
2400            lower = atof(pos);
2401            break;
2402            // rhs
2403          case 4:
2404            lower = atof(pos);
2405            upper = lower;
2406            break;
2407          }
2408          if (comma) {
2409            *comma = ',';
2410            pos = comma + 1;
2411          }
2412        }
2413        if (iRow >= 0) {
2414          if (rowLower_[iRow]>-1.0e20)
2415            lowerRowMove[iRow] = lower;
2416          else
2417            lowerRowMove[iRow]=0.0;
2418          if (rowUpper_[iRow]<1.0e20)
2419            upperRowMove[iRow] = upper;
2420          else
2421            upperRowMove[iRow] = lower;
2422        } else {
2423          nBadName++;
2424          if(saveLine[0]=='\0')
2425            strcpy(saveLine,line);
2426        }
2427      }
2428      sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2429      handler_->message(CLP_GENERAL,messages_)
2430        << line << CoinMessageEol;
2431      if (nBadName) {
2432        sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2433        handler_->message(CLP_GENERAL,messages_)
2434          << line << CoinMessageEol;
2435        returnCode=-1;
2436        good=false;
2437      }
2438      for (iRow = 0; iRow < numberRows_; iRow++) {
2439        free(rowNames[iRow]);
2440      }
2441      delete [] rowNames;
2442    } else {
2443      sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2444      handler_->message(CLP_GENERAL,messages_)
2445        << line << CoinMessageEol;
2446      returnCode=-1;
2447      good=false;
2448    }
2449  }
2450  if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
2451    if (!fgets(line, 200, fp)) {
2452      sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2453      handler_->message(CLP_GENERAL,messages_)
2454        << line << CoinMessageEol;
2455      fclose(fp);
2456      return -2;
2457    }
2458    std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2459    saveLine[0]='\0';
2460    int gotColumn[] = { -1, -1, -1, -1, -1};
2461    int orderColumn[5];
2462    assert(sizeof(gotColumn) == sizeof(orderColumn));
2463    nAcross = 0;
2464    pos = line;
2465    put = line;
2466    while (*pos >= ' ' && *pos != '\n') {
2467      if (*pos != ' ' && *pos != '\t') {
2468        *put = static_cast<char>(tolower(*pos));
2469        put++;
2470      }
2471      pos++;
2472    }
2473    *put = '\0';
2474    pos = line;
2475    int i;
2476    if (strncmp(line,"endata",6)&&good) {
2477      while (pos) {
2478        char * comma = strchr(pos, ',');
2479        if (comma)
2480          *comma = '\0';
2481        for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2482          if (headingsColumn[i] == pos) {
2483            if (gotColumn[i] < 0) {
2484              orderColumn[nAcross] = i;
2485              gotColumn[i] = nAcross++;
2486            } else {
2487              // duplicate
2488              good = false;
2489            }
2490            break;
2491          }
2492        }
2493        if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2494          good = false;
2495        if (comma) {
2496          *comma = ',';
2497          pos = comma + 1;
2498        } else {
2499          break;
2500        }
2501      }
2502      if (gotColumn[0] < 0 && gotColumn[1] < 0)
2503        good = false;
2504      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2505        good = false;
2506      if (gotColumn[0] >= 0 && !lengthNames())
2507        good = false;
2508      if (good) {
2509        char ** columnNames = new char * [numberColumns_];
2510        int iColumn;
2511        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2512          columnNames[iColumn] =
2513            CoinStrdup(columnName(iColumn).c_str());
2514        }
2515        lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2516        memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2517        upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2518        memset(upperColumnMove,0,numberColumns_*sizeof(double));
2519        objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2520        memset(objectiveMove,0,numberColumns_*sizeof(double));
2521        int nLine = 0;
2522        int nBadLine = 0;
2523        int nBadName = 0;
2524        bool goodLine=false;
2525        while (fgets(line, 200, fp)) {
2526          goodLine=true;
2527          if (!strncmp(line, "ENDATA", 6))
2528            break;
2529          goodLine=false;
2530          nLine++;
2531          iColumn = -1;
2532          double upper = 0.0;
2533          double lower = 0.0;
2534          double obj =0.0;
2535          char * pos = line;
2536          char * put = line;
2537          while (*pos >= ' ' && *pos != '\n') {
2538            if (*pos != ' ' && *pos != '\t') {
2539              *put = *pos;
2540              put++;
2541            }
2542            pos++;
2543          }
2544          *put = '\0';
2545          pos = line;
2546          for (int i = 0; i < nAcross; i++) {
2547            char * comma = strchr(pos, ',');
2548            if (comma) {
2549              *comma = '\0';
2550            } else if (i < nAcross - 1) {
2551              nBadLine++;
2552              break;
2553            }
2554            switch (orderColumn[i]) {
2555              // name
2556            case 0:
2557              // For large problems this could be slow
2558              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2559                if (!strcmp(columnNames[iColumn], pos))
2560                  break;
2561              }
2562              if (iColumn == numberColumns_)
2563                iColumn = -1;
2564              break;
2565              // number
2566            case 1:
2567              iColumn = atoi(pos);
2568              if (iColumn < 0 || iColumn >= numberColumns_)
2569                iColumn = -1;
2570              break;
2571              // lower
2572            case 2:
2573              upper = atof(pos);
2574              break;
2575              // upper
2576            case 3:
2577              lower = atof(pos);
2578              break;
2579              // objective
2580            case 4:
2581              obj = atof(pos);
2582              upper = lower;
2583              break;
2584            }
2585            if (comma) {
2586              *comma = ',';
2587              pos = comma + 1;
2588            }
2589          }
2590          if (iColumn >= 0) {
2591            if (columnLower_[iColumn]>-1.0e20)
2592              lowerColumnMove[iColumn] = lower;
2593            else
2594              lowerColumnMove[iColumn]=0.0;
2595            if (columnUpper_[iColumn]<1.0e20)
2596              upperColumnMove[iColumn] = upper;
2597            else
2598              upperColumnMove[iColumn] = lower;
2599            objectiveMove[iColumn] = obj;
2600          } else {
2601            nBadName++;
2602            if(saveLine[0]=='\0')
2603              strcpy(saveLine,line);
2604          }
2605        }
2606        sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2607        handler_->message(CLP_GENERAL,messages_)
2608          << line << CoinMessageEol;
2609        if (nBadName) {
2610          sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2611          handler_->message(CLP_GENERAL,messages_)
2612            << line << CoinMessageEol;
2613          returnCode=-1;
2614          good=false;
2615        }
2616        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2617          free(columnNames[iColumn]);
2618        }
2619        delete [] columnNames;
2620      } else {
2621        sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2622        handler_->message(CLP_GENERAL,messages_)
2623          << line << CoinMessageEol;
2624        returnCode=-1;
2625        good=false;
2626      }
2627    }
2628  }
2629  returnCode=-1;
2630  if (good) {
2631    // clean arrays
2632    if (lowerRowMove) {
2633      bool empty=true;
2634      for (int i=0;i<numberRows_;i++) {
2635        if (lowerRowMove[i]) {
2636          empty=false;
2637        break;
2638        }
2639      }
2640      if (empty) {
2641        delete [] lowerRowMove;
2642        lowerRowMove=NULL;
2643      }
2644    }
2645    if (upperRowMove) {
2646      bool empty=true;
2647      for (int i=0;i<numberRows_;i++) {
2648        if (upperRowMove[i]) {
2649          empty=false;
2650        break;
2651        }
2652      }
2653      if (empty) {
2654        delete [] upperRowMove;
2655        upperRowMove=NULL;
2656      }
2657    }
2658    if (lowerColumnMove) {
2659      bool empty=true;
2660      for (int i=0;i<numberColumns_;i++) {
2661        if (lowerColumnMove[i]) {
2662          empty=false;
2663        break;
2664        }
2665      }
2666      if (empty) {
2667        delete [] lowerColumnMove;
2668        lowerColumnMove=NULL;
2669      }
2670    }
2671    if (upperColumnMove) {
2672      bool empty=true;
2673      for (int i=0;i<numberColumns_;i++) {
2674        if (upperColumnMove[i]) {
2675          empty=false;
2676        break;
2677        }
2678      }
2679      if (empty) {
2680        delete [] upperColumnMove;
2681        upperColumnMove=NULL;
2682      }
2683    }
2684    if (objectiveMove) {
2685      bool empty=true;
2686      for (int i=0;i<numberColumns_;i++) {
2687        if (objectiveMove[i]) {
2688          empty=false;
2689        break;
2690        }
2691      }
2692      if (empty) {
2693        delete [] objectiveMove;
2694        objectiveMove=NULL;
2695      }
2696    }
2697    int saveScaling = scalingFlag_;
2698    scalingFlag_ = 0;
2699    int saveLogLevel = handler_->logLevel();
2700    if (detail>0&&!intervalTheta)
2701      handler_->setLogLevel(3);
2702    else
2703      handler_->setLogLevel(1);
2704    returnCode = parametrics(startTheta,endTheta,intervalTheta,
2705                             lowerColumnMove,upperColumnMove,
2706                             lowerRowMove,upperRowMove,
2707                             objectiveMove);
2708    scalingFlag_ = saveScaling;
2709    handler_->setLogLevel(saveLogLevel);
2710  }
2711  delete [] lowerRowMove;
2712  delete [] upperRowMove;
2713  delete [] lowerColumnMove;
2714  delete [] upperColumnMove;
2715  delete [] objectiveMove;
2716  fclose(fp);
2717  return returnCode;
2718}
2719int
2720ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
2721                                 const double * lowerChange, const double * upperChange,
2722                                 const double * changeObjective, ClpDataSave & data,
2723                                 bool canTryQuick)
2724{
2725  double startingTheta = paramData.startingTheta;
2726  double & endingTheta = paramData.endingTheta;
2727     // stuff is already at starting
2728     // For this crude version just try and go to end
2729     double change = 0.0;
2730     if (reportIncrement && canTryQuick) {
2731          endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2732          change = endingTheta - startingTheta;
2733     }
2734     int numberTotal = numberRows_ + numberColumns_;
2735     int i;
2736     for ( i = 0; i < numberTotal; i++) {
2737          lower_[i] += change * lowerChange[i];
2738          upper_[i] += change * upperChange[i];
2739          switch(getStatus(i)) {
2740
2741          case basic:
2742          case isFree:
2743          case superBasic:
2744               break;
2745          case isFixed:
2746          case atUpperBound:
2747               solution_[i] = upper_[i];
2748               break;
2749          case atLowerBound:
2750               solution_[i] = lower_[i];
2751               break;
2752          }
2753          cost_[i] += change * changeObjective[i];
2754     }
2755     problemStatus_ = -1;
2756
2757     // This says whether to restore things etc
2758     // startup will have factorized so can skip
2759     int factorType = 0;
2760     // Start check for cycles
2761     progress_.startCheck();
2762     // Say change made on first iteration
2763     changeMade_ = 1;
2764     /*
2765       Status of problem:
2766       0 - optimal
2767       1 - infeasible
2768       2 - unbounded
2769       -1 - iterating
2770       -2 - factorization wanted
2771       -3 - redo checking without factorization
2772       -4 - looks infeasible
2773     */
2774     while (problemStatus_ < 0) {
2775          int iRow, iColumn;
2776          // clear
2777          for (iRow = 0; iRow < 4; iRow++) {
2778               rowArray_[iRow]->clear();
2779          }
2780
2781          for (iColumn = 0; iColumn < 2; iColumn++) {
2782               columnArray_[iColumn]->clear();
2783          }
2784
2785          // give matrix (and model costs and bounds a chance to be
2786          // refreshed (normally null)
2787          matrix_->refresh(this);
2788          // may factorize, checks if problem finished
2789          statusOfProblemInParametrics(factorType, data);
2790          // Say good factorization
2791          factorType = 1;
2792          if (data.sparseThreshold_) {
2793               // use default at present
2794               factorization_->sparseThreshold(0);
2795               factorization_->goSparse();
2796          }
2797
2798          // exit if victory declared
2799          if (problemStatus_ >= 0 && 
2800              (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
2801               break;
2802
2803          // test for maximum iterations
2804          if (hitMaximumIterations()) {
2805               problemStatus_ = 3;
2806               break;
2807          }
2808          // Check event
2809          {
2810               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2811               if (status >= 0) {
2812                    problemStatus_ = 5;
2813                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
2814                    break;
2815               }
2816          }
2817          // Do iterations
2818          problemStatus_=-1;
2819          if (canTryQuick) {
2820               double * saveDuals = NULL;
2821               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2822          } else {
2823               whileIterating(paramData, reportIncrement,
2824                              changeObjective);
2825               startingTheta = endingTheta;
2826          }
2827     }
2828     if (!problemStatus_) {
2829          theta_ = change + startingTheta;
2830          eventHandler_->event(ClpEventHandler::theta);
2831          return 0;
2832     } else if (problemStatus_ == 10) {
2833          return -1;
2834     } else {
2835          return problemStatus_;
2836     }
2837}
2838/* Parametrics
2839   The code uses current bounds + theta * change (if change array not NULL)
2840   It starts at startingTheta and returns ending theta in endingTheta.
2841   If it can not reach input endingTheta return code will be 1 for infeasible,
2842   2 for unbounded, if error on ranges -1,  otherwise 0.
2843   Event handler may do more
2844   On exit endingTheta is maximum reached (can be used for next startingTheta)
2845*/
2846int
2847ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
2848                             const double * lowerChangeBound, const double * upperChangeBound,
2849                             const double * lowerChangeRhs, const double * upperChangeRhs)
2850{
2851  int savePerturbation = perturbation_;
2852  perturbation_ = 102; // switch off
2853  algorithm_ = -1;
2854  // extra region
2855  int maximumPivots = factorization_->maximumPivots();
2856  int numberDense = factorization_->numberDense();
2857  int length = numberRows_ + numberDense + maximumPivots;
2858  assert (!rowArray_[4]);
2859  rowArray_[4]=new CoinIndexedVector(length);
2860  assert (!rowArray_[5]);
2861  rowArray_[5]=new CoinIndexedVector(length);
2862
2863  // save data
2864  ClpDataSave data = saveData();
2865  int numberTotal = numberRows_ + numberColumns_;
2866  int ratio = (2*sizeof(int))/sizeof(double);
2867  assert (ratio==1||ratio==2);
2868  // allow for unscaled - even if not needed
2869  int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1;
2870  /*
2871    Save information and modify
2872  */
2873  double * saveLower = new double [lengthArrays];
2874  double * saveUpper = new double [lengthArrays];
2875  double * lowerCopy = saveLower+2*numberTotal;
2876  double * upperCopy = saveUpper+2*numberTotal;
2877  double * lowerChange = saveLower+numberTotal;
2878  double * upperChange = saveUpper+numberTotal;
2879  double * lowerGap = saveLower+4*numberTotal;
2880  double * lowerCoefficient = lowerGap+numberRows_;
2881  double * upperGap = saveUpper+4*numberTotal;
2882  double * upperCoefficient = upperGap+numberRows_;
2883  int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal+2*numberRows_))+2;
2884  int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal+2*numberRows_))+2;
2885  int * lowerActive = lowerList+numberTotal+1;
2886  int * upperActive = upperList+numberTotal+1;
2887  // To mark as odd
2888  char * markDone = reinterpret_cast<char *>(lowerActive+numberTotal);
2889  //memset(markDone,0,numberTotal);
2890  int * backwardBasic = upperActive+numberTotal;
2891  parametricsData paramData;
2892  paramData.lowerChange = lowerChange;
2893  paramData.lowerList=lowerList;
2894  paramData.upperChange = upperChange;
2895  paramData.upperList=upperList;
2896  paramData.markDone=markDone;
2897  paramData.backwardBasic=backwardBasic;
2898  paramData.lowerActive = lowerActive;
2899  paramData.lowerGap = lowerGap;
2900  paramData.lowerCoefficient = lowerCoefficient;
2901  paramData.upperActive = upperActive;
2902  paramData.upperGap = upperGap;
2903  paramData.upperCoefficient = upperCoefficient;
2904  // Find theta when bounds will cross over and create arrays
2905  memset(lowerChange, 0, numberTotal * sizeof(double));
2906  memset(upperChange, 0, numberTotal * sizeof(double));
2907  if (lowerChangeBound)
2908    memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
2909  if (upperChangeBound)
2910    memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
2911  if (lowerChangeRhs)
2912    memcpy(lowerChange+numberColumns_,
2913           lowerChangeRhs,numberRows_*sizeof(double));
2914  if (upperChangeRhs)
2915    memcpy(upperChange+numberColumns_,
2916           upperChangeRhs,numberRows_*sizeof(double));
2917  int nLowerChange=0;
2918  int nUpperChange=0;
2919  for (int i=0;i<numberColumns_;i++) {
2920    if (lowerChange[i]) { 
2921      lowerList[nLowerChange++]=i;
2922    }
2923    if (upperChange[i]) { 
2924      upperList[nUpperChange++]=i;
2925    }
2926  }
2927  lowerList[-2]=nLowerChange;
2928  upperList[-2]=nUpperChange;
2929  for (int i=numberColumns_;i<numberTotal;i++) {
2930    if (lowerChange[i]) { 
2931      lowerList[nLowerChange++]=i;
2932    }
2933    if (upperChange[i]) { 
2934      upperList[nUpperChange++]=i;
2935    }
2936  }
2937  lowerList[-1]=nLowerChange;
2938  upperList[-1]=nUpperChange;
2939  memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
2940  memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
2941  memcpy(lowerCopy+numberColumns_,
2942         rowLower_,numberRows_*sizeof(double));
2943  memcpy(upperCopy+numberColumns_,
2944         rowUpper_,numberRows_*sizeof(double));
2945  {
2946    //  extra for unscaled
2947    double * unscaledCopy;
2948    unscaledCopy = lowerCopy + numberTotal; 
2949    memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
2950    memcpy(unscaledCopy+numberColumns_,
2951           rowLower_,numberRows_*sizeof(double));
2952    unscaledCopy = upperCopy + numberTotal; 
2953    memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
2954    memcpy(unscaledCopy+numberColumns_,
2955           rowUpper_,numberRows_*sizeof(double));
2956  }
2957  double maxTheta = 1.0e50;
2958  double largestChange=0.0;
2959  for (int iRow = 0; iRow < numberRows_; iRow++) {
2960    double lower = rowLower_[iRow];
2961    double upper = rowUpper_[iRow];
2962    if (lower<-1.0e30)
2963      lowerChange[numberColumns_+iRow]=0.0;
2964    double chgLower = lowerChange[numberColumns_+iRow];
2965    largestChange=CoinMax(largestChange,fabs(chgLower));
2966    if (upper>1.0e30)
2967      upperChange[numberColumns_+iRow]=0.0;
2968    double chgUpper = upperChange[numberColumns_+iRow];
2969    largestChange=CoinMax(largestChange,fabs(chgUpper));
2970    if (lower > -1.0e30 && upper < 1.0e30) {
2971      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2972        maxTheta = (upper - lower) / (chgLower - chgUpper);
2973      }
2974    }
2975    lower+=startingTheta*chgLower;
2976    upper+=startingTheta*chgUpper;
2977#ifndef CLP_USER_DRIVEN
2978    if (lower > upper) {
2979      maxTheta = -1.0;
2980      break;
2981    }
2982#endif
2983    rowLower_[iRow]=lower;
2984    rowUpper_[iRow]=upper;
2985  }
2986  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2987    double lower = columnLower_[iColumn];
2988    double upper = columnUpper_[iColumn];
2989    if (lower<-1.0e30)
2990      lowerChange[iColumn]=0.0;
2991    double chgLower = lowerChange[iColumn];
2992    largestChange=CoinMax(largestChange,fabs(chgLower));
2993    if (upper>1.0e30)
2994      upperChange[iColumn]=0.0;
2995    double chgUpper = upperChange[iColumn];
2996    largestChange=CoinMax(largestChange,fabs(chgUpper));
2997    if (lower > -1.0e30 && upper < 1.0e30) {
2998      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2999        maxTheta = (upper - lower) / (chgLower - chgUpper);
3000      }
3001    }
3002    lower+=startingTheta*chgLower;
3003    upper+=startingTheta*chgUpper;
3004#ifndef CLP_USER_DRIVEN
3005    if (lower > upper) {
3006      maxTheta = -1.0;
3007      break;
3008    }
3009#endif
3010    columnLower_[iColumn]=lower;
3011    columnUpper_[iColumn]=upper;
3012  }
3013  if (maxTheta == 1.0e50)
3014    maxTheta = COIN_DBL_MAX;
3015  int returnCode=0;
3016  /*
3017    If user is in charge then they are sophisticated enough
3018    to know what they are doing.  They must be doing something clever
3019    with event handlers!
3020   */
3021#ifndef CLP_USER_DRIVEN
3022  if (maxTheta < 0.0) {
3023    // bad ranges or initial
3024    returnCode = -1;
3025  }
3026  if (maxTheta < endingTheta) {
3027    char line[100];
3028    sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
3029            endingTheta,maxTheta);
3030    handler_->message(CLP_GENERAL,messages_)
3031      << line << CoinMessageEol;
3032    endingTheta = maxTheta;
3033  }
3034  if (endingTheta < startingTheta) {
3035    // bad initial
3036    returnCode = -2;
3037  }
3038  paramData.maxTheta=CoinMin(maxTheta,endingTheta);
3039#else
3040  // accept user's version
3041  paramData.maxTheta=endingTheta;
3042#endif
3043  /* given largest change element choose acceptable end
3044     be safe and make sure difference < 0.1*tolerance */
3045  double acceptableDifference=0.1*primalTolerance_/
3046    CoinMax(largestChange,1.0);
3047  paramData.acceptableMaxTheta=maxTheta-acceptableDifference;
3048  bool swapped=false;
3049  // Dantzig
3050#define ALL_DANTZIG
3051#ifdef ALL_DANTZIG
3052  ClpDualRowPivot * savePivot = dualRowPivot_;
3053  dualRowPivot_ = new ClpDualRowDantzig();
3054  dualRowPivot_->setModel(this);
3055#else
3056  ClpDualRowPivot * savePivot = NULL;
3057#endif
3058  if (!returnCode) {
3059    assert (objective_->type()==1);
3060    objective_->setType(2); // in case matrix empty
3061    returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
3062    objective_->setType(1);
3063    if (!returnCode) {
3064      double saveDualBound=dualBound_;
3065      dualBound_=CoinMax(dualBound_,1.0e15);
3066      swapped=true;
3067      double * temp;
3068      memcpy(saveLower,lower_,numberTotal*sizeof(double));
3069      temp=saveLower;
3070      saveLower=lower_;
3071      lower_=temp;
3072      //columnLowerWork_ = lower_;
3073      //rowLowerWork_ = lower_ + numberColumns_;
3074      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
3075      temp=saveUpper;
3076      saveUpper=upper_;
3077      upper_=temp;
3078      //columnUpperWork_ = upper_;
3079      //rowUpperWork_ = upper_ + numberColumns_;
3080      if (rowScale_) {
3081        // scale saved and change arrays
3082        double * lowerChange = lower_+numberTotal;
3083        double * upperChange = upper_+numberTotal;
3084        double * lowerSave = lowerChange+numberTotal;
3085        double * upperSave = upperChange+numberTotal;
3086        for (int i=0;i<numberColumns_;i++) {
3087          double multiplier = inverseColumnScale_[i];
3088          if (lowerSave[i]>-1.0e20)
3089            lowerSave[i] *= multiplier;
3090          if (upperSave[i]<1.0e20)
3091            upperSave[i] *= multiplier;
3092          lowerChange[i] *= multiplier;
3093          upperChange[i] *= multiplier;
3094        }
3095        lowerChange += numberColumns_;
3096        upperChange += numberColumns_;
3097        lowerSave += numberColumns_;
3098        upperSave += numberColumns_;
3099        for (int i=0;i<numberRows_;i++) {
3100          double multiplier = rowScale_[i];
3101          if (lowerSave[i]>-1.0e20)
3102            lowerSave[i] *= multiplier;
3103          if (upperSave[i]<1.0e20)
3104            upperSave[i] *= multiplier;
3105          lowerChange[i] *= multiplier;
3106          upperChange[i] *= multiplier;
3107        }
3108      }
3109      //double saveEndingTheta = endingTheta;
3110      double * saveDuals = NULL;
3111      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3112      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3113        // probably can get rid of this if we adjust every change in theta
3114        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3115        //   sumPrimalInfeasibilities_);
3116        int pass=100;
3117        while(sumPrimalInfeasibilities_) {
3118          pass--;
3119          if (!pass)
3120            break;
3121          problemStatus_=-1;
3122          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3123            double value=solution_[iSequence];
3124            // remember scaling
3125            if (value<lower_[iSequence]-1.0e-9) {
3126              lower_[iSequence]=value;
3127              lowerCopy[iSequence]=value;
3128            } else if (value>upper_[iSequence]+1.0e-9) {
3129              upper_[iSequence]=value;
3130              upperCopy[iSequence]=value;
3131            }
3132          }
3133          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3134        }
3135      }
3136      if (!problemStatus_) {
3137        if (nLowerChange||nUpperChange) {
3138#ifndef ALL_DANTZIG
3139          // Dantzig
3140          savePivot = dualRowPivot_;
3141          dualRowPivot_ = new ClpDualRowDantzig();
3142          dualRowPivot_->setModel(this);
3143#endif
3144          //for (int i=0;i<numberRows_+numberColumns_;i++)
3145          //setFakeBound(i, noFake);
3146          // Now do parametrics
3147          handler_->message(CLP_PARAMETRICS_STATS, messages_)
3148            << startingTheta << objectiveValue() << CoinMessageEol;
3149          bool canSkipFactorization=true;
3150          while (!returnCode) {
3151            paramData.startingTheta=startingTheta;
3152            paramData.endingTheta=endingTheta;
3153            returnCode = parametricsLoop(paramData,
3154                                         data,canSkipFactorization);
3155            startingTheta=paramData.startingTheta;
3156            endingTheta=paramData.endingTheta;
3157            canSkipFactorization=false;
3158            if (!returnCode) {
3159              //startingTheta = endingTheta;
3160              //endingTheta = saveEndingTheta;
3161              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3162                << startingTheta << objectiveValue() << CoinMessageEol;
3163              if (startingTheta >= endingTheta-primalTolerance_
3164                  ||problemStatus_==2)
3165                break;
3166            } else if (returnCode == -1) {
3167              // trouble - do external solve
3168              abort(); //needToDoSomething = true;
3169            } else if (problemStatus_==1) {
3170              // can't move any further
3171              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3172                << endingTheta << objectiveValue() << CoinMessageEol;
3173              problemStatus_=0;
3174            }
3175          }
3176        }
3177        dualBound_ = saveDualBound;
3178        //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3179      } else {
3180        // check if empty
3181        //if (!numberRows_||!matrix_->getNumElements()) {
3182        // success
3183#ifdef CLP_USER_DRIVEN
3184        //theta_ = endingTheta;
3185        //eventHandler_->event(ClpEventHandler::theta);
3186#endif
3187        //}
3188      }
3189    }
3190    if (problemStatus_==2) {
3191      delete [] ray_;
3192      ray_ = new double [numberColumns_];
3193    }
3194    if (swapped&&lower_) {
3195      double * temp=saveLower;
3196      saveLower=lower_;
3197      lower_=temp;
3198      temp=saveUpper;
3199      saveUpper=upper_;
3200      upper_=temp;
3201    }
3202    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
3203  }   
3204  if (!scalingFlag_) {
3205    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3206    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3207    memcpy(rowLower_,lowerCopy+numberColumns_,
3208           numberRows_*sizeof(double));
3209    memcpy(rowUpper_,upperCopy+numberColumns_,
3210           numberRows_*sizeof(double));
3211  } else {
3212    //  extra for unscaled
3213    double * unscaledCopy;
3214    unscaledCopy = lowerCopy + numberTotal; 
3215    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3216    memcpy(rowLower_,unscaledCopy+numberColumns_,
3217           numberRows_*sizeof(double));
3218    unscaledCopy = upperCopy + numberTotal; 
3219    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3220    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3221           numberRows_*sizeof(double));
3222  }
3223  delete [] saveLower;
3224  delete [] saveUpper;
3225#ifdef ALL_DANTZIG
3226  if (savePivot) {
3227#endif
3228    delete dualRowPivot_;
3229    dualRowPivot_ = savePivot;
3230#ifdef ALL_DANTZIG
3231  }
3232#endif
3233  // Restore any saved stuff
3234  restoreData(data);
3235  perturbation_ = savePerturbation;
3236  delete rowArray_[4];
3237  rowArray_[4]=NULL;
3238  delete rowArray_[5];
3239  rowArray_[5]=NULL;
3240  char line[100];
3241  sprintf(line,"Ending theta %g\n", endingTheta);
3242  handler_->message(CLP_GENERAL,messages_)
3243    << line << CoinMessageEol;
3244  return problemStatus_;
3245}
3246int
3247ClpSimplexOther::parametricsLoop(parametricsData & paramData,
3248                                 ClpDataSave & data,bool canSkipFactorization)
3249{
3250  double & startingTheta = paramData.startingTheta;
3251  double & endingTheta = paramData.endingTheta;
3252  int numberTotal = numberRows_+numberColumns_;
3253  // stuff is already at starting
3254  const int * lowerList = paramData.lowerList;
3255  const int * upperList = paramData.upperList;
3256  problemStatus_ = -1;
3257  //double saveEndingTheta=endingTheta;
3258
3259  // This says whether to restore things etc
3260  // startup will have factorized so can skip
3261  int factorType = 0;
3262  // Start check for cycles
3263  progress_.startCheck();
3264  // Say change made on first iteration
3265  changeMade_ = 1;
3266  /*
3267    Status of problem:
3268    0 - optimal
3269    1 - infeasible
3270    2 - unbounded
3271    -1 - iterating
3272    -2 - factorization wanted
3273    -3 - redo checking without factorization
3274    -4 - looks infeasible
3275  */
3276  while (problemStatus_ < 0) {
3277    int iRow, iColumn;
3278    // clear
3279    for (iRow = 0; iRow < 6; iRow++) {
3280      rowArray_[iRow]->clear();
3281    }
3282   
3283    for (iColumn = 0; iColumn < 2; iColumn++) {
3284      columnArray_[iColumn]->clear();
3285    }
3286   
3287    // give matrix (and model costs and bounds a chance to be
3288    // refreshed (normally null)
3289    matrix_->refresh(this);
3290    // may factorize, checks if problem finished
3291    if (!canSkipFactorization)
3292      statusOfProblemInParametrics(factorType, data);
3293    canSkipFactorization=false;
3294    if (numberPrimalInfeasibilities_) {
3295      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3296        // treat as success
3297        problemStatus_=0;
3298        endingTheta=startingTheta;
3299        break;
3300      }
3301      // probably can get rid of this if we adjust every change in theta
3302      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3303      //     sumPrimalInfeasibilities_);
3304      const double * lowerChange = lower_+numberTotal;
3305      const double * upperChange = upper_+numberTotal;
3306      const double * startLower = lowerChange+numberTotal;
3307      const double * startUpper = upperChange+numberTotal;
3308      //startingTheta -= 1.0e-7;
3309      int nLowerChange = lowerList[-1];
3310      for (int i = 0; i < nLowerChange; i++) {
3311        int iSequence = lowerList[i];
3312        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3313      }
3314      int nUpperChange = upperList[-1];
3315      for (int i = 0; i < nUpperChange; i++) {
3316        int iSequence = upperList[i];
3317        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3318      }
3319      // adjust rhs in case dual uses
3320      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3321      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3322      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3323      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3324      if (rowScale_) {
3325        for (int i=0;i<numberColumns_;i++) {
3326          double multiplier = columnScale_[i];
3327          if (columnLower_[i]>-1.0e20)
3328            columnLower_[i] *= multiplier;
3329          if (columnUpper_[i]<1.0e20)
3330            columnUpper_[i] *= multiplier;
3331        }
3332        for (int i=0;i<numberRows_;i++) {
3333          double multiplier = inverseRowScale_[i];
3334          if (rowLower_[i]>-1.0e20)
3335            rowLower_[i] *= multiplier;
3336          if (rowUpper_[i]<1.0e20)
3337            rowUpper_[i] *= multiplier;
3338        }
3339      }
3340      double * saveDuals = NULL;
3341      problemStatus_=-1;
3342      ClpObjective * saveObjective = objective_;
3343      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3344      if (saveObjective!=objective_) {
3345        delete objective_;
3346        objective_=saveObjective;
3347      }
3348      int pass=100;
3349      double moved=0.0;
3350      while(sumPrimalInfeasibilities_) {
3351        //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_,
3352        //     sumPrimalInfeasibilities_);
3353        pass--;
3354        if (!pass)
3355          break;
3356        problemStatus_=-1;
3357        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3358          double value=solution_[iSequence];
3359          if (value<lower_[iSequence]-1.0e-9) {
3360            moved += lower_[iSequence]-value;
3361            lower_[iSequence]=value;
3362          } else if (value>upper_[iSequence]+1.0e-9) {
3363            moved += upper_[iSequence]-value;
3364            upper_[iSequence]=value;
3365          }
3366        }
3367        if (!moved) {
3368          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3369            double value=solution_[iSequence];
3370            if (value<lower_[iSequence]-1.0e-9) {
3371              moved += lower_[iSequence]-value;
3372              lower_[iSequence]=value;
3373            } else if (value>upper_[iSequence]+1.0e-9) {
3374              moved += upper_[iSequence]-value;
3375              upper_[iSequence]=value;
3376            }
3377          }
3378        }
3379        assert (moved);
3380        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3381      }
3382      // adjust
3383      //printf("Should adjust - moved %g\n",moved);
3384    }
3385    // Say good factorization
3386    factorType = 1;
3387    if (data.sparseThreshold_) {
3388      // use default at present
3389      factorization_->sparseThreshold(0);
3390      factorization_->goSparse();
3391    }
3392   
3393    // exit if victory declared
3394    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3395      break;
3396   
3397    // test for maximum iterations
3398    if (hitMaximumIterations()) {
3399      problemStatus_ = 3;
3400      break;
3401    }
3402#ifdef CLP_USER_DRIVEN
3403    // Check event
3404    {
3405      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3406      if (status >= 0) {
3407        problemStatus_ = 5;
3408        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3409        break;
3410      }
3411    }
3412#endif
3413    // Do iterations
3414    problemStatus_=-1;
3415    whileIterating(paramData, 0.0,
3416                   NULL);
3417    //startingTheta = endingTheta;
3418    //endingTheta = saveEndingTheta;
3419  }
3420  if (!problemStatus_/*||problemStatus_==2*/) {
3421    theta_ = endingTheta;
3422#ifdef CLP_USER_DRIVEN
3423    {
3424      double saveTheta=theta_;
3425      theta_ = endingTheta;
3426      int status=eventHandler_->event(ClpEventHandler::theta);
3427      if (status>=0&&status<10) {
3428        endingTheta=theta_;
3429        theta_=saveTheta;
3430        problemStatus_=-1;
3431      } else {
3432        if (status>=10) {
3433          problemStatus_=status-10;
3434          startingTheta=endingTheta;
3435        }
3436        theta_=saveTheta;
3437      }
3438    }
3439#endif
3440    return 0;
3441  } else if (problemStatus_ == 10) {
3442    return -1;
3443  } else {
3444    return problemStatus_;
3445  }
3446}
3447/* Checks if finished.  Updates status */
3448void
3449ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3450{
3451     if (type == 2) {
3452          // trouble - go to recovery
3453          problemStatus_ = 10;
3454          return;
3455     }
3456     if (problemStatus_ > -3 || factorization_->pivots()) {
3457          // factorize
3458          // later on we will need to recover from singularities
3459          // also we could skip if first time
3460          if (type) {
3461               // is factorization okay?
3462               if (internalFactorize(1)) {
3463                    // trouble - go to recovery
3464                    problemStatus_ = 10;
3465                    return;
3466               }
3467          }
3468          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3469               problemStatus_ = -3;
3470     }
3471     // at this stage status is -3 or -4 if looks infeasible
3472     // get primal and dual solutions
3473     gutsOfSolution(NULL, NULL);
3474     double realDualInfeasibilities = sumDualInfeasibilities_;
3475     // If bad accuracy treat as singular
3476     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3477          // trouble - go to recovery
3478          problemStatus_ = 10;
3479          return;
3480     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3481          // Can reduce tolerance
3482          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3483          factorization_->pivotTolerance(newTolerance);
3484     }
3485     // Check if looping
3486     int loop;
3487     if (type != 2)
3488          loop = progress_.looping();
3489     else
3490          loop = -1;
3491     if (loop >= 0) {
3492          problemStatus_ = loop; //exit if in loop
3493          if (!problemStatus_) {
3494               // declaring victory
3495               numberPrimalInfeasibilities_ = 0;
3496               sumPrimalInfeasibilities_ = 0.0;
3497          } else {
3498               problemStatus_ = 10; // instead - try other algorithm
3499          }
3500          return;
3501     } else if (loop < -1) {
3502          // something may have changed
3503          gutsOfSolution(NULL, NULL);
3504     }
3505     progressFlag_ = 0; //reset progress flag
3506     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3507          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3508                    << numberIterations_ << objectiveValue();
3509          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3510                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3511          handler_->printing(sumDualInfeasibilities_ > 0.0)
3512                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3513          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3514                             < numberDualInfeasibilities_)
3515                    << numberDualInfeasibilitiesWithoutFree_;
3516          handler_->message() << CoinMessageEol;
3517     }
3518#ifdef CLP_USER_DRIVEN
3519     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3520       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3521       if (status>=0) {
3522         // fix up
3523         for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3524           double value=solution_[iSequence];
3525           if (value<=lower_[iSequence]-primalTolerance_) {
3526             lower_[iSequence]=value;
3527           } else if (value>=upper_[iSequence]+primalTolerance_) {
3528             upper_[iSequence]=value;
3529           }
3530         }
3531         numberPrimalInfeasibilities_ = 0;
3532         sumPrimalInfeasibilities_ = 0.0;
3533       }
3534     }
3535#endif
3536     /* If we are primal feasible and any dual infeasibilities are on
3537        free variables then it is better to go to primal */
3538     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3539               numberDualInfeasibilities_) {
3540          problemStatus_ = 10;
3541          return;
3542     }
3543
3544     // check optimal
3545     // give code benefit of doubt
3546     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3547               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3548          // say optimal (with these bounds etc)
3549          numberDualInfeasibilities_ = 0;
3550          sumDualInfeasibilities_ = 0.0;
3551          numberPrimalInfeasibilities_ = 0;
3552          sumPrimalInfeasibilities_ = 0.0;
3553     }
3554     if (dualFeasible() || problemStatus_ == -4) {
3555          progress_.modifyObjective(objectiveValue_
3556                                    - sumDualInfeasibilities_ * dualBound_);
3557     }
3558     if (numberPrimalInfeasibilities_) {
3559          if (problemStatus_ == -4 || problemStatus_ == -5) {
3560               problemStatus_ = 1; // infeasible
3561          }
3562     } else if (numberDualInfeasibilities_) {
3563          // clean up
3564          problemStatus_ = 10;
3565     } else {
3566          problemStatus_ = 0;
3567     }
3568     lastGoodIteration_ = numberIterations_;
3569     if (problemStatus_ < 0) {
3570          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3571          if (sumDualInfeasibilities_)
3572               numberDualInfeasibilities_ = 1;
3573     }
3574     // Allow matrices to be sorted etc
3575     int fake = -999; // signal sort
3576     matrix_->correctSequence(this, fake, fake);
3577}
3578//static double lastThetaX=0.0;
3579/* This has the flow between re-factorizations
3580   Reasons to come out:
3581   -1 iterations etc
3582   -2 inaccuracy
3583   -3 slight inaccuracy (and done iterations)
3584   +0 looks optimal (might be unbounded - but we will investigate)
3585   +1 looks infeasible
3586   +3 max iterations
3587   +4 accuracy problems
3588*/
3589int
3590ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
3591                                const double * /*changeObjective*/)
3592{
3593  double & startingTheta = paramData.startingTheta;
3594  double & endingTheta = paramData.endingTheta;
3595  const double * lowerChange = paramData.lowerChange;
3596  const double * upperChange = paramData.upperChange;
3597  int numberTotal = numberColumns_ + numberRows_;
3598  const int * lowerList = paramData.lowerList;
3599  const int * upperList = paramData.upperList;
3600  //#define CLP_PARAMETRIC_DENSE_ARRAYS 2
3601#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3602  double * lowerGap = paramData.lowerGap;
3603  double * upperGap = paramData.upperGap;
3604  double * lowerCoefficient = paramData.lowerCoefficient;
3605  double * upperCoefficient = paramData.upperCoefficient;
3606#endif
3607  // do basic pointers
3608  int * backwardBasic = paramData.backwardBasic;
3609  for (int i=0;i<numberTotal;i++)
3610    backwardBasic[i]=-1;
3611  for (int i=0;i<numberRows_;i++) {
3612    int iPivot=pivotVariable_[i];
3613    backwardBasic[iPivot]=i;
3614  }
3615     {
3616          int i;
3617          for (i = 0; i < 4; i++) {
3618               rowArray_[i]->clear();
3619          }
3620          for (i = 0; i < 2; i++) {
3621               columnArray_[i]->clear();
3622          }
3623     }
3624     // if can't trust much and long way from optimal then relax
3625     if (largestPrimalError_ > 10.0)
3626          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3627     else
3628          factorization_->relaxAccuracyCheck(1.0);
3629     // status stays at -1 while iterating, >=0 finished, -2 to invert
3630     // status -3 to go to top without an invert
3631     int returnCode = -1;
3632     double lastTheta = startingTheta;
3633     double useTheta = startingTheta;
3634     while (problemStatus_ == -1) {
3635          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3636          // Get theta for bounds - we know can't crossover
3637          int pivotType = nextTheta(1, increaseTheta, paramData,
3638                                     NULL);
3639          useTheta += theta_;
3640          double change = useTheta - lastTheta;
3641          if (change>1.0e-14) {
3642            int n;
3643            n=lowerList[-1];
3644            for (int i=0;i<n;i++) {
3645              int iSequence = lowerList[i];
3646              double thisChange = change * lowerChange[iSequence];
3647              double newValue = lower_[iSequence] + thisChange;
3648              lower_[iSequence] = newValue;
3649#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3650              if (getStatus(iSequence)==basic) {
3651                int iRow=backwardBasic[iSequence];
3652                lowerGap[iRow] -= thisChange;
3653              } else if(getStatus(iSequence)==atLowerBound) {
3654                solution_[iSequence] = newValue;
3655              }
3656#else
3657              if(getStatus(iSequence)==atLowerBound) {
3658                solution_[iSequence] = newValue;
3659              }
3660#endif
3661#if 0
3662              // may have to adjust other bound
3663              double otherValue = upper_[iSequence];
3664              if (otherValue-newValue<dualBound_) {
3665                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3666                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3667                //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5);
3668              }
3669#endif
3670            }
3671            n=upperList[-1];
3672            for (int i=0;i<n;i++) {
3673              int iSequence = upperList[i];
3674              double thisChange = change * upperChange[iSequence];
3675              double newValue = upper_[iSequence] + thisChange;
3676              upper_[iSequence] = newValue;
3677              if(getStatus(iSequence)==atUpperBound||
3678                 getStatus(iSequence)==isFixed) {
3679                solution_[iSequence] = newValue;
3680#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3681              } else if (getStatus(iSequence)==basic) {
3682                int iRow=backwardBasic[iSequence];
3683                upperGap[iRow] += thisChange;
3684#endif
3685              }
3686              // may have to adjust other bound
3687              double otherValue = lower_[iSequence];
3688              if (newValue-otherValue<dualBound_) {
3689                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3690                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3691                //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5);
3692              }
3693            }
3694          }
3695          sequenceIn_=-1;
3696          if (pivotType) {
3697            if (useTheta>lastTheta+1.0e-9) {
3698              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3699                << useTheta << objectiveValue() << CoinMessageEol;
3700              lastTheta = useTheta;
3701            }
3702            problemStatus_ = -2;
3703            if (!factorization_->pivots()&&pivotRow_<0)
3704              problemStatus_=2;
3705#ifdef CLP_USER_DRIVEN
3706            {
3707              double saveTheta=theta_;
3708              theta_ = endingTheta;
3709              if (problemStatus_==2&&theta_>paramData.acceptableMaxTheta)
3710                theta_=COIN_DBL_MAX; // we have finished
3711              int status=eventHandler_->event(ClpEventHandler::theta);
3712              if (status>=0&&status<10) {
3713                endingTheta=theta_;
3714                problemStatus_=-1;
3715                continue;
3716              } else {
3717                if (status>=10)
3718                  problemStatus_=status-10;
3719                if (status<0)
3720                  startingTheta = useTheta;
3721              }
3722              theta_=saveTheta;
3723            }
3724#else
3725            startingTheta = useTheta;
3726#endif
3727            return 4;
3728          }
3729          // choose row to go out
3730          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3731          if (pivotRow_ >= 0) {
3732               // we found a pivot row
3733               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3734                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3735                              << pivotRow_
3736                              << CoinMessageEol;
3737               }
3738               // check accuracy of weights
3739               dualRowPivot_->checkAccuracy();
3740               // do ratio test for normal iteration
3741               double bestPossiblePivot = bestPivot();
3742               if (sequenceIn_ >= 0) {
3743                    // normal iteration
3744                    // update the incoming column
3745                    double btranAlpha = -alpha_ * directionOut_; // for check
3746#ifndef COIN_FAC_NEW
3747                    unpackPacked(rowArray_[1]);
3748#else
3749                    unpack(rowArray_[1]);
3750#endif
3751                    // and update dual weights (can do in parallel - with extra array)
3752                    rowArray_[2]->clear();
3753                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3754                                                          rowArray_[2],
3755                                                          rowArray_[3],
3756                                                          rowArray_[1]);
3757                    // see if update stable
3758#ifdef CLP_DEBUG
3759                    if ((handler_->logLevel() & 32))
3760                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3761#endif
3762                    double checkValue = 1.0e-7;
3763                    // if can't trust much and long way from optimal then relax
3764                    if (largestPrimalError_ > 10.0)
3765                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3766                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3767                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3768                         handler_->message(CLP_DUAL_CHECK, messages_)
3769                                   << btranAlpha
3770                                   << alpha_
3771                                   << CoinMessageEol;
3772                         // clear arrays
3773                         rowArray_[4]->clear();
3774                         if (factorization_->pivots()) {
3775                              dualRowPivot_->unrollWeights();
3776                              problemStatus_ = -2; // factorize now
3777                              rowArray_[0]->clear();
3778                              rowArray_[1]->clear();
3779                              columnArray_[0]->clear();
3780                              returnCode = -2;
3781                              break;
3782                         } else {
3783                              // take on more relaxed criterion
3784                              double test;
3785                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3786                                   test = 1.0e-1 * fabs(alpha_);
3787                              else
3788                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
3789                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3790                                        fabs(btranAlpha - alpha_) > test) {
3791                                   dualRowPivot_->unrollWeights();
3792                                   // need to reject something
3793                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
3794                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
3795                                             << x << sequenceWithin(sequenceOut_)
3796                                             << CoinMessageEol;
3797                                   setFlagged(sequenceOut_);
3798                                   progress_.clearBadTimes();
3799                                   lastBadIteration_ = numberIterations_; // say be more cautious
3800                                   rowArray_[0]->clear();
3801                                   rowArray_[1]->clear();
3802                                   columnArray_[0]->clear();
3803                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3804                                        //printf("I think should declare infeasible\n");
3805                                        problemStatus_ = 1;
3806                                        returnCode = 1;
3807                                        break;
3808                                   }
3809                                   continue;
3810                              }
3811                         }
3812                    }
3813                    // update duals BEFORE replaceColumn so can do updateColumn
3814                    double objectiveChange = 0.0;
3815                    // do duals first as variables may flip bounds
3816                    // rowArray_[0] and columnArray_[0] may have flips
3817                    // so use rowArray_[3] for work array from here on
3818                    int nswapped = 0;
3819                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3820                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3821#if CLP_CAN_HAVE_ZERO_OBJ
3822                    if ((specialOptions_&2097152)==0) {
3823#endif
3824                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3825                                                                                               rowArray_[2], theta_,
3826                                                                                               objectiveChange, false);
3827                      assert (!nswapped);
3828#if CLP_CAN_HAVE_ZERO_OBJ
3829                    } else {
3830                      rowArray_[0]->clear();
3831                      rowArray_[2]->clear();
3832                      columnArray_[0]->clear();
3833                    }
3834#endif
3835                    // which will change basic solution
3836                    if (nswapped) {
3837                      abort(); //needs testing
3838                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3839                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
3840                                                             1.0, objectiveChange);
3841                         // recompute dualOut_
3842                         valueOut_ = solution_[sequenceOut_];
3843                         if (directionOut_ < 0) {
3844                              dualOut_ = valueOut_ - upperOut_;
3845                         } else {
3846                              dualOut_ = lowerOut_ - valueOut_;
3847                         }
3848                    }
3849                    // amount primal will move
3850                    double movement = -dualOut_ * directionOut_ / alpha_;
3851                    // so objective should increase by fabs(dj)*movement
3852                    // but we already have objective change - so check will be good
3853                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3854#ifdef CLP_DEBUG
3855                         if (handler_->logLevel() & 32)
3856                              printf("movement %g, swap change %g, rest %g  * %g\n",
3857                                     objectiveChange + fabs(movement * dualIn_),
3858                                     objectiveChange, movement, dualIn_);
3859#endif
3860                         assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3861                         if(factorization_->pivots()) {
3862                              // going backwards - factorize
3863                              dualRowPivot_->unrollWeights();
3864                              problemStatus_ = -2; // factorize now
3865                              returnCode = -2;
3866                              break;
3867                         }
3868                    }
3869                    CoinAssert(fabs(dualOut_) < 1.0e50);
3870                    // if stable replace in basis
3871                    int updateStatus = factorization_->replaceColumn(this,
3872                                       rowArray_[2],
3873                                       rowArray_[1],
3874                                       pivotRow_,
3875                                       alpha_);
3876                    // if no pivots, bad update but reasonable alpha - take and invert
3877                    if (updateStatus == 2 &&
3878                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3879                         updateStatus = 4;
3880                    if (updateStatus == 1 || updateStatus == 4) {
3881                         // slight error
3882                         if (factorization_->pivots() > 5 || updateStatus == 4) {
3883                              problemStatus_ = -2; // factorize now
3884                              returnCode = -3;
3885                         }
3886                    } else if (updateStatus == 2) {
3887                         // major error
3888                         dualRowPivot_->unrollWeights();
3889                         // later we may need to unwind more e.g. fake bounds
3890                         if (factorization_->pivots()) {
3891                              problemStatus_ = -2; // factorize now
3892                              returnCode = -2;
3893                              break;
3894                         } else {
3895                              // need to reject something
3896                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3897                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3898                                        << x << sequenceWithin(sequenceOut_)
3899                                        << CoinMessageEol;
3900                              setFlagged(sequenceOut_);
3901                              progress_.clearBadTimes();
3902                              lastBadIteration_ = numberIterations_; // say be more cautious
3903                              rowArray_[0]->clear();
3904                              rowArray_[1]->clear();
3905                              columnArray_[0]->clear();
3906                              // make sure dual feasible
3907                              // look at all rows and columns
3908                              double objectiveChange = 0.0;
3909                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3910                                        0.0, objectiveChange, true);
3911                              continue;
3912                         }
3913                    } else if (updateStatus == 3) {
3914                         // out of memory
3915                         // increase space if not many iterations
3916                         if (factorization_->pivots() <
3917                                   0.5 * factorization_->maximumPivots() &&
3918                                   factorization_->pivots() < 200)
3919                              factorization_->areaFactor(
3920                                   factorization_->areaFactor() * 1.1);
3921                         problemStatus_ = -2; // factorize now
3922                    } else if (updateStatus == 5) {
3923                         problemStatus_ = -2; // factorize now
3924                    }
3925                    int * lowerActive = paramData.lowerActive;
3926                    int * upperActive = paramData.upperActive;
3927                    // update change vector
3928                    {
3929                      double * work = rowArray_[1]->denseVector();
3930                      int number = rowArray_[1]->getNumElements();
3931                      int * which = rowArray_[1]->getIndices();
3932                      assert (!rowArray_[4]->packedMode());
3933#ifndef COIN_FAC_NEW
3934                      assert (rowArray_[1]->packedMode());
3935#else
3936                      assert (!rowArray_[1]->packedMode());
3937#endif
3938                      double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3939                      double multiplier = -pivotValue/alpha_;
3940                      double * array=rowArray_[4]->denseVector();
3941#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3942                      int lowerN=lowerActive[-1];
3943                      int upperN=upperActive[-1];
3944#endif
3945                      if (multiplier) {
3946                        for (int i = 0; i < number; i++) {
3947                          int iRow = which[i];
3948#ifndef COIN_FAC_NEW
3949                          double alpha=multiplier*work[i];
3950#else
3951                          double alpha=multiplier*work[iRow];
3952#endif
3953#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3954                          double alpha3 = alpha+array[iRow];
3955                          int iSequence = pivotVariable_[iRow];
3956                          double oldLower = lowerCoefficient[iRow];
3957                          double oldUpper = upperCoefficient[iRow];
3958                          if (lower_[iSequence]>-1.0e30) {
3959                            //lowerGap[iRow]=value-lower_[iSequence];
3960                            double alpha2 = alpha3 + lowerChange[iSequence];
3961                            if (alpha2>1.0e-8)  {
3962                              lowerCoefficient[iRow]=alpha2;
3963                              if (!oldLower)
3964                                lowerActive[lowerN++]=iRow;
3965                            } else {
3966                              if (oldLower)
3967                                lowerCoefficient[iRow]=COIN_DBL_MIN;
3968                            }
3969                          } else {
3970                            if (oldLower)
3971                              lowerCoefficient[iRow]=COIN_DBL_MIN;
3972                          }
3973                          if (upper_[iSequence]<1.0e30) {
3974                            //upperGap[iRow]=-(value-upper_[iSequence]);
3975                            double alpha2 = -(alpha3+upperChange[iSequence]);
3976                            if (alpha2>1.0e-8) {
3977                              upperCoefficient[iRow]=alpha2;
3978                              if (!oldUpper)
3979                                upperActive[upperN++]=iRow;
3980                            } else {
3981                              if (oldUpper)
3982                                upperCoefficient[iRow]=COIN_DBL_MIN;
3983                            }
3984                          } else {
3985                            if (oldUpper)
3986                              upperCoefficient[iRow]=COIN_DBL_MIN;
3987                          }
3988#endif
3989                          rowArray_[4]->quickAdd(iRow,alpha);
3990                        }
3991                      }
3992                      pivotValue = array[pivotRow_];
3993                      // we want pivot to be -multiplier
3994                      rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3995#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3996                      assert (lowerN>=0&&lowerN<=numberRows_);
3997                      lowerActive[-1]=lowerN;
3998                      upperActive[-1]=upperN;
3999#endif
4000                    }
4001                    // update primal solution
4002#if CLP_CAN_HAVE_ZERO_OBJ
4003                    if ((specialOptions_&2097152)!=0) 
4004                      theta_=0.0;
4005#endif
4006                    if (theta_ < 0.0) {
4007#ifdef CLP_DEBUG
4008                         if (handler_->logLevel() & 32)
4009                              printf("negative theta %g\n", theta_);
4010#endif
4011                         theta_ = 0.0;
4012                    }
4013                    // do actual flips
4014                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
4015                    //rowArray_[1]->expand();
4016#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4017                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
4018                                                        movement,
4019                                                        objectiveChange);
4020#else
4021                    // do by hand
4022                    {
4023                      double * work = rowArray_[1]->denseVector();
4024                      int number = rowArray_[1]->getNumElements();
4025                      int * which = rowArray_[1]->getIndices();
4026                      int i;
4027                      if (rowArray_[1]->packedMode()) {
4028                        for (i = 0; i < number; i++) {
4029                          int iRow = which[i];
4030                          int iSequence = pivotVariable_[iRow];
4031                          double value = solution_[iSequence];
4032                          double change = movement * work[i];
4033                          value -= change;
4034                          if (lower_[iSequence]>-1.0e30)
4035                            lowerGap[iRow]=value-lower_[iSequence];
4036                          if (upper_[iSequence]<1.0e30)
4037                            upperGap[iRow]=-(value-upper_[iSequence]);
4038                          solution_[iSequence] = value;
4039                          objectiveChange -= change * cost_[iSequence];
4040                          work[i] = 0.0;
4041                        }
4042                      } else {
4043                        for (i = 0; i < number; i++) {
4044                          int iRow = which[i];
4045                          int iSequence = pivotVariable_[iRow];
4046                          double value = solution_[iSequence];
4047                          double change = movement * work[iRow];
4048                          value -= change;
4049                          solution_[iSequence] = value;
4050                          objectiveChange -= change * cost_[iSequence];
4051                          work[iRow] = 0.0;
4052                        }
4053                      }
4054                      rowArray_[1]->setNumElements(0);
4055                    }
4056#endif
4057                    // modify dualout
4058                    dualOut_ /= alpha_;
4059                    dualOut_ *= -directionOut_;
4060                    //setStatus(sequenceIn_,basic);
4061                    dj_[sequenceIn_] = 0.0;
4062                    //double oldValue = valueIn_;
4063                    if (directionIn_ == -1) {
4064                         // as if from upper bound
4065                         valueIn_ = upperIn_ + dualOut_;
4066                    } else {
4067                         // as if from lower bound
4068                         valueIn_ = lowerIn_ + dualOut_;
4069                    }
4070                    objectiveChange = 0.0;
4071#if CLP_CAN_HAVE_ZERO_OBJ
4072                    if ((specialOptions_&2097152)==0) {
4073#endif
4074                      for (int i=0;i<numberTotal;i++)
4075                        objectiveChange += solution_[i]*cost_[i];
4076                      objectiveChange -= objectiveValue_;
4077#if CLP_CAN_HAVE_ZERO_OBJ
4078                    }
4079#endif
4080                    // outgoing
4081                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
4082                    lowerOut_=lower_[sequenceOut_];
4083                    upperOut_=upper_[sequenceOut_];
4084                    // set dj to zero unless values pass
4085                    if (directionOut_ > 0) {
4086                         valueOut_ = lowerOut_;
4087                         dj_[sequenceOut_] = theta_;
4088#if CLP_CAN_HAVE_ZERO_OBJ>1
4089#ifdef COIN_REUSE_RANDOM
4090                         if ((specialOptions_&2097152)!=0) {
4091                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
4092                         }
4093#endif
4094#endif
4095                    } else {
4096                         valueOut_ = upperOut_;
4097                         dj_[sequenceOut_] = -theta_;
4098#if CLP_CAN_HAVE_ZERO_OBJ>1
4099#ifdef COIN_REUSE_RANDOM
4100                         if ((specialOptions_&2097152)!=0) {
4101                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
4102                         }
4103#endif
4104#endif
4105                    }
4106                    solution_[sequenceOut_] = valueOut_;
4107                    int whatNext = housekeeping(objectiveChange);
4108                    reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_);
4109                    assert (backwardBasic[sequenceOut_]==pivotRow_);
4110                    backwardBasic[sequenceOut_]=-1;
4111                    backwardBasic[sequenceIn_]=pivotRow_;
4112#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4113                    double value = solution_[sequenceIn_];
4114                    double alpha = rowArray_[4]->denseVector()[pivotRow_];
4115                    double oldLower = lowerCoefficient[pivotRow_];
4116                    double oldUpper = upperCoefficient[pivotRow_];
4117                    if (lower_[sequenceIn_]>-1.0e30) {
4118                      lowerGap[pivotRow_]=value-lower_[sequenceIn_];
4119                      double alpha2 = alpha + lowerChange[sequenceIn_];
4120                      if (alpha2>1.0e-8)  {
4121                        lowerCoefficient[pivotRow_]=alpha2;
4122                        if (!oldLower) {
4123                          int lowerN=lowerActive[-1];
4124                          assert (lowerN>=0&&lowerN<numberRows_);
4125                          lowerActive[lowerN]=pivotRow_;
4126                          lowerActive[-1]=lowerN+1;
4127                        }
4128                      } else {
4129                        if (oldLower)
4130                          lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4131                      }
4132                    } else {
4133                      if (oldLower)
4134                        lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4135                    }
4136                    if (upper_[sequenceIn_]<1.0e30) {
4137                      upperGap[pivotRow_]=-(value-upper_[sequenceIn_]);
4138                      double alpha2 = -(alpha+upperChange[sequenceIn_]);
4139                      if (alpha2>1.0e-8) {
4140                        upperCoefficient[pivotRow_]=alpha2;
4141                        if (!oldUpper) {
4142                          int upperN=upperActive[-1];
4143                          assert (upperN>=0&&upperN<numberRows_);
4144                          upperActive[upperN]=pivotRow_;
4145                          upperActive[-1]=upperN+1;
4146                        }
4147                      } else {
4148                        if (oldUpper)
4149                          upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4150                      }
4151                    } else {
4152                      if (oldUpper)
4153                        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4154                    }
4155#endif
4156                    {
4157                      char in[200],out[200];
4158                      int iSequence=sequenceIn_;
4159                      if (iSequence<numberColumns_) {
4160                        if (lengthNames_) 
4161                          strcpy(in,columnNames_[iSequence].c_str());
4162                         else 
4163                          sprintf(in,"C%7.7d",iSequence);
4164                      } else {
4165                        iSequence -= numberColumns_;
4166                        if (lengthNames_) 
4167                          strcpy(in,rowNames_[iSequence].c_str());
4168                         else 
4169                          sprintf(in,"R%7.7d",iSequence);
4170                      }
4171                      iSequence=sequenceOut_;
4172                      if (iSequence<numberColumns_) {
4173                        if (lengthNames_) 
4174                          strcpy(out,columnNames_[iSequence].c_str());
4175                         else 
4176                          sprintf(out,"C%7.7d",iSequence);
4177                      } else {
4178                        iSequence -= numberColumns_;
4179                        if (lengthNames_) 
4180                          strcpy(out,rowNames_[iSequence].c_str());
4181                         else 
4182                          sprintf(out,"R%7.7d",iSequence);
4183                      }
4184                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
4185                        << useTheta << objectiveValue() 
4186                        << in << out << CoinMessageEol;
4187                    }
4188                    if (useTheta>lastTheta+1.0e-9) {
4189                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
4190                        << useTheta << objectiveValue() << CoinMessageEol;
4191                      lastTheta = useTheta;
4192                    }
4193                    // and set bounds correctly
4194                    originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
4195                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
4196                    if (whatNext == 1) {
4197                         problemStatus_ = -2; // refactorize
4198                    } else if (whatNext == 2) {
4199                         // maximum iterations or equivalent
4200                         problemStatus_ = 3;
4201                         returnCode = 3;
4202                         break;
4203                    }
4204#ifdef CLP_USER_DRIVEN
4205                    // Check event
4206                    {
4207                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
4208                         if (status >= 0) {
4209                              problemStatus_ = 5;
4210                              secondaryStatus_ = ClpEventHandler::endOfIteration;
4211                              returnCode = 4;
4212                              break;
4213                         }
4214                    }
4215#endif
4216               } else {
4217                    // no incoming column is valid
4218#ifdef CLP_USER_DRIVEN
4219                 rowArray_[0]->clear();
4220                 columnArray_[0]->clear();
4221                 theta_ = useTheta;
4222                 lastTheta = useTheta;
4223                 int action = eventHandler_->event(ClpEventHandler::noTheta);
4224                 if (action>=0) {
4225                   endingTheta=theta_;
4226                   theta_ = 0.0;
4227                   //adjust [4] from handler - but
4228                   //rowArray_[4]->clear(); // temp
4229                   if (action>=0&&action<10)
4230                     problemStatus_=-1; // carry on
4231                   else if (action==15)
4232                     problemStatus_ =5; // say stopped
4233                   returnCode = 1;
4234                   if (action==0||action>=10) 
4235                     break;
4236                   else
4237                     continue;
4238                 } else {
4239                 theta_ = 0.0;
4240                 }
4241#endif
4242                    pivotRow_ = -1;
4243#ifdef CLP_DEBUG
4244                    if (handler_->logLevel() & 32)
4245                         printf("** no column pivot\n");
4246#endif
4247                    if (factorization_->pivots() < 10) { 
4248                         // If we have just factorized and infeasibility reasonable say infeas
4249                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
4250                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
4251                                        || (specialOptions_ & 64) == 0) {
4252                                   // say infeasible
4253                                   problemStatus_ = 1;
4254                                   // unless primal feasible!!!!
4255                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
4256                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
4257                                   if (numberDualInfeasibilities_)
4258                                        problemStatus_ = 10;
4259                                   rowArray_[0]->clear();
4260                                   columnArray_[0]->clear();
4261                              }
4262                         }
4263                         // If special option set - put off as long as possible
4264                         if ((specialOptions_ & 64) == 0) {
4265                              problemStatus_ = -4; //say looks infeasible
4266                         } else {
4267                              // flag
4268                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
4269                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
4270                                        << x << sequenceWithin(sequenceOut_)
4271                                        << CoinMessageEol;
4272                              setFlagged(sequenceOut_);
4273                              if (!factorization_->pivots()) {
4274                                   rowArray_[0]->clear();
4275                                   columnArray_[0]->clear();
4276                                   continue;
4277                              }
4278                         }
4279                    }
4280                    rowArray_[0]->clear();
4281                    columnArray_[0]->clear();
4282                    returnCode = 1;
4283                    break;
4284               }
4285          } else {
4286               // no pivot row
4287#ifdef CLP_USER_DRIVEN
4288            {
4289              double saveTheta=theta_;
4290              theta_ = endingTheta;
4291              int status=eventHandler_->event(ClpEventHandler::theta);
4292              if (status>=0&&status<10) {
4293                endingTheta=theta_;
4294                theta_=saveTheta;
4295                continue;
4296              } else {
4297                theta_=saveTheta;
4298              }
4299            }
4300#endif
4301#ifdef CLP_DEBUG
4302               if (handler_->logLevel() & 32)
4303                    printf("** no row pivot\n");
4304#endif
4305               int numberPivots = factorization_->pivots();
4306               bool specialCase;
4307               int useNumberFake;
4308               returnCode = 0;
4309               if (numberPivots < 20 &&
4310                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
4311                         && dualBound_ > 1.0e8) {
4312                    specialCase = true;
4313                    // as dual bound high - should be okay
4314                    useNumberFake = 0;
4315               } else {
4316                    specialCase = false;
4317                    useNumberFake = numberFake_;
4318               }
4319               if (!numberPivots || specialCase) {
4320                    // may have crept through - so may be optimal
4321                    // check any flagged variables
4322                    int iRow;
4323                    for (iRow = 0; iRow < numberRows_; iRow++) {
4324                         int iPivot = pivotVariable_[iRow];
4325                         if (flagged(iPivot))
4326                              break;
4327                    }
4328                    if (iRow < numberRows_ && numberPivots) {
4329                         // try factorization
4330                         returnCode = -2;
4331                    }
4332
4333                    if (useNumberFake || numberDualInfeasibilities_) {
4334                         // may be dual infeasible
4335                         problemStatus_ = -5;
4336                    } else {
4337                         if (iRow < numberRows_) {
4338                              problemStatus_ = -5;
4339                         } else {
4340                              if (numberPivots) {
4341                                   // objective may be wrong
4342                                   objectiveValue_ = innerProduct(cost_,
4343                                                                  numberColumns_ + numberRows_,
4344                                                                  solution_);
4345                                   objectiveValue_ += objective_->nonlinearOffset();
4346                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
4347                                   if ((specialOptions_ & 16384) == 0) {
4348                                        // and dual_ may be wrong (i.e. for fixed or basic)
4349                                        CoinIndexedVector * arrayVector = rowArray_[1];
4350                                        arrayVector->clear();
4351                                        int iRow;
4352                                        double * array = arrayVector->denseVector();
4353                                        /* Use dual_ instead of array
4354                                           Even though dual_ is only numberRows_ long this is
4355                                           okay as gets permuted to longer rowArray_[2]
4356                                        */
4357                                        arrayVector->setDenseVector(dual_);
4358                                        int * index = arrayVector->getIndices();
4359                                        int number = 0;
4360                                        for (iRow = 0; iRow < numberRows_; iRow++) {
4361                                             int iPivot = pivotVariable_[iRow];
4362                                             double value = cost_[iPivot];
4363                                             dual_[iRow] = value;
4364                                             if (value) {
4365                                                  index[number++] = iRow;
4366                                             }
4367                                        }
4368                                        arrayVector->setNumElements(number);
4369                                        // Extended duals before "updateTranspose"
4370                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
4371                                        // Btran basic costs
4372                                        rowArray_[2]->clear();
4373                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
4374                                        // and return vector
4375                                        arrayVector->setDenseVector(array);
4376                                   }
4377                              }
4378                              problemStatus_ = 0;
4379                              sumPrimalInfeasibilities_ = 0.0;
4380                              if ((specialOptions_&(1024 + 16384)) != 0) {
4381                                   CoinIndexedVector * arrayVector = rowArray_[1];
4382                                   arrayVector->clear();
4383                                   double * rhs = arrayVector->denseVector();
4384                                   times(1.0, solution_, rhs);
4385                                   bool bad2 = false;
4386                                   int i;
4387                                   for ( i = 0; i < numberRows_; i++) {
4388                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
4389                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
4390                                             bad2 = true;
4391                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
4392                                        }
4393                                        rhs[i] = 0.0;
4394                                   }
4395                                   for ( i = 0; i < numberColumns_; i++) {
4396                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
4397                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4398                                             bad2 = true;
4399                                        }
4400                                   }
4401                                   if (bad2) {
4402                                        problemStatus_ = -3;
4403                                        returnCode = -2;
4404                                        // Force to re-factorize early next time
4405                                        int numberPivots = factorization_->pivots();
4406                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4407                                   }
4408                              }
4409                         }
4410                    }
4411               } else {
4412                    problemStatus_ = -3;
4413                    returnCode = -2;
4414                    // Force to re-factorize early next time
4415                    int numberPivots = factorization_->pivots();
4416                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4417               }
4418               break;
4419          }
4420     }
4421     startingTheta = lastTheta+theta_;
4422     return returnCode;
4423}
4424#if 0
4425static int zzzzzz=0;
4426int zzzzzzOther=0;
4427#endif
4428// Finds best possible pivot
4429double 
4430ClpSimplexOther::bestPivot(bool justColumns)
4431{
4432  // Get good size for pivot
4433  // Allow first few iterations to take tiny
4434  double acceptablePivot = 1.0e-9;
4435  if (numberIterations_ > 100)
4436    acceptablePivot = 1.0e-8;
4437  if (factorization_->pivots() > 10 ||
4438      (factorization_->pivots() && sumDualInfeasibilities_))
4439    acceptablePivot = 1.0e-5; // if we have iterated be more strict
4440  else if (factorization_->pivots() > 5)
4441    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4442  else if (factorization_->pivots())
4443    acceptablePivot = 1.0e-8; // relax
4444  double bestPossiblePivot = 1.0;
4445  // get sign for finding row of tableau
4446  // normal iteration
4447  // create as packed
4448  double direction = directionOut_;
4449#ifndef COIN_FAC_NEW
4450  rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4451#else
4452  rowArray_[0]->createOneUnpackedElement(pivotRow_, direction);
4453#endif
4454  factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4455  // put row of tableau in rowArray[0] and columnArray[0]
4456  matrix_->transposeTimes(this, -1.0,
4457                          rowArray_[0], rowArray_[3], columnArray_[0]);
4458  sequenceIn_=-1;
4459  if (justColumns)
4460    rowArray_[0]->clear();
4461  // do ratio test for normal iteration
4462  bestPossiblePivot = 
4463    reinterpret_cast<ClpSimplexDual *> 
4464    ( this)->dualColumn(rowArray_[0],
4465                        columnArray_[0], columnArray_[1],
4466                        rowArray_[3], acceptablePivot, NULL);
4467  return bestPossiblePivot;
4468}
4469// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4470int
4471ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData,
4472                           const double * /*changeObjective*/)
4473{
4474  const double * lowerChange = paramData.lowerChange;
4475  const double * upperChange = paramData.upperChange;
4476  const int * lowerList = paramData.lowerList;
4477  const int * upperList = paramData.upperList;
4478  int iSequence;
4479  bool toLower = false;
4480  //assert (type==1);
4481  // may need to decide based on model?
4482  bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4483  double * array = rowArray_[4]->denseVector();
4484  //rowArray_[4]->checkClean();
4485  const int * row = matrix_->getIndices();
4486  const int * columnLength = matrix_->getVectorLengths();
4487  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4488  const double * elementByColumn = matrix_->getElements();
4489#if 0
4490  double tempArray[5000];
4491  bool checkIt=false;
4492  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
4493    memcpy(tempArray,array,numberRows_*sizeof(double));
4494    checkIt=true;
4495    needFullUpdate=true;
4496  }
4497#endif
4498#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4499  double * lowerGap = paramData.lowerGap;
4500  double * upperGap = paramData.upperGap;
4501  double * lowerCoefficient = paramData.lowerCoefficient;
4502  double * upperCoefficient = paramData.upperCoefficient;
4503  int * lowerActive=paramData.lowerActive;
4504  int * upperActive=paramData.upperActive;
4505#endif
4506  if (!factorization_->pivots()||needFullUpdate) {
4507    //zzzzzz=0;
4508    rowArray_[4]->clear();
4509    // get change
4510    if (!rowScale_) {
4511      int n;
4512      n=lowerList[-2];
4513      int i;
4514      for (i=0;i<n;i++) {
4515        int iSequence = lowerList[i];
4516        assert (iSequence<numberColumns_);
4517        if (getColumnStatus(iSequence)==atLowerBound) {
4518          double value=lowerChange[iSequence];
4519          for (CoinBigIndex j = columnStart[iSequence];
4520               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4521            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4522          }
4523        }
4524      }
4525      n=lowerList[-1];
4526      const double * change = lowerChange+numberColumns_;
4527      for (;i<n;i++) {
4528        int iSequence = lowerList[i]-numberColumns_;
4529        assert (iSequence>=0);
4530        if (getRowStatus(iSequence)==atLowerBound) {
4531          double value=change[iSequence];
4532          rowArray_[4]->quickAdd(iSequence, -value);
4533        }
4534      }
4535      n=upperList[-2];
4536      for (i=0;i<n;i++) {
4537        int iSequence = upperList[i];
4538        assert (iSequence<numberColumns_);
4539        if (getColumnStatus(iSequence)==atUpperBound) {
4540          double value=upperChange[iSequence];
4541          for (CoinBigIndex j = columnStart[iSequence];
4542               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4543            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4544          }
4545        }
4546      }
4547      n=upperList[-1];
4548      change = upperChange+numberColumns_;
4549      for (;i<n;i++) {
4550        int iSequence = upperList[i]-numberColumns_;
4551        assert (iSequence>=0);
4552        if (getRowStatus(iSequence)==atUpperBound) {
4553          double value=change[iSequence];
4554          rowArray_[4]->quickAdd(iSequence, -value);
4555        }
4556      }
4557    } else {
4558      int n;
4559      n=lowerList[-2];
4560      int i;
4561      for (i=0;i<n;i++) {
4562        int iSequence = lowerList[i];
4563        assert (iSequence<numberColumns_);
4564        if (getColumnStatus(iSequence)==atLowerBound) {
4565          double value=lowerChange[iSequence];
4566          // apply scaling
4567          double scale = columnScale_[iSequence];
4568          for (CoinBigIndex j = columnStart[iSequence];
4569               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4570            int iRow = row[j];
4571            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4572          }
4573        }
4574      }
4575      n=lowerList[-1];
4576      const double * change = lowerChange+numberColumns_;
4577      for (;i<n;i++) {
4578        int iSequence = lowerList[i]-numberColumns_;
4579        assert (iSequence>=0);
4580        if (getRowStatus(iSequence)==atLowerBound) {
4581          double value=change[iSequence];
4582          rowArray_[4]->quickAdd(iSequence, -value);
4583        }
4584      }
4585      n=upperList[-2];
4586      for (i=0;i<n;i++) {
4587        int iSequence = upperList[i];
4588        assert (iSequence<numberColumns_);
4589        if (getColumnStatus(iSequence)==atUpperBound) {
4590          double value=upperChange[iSequence];
4591          // apply scaling
4592          double scale = columnScale_[iSequence];
4593          for (CoinBigIndex j = columnStart[iSequence];
4594               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4595            int iRow = row[j];
4596            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4597          }
4598        }
4599      }
4600      n=upperList[-1];
4601      change = upperChange+numberColumns_;
4602      for (;i<n;i++) {
4603        int iSequence = upperList[i]-numberColumns_;
4604        assert (iSequence>=0);
4605        if (getRowStatus(iSequence)==atUpperBound) {
4606          double value=change[iSequence];
4607          rowArray_[4]->quickAdd(iSequence, -value);
4608        }
4609      }
4610    }
4611    // ftran it
4612    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4613#if 0
4614    if (checkIt) {
4615      for (int i=0;i<numberRows_;i++) {
4616        assert (fabs(tempArray[i]-array[i])<1.0e-8);
4617      }
4618    }
4619#endif
4620#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4621    /* later for sparse - keep like CoinIndexedvector
4622       and just redo here */
4623    int lowerN=0;
4624    int upperN=0;
4625    memset(lowerCoefficient,0,numberRows_*sizeof(double));
4626    memset(upperCoefficient,0,numberRows_*sizeof(double));
4627    for (int iRow=0;iRow<numberRows_;iRow++) {
4628      iSequence = pivotVariable_[iRow];
4629      double currentSolution = solution_[iSequence];
4630      double alpha = array[iRow];
4631      double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4632      double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4633      if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4634        double currentLower = lower_[iSequence];
4635        ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4636        double gap=currentSolution-currentLower;
4637        lowerGap[iRow]=gap;
4638        lowerCoefficient[iRow]=thetaCoefficientLower;
4639        lowerActive[lowerN++]=iRow;
4640        //} else {
4641        //lowerCoefficient[iRow]=0.0;
4642      }
4643      if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4644        double currentUpper = upper_[iSequence];
4645        ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4646        double gap2=-(currentSolution-currentUpper); //positive
4647        upperGap[iRow]=gap2;
4648        upperCoefficient[iRow]=-thetaCoefficientUpper;
4649        upperActive[upperN++]=iRow;
4650      }
4651    }
4652    assert (lowerN>=0&&lowerN<=numberRows_);
4653    lowerActive[-1]=lowerN;
4654    upperActive[-1]=upperN;
4655#endif
4656  } else if (sequenceIn_>=0) {
4657    //assert (sequenceIn_>=0);
4658    assert (sequenceOut_>=0);
4659    assert (sequenceIn_!=sequenceOut_);
4660    double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4661    int needed=0;
4662    assert (!rowArray_[5]->getNumElements());
4663    if (change) {
4664      if (sequenceIn_<numberColumns_) {
4665        if (!rowScale_) {
4666          for (CoinBigIndex i = columnStart[sequenceIn_];
4667               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4668            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4669          }
4670        } else {
4671          // apply scaling
4672          double scale = columnScale_[sequenceIn_];
4673          for (CoinBigIndex i = columnStart[sequenceIn_];
4674               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4675            int iRow = row[i];
4676            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4677          }
4678        }
4679      } else {
4680        rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4681      }
4682      needed++;
4683    }
4684    if (getStatus(sequenceOut_)==atLowerBound)
4685      change=lowerChange[sequenceOut_];
4686    else
4687      change=upperChange[sequenceOut_];
4688    if (change) {
4689      if (sequenceOut_<numberColumns_) {
4690        if (!rowScale_) {
4691          for (CoinBigIndex i = columnStart[sequenceOut_];
4692               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4693            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4694          }
4695        } else {
4696          // apply scaling
4697          double scale = columnScale_[sequenceOut_];
4698          for (CoinBigIndex i = columnStart[sequenceOut_];
4699               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4700            int iRow = row[i];
4701            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4702          }
4703        }
4704      } else {
4705        rowArray_[5]->quickAdd(sequenceOut_-numberColumns_,-change);
4706      }
4707      needed++;
4708    }
4709    //printf("seqin %d seqout %d needed %d\n",
4710    //     sequenceIn_,sequenceOut_,needed);
4711    if (needed) {
4712      // ftran it
4713      factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4714      // add
4715      double * array5 = rowArray_[5]->denseVector();
4716      int * index5 = rowArray_[5]->getIndices();
4717      int number5 = rowArray_[5]->getNumElements();
4718#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4719      int lowerN=lowerActive[-1];
4720      int upperN=upperActive[-1];
4721      int nIn4=rowArray_[4]->getNumElements();
4722      int * index4 = rowArray_[4]->getIndices();
4723#endif
4724      for (int i = 0; i < number5; i++) {
4725        int iPivot = index5[i];
4726#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4727        rowArray_[4]->quickAdd(iPivot,array5[iPivot]);
4728#else
4729        /* later for sparse - modify here */
4730        int iSequence = pivotVariable_[iPivot];
4731        double currentSolution = solution_[iSequence];
4732        double currentAlpha = array[iPivot];
4733        double alpha5 = array5[iPivot];
4734        double alpha = currentAlpha+alpha5;
4735        if (currentAlpha) {
4736          if (alpha) {
4737            array[iPivot] = alpha;
4738          } else {
4739            array[iPivot] = COIN_DBL_MIN;
4740          }
4741        } else {
4742          index4[nIn4++] = iPivot;
4743          array[iPivot] = alpha;
4744        }
4745        double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4746        double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4747        double oldLower = lowerCoefficient[iPivot];
4748        double oldUpper = upperCoefficient[iPivot];
4749        if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4750          double currentLower = lower_[iSequence];
4751          ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4752          double gap=currentSolution-currentLower;
4753          lowerGap[iPivot]=gap;
4754          lowerCoefficient[iPivot]=thetaCoefficientLower;
4755          if (!oldLower)
4756            lowerActive[lowerN++]=iPivot;
4757        } else {
4758          if (oldLower)
4759            lowerCoefficient[iPivot]=COIN_DBL_MIN;
4760        }
4761        if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4762          double currentUpper = upper_[iSequence];
4763          ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4764          double gap2=-(currentSolution-currentUpper); //positive
4765          upperGap[iPivot]=gap2;
4766          upperCoefficient[iPivot]=-thetaCoefficientUpper;
4767          if (!oldUpper)
4768            upperActive[upperN++]=iPivot;
4769        } else {
4770          if (oldUpper)
4771            upperCoefficient[iPivot]=COIN_DBL_MIN;
4772        }
4773#endif
4774        array5[iPivot]=0.0;
4775      }
4776      rowArray_[5]->setNumElements(0);
4777#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4778      rowArray_[4]->setNumElements(nIn4);
4779      assert (lowerN>=0&&lowerN<=numberRows_);
4780      lowerActive[-1]=lowerN;
4781      upperActive[-1]=upperN;
4782#endif
4783    }
4784  }
4785  const int * index = rowArray_[4]->getIndices();
4786  int number = rowArray_[4]->getNumElements();
4787#define TESTXX 0
4788#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4789  int * markDone = reinterpret_cast<int *>(paramData.markDone);
4790  int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT;
4791  memset(markDone,0,nToZero*sizeof(int));
4792  const int * backwardBasic = paramData.backwardBasic;
4793#endif
4794  // first ones with alpha
4795  double theta1=maxTheta;
4796  int pivotRow1=-1;
4797#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4798  int pivotRow2=-1;
4799  double theta2=maxTheta;
4800#endif
4801#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4802  for (int i=0;i<number;i++) {
4803    int iPivot=index[i];
4804    iSequence = pivotVariable_[iPivot];
4805    //assert(!markDone[iSequence]);
4806    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4807    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4808    markDone[word] |= ( 1 << bit );
4809    // solution value will be sol - theta*alpha
4810    // bounds will be bounds + change *theta
4811    double currentSolution = solution_[iSequence];
4812    double alpha = array[iPivot];
4813    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4814    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4815    if (thetaCoefficientLower > 1.0e-8) {
4816      double currentLower = lower_[iSequence];
4817      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4818      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4819      double gap=currentSolution-currentLower;
4820      if (thetaCoefficientLower*theta1>gap) {
4821        theta1 = gap/thetaCoefficientLower;
4822        //toLower=true;
4823        pivotRow1=iPivot;
4824      }
4825    }
4826    if (thetaCoefficientUpper < -1.0e-8) {
4827      double currentUpper = upper_[iSequence];
4828      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4829      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4830      double gap2=currentSolution-currentUpper; //negative
4831      if (thetaCoefficientUpper*theta2<gap2) {
4832        theta2 = gap2/thetaCoefficientUpper;
4833        //toLower=false;
4834        pivotRow2=iPivot;
4835      }
4836    }
4837  }
4838  // now others
4839  int nLook=lowerList[-1];
4840  for (int i=0;i<nLook;i++) {
4841    int iSequence = lowerList[i];
4842    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4843    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4844    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4845      double currentSolution = solution_[iSequence];
4846      double currentLower = lower_[iSequence];
4847      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4848      double thetaCoefficient = lowerChange[iSequence];
4849      if (thetaCoefficient > 0.0) {
4850        double gap=currentSolution-currentLower;
4851        if (thetaCoefficient*theta1>gap) {
4852          theta1 = gap/thetaCoefficient;
4853          //toLower=true;
4854          pivotRow1 = backwardBasic[iSequence];
4855        }
4856      }
4857    }
4858  }
4859  nLook=upperList[-1];
4860  for (int i=0;i<nLook;i++) {
4861    int iSequence = upperList[i];
4862    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4863    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4864    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4865      double currentSolution = solution_[iSequence];
4866      double currentUpper = upper_[iSequence];
4867      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4868      double thetaCoefficient = upperChange[iSequence];
4869      if (thetaCoefficient < 0) {
4870        double gap=currentSolution-currentUpper; //negative
4871        if (thetaCoefficient*theta2<gap) {
4872          theta2 = gap/thetaCoefficient;
4873          //toLower=false;
4874          pivotRow2 = backwardBasic[iSequence];
4875        }
4876      }
4877    }
4878  }
4879  if (theta2<theta1) {
4880    theta_=theta2;
4881    toLower=false;
4882    pivotRow_=pivotRow2;
4883  } else {
4884    theta_=theta1;
4885    toLower=true;
4886    pivotRow_=pivotRow1;
4887  }
4888#if 0 //TESTXX
4889#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4890  {
4891    double * checkArray = new double[numberRows_];
4892    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4893    int lowerN=lowerActive[-1];
4894    for (int i=0;i<lowerN;i++) {
4895      int iRow=lowerActive[i];
4896      int iSequence = pivotVariable_[iRow];
4897      double alpha = array[iRow];
4898      double thetaCoefficient = lowerChange[iSequence] + alpha;
4899      if (thetaCoefficient > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4900        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4901        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4902          abort();
4903        }
4904      } else {
4905        assert (fabs(checkArray[iRow])<1.0e-12);
4906        if (fabs(checkArray[iRow])>1.0e-12) {
4907          abort();
4908        }
4909      }
4910      checkArray[iRow]=0.0;
4911    }
4912    for (int i=0;i<numberRows_;i++) {
4913      assert (!checkArray[i]);
4914      if (checkArray[i])
4915        abort();
4916    }
4917    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
4918    int upperN=upperActive[-1];
4919    for (int i=0;i<upperN;i++) {
4920      int iRow=upperActive[i];
4921      int iSequence = pivotVariable_[iRow];
4922      double alpha = array[iRow];
4923      double thetaCoefficient = -(upperChange[iSequence] + alpha);
4924      if (thetaCoefficient > 1.0e-8&&upper_[iSequence]<1.0e30) {
4925        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4926        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4927          abort();
4928        }
4929      } else {
4930        assert (fabs(checkArray[iRow])<1.0e-12);
4931        if (fabs(checkArray[iRow])>1.0e-12) {
4932          abort();
4933        }
4934      }
4935      checkArray[iRow]=0.0;
4936    }
4937    for (int i=0;i<numberRows_;i++) {
4938      assert (!checkArray[i]);
4939      if (checkArray[i])
4940        abort();
4941    }
4942    delete [] checkArray;
4943  }
4944  double theta3=maxTheta;
4945  int pivotRow3=-1;
4946  int lowerN=lowerActive[-1];
4947  for (int i=0;i<lowerN;i++) {
4948    int iRow=lowerActive[i];
4949    double lowerC = lowerCoefficient[iRow];
4950    double gap=lowerGap[iRow];
4951    if (toLower&&iRow==pivotRow_) {
4952      assert (lowerC*theta3>gap-1.0e-8);
4953      if (lowerC*theta3<gap-1.0e-8)
4954        abort();
4955    }
4956    if (lowerC*theta3>gap&&lowerC!=COIN_DBL_MIN) {
4957      theta3 = gap/lowerC;
4958      pivotRow3=iRow;
4959    }
4960  }
4961  int pivotRow4=pivotRow3;
4962  double theta4=theta3;
4963  int upperN=upperActive[-1];
4964  for (int i=0;i<upperN;i++) {
4965    int iRow=upperActive[i];
4966    double upperC = upperCoefficient[iRow];
4967    double gap=upperGap[iRow];
4968    if (!toLower&&iRow==pivotRow_) {
4969      assert (upperC*theta3>gap-1.0e-8);
4970      if (upperC*theta3<gap-1.0e-8)
4971        abort();
4972    }
4973    if (upperC*theta4>gap&&upperC!=COIN_DBL_MIN) {
4974      theta4 = gap/upperC;
4975      pivotRow4=iRow;
4976    }
4977  }
4978  bool toLower3;
4979  if (theta4<theta3) {
4980    theta3=theta4;
4981    toLower3=false;
4982    pivotRow3=pivotRow4;
4983  } else {
4984    toLower3=true;
4985  }
4986  if (fabs(theta3-theta_)>1.0e-8)
4987    abort();
4988  if (toLower!=toLower3||pivotRow_!=pivotRow3) {
4989    printf("bad piv - good %d %g %s, bad %d %g %s\n",pivotRow_,theta_,toLower ? "toLower" : "toUpper",
4990           pivotRow3,theta3,toLower3 ? "toLower" : "toUpper");
4991    //zzzzzz++;
4992    if (true/*zzzzzz>zzzzzzOther*/) {
4993      printf("Swapping\n");
4994      pivotRow_=pivotRow3;
4995      theta_=theta3;
4996      toLower=toLower3;
4997    }
4998  }
4999#endif
5000#endif
5001#else
5002#if 0 //CLP_PARAMETRIC_DENSE_ARRAYS==2
5003  {
5004    double * checkArray = new double[numberRows_];
5005    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
5006    int lowerN=lowerActive[-1];
5007    for (int i=0;i<lowerN;i++) {
5008      int iRow=lowerActive[i];
5009      checkArray[iRow]=0.0;
5010    }
5011    for (int i=0;i<numberRows_;i++) {
5012      assert (!checkArray[i]);
5013      if (checkArray[i])
5014        abort();
5015    }
5016    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
5017    int upperN=upperActive[-1];
5018    for (int i=0;i<upperN;i++) {
5019      int iRow=upperActive[i];
5020      checkArray[iRow]=0.0;
5021    }
5022    for (int i=0;i<numberRows_;i++) {
5023      assert (!checkArray[i]);
5024      if (checkArray[i])
5025        abort();
5026    }
5027    delete [] checkArray;
5028  }
5029#endif
5030  int lowerN=lowerActive[-1];
5031  for (int i=0;i<lowerN;i++) {
5032    int iRow=lowerActive[i];
5033    double lowerC = lowerCoefficient[iRow];
5034    double gap=lowerGap[iRow];
5035    if (lowerC*theta1>gap&&lowerC!=COIN_DBL_MIN) {
5036      theta1 = gap/lowerC;
5037      pivotRow1=iRow;
5038    }
5039  }
5040  pivotRow_=pivotRow1;
5041  theta_=theta1;
5042  int upperN=upperActive[-1];
5043  for (int i=0;i<upperN;i++) {
5044    int iRow=upperActive[i];
5045    double upperC = upperCoefficient[iRow];
5046    double gap=upperGap[iRow];
5047    if (upperC*theta1>gap&&upperC!=COIN_DBL_MIN) {
5048      theta1 = gap/upperC;
5049      pivotRow1=iRow;
5050    }
5051  }
5052  if (theta1<theta_) {
5053    theta_=theta1;
5054    toLower=false;
5055    pivotRow_=pivotRow1;
5056  } else {
5057    toLower=true;
5058  }
5059#endif
5060  theta_ = CoinMax(theta_,0.0);
5061  if (theta_>1.0e-15) {
5062    // update solution
5063    for (int iRow = 0; iRow < number; iRow++) {
5064      int iPivot = index[iRow];
5065      iSequence = pivotVariable_[iPivot];
5066      // solution value will be sol - theta*alpha
5067      double alpha = array[iPivot];
5068      double currentSolution = solution_[iSequence] - theta_ * alpha;
5069      solution_[iSequence] =currentSolution;
5070#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5071      if (lower_[iSequence]>-1.0e30)
5072        lowerGap[iPivot]=currentSolution-lower_[iSequence];
5073      if (upper_[iSequence]<1.0e30)
5074        upperGap[iPivot]=-(currentSolution-upper_[iSequence]);
5075#endif
5076    }
5077  }
5078#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5079  if (pivotRow_>=0&&false) {
5080    double oldValue = upperCoefficient[pivotRow_];
5081    double value = array[pivotRow_];
5082    if (value) {
5083      if (!oldValue) {
5084        int upperN=upperActive[-1];
5085        assert (upperN>=0&&upperN<numberRows_);
5086        upperActive[upperN]=pivotRow_;
5087        upperActive[-1]=upperN+1;
5088      }
5089    } else {
5090      if (oldValue)
5091        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
5092    }
5093  }
5094#endif
5095#if 0
5096  for (int i=0;i<numberTotal;i++)
5097    assert(!markDone[i]);
5098#endif
5099  if (pivotRow_ >= 0) {
5100    sequenceOut_ = pivotVariable_[pivotRow_];
5101    valueOut_ = solution_[sequenceOut_];
5102    lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
5103    upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
5104    if (!toLower) {
5105      directionOut_ = -1;
5106      dualOut_ = valueOut_ - upperOut_;
5107    } else {
5108      directionOut_ = 1;
5109      dualOut_ = lowerOut_ - valueOut_;
5110    }
5111    return 0;
5112  } else {
5113    //theta_=0.0;
5114    return -1;
5115  }
5116}
5117// Restores bound to original bound
5118void 
5119ClpSimplexOther::originalBound(int iSequence, double theta, 
5120                               const double * lowerChange,
5121                               const double * upperChange)
5122{
5123     if (getFakeBound(iSequence) != noFake) {
5124          numberFake_--;
5125          setFakeBound(iSequence, noFake);
5126          if (iSequence >= numberColumns_) {
5127               // rows
5128               int iRow = iSequence - numberColumns_;
5129               rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
5130               rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
5131               if (rowScale_) {
5132                    if (rowLowerWork_[iRow] > -1.0e50)
5133                         rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5134                    if (rowUpperWork_[iRow] < 1.0e50)
5135                         rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5136               } else if (rhsScale_ != 1.0) {
5137                    if (rowLowerWork_[iRow] > -1.0e50)
5138                         rowLowerWork_[iRow] *= rhsScale_;
5139                    if (rowUpperWork_[iRow] < 1.0e50)
5140                         rowUpperWork_[iRow] *= rhsScale_;
5141               }
5142          } else {
5143               // columns
5144               columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
5145               columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
5146               if (rowScale_) {
5147                    double multiplier = 1.0 * inverseColumnScale_[iSequence];
5148                    if (columnLowerWork_[iSequence] > -1.0e50)
5149                         columnLowerWork_[iSequence] *= multiplier * rhsScale_;
5150                    if (columnUpperWork_[iSequence] < 1.0e50)
5151                         columnUpperWork_[iSequence] *= multiplier * rhsScale_;
5152               } else if (rhsScale_ != 1.0) {
5153                    if (columnLowerWork_[iSequence] > -1.0e50)
5154                         columnLowerWork_[iSequence] *= rhsScale_;
5155                    if (columnUpperWork_[iSequence] < 1.0e50)
5156                         columnUpperWork_[iSequence] *= rhsScale_;
5157               }
5158          }
5159     }
5160}
5161/* Expands out all possible combinations for a knapsack
5162   If buildObj NULL then just computes space needed - returns number elements
5163   On entry numberOutput is maximum allowed, on exit it is number needed or
5164   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
5165   least space to return values which reconstruct input.
5166   Rows returned will be original rows but no entries will be returned for
5167   any rows all of whose entries are in knapsack.  So up to user to allow for this.
5168   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
5169   in expanded knapsack.  Values in buildRow and buildElement;
5170*/
5171int
5172ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
5173                                double * buildObj, CoinBigIndex * buildStart,
5174                                int * buildRow, double * buildElement, int reConstruct) const
5175{
5176     int iRow;
5177     int iColumn;
5178     // Get column copy
5179     CoinPackedMatrix * columnCopy = matrix();
5180     // Get a row copy in standard format
5181     CoinPackedMatrix matrixByRow;
5182     matrixByRow.reverseOrderedCopyOf(*columnCopy);
5183     const double * elementByRow = matrixByRow.getElements();
5184     const int * column = matrixByRow.getIndices();
5185     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5186     const int * rowLength = matrixByRow.getVectorLengths();
5187     CoinBigIndex j;
5188     int * whichColumn = new int [numberColumns_];
5189     int * whichRow = new int [numberRows_];
5190     int numJ = 0;
5191     // Get what other columns can compensate for
5192     double * lo = new double [numberRows_];
5193     double * high = new double [numberRows_];
5194     {
5195          // Use to get tight column bounds
5196          ClpSimplex tempModel(*this);
5197          tempModel.tightenPrimalBounds(0.0, 0, true);
5198          // Now another model without knapsacks
5199          int nCol = 0;
5200          for (iRow = 0; iRow < numberRows_; iRow++) {
5201               whichRow[iRow] = iRow;
5202          }
5203          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
5204               whichColumn[iColumn] = -1;
5205          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
5206               int iColumn = column[j];
5207               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
5208                    whichColumn[iColumn] = 0;
5209               } else {
5210                    assert (!columnLower_[iColumn]); // fix later
5211               }
5212          }
5213          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5214               if (whichColumn[iColumn] < 0)
5215                    whichColumn[nCol++] = iColumn;
5216          }
5217          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
5218          // Row copy
5219          CoinPackedMatrix matrixByRow;
5220          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
5221          const double * elementByRow = matrixByRow.getElements();
5222          const int * column = matrixByRow.getIndices();
5223          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5224          const int * rowLength = matrixByRow.getVectorLengths();
5225          const double * columnLower = tempModel2.getColLower();
5226          const double * columnUpper = tempModel2.getColUpper();
5227          for (iRow = 0; iRow < numberRows_; iRow++) {
5228               lo[iRow] = -COIN_DBL_MAX;
5229               high[iRow] = COIN_DBL_MAX;
5230               if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
5231
5232                    // possible row
5233                    int infiniteUpper = 0;
5234