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

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

clean up code a bit

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 352.9 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1833 2011-12-06 14:23:30Z 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.lowerChange = chgLower;
2052                 paramData.upperChange = chgUpper;
2053                    returnCode = parametricsLoop(paramData, reportIncrement,
2054                                                 chgLower, chgUpper, chgObjective, data,
2055                                                 canTryQuick);
2056                 startingTheta=paramData.startingTheta;
2057                 endingTheta=paramData.endingTheta;
2058                    if (!returnCode) {
2059                         //double change = endingTheta-startingTheta;
2060                         startingTheta = endingTheta;
2061                         endingTheta = saveEndingTheta;
2062                         //for (int i=0;i<numberTotal;i++) {
2063                         //lower_[i] += change*chgLower[i];
2064                         //upper_[i] += change*chgUpper[i];
2065                         //cost_[i] += change*chgObjective[i];
2066                         //}
2067                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
2068                           << startingTheta << objectiveValue() << CoinMessageEol;
2069                         if (startingTheta >= endingTheta)
2070                              break;
2071                    } else if (returnCode == -1) {
2072                         // trouble - do external solve
2073                         needToDoSomething = true;
2074                    } else if (problemStatus_==1) {
2075                      // can't move any further
2076                      if (!canTryQuick) {
2077                        handler_->message(CLP_PARAMETRICS_STATS, messages_)
2078                          << endingTheta << objectiveValue() << CoinMessageEol;
2079                        problemStatus_=0;
2080                      }
2081                    } else {
2082                         abort();
2083                    }
2084               }
2085          }
2086          reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2087
2088          delete dualRowPivot_;
2089          dualRowPivot_ = savePivot;
2090          // Restore any saved stuff
2091          restoreData(data);
2092          if (needToDoSomething) {
2093               double saveStartingTheta = startingTheta; // known to be feasible
2094               int cleanedUp = 1;
2095               while (cleanedUp) {
2096                    // tweak
2097                    if (cleanedUp == 1) {
2098                         if (!reportIncrement)
2099                              startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2100                         else
2101                              startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2102                    } else {
2103                         // restoring to go slowly
2104                         startingTheta = saveStartingTheta;
2105                    }
2106                    // only works if not scaled
2107                    int i;
2108                    const double * obj1 = objective();
2109                    double * obj2 = copyModel.objective();
2110                    const double * lower1 = columnLower_;
2111                    double * lower2 = copyModel.columnLower();
2112                    const double * upper1 = columnUpper_;
2113                    double * upper2 = copyModel.columnUpper();
2114                    for (i = 0; i < numberColumns_; i++) {
2115                         obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2116                         lower2[i] = lower1[i] + startingTheta * chgLower[i];
2117                         upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2118                    }
2119                    lower1 = rowLower_;
2120                    lower2 = copyModel.rowLower();
2121                    upper1 = rowUpper_;
2122                    upper2 = copyModel.rowUpper();
2123                    for (i = 0; i < numberRows_; i++) {
2124                         lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2125                         upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2126                    }
2127                    copyModel.dual();
2128                    if (copyModel.problemStatus()) {
2129                      char line[100];
2130                      sprintf(line,"Can not get to theta of %g\n", startingTheta);
2131                      handler_->message(CLP_GENERAL,messages_)
2132                        << line << CoinMessageEol;
2133                         canTryQuick = false; // do slowly to get exact amount
2134                         // back to last known good
2135                         if (cleanedUp == 1)
2136                              cleanedUp = 2;
2137                         else
2138                              abort();
2139                    } else {
2140                         // and move stuff back
2141                         int numberTotal = numberRows_ + numberColumns_;
2142                         CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2143                         CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2144                         CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2145                         cleanedUp = 0;
2146                    }
2147               }
2148          }
2149          delete [] chgLower;
2150          delete [] chgUpper;
2151          delete [] chgObjective;
2152     }
2153     perturbation_ = savePerturbation;
2154     char line[100];
2155     sprintf(line,"Ending theta %g\n", endingTheta);
2156     handler_->message(CLP_GENERAL,messages_)
2157       << line << CoinMessageEol;
2158     return problemStatus_;
2159}
2160/* Version of parametrics which reads from file
2161   See CbcClpParam.cpp for details of format
2162   Returns -2 if unable to open file */
2163int 
2164ClpSimplexOther::parametrics(const char * dataFile)
2165{
2166  int returnCode=-2;
2167  FILE *fp = fopen(dataFile, "r");
2168  char line[200];
2169  if (!fp) {
2170    handler_->message(CLP_UNABLE_OPEN, messages_)
2171      << dataFile << CoinMessageEol;
2172    return -2;
2173  }
2174
2175  if (!fgets(line, 200, fp)) {
2176    sprintf(line,"Empty parametrics file %s?",dataFile);
2177    handler_->message(CLP_GENERAL,messages_)
2178      << line << CoinMessageEol;
2179    fclose(fp);
2180    return -2;
2181  }
2182  char * pos = line;
2183  char * put = line;
2184  while (*pos >= ' ' && *pos != '\n') {
2185    if (*pos != ' ' && *pos != '\t') {
2186      *put = static_cast<char>(tolower(*pos));
2187      put++;
2188    }
2189    pos++;
2190  }
2191  *put = '\0';
2192  pos = line;
2193  double startTheta=0.0;
2194  double endTheta=0.0;
2195  double intervalTheta=COIN_DBL_MAX;
2196  int detail=0;
2197  bool good = true;
2198  while (good) {
2199    good=false;
2200    // check ROWS
2201    char * comma = strchr(pos, ',');
2202    if (!comma)
2203      break;
2204    *comma = '\0';
2205    if (strcmp(pos,"rows"))
2206      break;
2207    *comma = ',';
2208    pos = comma+1;
2209    // check lower theta
2210    comma = strchr(pos, ',');
2211    if (!comma)
2212      break;
2213    *comma = '\0';
2214    startTheta = atof(pos);
2215    *comma = ',';
2216    pos = comma+1;
2217    // check upper theta
2218    comma = strchr(pos, ',');
2219    good=true;
2220    if (comma)
2221      *comma = '\0';
2222    endTheta = atof(pos);
2223    if (comma) {
2224      *comma = ',';
2225      pos = comma+1;
2226      comma = strchr(pos, ',');
2227      if (comma)
2228        *comma = '\0';
2229      intervalTheta = atof(pos);
2230      if (comma) {
2231        *comma = ',';
2232        pos = comma+1;
2233        comma = strchr(pos, ',');
2234        if (comma)
2235          *comma = '\0';
2236        detail = atoi(pos);
2237        if (comma) 
2238        *comma = ',';
2239      }
2240    }
2241    break;
2242  }
2243  if (good) {
2244    if (startTheta<0.0||
2245        startTheta>endTheta||
2246        intervalTheta<0.0)
2247      good=false;
2248    if (detail<0||detail>1)
2249      good=false;
2250  }
2251  if (intervalTheta>=endTheta)
2252    intervalTheta=0.0;
2253  if (!good) {
2254    sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2255    handler_->message(CLP_GENERAL,messages_)
2256      << line << CoinMessageEol;
2257    fclose(fp);
2258    return -2;
2259  }
2260  if (!fgets(line, 200, fp)) {
2261    sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2262    handler_->message(CLP_GENERAL,messages_)
2263      << line << CoinMessageEol;
2264    fclose(fp);
2265    return -2;
2266  }
2267  double * lowerRowMove = NULL;
2268  double * upperRowMove = NULL;
2269  double * lowerColumnMove = NULL;
2270  double * upperColumnMove = NULL;
2271  double * objectiveMove = NULL;
2272  char saveLine[200];
2273  saveLine[0]='\0';
2274  std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2275  int gotRow[] = { -1, -1, -1, -1, -1};
2276  int orderRow[5];
2277  assert(sizeof(gotRow) == sizeof(orderRow));
2278  int nAcross = 0;
2279  pos = line;
2280  put = line;
2281  while (*pos >= ' ' && *pos != '\n') {
2282    if (*pos != ' ' && *pos != '\t') {
2283      *put = static_cast<char>(tolower(*pos));
2284      put++;
2285    }
2286    pos++;
2287  }
2288  *put = '\0';
2289  pos = line;
2290  int i;
2291  good = true;
2292  if (strncmp(line,"column",6)) {
2293    while (pos) {
2294      char * comma = strchr(pos, ',');
2295      if (comma)
2296        *comma = '\0';
2297      for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2298        if (headingsRow[i] == pos) {
2299          if (gotRow[i] < 0) {
2300            orderRow[nAcross] = i;
2301            gotRow[i] = nAcross++;
2302          } else {
2303            // duplicate
2304            good = false;
2305          }
2306          break;
2307        }
2308      }
2309      if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2310        good = false;
2311      if (comma) {
2312        *comma = ',';
2313        pos = comma + 1;
2314      } else {
2315        break;
2316      }
2317    }
2318    if (gotRow[0] < 0 && gotRow[1] < 0)
2319      good = false;
2320    if (gotRow[0] >= 0 && gotRow[1] >= 0)
2321      good = false;
2322    if (gotRow[0] >= 0 && !lengthNames())
2323      good = false;
2324    if (gotRow[4]<0) {
2325      if (gotRow[2] < 0 && gotRow[3] >= 0)
2326        good = false;
2327      else if (gotRow[3] < 0 && gotRow[2] >= 0)
2328        good = false;
2329    } else if (gotRow[2]>=0||gotRow[3]>=0) {
2330      good = false;
2331    }
2332    if (good) {
2333      char ** rowNames = new char * [numberRows_];
2334      int iRow;
2335      for (iRow = 0; iRow < numberRows_; iRow++) {
2336        rowNames[iRow] =
2337          CoinStrdup(rowName(iRow).c_str());
2338      }
2339      lowerRowMove = new double [numberRows_];
2340      memset(lowerRowMove,0,numberRows_*sizeof(double));
2341      upperRowMove = new double [numberRows_];
2342      memset(upperRowMove,0,numberRows_*sizeof(double));
2343      int nLine = 0;
2344      int nBadLine = 0;
2345      int nBadName = 0;
2346      bool goodLine=false;
2347      while (fgets(line, 200, fp)) {
2348        goodLine=true;
2349        if (!strncmp(line, "ENDATA", 6)||
2350            !strncmp(line, "COLUMN",6))
2351          break;
2352        goodLine=false;
2353        nLine++;
2354        iRow = -1;
2355        double upper = 0.0;
2356        double lower = 0.0;
2357        char * pos = line;
2358        char * put = line;
2359        while (*pos >= ' ' && *pos != '\n') {
2360          if (*pos != ' ' && *pos != '\t') {
2361            *put = *pos;
2362            put++;
2363          }
2364          pos++;
2365        }
2366        *put = '\0';
2367        pos = line;
2368        for (int i = 0; i < nAcross; i++) {
2369          char * comma = strchr(pos, ',');
2370          if (comma) {
2371            *comma = '\0';
2372          } else if (i < nAcross - 1) {
2373            nBadLine++;
2374            break;
2375          }
2376          switch (orderRow[i]) {
2377            // name
2378          case 0:
2379            // For large problems this could be slow
2380            for (iRow = 0; iRow < numberRows_; iRow++) {
2381              if (!strcmp(rowNames[iRow], pos))
2382                break;
2383            }
2384            if (iRow == numberRows_)
2385              iRow = -1;
2386            break;
2387            // number
2388          case 1:
2389            iRow = atoi(pos);
2390            if (iRow < 0 || iRow >= numberRows_)
2391              iRow = -1;
2392            break;
2393            // lower
2394          case 2:
2395            upper = atof(pos);
2396            break;
2397            // upper
2398          case 3:
2399            lower = atof(pos);
2400            break;
2401            // rhs
2402          case 4:
2403            lower = atof(pos);
2404            upper = lower;
2405            break;
2406          }
2407          if (comma) {
2408            *comma = ',';
2409            pos = comma + 1;
2410          }
2411        }
2412        if (iRow >= 0) {
2413          if (rowLower_[iRow]>-1.0e20)
2414            lowerRowMove[iRow] = lower;
2415          else
2416            lowerRowMove[iRow]=0.0;
2417          if (rowUpper_[iRow]<1.0e20)
2418            upperRowMove[iRow] = upper;
2419          else
2420            upperRowMove[iRow] = lower;
2421        } else {
2422          nBadName++;
2423          if(saveLine[0]=='\0')
2424            strcpy(saveLine,line);
2425        }
2426      }
2427      sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2428      handler_->message(CLP_GENERAL,messages_)
2429        << line << CoinMessageEol;
2430      if (nBadName) {
2431        sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2432        handler_->message(CLP_GENERAL,messages_)
2433          << line << CoinMessageEol;
2434        returnCode=-1;
2435        good=false;
2436      }
2437      for (iRow = 0; iRow < numberRows_; iRow++) {
2438        free(rowNames[iRow]);
2439      }
2440      delete [] rowNames;
2441    } else {
2442      sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2443      handler_->message(CLP_GENERAL,messages_)
2444        << line << CoinMessageEol;
2445      returnCode=-1;
2446      good=false;
2447    }
2448  }
2449  if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
2450    if (!fgets(line, 200, fp)) {
2451      sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2452      handler_->message(CLP_GENERAL,messages_)
2453        << line << CoinMessageEol;
2454      fclose(fp);
2455      return -2;
2456    }
2457    std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2458    saveLine[0]='\0';
2459    int gotColumn[] = { -1, -1, -1, -1, -1};
2460    int orderColumn[5];
2461    assert(sizeof(gotColumn) == sizeof(orderColumn));
2462    nAcross = 0;
2463    pos = line;
2464    put = line;
2465    while (*pos >= ' ' && *pos != '\n') {
2466      if (*pos != ' ' && *pos != '\t') {
2467        *put = static_cast<char>(tolower(*pos));
2468        put++;
2469      }
2470      pos++;
2471    }
2472    *put = '\0';
2473    pos = line;
2474    int i;
2475    if (strncmp(line,"endata",6)&&good) {
2476      while (pos) {
2477        char * comma = strchr(pos, ',');
2478        if (comma)
2479          *comma = '\0';
2480        for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2481          if (headingsColumn[i] == pos) {
2482            if (gotColumn[i] < 0) {
2483              orderColumn[nAcross] = i;
2484              gotColumn[i] = nAcross++;
2485            } else {
2486              // duplicate
2487              good = false;
2488            }
2489            break;
2490          }
2491        }
2492        if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2493          good = false;
2494        if (comma) {
2495          *comma = ',';
2496          pos = comma + 1;
2497        } else {
2498          break;
2499        }
2500      }
2501      if (gotColumn[0] < 0 && gotColumn[1] < 0)
2502        good = false;
2503      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2504        good = false;
2505      if (gotColumn[0] >= 0 && !lengthNames())
2506        good = false;
2507      if (good) {
2508        char ** columnNames = new char * [numberColumns_];
2509        int iColumn;
2510        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2511          columnNames[iColumn] =
2512            CoinStrdup(columnName(iColumn).c_str());
2513        }
2514        lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2515        memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2516        upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2517        memset(upperColumnMove,0,numberColumns_*sizeof(double));
2518        objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2519        memset(objectiveMove,0,numberColumns_*sizeof(double));
2520        int nLine = 0;
2521        int nBadLine = 0;
2522        int nBadName = 0;
2523        bool goodLine=false;
2524        while (fgets(line, 200, fp)) {
2525          goodLine=true;
2526          if (!strncmp(line, "ENDATA", 6))
2527            break;
2528          goodLine=false;
2529          nLine++;
2530          iColumn = -1;
2531          double upper = 0.0;
2532          double lower = 0.0;
2533          double obj =0.0;
2534          char * pos = line;
2535          char * put = line;
2536          while (*pos >= ' ' && *pos != '\n') {
2537            if (*pos != ' ' && *pos != '\t') {
2538              *put = *pos;
2539              put++;
2540            }
2541            pos++;
2542          }
2543          *put = '\0';
2544          pos = line;
2545          for (int i = 0; i < nAcross; i++) {
2546            char * comma = strchr(pos, ',');
2547            if (comma) {
2548              *comma = '\0';
2549            } else if (i < nAcross - 1) {
2550              nBadLine++;
2551              break;
2552            }
2553            switch (orderColumn[i]) {
2554              // name
2555            case 0:
2556              // For large problems this could be slow
2557              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2558                if (!strcmp(columnNames[iColumn], pos))
2559                  break;
2560              }
2561              if (iColumn == numberColumns_)
2562                iColumn = -1;
2563              break;
2564              // number
2565            case 1:
2566              iColumn = atoi(pos);
2567              if (iColumn < 0 || iColumn >= numberColumns_)
2568                iColumn = -1;
2569              break;
2570              // lower
2571            case 2:
2572              upper = atof(pos);
2573              break;
2574              // upper
2575            case 3:
2576              lower = atof(pos);
2577              break;
2578              // objective
2579            case 4:
2580              obj = atof(pos);
2581              upper = lower;
2582              break;
2583            }
2584            if (comma) {
2585              *comma = ',';
2586              pos = comma + 1;
2587            }
2588          }
2589          if (iColumn >= 0) {
2590            if (columnLower_[iColumn]>-1.0e20)
2591              lowerColumnMove[iColumn] = lower;
2592            else
2593              lowerColumnMove[iColumn]=0.0;
2594            if (columnUpper_[iColumn]<1.0e20)
2595              upperColumnMove[iColumn] = upper;
2596            else
2597              upperColumnMove[iColumn] = lower;
2598            objectiveMove[iColumn] = obj;
2599          } else {
2600            nBadName++;
2601            if(saveLine[0]=='\0')
2602              strcpy(saveLine,line);
2603          }
2604        }
2605        sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2606        handler_->message(CLP_GENERAL,messages_)
2607          << line << CoinMessageEol;
2608        if (nBadName) {
2609          sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2610          handler_->message(CLP_GENERAL,messages_)
2611            << line << CoinMessageEol;
2612          returnCode=-1;
2613          good=false;
2614        }
2615        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2616          free(columnNames[iColumn]);
2617        }
2618        delete [] columnNames;
2619      } else {
2620        sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2621        handler_->message(CLP_GENERAL,messages_)
2622          << line << CoinMessageEol;
2623        returnCode=-1;
2624        good=false;
2625      }
2626    }
2627  }
2628  returnCode=-1;
2629  if (good) {
2630    // clean arrays
2631    if (lowerRowMove) {
2632      bool empty=true;
2633      for (int i=0;i<numberRows_;i++) {
2634        if (lowerRowMove[i]) {
2635          empty=false;
2636        break;
2637        }
2638      }
2639      if (empty) {
2640        delete [] lowerRowMove;
2641        lowerRowMove=NULL;
2642      }
2643    }
2644    if (upperRowMove) {
2645      bool empty=true;
2646      for (int i=0;i<numberRows_;i++) {
2647        if (upperRowMove[i]) {
2648          empty=false;
2649        break;
2650        }
2651      }
2652      if (empty) {
2653        delete [] upperRowMove;
2654        upperRowMove=NULL;
2655      }
2656    }
2657    if (lowerColumnMove) {
2658      bool empty=true;
2659      for (int i=0;i<numberColumns_;i++) {
2660        if (lowerColumnMove[i]) {
2661          empty=false;
2662        break;
2663        }
2664      }
2665      if (empty) {
2666        delete [] lowerColumnMove;
2667        lowerColumnMove=NULL;
2668      }
2669    }
2670    if (upperColumnMove) {
2671      bool empty=true;
2672      for (int i=0;i<numberColumns_;i++) {
2673        if (upperColumnMove[i]) {
2674          empty=false;
2675        break;
2676        }
2677      }
2678      if (empty) {
2679        delete [] upperColumnMove;
2680        upperColumnMove=NULL;
2681      }
2682    }
2683    if (objectiveMove) {
2684      bool empty=true;
2685      for (int i=0;i<numberColumns_;i++) {
2686        if (objectiveMove[i]) {
2687          empty=false;
2688        break;
2689        }
2690      }
2691      if (empty) {
2692        delete [] objectiveMove;
2693        objectiveMove=NULL;
2694      }
2695    }
2696    int saveScaling = scalingFlag_;
2697    scalingFlag_ = 0;
2698    int saveLogLevel = handler_->logLevel();
2699    if (detail>0&&!intervalTheta)
2700      handler_->setLogLevel(3);
2701    else
2702      handler_->setLogLevel(1);
2703    returnCode = parametrics(startTheta,endTheta,intervalTheta,
2704                             lowerColumnMove,upperColumnMove,
2705                             lowerRowMove,upperRowMove,
2706                             objectiveMove);
2707    scalingFlag_ = saveScaling;
2708    handler_->setLogLevel(saveLogLevel);
2709  }
2710  delete [] lowerRowMove;
2711  delete [] upperRowMove;
2712  delete [] lowerColumnMove;
2713  delete [] upperColumnMove;
2714  delete [] objectiveMove;
2715  fclose(fp);
2716  return returnCode;
2717}
2718int
2719ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
2720                                 const double * lowerChange, const double * upperChange,
2721                                 const double * changeObjective, ClpDataSave & data,
2722                                 bool canTryQuick)
2723{
2724  double startingTheta = paramData.startingTheta;
2725  double & endingTheta = paramData.endingTheta;
2726     // stuff is already at starting
2727     // For this crude version just try and go to end
2728     double change = 0.0;
2729     if (reportIncrement && canTryQuick) {
2730          endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2731          change = endingTheta - startingTheta;
2732     }
2733     int numberTotal = numberRows_ + numberColumns_;
2734     int i;
2735     for ( i = 0; i < numberTotal; i++) {
2736          lower_[i] += change * lowerChange[i];
2737          upper_[i] += change * upperChange[i];
2738          switch(getStatus(i)) {
2739
2740          case basic:
2741          case isFree:
2742          case superBasic:
2743               break;
2744          case isFixed:
2745          case atUpperBound:
2746               solution_[i] = upper_[i];
2747               break;
2748          case atLowerBound:
2749               solution_[i] = lower_[i];
2750               break;
2751          }
2752          cost_[i] += change * changeObjective[i];
2753     }
2754     problemStatus_ = -1;
2755
2756     // This says whether to restore things etc
2757     // startup will have factorized so can skip
2758     int factorType = 0;
2759     // Start check for cycles
2760     progress_.startCheck();
2761     // Say change made on first iteration
2762     changeMade_ = 1;
2763     /*
2764       Status of problem:
2765       0 - optimal
2766       1 - infeasible
2767       2 - unbounded
2768       -1 - iterating
2769       -2 - factorization wanted
2770       -3 - redo checking without factorization
2771       -4 - looks infeasible
2772     */
2773     while (problemStatus_ < 0) {
2774          int iRow, iColumn;
2775          // clear
2776          for (iRow = 0; iRow < 4; iRow++) {
2777               rowArray_[iRow]->clear();
2778          }
2779
2780          for (iColumn = 0; iColumn < 2; iColumn++) {
2781               columnArray_[iColumn]->clear();
2782          }
2783
2784          // give matrix (and model costs and bounds a chance to be
2785          // refreshed (normally null)
2786          matrix_->refresh(this);
2787          // may factorize, checks if problem finished
2788          statusOfProblemInParametrics(factorType, data);
2789          // Say good factorization
2790          factorType = 1;
2791          if (data.sparseThreshold_) {
2792               // use default at present
2793               factorization_->sparseThreshold(0);
2794               factorization_->goSparse();
2795          }
2796
2797          // exit if victory declared
2798          if (problemStatus_ >= 0 && 
2799              (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
2800               break;
2801
2802          // test for maximum iterations
2803          if (hitMaximumIterations()) {
2804               problemStatus_ = 3;
2805               break;
2806          }
2807          // Check event
2808          {
2809               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2810               if (status >= 0) {
2811                    problemStatus_ = 5;
2812                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
2813                    break;
2814               }
2815          }
2816          // Do iterations
2817          problemStatus_=-1;
2818          if (canTryQuick) {
2819               double * saveDuals = NULL;
2820               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2821          } else {
2822               whileIterating(paramData, reportIncrement,
2823                              changeObjective);
2824               startingTheta = endingTheta;
2825          }
2826     }
2827     if (!problemStatus_) {
2828          theta_ = change + startingTheta;
2829          eventHandler_->event(ClpEventHandler::theta);
2830          return 0;
2831     } else if (problemStatus_ == 10) {
2832          return -1;
2833     } else {
2834          return problemStatus_;
2835     }
2836}
2837/* Parametrics
2838   The code uses current bounds + theta * change (if change array not NULL)
2839   It starts at startingTheta and returns ending theta in endingTheta.
2840   If it can not reach input endingTheta return code will be 1 for infeasible,
2841   2 for unbounded, if error on ranges -1,  otherwise 0.
2842   Event handler may do more
2843   On exit endingTheta is maximum reached (can be used for next startingTheta)
2844*/
2845int
2846ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
2847                             const double * lowerChangeBound, const double * upperChangeBound,
2848                             const double * lowerChangeRhs, const double * upperChangeRhs)
2849{
2850  int savePerturbation = perturbation_;
2851  perturbation_ = 102; // switch off
2852  algorithm_ = -1;
2853  // extra region
2854  int maximumPivots = factorization_->maximumPivots();
2855  int numberDense = factorization_->numberDense();
2856  int length = numberRows_ + numberDense + maximumPivots;
2857  assert (!rowArray_[4]);
2858  rowArray_[4]=new CoinIndexedVector(length);
2859  assert (!rowArray_[5]);
2860  rowArray_[5]=new CoinIndexedVector(length);
2861
2862  // save data
2863  ClpDataSave data = saveData();
2864  int numberTotal = numberRows_ + numberColumns_;
2865  int ratio = (2*sizeof(int))/sizeof(double);
2866  assert (ratio==1||ratio==2);
2867  // allow for unscaled - even if not needed
2868  int lengthArrays = 4*numberTotal+(3*numberTotal+2)*ratio+2*numberRows_+1;
2869  /*
2870    Save information and modify
2871  */
2872  double * saveLower = new double [lengthArrays];
2873  double * saveUpper = new double [lengthArrays];
2874  double * lowerCopy = saveLower+2*numberTotal;
2875  double * upperCopy = saveUpper+2*numberTotal;
2876  double * lowerChange = saveLower+numberTotal;
2877  double * upperChange = saveUpper+numberTotal;
2878  double * lowerGap = saveLower+4*numberTotal;
2879  double * lowerCoefficient = lowerGap+numberRows_;
2880  double * upperGap = saveUpper+4*numberTotal;
2881  double * upperCoefficient = upperGap+numberRows_;
2882  int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal+2*numberRows_))+2;
2883  int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal+2*numberRows_))+2;
2884  int * lowerActive = lowerList+numberTotal+1;
2885  int * upperActive = upperList+numberTotal+1;
2886  // To mark as odd
2887  char * markDone = reinterpret_cast<char *>(lowerActive+numberTotal);
2888  //memset(markDone,0,numberTotal);
2889  int * backwardBasic = upperActive+numberTotal;
2890  parametricsData paramData;
2891  paramData.lowerChange = lowerChange;
2892  paramData.lowerList=lowerList;
2893  paramData.upperChange = upperChange;
2894  paramData.upperList=upperList;
2895  paramData.markDone=markDone;
2896  paramData.backwardBasic=backwardBasic;
2897  paramData.lowerActive = lowerActive;
2898  paramData.lowerGap = lowerGap;
2899  paramData.lowerCoefficient = lowerCoefficient;
2900  paramData.upperActive = upperActive;
2901  paramData.upperGap = upperGap;
2902  paramData.upperCoefficient = upperCoefficient;
2903  // Find theta when bounds will cross over and create arrays
2904  memset(lowerChange, 0, numberTotal * sizeof(double));
2905  memset(upperChange, 0, numberTotal * sizeof(double));
2906  if (lowerChangeBound)
2907    memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
2908  if (upperChangeBound)
2909    memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
2910  if (lowerChangeRhs)
2911    memcpy(lowerChange+numberColumns_,
2912           lowerChangeRhs,numberRows_*sizeof(double));
2913  if (upperChangeRhs)
2914    memcpy(upperChange+numberColumns_,
2915           upperChangeRhs,numberRows_*sizeof(double));
2916  int nLowerChange=0;
2917  int nUpperChange=0;
2918  for (int i=0;i<numberColumns_;i++) {
2919    if (lowerChange[i]) { 
2920      lowerList[nLowerChange++]=i;
2921    }
2922    if (upperChange[i]) { 
2923      upperList[nUpperChange++]=i;
2924    }
2925  }
2926  lowerList[-2]=nLowerChange;
2927  upperList[-2]=nUpperChange;
2928  for (int i=numberColumns_;i<numberTotal;i++) {
2929    if (lowerChange[i]) { 
2930      lowerList[nLowerChange++]=i;
2931    }
2932    if (upperChange[i]) { 
2933      upperList[nUpperChange++]=i;
2934    }
2935  }
2936  lowerList[-1]=nLowerChange;
2937  upperList[-1]=nUpperChange;
2938  memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
2939  memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
2940  memcpy(lowerCopy+numberColumns_,
2941         rowLower_,numberRows_*sizeof(double));
2942  memcpy(upperCopy+numberColumns_,
2943         rowUpper_,numberRows_*sizeof(double));
2944  {
2945    //  extra for unscaled
2946    double * unscaledCopy;
2947    unscaledCopy = lowerCopy + numberTotal; 
2948    memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
2949    memcpy(unscaledCopy+numberColumns_,
2950           rowLower_,numberRows_*sizeof(double));
2951    unscaledCopy = upperCopy + numberTotal; 
2952    memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
2953    memcpy(unscaledCopy+numberColumns_,
2954           rowUpper_,numberRows_*sizeof(double));
2955  }
2956  double maxTheta = 1.0e50;
2957  for (int iRow = 0; iRow < numberRows_; iRow++) {
2958    double lower = rowLower_[iRow];
2959    double upper = rowUpper_[iRow];
2960    if (lower<-1.0e30)
2961      lowerChange[numberColumns_+iRow]=0.0;
2962    double chgLower = lowerChange[numberColumns_+iRow];
2963    if (upper>1.0e30)
2964      upperChange[numberColumns_+iRow]=0.0;
2965    double chgUpper = upperChange[numberColumns_+iRow];
2966    if (lower > -1.0e30 && upper < 1.0e30) {
2967      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2968        maxTheta = (upper - lower) / (chgLower - chgUpper);
2969      }
2970    }
2971    lower+=startingTheta*chgLower;
2972    upper+=startingTheta*chgUpper;
2973    if (lower > upper) {
2974      maxTheta = -1.0;
2975      break;
2976    }
2977    rowLower_[iRow]=lower;
2978    rowUpper_[iRow]=upper;
2979  }
2980  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2981    double lower = columnLower_[iColumn];
2982    double upper = columnUpper_[iColumn];
2983    if (lower<-1.0e30)
2984      lowerChange[iColumn]=0.0;
2985    double chgLower = lowerChange[iColumn];
2986    if (upper>1.0e30)
2987      upperChange[iColumn]=0.0;
2988    double chgUpper = upperChange[iColumn];
2989    if (lower > -1.0e30 && upper < 1.0e30) {
2990      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2991        maxTheta = (upper - lower) / (chgLower - chgUpper);
2992      }
2993    }
2994    lower+=startingTheta*chgLower;
2995    upper+=startingTheta*chgUpper;
2996    if (lower > upper) {
2997      maxTheta = -1.0;
2998      break;
2999    }
3000    columnLower_[iColumn]=lower;
3001    columnUpper_[iColumn]=upper;
3002  }
3003  if (maxTheta == 1.0e50)
3004    maxTheta = COIN_DBL_MAX;
3005  int returnCode=0;
3006  if (maxTheta < 0.0) {
3007    // bad ranges or initial
3008    returnCode = -1;
3009  }
3010  if (maxTheta < endingTheta) {
3011    char line[100];
3012    sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
3013            endingTheta,maxTheta);
3014    handler_->message(CLP_GENERAL,messages_)
3015      << line << CoinMessageEol;
3016    endingTheta = maxTheta;
3017  }
3018  if (endingTheta < startingTheta) {
3019    // bad initial
3020    returnCode = -2;
3021  }
3022  bool swapped=false;
3023  // Dantzig
3024#define ALL_DANTZIG
3025#ifdef ALL_DANTZIG
3026  ClpDualRowPivot * savePivot = dualRowPivot_;
3027  dualRowPivot_ = new ClpDualRowDantzig();
3028  dualRowPivot_->setModel(this);
3029#else
3030  ClpDualRowPivot * savePivot = NULL;
3031#endif
3032  if (!returnCode) {
3033    assert (objective_->type()==1);
3034    objective_->setType(2); // in case matrix empty
3035    returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
3036    objective_->setType(1);
3037    if (!returnCode) {
3038      double saveDualBound=dualBound_;
3039      dualBound_=CoinMax(dualBound_,1.0e15);
3040      swapped=true;
3041      double * temp;
3042      memcpy(saveLower,lower_,numberTotal*sizeof(double));
3043      temp=saveLower;
3044      saveLower=lower_;
3045      lower_=temp;
3046      //columnLowerWork_ = lower_;
3047      //rowLowerWork_ = lower_ + numberColumns_;
3048      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
3049      temp=saveUpper;
3050      saveUpper=upper_;
3051      upper_=temp;
3052      //columnUpperWork_ = upper_;
3053      //rowUpperWork_ = upper_ + numberColumns_;
3054      if (rowScale_) {
3055        // scale saved and change arrays
3056        double * lowerChange = lower_+numberTotal;
3057        double * upperChange = upper_+numberTotal;
3058        double * lowerSave = lowerChange+numberTotal;
3059        double * upperSave = upperChange+numberTotal;
3060        for (int i=0;i<numberColumns_;i++) {
3061          double multiplier = inverseColumnScale_[i];
3062          if (lowerSave[i]>-1.0e20)
3063            lowerSave[i] *= multiplier;
3064          if (upperSave[i]<1.0e20)
3065            upperSave[i] *= multiplier;
3066          lowerChange[i] *= multiplier;
3067          upperChange[i] *= multiplier;
3068        }
3069        lowerChange += numberColumns_;
3070        upperChange += numberColumns_;
3071        lowerSave += numberColumns_;
3072        upperSave += numberColumns_;
3073        for (int i=0;i<numberRows_;i++) {
3074          double multiplier = rowScale_[i];
3075          if (lowerSave[i]>-1.0e20)
3076            lowerSave[i] *= multiplier;
3077          if (upperSave[i]<1.0e20)
3078            upperSave[i] *= multiplier;
3079          lowerChange[i] *= multiplier;
3080          upperChange[i] *= multiplier;
3081        }
3082      }
3083      //double saveEndingTheta = endingTheta;
3084      double * saveDuals = NULL;
3085      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3086      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3087        // probably can get rid of this if we adjust every change in theta
3088        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3089        //   sumPrimalInfeasibilities_);
3090        int pass=100;
3091        while(sumPrimalInfeasibilities_) {
3092          pass--;
3093          if (!pass)
3094            break;
3095          problemStatus_=-1;
3096          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3097            double value=solution_[iSequence];
3098            // remember scaling
3099            if (value<lower_[iSequence]-1.0e-9) {
3100              lower_[iSequence]=value;
3101              lowerCopy[iSequence]=value;
3102            } else if (value>upper_[iSequence]+1.0e-9) {
3103              upper_[iSequence]=value;
3104              upperCopy[iSequence]=value;
3105            }
3106          }
3107          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3108        }
3109      }
3110      assert (!problemStatus_);
3111      if (nLowerChange||nUpperChange) {
3112#ifndef ALL_DANTZIG
3113        // Dantzig
3114        savePivot = dualRowPivot_;
3115        dualRowPivot_ = new ClpDualRowDantzig();
3116        dualRowPivot_->setModel(this);
3117#endif
3118        //for (int i=0;i<numberRows_+numberColumns_;i++)
3119        //setFakeBound(i, noFake);
3120        // Now do parametrics
3121        handler_->message(CLP_PARAMETRICS_STATS, messages_)
3122          << startingTheta << objectiveValue() << CoinMessageEol;
3123        bool canSkipFactorization=true;
3124        while (!returnCode) {
3125                 paramData.startingTheta=startingTheta;
3126                 paramData.endingTheta=endingTheta;
3127                 returnCode = parametricsLoop(paramData,
3128                                       data,canSkipFactorization);
3129                 startingTheta=paramData.startingTheta;
3130                 endingTheta=paramData.endingTheta;
3131          canSkipFactorization=false;
3132          if (!returnCode) {
3133            //startingTheta = endingTheta;
3134            //endingTheta = saveEndingTheta;
3135            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3136              << startingTheta << objectiveValue() << CoinMessageEol;
3137            if (startingTheta >= endingTheta-primalTolerance_
3138                ||problemStatus_==2)
3139              break;
3140          } else if (returnCode == -1) {
3141            // trouble - do external solve
3142            abort(); //needToDoSomething = true;
3143          } else if (problemStatus_==1) {
3144            // can't move any further
3145            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3146              << endingTheta << objectiveValue() << CoinMessageEol;
3147            problemStatus_=0;
3148          }
3149        }
3150      }
3151      dualBound_ = saveDualBound;
3152      //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3153    } else {
3154      // check if empty
3155      //if (!numberRows_||!matrix_->getNumElements()) {
3156        // success
3157#ifdef CLP_USER_DRIVEN
3158      //theta_ = endingTheta;
3159      //eventHandler_->event(ClpEventHandler::theta);
3160#endif
3161      //}
3162    }
3163    if (problemStatus_==2) {
3164      delete [] ray_;
3165      ray_ = new double [numberColumns_];
3166    }
3167    if (swapped&&lower_) {
3168      double * temp=saveLower;
3169      saveLower=lower_;
3170      lower_=temp;
3171      temp=saveUpper;
3172      saveUpper=upper_;
3173      upper_=temp;
3174    }
3175    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
3176  }   
3177  if (!scalingFlag_) {
3178    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3179    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3180    memcpy(rowLower_,lowerCopy+numberColumns_,
3181           numberRows_*sizeof(double));
3182    memcpy(rowUpper_,upperCopy+numberColumns_,
3183           numberRows_*sizeof(double));
3184  } else {
3185    //  extra for unscaled
3186    double * unscaledCopy;
3187    unscaledCopy = lowerCopy + numberTotal; 
3188    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3189    memcpy(rowLower_,unscaledCopy+numberColumns_,
3190           numberRows_*sizeof(double));
3191    unscaledCopy = upperCopy + numberTotal; 
3192    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3193    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3194           numberRows_*sizeof(double));
3195  }
3196  delete [] saveLower;
3197  delete [] saveUpper;
3198#ifdef ALL_DANTZIG
3199  if (savePivot) {
3200#endif
3201    delete dualRowPivot_;
3202    dualRowPivot_ = savePivot;
3203#ifdef ALL_DANTZIG
3204  }
3205#endif
3206  // Restore any saved stuff
3207  restoreData(data);
3208  perturbation_ = savePerturbation;
3209  delete rowArray_[4];
3210  rowArray_[4]=NULL;
3211  delete rowArray_[5];
3212  rowArray_[5]=NULL;
3213  char line[100];
3214  sprintf(line,"Ending theta %g\n", endingTheta);
3215  handler_->message(CLP_GENERAL,messages_)
3216    << line << CoinMessageEol;
3217  return problemStatus_;
3218}
3219int
3220ClpSimplexOther::parametricsLoop(parametricsData & paramData,
3221                                 ClpDataSave & data,bool canSkipFactorization)
3222{
3223  double & startingTheta = paramData.startingTheta;
3224  double & endingTheta = paramData.endingTheta;
3225  int numberTotal = numberRows_+numberColumns_;
3226  // stuff is already at starting
3227  const int * lowerList = paramData.lowerList;
3228  const int * upperList = paramData.upperList;
3229  problemStatus_ = -1;
3230  //double saveEndingTheta=endingTheta;
3231
3232  // This says whether to restore things etc
3233  // startup will have factorized so can skip
3234  int factorType = 0;
3235  // Start check for cycles
3236  progress_.startCheck();
3237  // Say change made on first iteration
3238  changeMade_ = 1;
3239  /*
3240    Status of problem:
3241    0 - optimal
3242    1 - infeasible
3243    2 - unbounded
3244    -1 - iterating
3245    -2 - factorization wanted
3246    -3 - redo checking without factorization
3247    -4 - looks infeasible
3248  */
3249  while (problemStatus_ < 0) {
3250    int iRow, iColumn;
3251    // clear
3252    for (iRow = 0; iRow < 6; iRow++) {
3253      rowArray_[iRow]->clear();
3254    }
3255   
3256    for (iColumn = 0; iColumn < 2; iColumn++) {
3257      columnArray_[iColumn]->clear();
3258    }
3259   
3260    // give matrix (and model costs and bounds a chance to be
3261    // refreshed (normally null)
3262    matrix_->refresh(this);
3263    // may factorize, checks if problem finished
3264    if (!canSkipFactorization)
3265      statusOfProblemInParametrics(factorType, data);
3266    canSkipFactorization=false;
3267    if (numberPrimalInfeasibilities_) {
3268      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3269        // treat as success
3270        problemStatus_=0;
3271        endingTheta=startingTheta;
3272        break;
3273      }
3274      // probably can get rid of this if we adjust every change in theta
3275      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3276      //     sumPrimalInfeasibilities_);
3277      const double * lowerChange = lower_+numberTotal;
3278      const double * upperChange = upper_+numberTotal;
3279      const double * startLower = lowerChange+numberTotal;
3280      const double * startUpper = upperChange+numberTotal;
3281      //startingTheta -= 1.0e-7;
3282      int nLowerChange = lowerList[-1];
3283      for (int i = 0; i < nLowerChange; i++) {
3284        int iSequence = lowerList[i];
3285        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3286      }
3287      int nUpperChange = upperList[-1];
3288      for (int i = 0; i < nUpperChange; i++) {
3289        int iSequence = upperList[i];
3290        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3291      }
3292      // adjust rhs in case dual uses
3293      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3294      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3295      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3296      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3297      if (rowScale_) {
3298        for (int i=0;i<numberColumns_;i++) {
3299          double multiplier = columnScale_[i];
3300          if (columnLower_[i]>-1.0e20)
3301            columnLower_[i] *= multiplier;
3302          if (columnUpper_[i]<1.0e20)
3303            columnUpper_[i] *= multiplier;
3304        }
3305        for (int i=0;i<numberRows_;i++) {
3306          double multiplier = inverseRowScale_[i];
3307          if (rowLower_[i]>-1.0e20)
3308            rowLower_[i] *= multiplier;
3309          if (rowUpper_[i]<1.0e20)
3310            rowUpper_[i] *= multiplier;
3311        }
3312      }
3313      double * saveDuals = NULL;
3314      problemStatus_=-1;
3315      ClpObjective * saveObjective = objective_;
3316      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3317      if (saveObjective!=objective_) {
3318        delete objective_;
3319        objective_=saveObjective;
3320      }
3321      int pass=100;
3322      double moved=0.0;
3323      while(sumPrimalInfeasibilities_) {
3324        //printf("INFEAS pass %d %d %g\n",100-pass,numberPrimalInfeasibilities_,
3325        //     sumPrimalInfeasibilities_);
3326        pass--;
3327        if (!pass)
3328          break;
3329        problemStatus_=-1;
3330        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3331          double value=solution_[iSequence];
3332          if (value<lower_[iSequence]-1.0e-9) {
3333            moved += lower_[iSequence]-value;
3334            lower_[iSequence]=value;
3335          } else if (value>upper_[iSequence]+1.0e-9) {
3336            moved += upper_[iSequence]-value;
3337            upper_[iSequence]=value;
3338          }
3339        }
3340        if (!moved) {
3341          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3342            double value=solution_[iSequence];
3343            if (value<lower_[iSequence]-1.0e-9) {
3344              moved += lower_[iSequence]-value;
3345              lower_[iSequence]=value;
3346            } else if (value>upper_[iSequence]+1.0e-9) {
3347              moved += upper_[iSequence]-value;
3348              upper_[iSequence]=value;
3349            }
3350          }
3351        }
3352        assert (moved);
3353        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3354      }
3355      // adjust
3356      //printf("Should adjust - moved %g\n",moved);
3357    }
3358    // Say good factorization
3359    factorType = 1;
3360    if (data.sparseThreshold_) {
3361      // use default at present
3362      factorization_->sparseThreshold(0);
3363      factorization_->goSparse();
3364    }
3365   
3366    // exit if victory declared
3367    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3368      break;
3369   
3370    // test for maximum iterations
3371    if (hitMaximumIterations()) {
3372      problemStatus_ = 3;
3373      break;
3374    }
3375#ifdef CLP_USER_DRIVEN
3376    // Check event
3377    {
3378      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3379      if (status >= 0) {
3380        problemStatus_ = 5;
3381        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3382        break;
3383      }
3384    }
3385#endif
3386    // Do iterations
3387    problemStatus_=-1;
3388    whileIterating(paramData, 0.0,
3389                   NULL);
3390    //startingTheta = endingTheta;
3391    //endingTheta = saveEndingTheta;
3392  }
3393  if (!problemStatus_/*||problemStatus_==2*/) {
3394    theta_ = endingTheta;
3395#ifdef CLP_USER_DRIVEN
3396    {
3397      double saveTheta=theta_;
3398      theta_ = endingTheta;
3399      int status=eventHandler_->event(ClpEventHandler::theta);
3400      if (status>=0&&status<10) {
3401        endingTheta=theta_;
3402        theta_=saveTheta;
3403        problemStatus_=-1;
3404      } else {
3405        if (status>=10) {
3406          problemStatus_=status-10;
3407          startingTheta=endingTheta;
3408        }
3409        theta_=saveTheta;
3410      }
3411    }
3412#endif
3413    return 0;
3414  } else if (problemStatus_ == 10) {
3415    return -1;
3416  } else {
3417    return problemStatus_;
3418  }
3419}
3420/* Checks if finished.  Updates status */
3421void
3422ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3423{
3424     if (type == 2) {
3425          // trouble - go to recovery
3426          problemStatus_ = 10;
3427          return;
3428     }
3429     if (problemStatus_ > -3 || factorization_->pivots()) {
3430          // factorize
3431          // later on we will need to recover from singularities
3432          // also we could skip if first time
3433          if (type) {
3434               // is factorization okay?
3435               if (internalFactorize(1)) {
3436                    // trouble - go to recovery
3437                    problemStatus_ = 10;
3438                    return;
3439               }
3440          }
3441          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3442               problemStatus_ = -3;
3443     }
3444     // at this stage status is -3 or -4 if looks infeasible
3445     // get primal and dual solutions
3446     gutsOfSolution(NULL, NULL);
3447     double realDualInfeasibilities = sumDualInfeasibilities_;
3448     // If bad accuracy treat as singular
3449     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3450          // trouble - go to recovery
3451          problemStatus_ = 10;
3452          return;
3453     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3454          // Can reduce tolerance
3455          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3456          factorization_->pivotTolerance(newTolerance);
3457     }
3458     // Check if looping
3459     int loop;
3460     if (type != 2)
3461          loop = progress_.looping();
3462     else
3463          loop = -1;
3464     if (loop >= 0) {
3465          problemStatus_ = loop; //exit if in loop
3466          if (!problemStatus_) {
3467               // declaring victory
3468               numberPrimalInfeasibilities_ = 0;
3469               sumPrimalInfeasibilities_ = 0.0;
3470          } else {
3471               problemStatus_ = 10; // instead - try other algorithm
3472          }
3473          return;
3474     } else if (loop < -1) {
3475          // something may have changed
3476          gutsOfSolution(NULL, NULL);
3477     }
3478     progressFlag_ = 0; //reset progress flag
3479     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3480          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3481                    << numberIterations_ << objectiveValue();
3482          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3483                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3484          handler_->printing(sumDualInfeasibilities_ > 0.0)
3485                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3486          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3487                             < numberDualInfeasibilities_)
3488                    << numberDualInfeasibilitiesWithoutFree_;
3489          handler_->message() << CoinMessageEol;
3490     }
3491#ifdef CLP_USER_DRIVEN
3492     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3493       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3494       if (status>=0) {
3495         // fix up
3496         for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3497           double value=solution_[iSequence];
3498           if (value<=lower_[iSequence]-primalTolerance_) {
3499             lower_[iSequence]=value;
3500           } else if (value>=upper_[iSequence]+primalTolerance_) {
3501             upper_[iSequence]=value;
3502           }
3503         }
3504         numberPrimalInfeasibilities_ = 0;
3505         sumPrimalInfeasibilities_ = 0.0;
3506       }
3507     }
3508#endif
3509     /* If we are primal feasible and any dual infeasibilities are on
3510        free variables then it is better to go to primal */
3511     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3512               numberDualInfeasibilities_) {
3513          problemStatus_ = 10;
3514          return;
3515     }
3516
3517     // check optimal
3518     // give code benefit of doubt
3519     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3520               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3521          // say optimal (with these bounds etc)
3522          numberDualInfeasibilities_ = 0;
3523          sumDualInfeasibilities_ = 0.0;
3524          numberPrimalInfeasibilities_ = 0;
3525          sumPrimalInfeasibilities_ = 0.0;
3526     }
3527     if (dualFeasible() || problemStatus_ == -4) {
3528          progress_.modifyObjective(objectiveValue_
3529                                    - sumDualInfeasibilities_ * dualBound_);
3530     }
3531     if (numberPrimalInfeasibilities_) {
3532          if (problemStatus_ == -4 || problemStatus_ == -5) {
3533               problemStatus_ = 1; // infeasible
3534          }
3535     } else if (numberDualInfeasibilities_) {
3536          // clean up
3537          problemStatus_ = 10;
3538     } else {
3539          problemStatus_ = 0;
3540     }
3541     lastGoodIteration_ = numberIterations_;
3542     if (problemStatus_ < 0) {
3543          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3544          if (sumDualInfeasibilities_)
3545               numberDualInfeasibilities_ = 1;
3546     }
3547     // Allow matrices to be sorted etc
3548     int fake = -999; // signal sort
3549     matrix_->correctSequence(this, fake, fake);
3550}
3551//static double lastThetaX=0.0;
3552/* This has the flow between re-factorizations
3553   Reasons to come out:
3554   -1 iterations etc
3555   -2 inaccuracy
3556   -3 slight inaccuracy (and done iterations)
3557   +0 looks optimal (might be unbounded - but we will investigate)
3558   +1 looks infeasible
3559   +3 max iterations
3560   +4 accuracy problems
3561*/
3562int
3563ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
3564                                const double * /*changeObjective*/)
3565{
3566  double & startingTheta = paramData.startingTheta;
3567  double & endingTheta = paramData.endingTheta;
3568  const double * lowerChange = paramData.lowerChange;
3569  const double * upperChange = paramData.upperChange;
3570  int numberTotal = numberColumns_ + numberRows_;
3571  const int * lowerList = paramData.lowerList;
3572  const int * upperList = paramData.upperList;
3573  //#define CLP_PARAMETRIC_DENSE_ARRAYS 2
3574#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3575  double * lowerGap = paramData.lowerGap;
3576  double * upperGap = paramData.upperGap;
3577  double * lowerCoefficient = paramData.lowerCoefficient;
3578  double * upperCoefficient = paramData.upperCoefficient;
3579#endif
3580  // do basic pointers
3581  int * backwardBasic = paramData.backwardBasic;
3582  for (int i=0;i<numberTotal;i++)
3583    backwardBasic[i]=-1;
3584  for (int i=0;i<numberRows_;i++) {
3585    int iPivot=pivotVariable_[i];
3586    backwardBasic[iPivot]=i;
3587  }
3588     {
3589          int i;
3590          for (i = 0; i < 4; i++) {
3591               rowArray_[i]->clear();
3592          }
3593          for (i = 0; i < 2; i++) {
3594               columnArray_[i]->clear();
3595          }
3596     }
3597     // if can't trust much and long way from optimal then relax
3598     if (largestPrimalError_ > 10.0)
3599          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3600     else
3601          factorization_->relaxAccuracyCheck(1.0);
3602     // status stays at -1 while iterating, >=0 finished, -2 to invert
3603     // status -3 to go to top without an invert
3604     int returnCode = -1;
3605     double lastTheta = startingTheta;
3606     double useTheta = startingTheta;
3607     while (problemStatus_ == -1) {
3608          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3609          // Get theta for bounds - we know can't crossover
3610          int pivotType = nextTheta(1, increaseTheta, paramData,
3611                                     NULL);
3612          useTheta += theta_;
3613          double change = useTheta - lastTheta;
3614          if (change>1.0e-14) {
3615            int n;
3616            n=lowerList[-1];
3617            for (int i=0;i<n;i++) {
3618              int iSequence = lowerList[i];
3619              double thisChange = change * lowerChange[iSequence];
3620              double newValue = lower_[iSequence] + thisChange;
3621              lower_[iSequence] = newValue;
3622#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3623              if (getStatus(iSequence)==basic) {
3624                int iRow=backwardBasic[iSequence];
3625                lowerGap[iRow] -= thisChange;
3626              } else if(getStatus(iSequence)==atLowerBound) {
3627                solution_[iSequence] = newValue;
3628              }
3629#else
3630              if(getStatus(iSequence)==atLowerBound) {
3631                solution_[iSequence] = newValue;
3632              }
3633#endif
3634#if 0
3635              // may have to adjust other bound
3636              double otherValue = upper_[iSequence];
3637              if (otherValue-newValue<dualBound_) {
3638                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3639                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3640                //ClpTraceDebug (fabs(lower_[iSequence]-newValue)<1.0e-5);
3641              }
3642#endif
3643            }
3644            n=upperList[-1];
3645            for (int i=0;i<n;i++) {
3646              int iSequence = upperList[i];
3647              double thisChange = change * upperChange[iSequence];
3648              double newValue = upper_[iSequence] + thisChange;
3649              upper_[iSequence] = newValue;
3650              if(getStatus(iSequence)==atUpperBound||
3651                 getStatus(iSequence)==isFixed) {
3652                solution_[iSequence] = newValue;
3653#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3654              } else if (getStatus(iSequence)==basic) {
3655                int iRow=backwardBasic[iSequence];
3656                upperGap[iRow] += thisChange;
3657#endif
3658              }
3659              // may have to adjust other bound
3660              double otherValue = lower_[iSequence];
3661              if (newValue-otherValue<dualBound_) {
3662                //originalBound(iSequence,useTheta,lowerChange,upperChange);
3663                //reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(iSequence);
3664                //ClpTraceDebug (fabs(upper_[iSequence]-newValue)<1.0e-5);
3665              }
3666            }
3667          }
3668          sequenceIn_=-1;
3669          if (pivotType) {
3670            if (useTheta>lastTheta+1.0e-9) {
3671              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3672                << useTheta << objectiveValue() << CoinMessageEol;
3673              lastTheta = useTheta;
3674            }
3675            problemStatus_ = -2;
3676            if (!factorization_->pivots()&&pivotRow_<0)
3677              problemStatus_=2;
3678#ifdef CLP_USER_DRIVEN
3679            {
3680              double saveTheta=theta_;
3681              theta_ = endingTheta;
3682              int status=eventHandler_->event(ClpEventHandler::theta);
3683              if (status>=0&&status<10) {
3684                endingTheta=theta_;
3685                theta_=saveTheta;
3686                problemStatus_=-1;
3687                continue;
3688              } else {
3689                if (status>=10)
3690                  problemStatus_=status-10;
3691                if (status<0)
3692                  startingTheta = useTheta;
3693                theta_=saveTheta;
3694              }
3695            }
3696#else
3697            startingTheta = useTheta;
3698#endif
3699            return 4;
3700          }
3701          // choose row to go out
3702          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3703          if (pivotRow_ >= 0) {
3704               // we found a pivot row
3705               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3706                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3707                              << pivotRow_
3708                              << CoinMessageEol;
3709               }
3710               // check accuracy of weights
3711               dualRowPivot_->checkAccuracy();
3712               // do ratio test for normal iteration
3713               double bestPossiblePivot = bestPivot();
3714               if (sequenceIn_ >= 0) {
3715                    // normal iteration
3716                    // update the incoming column
3717                    double btranAlpha = -alpha_ * directionOut_; // for check
3718#ifndef COIN_FAC_NEW
3719                    unpackPacked(rowArray_[1]);
3720#else
3721                    unpack(rowArray_[1]);
3722#endif
3723                    // and update dual weights (can do in parallel - with extra array)
3724                    rowArray_[2]->clear();
3725                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3726                                                          rowArray_[2],
3727                                                          rowArray_[3],
3728                                                          rowArray_[1]);
3729                    // see if update stable
3730#ifdef CLP_DEBUG
3731                    if ((handler_->logLevel() & 32))
3732                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3733#endif
3734                    double checkValue = 1.0e-7;
3735                    // if can't trust much and long way from optimal then relax
3736                    if (largestPrimalError_ > 10.0)
3737                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3738                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3739                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3740                         handler_->message(CLP_DUAL_CHECK, messages_)
3741                                   << btranAlpha
3742                                   << alpha_
3743                                   << CoinMessageEol;
3744                         // clear arrays
3745                         rowArray_[4]->clear();
3746                         if (factorization_->pivots()) {
3747                              dualRowPivot_->unrollWeights();
3748                              problemStatus_ = -2; // factorize now
3749                              rowArray_[0]->clear();
3750                              rowArray_[1]->clear();
3751                              columnArray_[0]->clear();
3752                              returnCode = -2;
3753                              break;
3754                         } else {
3755                              // take on more relaxed criterion
3756                              double test;
3757                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3758                                   test = 1.0e-1 * fabs(alpha_);
3759                              else
3760                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
3761                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3762                                        fabs(btranAlpha - alpha_) > test) {
3763                                   dualRowPivot_->unrollWeights();
3764                                   // need to reject something
3765                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
3766                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
3767                                             << x << sequenceWithin(sequenceOut_)
3768                                             << CoinMessageEol;
3769                                   setFlagged(sequenceOut_);
3770                                   progress_.clearBadTimes();
3771                                   lastBadIteration_ = numberIterations_; // say be more cautious
3772                                   rowArray_[0]->clear();
3773                                   rowArray_[1]->clear();
3774                                   columnArray_[0]->clear();
3775                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3776                                        //printf("I think should declare infeasible\n");
3777                                        problemStatus_ = 1;
3778                                        returnCode = 1;
3779                                        break;
3780                                   }
3781                                   continue;
3782                              }
3783                         }
3784                    }
3785                    // update duals BEFORE replaceColumn so can do updateColumn
3786                    double objectiveChange = 0.0;
3787                    // do duals first as variables may flip bounds
3788                    // rowArray_[0] and columnArray_[0] may have flips
3789                    // so use rowArray_[3] for work array from here on
3790                    int nswapped = 0;
3791                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3792                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3793#if CLP_CAN_HAVE_ZERO_OBJ
3794                    if ((specialOptions_&2097152)==0) {
3795#endif
3796                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3797                                                                                               rowArray_[2], theta_,
3798                                                                                               objectiveChange, false);
3799                      assert (!nswapped);
3800#if CLP_CAN_HAVE_ZERO_OBJ
3801                    } else {
3802                      rowArray_[0]->clear();
3803                      rowArray_[2]->clear();
3804                      columnArray_[0]->clear();
3805                    }
3806#endif
3807                    // which will change basic solution
3808                    if (nswapped) {
3809                      abort(); //needs testing
3810                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3811                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
3812                                                             1.0, objectiveChange);
3813                         // recompute dualOut_
3814                         valueOut_ = solution_[sequenceOut_];
3815                         if (directionOut_ < 0) {
3816                              dualOut_ = valueOut_ - upperOut_;
3817                         } else {
3818                              dualOut_ = lowerOut_ - valueOut_;
3819                         }
3820                    }
3821                    // amount primal will move
3822                    double movement = -dualOut_ * directionOut_ / alpha_;
3823                    // so objective should increase by fabs(dj)*movement
3824                    // but we already have objective change - so check will be good
3825                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3826#ifdef CLP_DEBUG
3827                         if (handler_->logLevel() & 32)
3828                              printf("movement %g, swap change %g, rest %g  * %g\n",
3829                                     objectiveChange + fabs(movement * dualIn_),
3830                                     objectiveChange, movement, dualIn_);
3831#endif
3832                         assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3833                         if(factorization_->pivots()) {
3834                              // going backwards - factorize
3835                              dualRowPivot_->unrollWeights();
3836                              problemStatus_ = -2; // factorize now
3837                              returnCode = -2;
3838                              break;
3839                         }
3840                    }
3841                    CoinAssert(fabs(dualOut_) < 1.0e50);
3842                    // if stable replace in basis
3843                    int updateStatus = factorization_->replaceColumn(this,
3844                                       rowArray_[2],
3845                                       rowArray_[1],
3846                                       pivotRow_,
3847                                       alpha_);
3848                    // if no pivots, bad update but reasonable alpha - take and invert
3849                    if (updateStatus == 2 &&
3850                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3851                         updateStatus = 4;
3852                    if (updateStatus == 1 || updateStatus == 4) {
3853                         // slight error
3854                         if (factorization_->pivots() > 5 || updateStatus == 4) {
3855                              problemStatus_ = -2; // factorize now
3856                              returnCode = -3;
3857                         }
3858                    } else if (updateStatus == 2) {
3859                         // major error
3860                         dualRowPivot_->unrollWeights();
3861                         // later we may need to unwind more e.g. fake bounds
3862                         if (factorization_->pivots()) {
3863                              problemStatus_ = -2; // factorize now
3864                              returnCode = -2;
3865                              break;
3866                         } else {
3867                              // need to reject something
3868                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3869                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3870                                        << x << sequenceWithin(sequenceOut_)
3871                                        << CoinMessageEol;
3872                              setFlagged(sequenceOut_);
3873                              progress_.clearBadTimes();
3874                              lastBadIteration_ = numberIterations_; // say be more cautious
3875                              rowArray_[0]->clear();
3876                              rowArray_[1]->clear();
3877                              columnArray_[0]->clear();
3878                              // make sure dual feasible
3879                              // look at all rows and columns
3880                              double objectiveChange = 0.0;
3881                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3882                                        0.0, objectiveChange, true);
3883                              continue;
3884                         }
3885                    } else if (updateStatus == 3) {
3886                         // out of memory
3887                         // increase space if not many iterations
3888                         if (factorization_->pivots() <
3889                                   0.5 * factorization_->maximumPivots() &&
3890                                   factorization_->pivots() < 200)
3891                              factorization_->areaFactor(
3892                                   factorization_->areaFactor() * 1.1);
3893                         problemStatus_ = -2; // factorize now
3894                    } else if (updateStatus == 5) {
3895                         problemStatus_ = -2; // factorize now
3896                    }
3897                    int * lowerActive = paramData.lowerActive;
3898                    int * upperActive = paramData.upperActive;
3899                    // update change vector
3900                    {
3901                      double * work = rowArray_[1]->denseVector();
3902                      int number = rowArray_[1]->getNumElements();
3903                      int * which = rowArray_[1]->getIndices();
3904                      assert (!rowArray_[4]->packedMode());
3905#ifndef COIN_FAC_NEW
3906                      assert (rowArray_[1]->packedMode());
3907#else
3908                      assert (!rowArray_[1]->packedMode());
3909#endif
3910                      double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3911                      double multiplier = -pivotValue/alpha_;
3912                      double * array=rowArray_[4]->denseVector();
3913#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3914                      int lowerN=lowerActive[-1];
3915                      int upperN=upperActive[-1];
3916#endif
3917                      if (multiplier) {
3918                        for (int i = 0; i < number; i++) {
3919                          int iRow = which[i];
3920#ifndef COIN_FAC_NEW
3921                          double alpha=multiplier*work[i];
3922#else
3923                          double alpha=multiplier*work[iRow];
3924#endif
3925#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3926                          double alpha3 = alpha+array[iRow];
3927                          int iSequence = pivotVariable_[iRow];
3928                          double oldLower = lowerCoefficient[iRow];
3929                          double oldUpper = upperCoefficient[iRow];
3930                          if (lower_[iSequence]>-1.0e30) {
3931                            //lowerGap[iRow]=value-lower_[iSequence];
3932                            double alpha2 = alpha3 + lowerChange[iSequence];
3933                            if (alpha2>1.0e-8)  {
3934                              lowerCoefficient[iRow]=alpha2;
3935                              if (!oldLower)
3936                                lowerActive[lowerN++]=iRow;
3937                            } else {
3938                              if (oldLower)
3939                                lowerCoefficient[iRow]=COIN_DBL_MIN;
3940                            }
3941                          } else {
3942                            if (oldLower)
3943                              lowerCoefficient[iRow]=COIN_DBL_MIN;
3944                          }
3945                          if (upper_[iSequence]<1.0e30) {
3946                            //upperGap[iRow]=-(value-upper_[iSequence]);
3947                            double alpha2 = -(alpha3+upperChange[iSequence]);
3948                            if (alpha2>1.0e-8) {
3949                              upperCoefficient[iRow]=alpha2;
3950                              if (!oldUpper)
3951                                upperActive[upperN++]=iRow;
3952                            } else {
3953                              if (oldUpper)
3954                                upperCoefficient[iRow]=COIN_DBL_MIN;
3955                            }
3956                          } else {
3957                            if (oldUpper)
3958                              upperCoefficient[iRow]=COIN_DBL_MIN;
3959                          }
3960#endif
3961                          rowArray_[4]->quickAdd(iRow,alpha);
3962                        }
3963                      }
3964                      pivotValue = array[pivotRow_];
3965                      // we want pivot to be -multiplier
3966                      rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3967#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
3968                      assert (lowerN>=0&&lowerN<=numberRows_);
3969                      lowerActive[-1]=lowerN;
3970                      upperActive[-1]=upperN;
3971#endif
3972                    }
3973                    // update primal solution
3974#if CLP_CAN_HAVE_ZERO_OBJ
3975                    if ((specialOptions_&2097152)!=0) 
3976                      theta_=0.0;
3977#endif
3978                    if (theta_ < 0.0) {
3979#ifdef CLP_DEBUG
3980                         if (handler_->logLevel() & 32)
3981                              printf("negative theta %g\n", theta_);
3982#endif
3983                         theta_ = 0.0;
3984                    }
3985                    // do actual flips
3986                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
3987                    //rowArray_[1]->expand();
3988#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
3989                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
3990                                                        movement,
3991                                                        objectiveChange);
3992#else
3993                    // do by hand
3994                    {
3995                      double * work = rowArray_[1]->denseVector();
3996                      int number = rowArray_[1]->getNumElements();
3997                      int * which = rowArray_[1]->getIndices();
3998                      int i;
3999                      if (rowArray_[1]->packedMode()) {
4000                        for (i = 0; i < number; i++) {
4001                          int iRow = which[i];
4002                          int iSequence = pivotVariable_[iRow];
4003                          double value = solution_[iSequence];
4004                          double change = movement * work[i];
4005                          value -= change;
4006                          if (lower_[iSequence]>-1.0e30)
4007                            lowerGap[iRow]=value-lower_[iSequence];
4008                          if (upper_[iSequence]<1.0e30)
4009                            upperGap[iRow]=-(value-upper_[iSequence]);
4010                          solution_[iSequence] = value;
4011                          objectiveChange -= change * cost_[iSequence];
4012                          work[i] = 0.0;
4013                        }
4014                      } else {
4015                        for (i = 0; i < number; i++) {
4016                          int iRow = which[i];
4017                          int iSequence = pivotVariable_[iRow];
4018                          double value = solution_[iSequence];
4019                          double change = movement * work[iRow];
4020                          value -= change;
4021                          solution_[iSequence] = value;
4022                          objectiveChange -= change * cost_[iSequence];
4023                          work[iRow] = 0.0;
4024                        }
4025                      }
4026                      rowArray_[1]->setNumElements(0);
4027                    }
4028#endif
4029                    // modify dualout
4030                    dualOut_ /= alpha_;
4031                    dualOut_ *= -directionOut_;
4032                    //setStatus(sequenceIn_,basic);
4033                    dj_[sequenceIn_] = 0.0;
4034                    //double oldValue = valueIn_;
4035                    if (directionIn_ == -1) {
4036                         // as if from upper bound
4037                         valueIn_ = upperIn_ + dualOut_;
4038                    } else {
4039                         // as if from lower bound
4040                         valueIn_ = lowerIn_ + dualOut_;
4041                    }
4042                    objectiveChange = 0.0;
4043#if CLP_CAN_HAVE_ZERO_OBJ
4044                    if ((specialOptions_&2097152)==0) {
4045#endif
4046                      for (int i=0;i<numberTotal;i++)
4047                        objectiveChange += solution_[i]*cost_[i];
4048                      objectiveChange -= objectiveValue_;
4049#if CLP_CAN_HAVE_ZERO_OBJ
4050                    }
4051#endif
4052                    // outgoing
4053                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
4054                    lowerOut_=lower_[sequenceOut_];
4055                    upperOut_=upper_[sequenceOut_];
4056                    // set dj to zero unless values pass
4057                    if (directionOut_ > 0) {
4058                         valueOut_ = lowerOut_;
4059                         dj_[sequenceOut_] = theta_;
4060#if CLP_CAN_HAVE_ZERO_OBJ>1
4061#ifdef COIN_REUSE_RANDOM
4062                         if ((specialOptions_&2097152)!=0) {
4063                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
4064                         }
4065#endif
4066#endif
4067                    } else {
4068                         valueOut_ = upperOut_;
4069                         dj_[sequenceOut_] = -theta_;
4070#if CLP_CAN_HAVE_ZERO_OBJ>1
4071#ifdef COIN_REUSE_RANDOM
4072                         if ((specialOptions_&2097152)!=0) {
4073                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
4074                         }
4075#endif
4076#endif
4077                    }
4078                    solution_[sequenceOut_] = valueOut_;
4079                    int whatNext = housekeeping(objectiveChange);
4080                    reinterpret_cast<ClpSimplexDual *>(this)->originalBound(sequenceIn_);
4081                    assert (backwardBasic[sequenceOut_]==pivotRow_);
4082                    backwardBasic[sequenceOut_]=-1;
4083                    backwardBasic[sequenceIn_]=pivotRow_;
4084#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4085                    double value = solution_[sequenceIn_];
4086                    double alpha = rowArray_[4]->denseVector()[pivotRow_];
4087                    double oldLower = lowerCoefficient[pivotRow_];
4088                    double oldUpper = upperCoefficient[pivotRow_];
4089                    if (lower_[sequenceIn_]>-1.0e30) {
4090                      lowerGap[pivotRow_]=value-lower_[sequenceIn_];
4091                      double alpha2 = alpha + lowerChange[sequenceIn_];
4092                      if (alpha2>1.0e-8)  {
4093                        lowerCoefficient[pivotRow_]=alpha2;
4094                        if (!oldLower) {
4095                          int lowerN=lowerActive[-1];
4096                          assert (lowerN>=0&&lowerN<numberRows_);
4097                          lowerActive[lowerN]=pivotRow_;
4098                          lowerActive[-1]=lowerN+1;
4099                        }
4100                      } else {
4101                        if (oldLower)
4102                          lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4103                      }
4104                    } else {
4105                      if (oldLower)
4106                        lowerCoefficient[pivotRow_]=COIN_DBL_MIN;
4107                    }
4108                    if (upper_[sequenceIn_]<1.0e30) {
4109                      upperGap[pivotRow_]=-(value-upper_[sequenceIn_]);
4110                      double alpha2 = -(alpha+upperChange[sequenceIn_]);
4111                      if (alpha2>1.0e-8) {
4112                        upperCoefficient[pivotRow_]=alpha2;
4113                        if (!oldUpper) {
4114                          int upperN=upperActive[-1];
4115                          assert (upperN>=0&&upperN<numberRows_);
4116                          upperActive[upperN]=pivotRow_;
4117                          upperActive[-1]=upperN+1;
4118                        }
4119                      } else {
4120                        if (oldUpper)
4121                          upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4122                      }
4123                    } else {
4124                      if (oldUpper)
4125                        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
4126                    }
4127#endif
4128                    {
4129                      char in[200],out[200];
4130                      int iSequence=sequenceIn_;
4131                      if (iSequence<numberColumns_) {
4132                        if (lengthNames_) 
4133                          strcpy(in,columnNames_[iSequence].c_str());
4134                         else 
4135                          sprintf(in,"C%7.7d",iSequence);
4136                      } else {
4137                        iSequence -= numberColumns_;
4138                        if (lengthNames_) 
4139                          strcpy(in,rowNames_[iSequence].c_str());
4140                         else 
4141                          sprintf(in,"R%7.7d",iSequence);
4142                      }
4143                      iSequence=sequenceOut_;
4144                      if (iSequence<numberColumns_) {
4145                        if (lengthNames_) 
4146                          strcpy(out,columnNames_[iSequence].c_str());
4147                         else 
4148                          sprintf(out,"C%7.7d",iSequence);
4149                      } else {
4150                        iSequence -= numberColumns_;
4151                        if (lengthNames_) 
4152                          strcpy(out,rowNames_[iSequence].c_str());
4153                         else 
4154                          sprintf(out,"R%7.7d",iSequence);
4155                      }
4156                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
4157                        << useTheta << objectiveValue() 
4158                        << in << out << CoinMessageEol;
4159                    }
4160                    if (useTheta>lastTheta+1.0e-9) {
4161                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
4162                        << useTheta << objectiveValue() << CoinMessageEol;
4163                      lastTheta = useTheta;
4164                    }
4165                    // and set bounds correctly
4166                    originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
4167                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
4168                    if (whatNext == 1) {
4169                         problemStatus_ = -2; // refactorize
4170                    } else if (whatNext == 2) {
4171                         // maximum iterations or equivalent
4172                         problemStatus_ = 3;
4173                         returnCode = 3;
4174                         break;
4175                    }
4176#ifdef CLP_USER_DRIVEN
4177                    // Check event
4178                    {
4179                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
4180                         if (status >= 0) {
4181                              problemStatus_ = 5;
4182                              secondaryStatus_ = ClpEventHandler::endOfIteration;
4183                              returnCode = 4;
4184                              break;
4185                         }
4186                    }
4187#endif
4188               } else {
4189                    // no incoming column is valid
4190#ifdef CLP_USER_DRIVEN
4191                 rowArray_[0]->clear();
4192                 columnArray_[0]->clear();
4193                 theta_ = useTheta;
4194                 lastTheta = useTheta;
4195                 int action = eventHandler_->event(ClpEventHandler::noTheta);
4196                 if (action>=0) {
4197                   endingTheta=theta_;
4198                   theta_ = 0.0;
4199                   //adjust [4] from handler - but
4200                   //rowArray_[4]->clear(); // temp
4201                   if (action>=0&&action<10)
4202                     problemStatus_=-1; // carry on
4203                   else if (action==15)
4204                     problemStatus_ =5; // say stopped
4205                   returnCode = 1;
4206                   if (action==0||action>=10) 
4207                     break;
4208                   else
4209                     continue;
4210                 } else {
4211                 theta_ = 0.0;
4212                 }
4213#endif
4214                    pivotRow_ = -1;
4215#ifdef CLP_DEBUG
4216                    if (handler_->logLevel() & 32)
4217                         printf("** no column pivot\n");
4218#endif
4219                    if (factorization_->pivots() < 10) { 
4220                         // If we have just factorized and infeasibility reasonable say infeas
4221                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
4222                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
4223                                        || (specialOptions_ & 64) == 0) {
4224                                   // say infeasible
4225                                   problemStatus_ = 1;
4226                                   // unless primal feasible!!!!
4227                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
4228                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
4229                                   if (numberDualInfeasibilities_)
4230                                        problemStatus_ = 10;
4231                                   rowArray_[0]->clear();
4232                                   columnArray_[0]->clear();
4233                              }
4234                         }
4235                         // If special option set - put off as long as possible
4236                         if ((specialOptions_ & 64) == 0) {
4237                              problemStatus_ = -4; //say looks infeasible
4238                         } else {
4239                              // flag
4240                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
4241                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
4242                                        << x << sequenceWithin(sequenceOut_)
4243                                        << CoinMessageEol;
4244                              setFlagged(sequenceOut_);
4245                              if (!factorization_->pivots()) {
4246                                   rowArray_[0]->clear();
4247                                   columnArray_[0]->clear();
4248                                   continue;
4249                              }
4250                         }
4251                    }
4252                    rowArray_[0]->clear();
4253                    columnArray_[0]->clear();
4254                    returnCode = 1;
4255                    break;
4256               }
4257          } else {
4258               // no pivot row
4259#ifdef CLP_USER_DRIVEN
4260            {
4261              double saveTheta=theta_;
4262              theta_ = endingTheta;
4263              int status=eventHandler_->event(ClpEventHandler::theta);
4264              if (status>=0&&status<10) {
4265                endingTheta=theta_;
4266                theta_=saveTheta;
4267                continue;
4268              } else {
4269                theta_=saveTheta;
4270              }
4271            }
4272#endif
4273#ifdef CLP_DEBUG
4274               if (handler_->logLevel() & 32)
4275                    printf("** no row pivot\n");
4276#endif
4277               int numberPivots = factorization_->pivots();
4278               bool specialCase;
4279               int useNumberFake;
4280               returnCode = 0;
4281               if (numberPivots < 20 &&
4282                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
4283                         && dualBound_ > 1.0e8) {
4284                    specialCase = true;
4285                    // as dual bound high - should be okay
4286                    useNumberFake = 0;
4287               } else {
4288                    specialCase = false;
4289                    useNumberFake = numberFake_;
4290               }
4291               if (!numberPivots || specialCase) {
4292                    // may have crept through - so may be optimal
4293                    // check any flagged variables
4294                    int iRow;
4295                    for (iRow = 0; iRow < numberRows_; iRow++) {
4296                         int iPivot = pivotVariable_[iRow];
4297                         if (flagged(iPivot))
4298                              break;
4299                    }
4300                    if (iRow < numberRows_ && numberPivots) {
4301                         // try factorization
4302                         returnCode = -2;
4303                    }
4304
4305                    if (useNumberFake || numberDualInfeasibilities_) {
4306                         // may be dual infeasible
4307                         problemStatus_ = -5;
4308                    } else {
4309                         if (iRow < numberRows_) {
4310                              problemStatus_ = -5;
4311                         } else {
4312                              if (numberPivots) {
4313                                   // objective may be wrong
4314                                   objectiveValue_ = innerProduct(cost_,
4315                                                                  numberColumns_ + numberRows_,
4316                                                                  solution_);
4317                                   objectiveValue_ += objective_->nonlinearOffset();
4318                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
4319                                   if ((specialOptions_ & 16384) == 0) {
4320                                        // and dual_ may be wrong (i.e. for fixed or basic)
4321                                        CoinIndexedVector * arrayVector = rowArray_[1];
4322                                        arrayVector->clear();
4323                                        int iRow;
4324                                        double * array = arrayVector->denseVector();
4325                                        /* Use dual_ instead of array
4326                                           Even though dual_ is only numberRows_ long this is
4327                                           okay as gets permuted to longer rowArray_[2]
4328                                        */
4329                                        arrayVector->setDenseVector(dual_);
4330                                        int * index = arrayVector->getIndices();
4331                                        int number = 0;
4332                                        for (iRow = 0; iRow < numberRows_; iRow++) {
4333                                             int iPivot = pivotVariable_[iRow];
4334                                             double value = cost_[iPivot];
4335                                             dual_[iRow] = value;
4336                                             if (value) {
4337                                                  index[number++] = iRow;
4338                                             }
4339                                        }
4340                                        arrayVector->setNumElements(number);
4341                                        // Extended duals before "updateTranspose"
4342                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
4343                                        // Btran basic costs
4344                                        rowArray_[2]->clear();
4345                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
4346                                        // and return vector
4347                                        arrayVector->setDenseVector(array);
4348                                   }
4349                              }
4350                              problemStatus_ = 0;
4351                              sumPrimalInfeasibilities_ = 0.0;
4352                              if ((specialOptions_&(1024 + 16384)) != 0) {
4353                                   CoinIndexedVector * arrayVector = rowArray_[1];
4354                                   arrayVector->clear();
4355                                   double * rhs = arrayVector->denseVector();
4356                                   times(1.0, solution_, rhs);
4357                                   bool bad2 = false;
4358                                   int i;
4359                                   for ( i = 0; i < numberRows_; i++) {
4360                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
4361                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
4362                                             bad2 = true;
4363                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
4364                                        }
4365                                        rhs[i] = 0.0;
4366                                   }
4367                                   for ( i = 0; i < numberColumns_; i++) {
4368                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
4369                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4370                                             bad2 = true;
4371                                        }
4372                                   }
4373                                   if (bad2) {
4374                                        problemStatus_ = -3;
4375                                        returnCode = -2;
4376                                        // Force to re-factorize early next time
4377                                        int numberPivots = factorization_->pivots();
4378                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4379                                   }
4380                              }
4381                         }
4382                    }
4383               } else {
4384                    problemStatus_ = -3;
4385                    returnCode = -2;
4386                    // Force to re-factorize early next time
4387                    int numberPivots = factorization_->pivots();
4388                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4389               }
4390               break;
4391          }
4392     }
4393     startingTheta = lastTheta+theta_;
4394     return returnCode;
4395}
4396#if 0
4397static int zzzzzz=0;
4398int zzzzzzOther=0;
4399#endif
4400// Finds best possible pivot
4401double 
4402ClpSimplexOther::bestPivot(bool justColumns)
4403{
4404  // Get good size for pivot
4405  // Allow first few iterations to take tiny
4406  double acceptablePivot = 1.0e-9;
4407  if (numberIterations_ > 100)
4408    acceptablePivot = 1.0e-8;
4409  if (factorization_->pivots() > 10 ||
4410      (factorization_->pivots() && sumDualInfeasibilities_))
4411    acceptablePivot = 1.0e-5; // if we have iterated be more strict
4412  else if (factorization_->pivots() > 5)
4413    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4414  else if (factorization_->pivots())
4415    acceptablePivot = 1.0e-8; // relax
4416  double bestPossiblePivot = 1.0;
4417  // get sign for finding row of tableau
4418  // normal iteration
4419  // create as packed
4420  double direction = directionOut_;
4421#ifndef COIN_FAC_NEW
4422  rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4423#else
4424  rowArray_[0]->createOneUnpackedElement(pivotRow_, direction);
4425#endif
4426  factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4427  // put row of tableau in rowArray[0] and columnArray[0]
4428  matrix_->transposeTimes(this, -1.0,
4429                          rowArray_[0], rowArray_[3], columnArray_[0]);
4430  sequenceIn_=-1;
4431  if (justColumns)
4432    rowArray_[0]->clear();
4433  // do ratio test for normal iteration
4434  bestPossiblePivot = 
4435    reinterpret_cast<ClpSimplexDual *> 
4436    ( this)->dualColumn(rowArray_[0],
4437                        columnArray_[0], columnArray_[1],
4438                        rowArray_[3], acceptablePivot, NULL);
4439  return bestPossiblePivot;
4440}
4441// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4442int
4443ClpSimplexOther::nextTheta(int /*type*/, double maxTheta, parametricsData & paramData,
4444                           const double * /*changeObjective*/)
4445{
4446  const double * lowerChange = paramData.lowerChange;
4447  const double * upperChange = paramData.upperChange;
4448  const int * lowerList = paramData.lowerList;
4449  const int * upperList = paramData.upperList;
4450  int iSequence;
4451  bool toLower = false;
4452  //assert (type==1);
4453  // may need to decide based on model?
4454  bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4455  double * array = rowArray_[4]->denseVector();
4456  //rowArray_[4]->checkClean();
4457  const int * row = matrix_->getIndices();
4458  const int * columnLength = matrix_->getVectorLengths();
4459  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4460  const double * elementByColumn = matrix_->getElements();
4461#if 0
4462  double tempArray[5000];
4463  bool checkIt=false;
4464  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
4465    memcpy(tempArray,array,numberRows_*sizeof(double));
4466    checkIt=true;
4467    needFullUpdate=true;
4468  }
4469#endif
4470#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4471  double * lowerGap = paramData.lowerGap;
4472  double * upperGap = paramData.upperGap;
4473  double * lowerCoefficient = paramData.lowerCoefficient;
4474  double * upperCoefficient = paramData.upperCoefficient;
4475  int * lowerActive=paramData.lowerActive;
4476  int * upperActive=paramData.upperActive;
4477#endif
4478  if (!factorization_->pivots()||needFullUpdate) {
4479    //zzzzzz=0;
4480    rowArray_[4]->clear();
4481    // get change
4482    if (!rowScale_) {
4483      int n;
4484      n=lowerList[-2];
4485      int i;
4486      for (i=0;i<n;i++) {
4487        int iSequence = lowerList[i];
4488        assert (iSequence<numberColumns_);
4489        if (getColumnStatus(iSequence)==atLowerBound) {
4490          double value=lowerChange[iSequence];
4491          for (CoinBigIndex j = columnStart[iSequence];
4492               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4493            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4494          }
4495        }
4496      }
4497      n=lowerList[-1];
4498      const double * change = lowerChange+numberColumns_;
4499      for (;i<n;i++) {
4500        int iSequence = lowerList[i]-numberColumns_;
4501        assert (iSequence>=0);
4502        if (getRowStatus(iSequence)==atLowerBound) {
4503          double value=change[iSequence];
4504          rowArray_[4]->quickAdd(iSequence, -value);
4505        }
4506      }
4507      n=upperList[-2];
4508      for (i=0;i<n;i++) {
4509        int iSequence = upperList[i];
4510        assert (iSequence<numberColumns_);
4511        if (getColumnStatus(iSequence)==atUpperBound) {
4512          double value=upperChange[iSequence];
4513          for (CoinBigIndex j = columnStart[iSequence];
4514               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4515            rowArray_[4]->quickAdd(row[j], elementByColumn[j]*value);
4516          }
4517        }
4518      }
4519      n=upperList[-1];
4520      change = upperChange+numberColumns_;
4521      for (;i<n;i++) {
4522        int iSequence = upperList[i]-numberColumns_;
4523        assert (iSequence>=0);
4524        if (getRowStatus(iSequence)==atUpperBound) {
4525          double value=change[iSequence];
4526          rowArray_[4]->quickAdd(iSequence, -value);
4527        }
4528      }
4529    } else {
4530      int n;
4531      n=lowerList[-2];
4532      int i;
4533      for (i=0;i<n;i++) {
4534        int iSequence = lowerList[i];
4535        assert (iSequence<numberColumns_);
4536        if (getColumnStatus(iSequence)==atLowerBound) {
4537          double value=lowerChange[iSequence];
4538          // apply scaling
4539          double scale = columnScale_[iSequence];
4540          for (CoinBigIndex j = columnStart[iSequence];
4541               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4542            int iRow = row[j];
4543            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4544          }
4545        }
4546      }
4547      n=lowerList[-1];
4548      const double * change = lowerChange+numberColumns_;
4549      for (;i<n;i++) {
4550        int iSequence = lowerList[i]-numberColumns_;
4551        assert (iSequence>=0);
4552        if (getRowStatus(iSequence)==atLowerBound) {
4553          double value=change[iSequence];
4554          rowArray_[4]->quickAdd(iSequence, -value);
4555        }
4556      }
4557      n=upperList[-2];
4558      for (i=0;i<n;i++) {
4559        int iSequence = upperList[i];
4560        assert (iSequence<numberColumns_);
4561        if (getColumnStatus(iSequence)==atUpperBound) {
4562          double value=upperChange[iSequence];
4563          // apply scaling
4564          double scale = columnScale_[iSequence];
4565          for (CoinBigIndex j = columnStart[iSequence];
4566               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4567            int iRow = row[j];
4568            rowArray_[4]->quickAdd(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4569          }
4570        }
4571      }
4572      n=upperList[-1];
4573      change = upperChange+numberColumns_;
4574      for (;i<n;i++) {
4575        int iSequence = upperList[i]-numberColumns_;
4576        assert (iSequence>=0);
4577        if (getRowStatus(iSequence)==atUpperBound) {
4578          double value=change[iSequence];
4579          rowArray_[4]->quickAdd(iSequence, -value);
4580        }
4581      }
4582    }
4583    // ftran it
4584    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4585#if 0
4586    if (checkIt) {
4587      for (int i=0;i<numberRows_;i++) {
4588        assert (fabs(tempArray[i]-array[i])<1.0e-8);
4589      }
4590    }
4591#endif
4592#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4593    /* later for sparse - keep like CoinIndexedvector
4594       and just redo here */
4595    int lowerN=0;
4596    int upperN=0;
4597    memset(lowerCoefficient,0,numberRows_*sizeof(double));
4598    memset(upperCoefficient,0,numberRows_*sizeof(double));
4599    for (int iRow=0;iRow<numberRows_;iRow++) {
4600      iSequence = pivotVariable_[iRow];
4601      double currentSolution = solution_[iSequence];
4602      double alpha = array[iRow];
4603      double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4604      double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4605      if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4606        double currentLower = lower_[iSequence];
4607        ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4608        double gap=currentSolution-currentLower;
4609        lowerGap[iRow]=gap;
4610        lowerCoefficient[iRow]=thetaCoefficientLower;
4611        lowerActive[lowerN++]=iRow;
4612        //} else {
4613        //lowerCoefficient[iRow]=0.0;
4614      }
4615      if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4616        double currentUpper = upper_[iSequence];
4617        ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4618        double gap2=-(currentSolution-currentUpper); //positive
4619        upperGap[iRow]=gap2;
4620        upperCoefficient[iRow]=-thetaCoefficientUpper;
4621        upperActive[upperN++]=iRow;
4622      }
4623    }
4624    assert (lowerN>=0&&lowerN<=numberRows_);
4625    lowerActive[-1]=lowerN;
4626    upperActive[-1]=upperN;
4627#endif
4628  } else if (sequenceIn_>=0) {
4629    //assert (sequenceIn_>=0);
4630    assert (sequenceOut_>=0);
4631    assert (sequenceIn_!=sequenceOut_);
4632    double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4633    int needed=0;
4634    assert (!rowArray_[5]->getNumElements());
4635    if (change) {
4636      if (sequenceIn_<numberColumns_) {
4637        if (!rowScale_) {
4638          for (CoinBigIndex i = columnStart[sequenceIn_];
4639               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4640            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4641          }
4642        } else {
4643          // apply scaling
4644          double scale = columnScale_[sequenceIn_];
4645          for (CoinBigIndex i = columnStart[sequenceIn_];
4646               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4647            int iRow = row[i];
4648            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4649          }
4650        }
4651      } else {
4652        rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4653      }
4654      needed++;
4655    }
4656    if (getStatus(sequenceOut_)==atLowerBound)
4657      change=lowerChange[sequenceOut_];
4658    else
4659      change=upperChange[sequenceOut_];
4660    if (change) {
4661      if (sequenceOut_<numberColumns_) {
4662        if (!rowScale_) {
4663          for (CoinBigIndex i = columnStart[sequenceOut_];
4664               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4665            rowArray_[5]->quickAdd(row[i], elementByColumn[i]*change);
4666          }
4667        } else {
4668          // apply scaling
4669          double scale = columnScale_[sequenceOut_];
4670          for (CoinBigIndex i = columnStart[sequenceOut_];
4671               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4672            int iRow = row[i];
4673            rowArray_[5]->quickAdd(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4674          }
4675        }
4676      } else {
4677        rowArray_[5]->quickAdd(sequenceOut_-numberColumns_,-change);
4678      }
4679      needed++;
4680    }
4681    //printf("seqin %d seqout %d needed %d\n",
4682    //     sequenceIn_,sequenceOut_,needed);
4683    if (needed) {
4684      // ftran it
4685      factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4686      // add
4687      double * array5 = rowArray_[5]->denseVector();
4688      int * index5 = rowArray_[5]->getIndices();
4689      int number5 = rowArray_[5]->getNumElements();
4690#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4691      int lowerN=lowerActive[-1];
4692      int upperN=upperActive[-1];
4693      int nIn4=rowArray_[4]->getNumElements();
4694      int * index4 = rowArray_[4]->getIndices();
4695#endif
4696      for (int i = 0; i < number5; i++) {
4697        int iPivot = index5[i];
4698#ifndef CLP_PARAMETRIC_DENSE_ARRAYS
4699        rowArray_[4]->quickAdd(iPivot,array5[iPivot]);
4700#else
4701        /* later for sparse - modify here */
4702        int iSequence = pivotVariable_[iPivot];
4703        double currentSolution = solution_[iSequence];
4704        double currentAlpha = array[iPivot];
4705        double alpha5 = array5[iPivot];
4706        double alpha = currentAlpha+alpha5;
4707        if (currentAlpha) {
4708          if (alpha) {
4709            array[iPivot] = alpha;
4710          } else {
4711            array[iPivot] = COIN_DBL_MIN;
4712          }
4713        } else {
4714          index4[nIn4++] = iPivot;
4715          array[iPivot] = alpha;
4716        }
4717        double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4718        double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4719        double oldLower = lowerCoefficient[iPivot];
4720        double oldUpper = upperCoefficient[iPivot];
4721        if (thetaCoefficientLower > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4722          double currentLower = lower_[iSequence];
4723          ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4724          double gap=currentSolution-currentLower;
4725          lowerGap[iPivot]=gap;
4726          lowerCoefficient[iPivot]=thetaCoefficientLower;
4727          if (!oldLower)
4728            lowerActive[lowerN++]=iPivot;
4729        } else {
4730          if (oldLower)
4731            lowerCoefficient[iPivot]=COIN_DBL_MIN;
4732        }
4733        if (thetaCoefficientUpper < -1.0e-8&&upper_[iSequence]<1.0e30) {
4734          double currentUpper = upper_[iSequence];
4735          ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4736          double gap2=-(currentSolution-currentUpper); //positive
4737          upperGap[iPivot]=gap2;
4738          upperCoefficient[iPivot]=-thetaCoefficientUpper;
4739          if (!oldUpper)
4740            upperActive[upperN++]=iPivot;
4741        } else {
4742          if (oldUpper)
4743            upperCoefficient[iPivot]=COIN_DBL_MIN;
4744        }
4745#endif
4746        array5[iPivot]=0.0;
4747      }
4748      rowArray_[5]->setNumElements(0);
4749#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4750      rowArray_[4]->setNumElements(nIn4);
4751      assert (lowerN>=0&&lowerN<=numberRows_);
4752      lowerActive[-1]=lowerN;
4753      upperActive[-1]=upperN;
4754#endif
4755    }
4756  }
4757  const int * index = rowArray_[4]->getIndices();
4758  int number = rowArray_[4]->getNumElements();
4759#define TESTXX 0
4760#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4761  int * markDone = reinterpret_cast<int *>(paramData.markDone);
4762  int nToZero=(numberRows_+numberColumns_+COIN_ANY_BITS_PER_INT-1)>>COIN_ANY_SHIFT_PER_INT;
4763  memset(markDone,0,nToZero*sizeof(int));
4764  const int * backwardBasic = paramData.backwardBasic;
4765#endif
4766  // first ones with alpha
4767  double theta1=maxTheta;
4768  int pivotRow1=-1;
4769#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4770  int pivotRow2=-1;
4771  double theta2=maxTheta;
4772#endif
4773#ifndef CLP_PARAMETRIC_DENSE_ARRAYS //TESTXX
4774  for (int i=0;i<number;i++) {
4775    int iPivot=index[i];
4776    iSequence = pivotVariable_[iPivot];
4777    //assert(!markDone[iSequence]);
4778    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4779    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4780    markDone[word] |= ( 1 << bit );
4781    // solution value will be sol - theta*alpha
4782    // bounds will be bounds + change *theta
4783    double currentSolution = solution_[iSequence];
4784    double alpha = array[iPivot];
4785    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4786    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4787    if (thetaCoefficientLower > 1.0e-8) {
4788      double currentLower = lower_[iSequence];
4789      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4790      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4791      double gap=currentSolution-currentLower;
4792      if (thetaCoefficientLower*theta1>gap) {
4793        theta1 = gap/thetaCoefficientLower;
4794        //toLower=true;
4795        pivotRow1=iPivot;
4796      }
4797    }
4798    if (thetaCoefficientUpper < -1.0e-8) {
4799      double currentUpper = upper_[iSequence];
4800      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4801      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4802      double gap2=currentSolution-currentUpper; //negative
4803      if (thetaCoefficientUpper*theta2<gap2) {
4804        theta2 = gap2/thetaCoefficientUpper;
4805        //toLower=false;
4806        pivotRow2=iPivot;
4807      }
4808    }
4809  }
4810  // now others
4811  int nLook=lowerList[-1];
4812  for (int i=0;i<nLook;i++) {
4813    int iSequence = lowerList[i];
4814    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4815    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4816    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4817      double currentSolution = solution_[iSequence];
4818      double currentLower = lower_[iSequence];
4819      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4820      double thetaCoefficient = lowerChange[iSequence];
4821      if (thetaCoefficient > 0.0) {
4822        double gap=currentSolution-currentLower;
4823        if (thetaCoefficient*theta1>gap) {
4824          theta1 = gap/thetaCoefficient;
4825          //toLower=true;
4826          pivotRow1 = backwardBasic[iSequence];
4827        }
4828      }
4829    }
4830  }
4831  nLook=upperList[-1];
4832  for (int i=0;i<nLook;i++) {
4833    int iSequence = upperList[i];
4834    int word = iSequence >> COIN_ANY_SHIFT_PER_INT;
4835    int bit = iSequence & COIN_ANY_MASK_PER_INT;
4836    if (getColumnStatus(iSequence)==basic&&(markDone[word]&(1<<bit))==0) {
4837      double currentSolution = solution_[iSequence];
4838      double currentUpper = upper_[iSequence];
4839      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4840      double thetaCoefficient = upperChange[iSequence];
4841      if (thetaCoefficient < 0) {
4842        double gap=currentSolution-currentUpper; //negative
4843        if (thetaCoefficient*theta2<gap) {
4844          theta2 = gap/thetaCoefficient;
4845          //toLower=false;
4846          pivotRow2 = backwardBasic[iSequence];
4847        }
4848      }
4849    }
4850  }
4851  if (theta2<theta1) {
4852    theta_=theta2;
4853    toLower=false;
4854    pivotRow_=pivotRow2;
4855  } else {
4856    theta_=theta1;
4857    toLower=true;
4858    pivotRow_=pivotRow1;
4859  }
4860#if 0 //TESTXX
4861#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
4862  {
4863    double * checkArray = new double[numberRows_];
4864    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4865    int lowerN=lowerActive[-1];
4866    for (int i=0;i<lowerN;i++) {
4867      int iRow=lowerActive[i];
4868      int iSequence = pivotVariable_[iRow];
4869      double alpha = array[iRow];
4870      double thetaCoefficient = lowerChange[iSequence] + alpha;
4871      if (thetaCoefficient > 1.0e-8&&lower_[iSequence]>-1.0e30) {
4872        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4873        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4874          abort();
4875        }
4876      } else {
4877        assert (fabs(checkArray[iRow])<1.0e-12);
4878        if (fabs(checkArray[iRow])>1.0e-12) {
4879          abort();
4880        }
4881      }
4882      checkArray[iRow]=0.0;
4883    }
4884    for (int i=0;i<numberRows_;i++) {
4885      assert (!checkArray[i]);
4886      if (checkArray[i])
4887        abort();
4888    }
4889    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
4890    int upperN=upperActive[-1];
4891    for (int i=0;i<upperN;i++) {
4892      int iRow=upperActive[i];
4893      int iSequence = pivotVariable_[iRow];
4894      double alpha = array[iRow];
4895      double thetaCoefficient = -(upperChange[iSequence] + alpha);
4896      if (thetaCoefficient > 1.0e-8&&upper_[iSequence]<1.0e30) {
4897        assert(fabs(checkArray[iRow]-thetaCoefficient)<1.0e-5);
4898        if(fabs(checkArray[iRow]-thetaCoefficient)>1.0e-5) {
4899          abort();
4900        }
4901      } else {
4902        assert (fabs(checkArray[iRow])<1.0e-12);
4903        if (fabs(checkArray[iRow])>1.0e-12) {
4904          abort();
4905        }
4906      }
4907      checkArray[iRow]=0.0;
4908    }
4909    for (int i=0;i<numberRows_;i++) {
4910      assert (!checkArray[i]);
4911      if (checkArray[i])
4912        abort();
4913    }
4914    delete [] checkArray;
4915  }
4916  double theta3=maxTheta;
4917  int pivotRow3=-1;
4918  int lowerN=lowerActive[-1];
4919  for (int i=0;i<lowerN;i++) {
4920    int iRow=lowerActive[i];
4921    double lowerC = lowerCoefficient[iRow];
4922    double gap=lowerGap[iRow];
4923    if (toLower&&iRow==pivotRow_) {
4924      assert (lowerC*theta3>gap-1.0e-8);
4925      if (lowerC*theta3<gap-1.0e-8)
4926        abort();
4927    }
4928    if (lowerC*theta3>gap&&lowerC!=COIN_DBL_MIN) {
4929      theta3 = gap/lowerC;
4930      pivotRow3=iRow;
4931    }
4932  }
4933  int pivotRow4=pivotRow3;
4934  double theta4=theta3;
4935  int upperN=upperActive[-1];
4936  for (int i=0;i<upperN;i++) {
4937    int iRow=upperActive[i];
4938    double upperC = upperCoefficient[iRow];
4939    double gap=upperGap[iRow];
4940    if (!toLower&&iRow==pivotRow_) {
4941      assert (upperC*theta3>gap-1.0e-8);
4942      if (upperC*theta3<gap-1.0e-8)
4943        abort();
4944    }
4945    if (upperC*theta4>gap&&upperC!=COIN_DBL_MIN) {
4946      theta4 = gap/upperC;
4947      pivotRow4=iRow;
4948    }
4949  }
4950  bool toLower3;
4951  if (theta4<theta3) {
4952    theta3=theta4;
4953    toLower3=false;
4954    pivotRow3=pivotRow4;
4955  } else {
4956    toLower3=true;
4957  }
4958  if (fabs(theta3-theta_)>1.0e-8)
4959    abort();
4960  if (toLower!=toLower3||pivotRow_!=pivotRow3) {
4961    printf("bad piv - good %d %g %s, bad %d %g %s\n",pivotRow_,theta_,toLower ? "toLower" : "toUpper",
4962           pivotRow3,theta3,toLower3 ? "toLower" : "toUpper");
4963    //zzzzzz++;
4964    if (true/*zzzzzz>zzzzzzOther*/) {
4965      printf("Swapping\n");
4966      pivotRow_=pivotRow3;
4967      theta_=theta3;
4968      toLower=toLower3;
4969    }
4970  }
4971#endif
4972#endif
4973#else
4974#if 0 //CLP_PARAMETRIC_DENSE_ARRAYS==2
4975  {
4976    double * checkArray = new double[numberRows_];
4977    memcpy(checkArray,lowerCoefficient,numberRows_*sizeof(double));
4978    int lowerN=lowerActive[-1];
4979    for (int i=0;i<lowerN;i++) {
4980      int iRow=lowerActive[i];
4981      checkArray[iRow]=0.0;
4982    }
4983    for (int i=0;i<numberRows_;i++) {
4984      assert (!checkArray[i]);
4985      if (checkArray[i])
4986        abort();
4987    }
4988    memcpy(checkArray,upperCoefficient,numberRows_*sizeof(double));
4989    int upperN=upperActive[-1];
4990    for (int i=0;i<upperN;i++) {
4991      int iRow=upperActive[i];
4992      checkArray[iRow]=0.0;
4993    }
4994    for (int i=0;i<numberRows_;i++) {
4995      assert (!checkArray[i]);
4996      if (checkArray[i])
4997        abort();
4998    }
4999    delete [] checkArray;
5000  }
5001#endif
5002  int lowerN=lowerActive[-1];
5003  for (int i=0;i<lowerN;i++) {
5004    int iRow=lowerActive[i];
5005    double lowerC = lowerCoefficient[iRow];
5006    double gap=lowerGap[iRow];
5007    if (lowerC*theta1>gap&&lowerC!=COIN_DBL_MIN) {
5008      theta1 = gap/lowerC;
5009      pivotRow1=iRow;
5010    }
5011  }
5012  pivotRow_=pivotRow1;
5013  theta_=theta1;
5014  int upperN=upperActive[-1];
5015  for (int i=0;i<upperN;i++) {
5016    int iRow=upperActive[i];
5017    double upperC = upperCoefficient[iRow];
5018    double gap=upperGap[iRow];
5019    if (upperC*theta1>gap&&upperC!=COIN_DBL_MIN) {
5020      theta1 = gap/upperC;
5021      pivotRow1=iRow;
5022    }
5023  }
5024  if (theta1<theta_) {
5025    theta_=theta1;
5026    toLower=false;
5027    pivotRow_=pivotRow1;
5028  } else {
5029    toLower=true;
5030  }
5031#endif
5032  theta_ = CoinMax(theta_,0.0);
5033  if (theta_>1.0e-15) {
5034    // update solution
5035    for (int iRow = 0; iRow < number; iRow++) {
5036      int iPivot = index[iRow];
5037      iSequence = pivotVariable_[iPivot];
5038      // solution value will be sol - theta*alpha
5039      double alpha = array[iPivot];
5040      double currentSolution = solution_[iSequence] - theta_ * alpha;
5041      solution_[iSequence] =currentSolution;
5042#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5043      if (lower_[iSequence]>-1.0e30)
5044        lowerGap[iPivot]=currentSolution-lower_[iSequence];
5045      if (upper_[iSequence]<1.0e30)
5046        upperGap[iPivot]=-(currentSolution-upper_[iSequence]);
5047#endif
5048    }
5049  }
5050#ifdef CLP_PARAMETRIC_DENSE_ARRAYS
5051  if (pivotRow_>=0&&false) {
5052    double oldValue = upperCoefficient[pivotRow_];
5053    double value = array[pivotRow_];
5054    if (value) {
5055      if (!oldValue) {
5056        int upperN=upperActive[-1];
5057        assert (upperN>=0&&upperN<numberRows_);
5058        upperActive[upperN]=pivotRow_;
5059        upperActive[-1]=upperN+1;
5060      }
5061    } else {
5062      if (oldValue)
5063        upperCoefficient[pivotRow_]=COIN_DBL_MIN;
5064    }
5065  }
5066#endif
5067#if 0
5068  for (int i=0;i<numberTotal;i++)
5069    assert(!markDone[i]);
5070#endif
5071  if (pivotRow_ >= 0) {
5072    sequenceOut_ = pivotVariable_[pivotRow_];
5073    valueOut_ = solution_[sequenceOut_];
5074    lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
5075    upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
5076    if (!toLower) {
5077      directionOut_ = -1;
5078      dualOut_ = valueOut_ - upperOut_;
5079    } else {
5080      directionOut_ = 1;
5081      dualOut_ = lowerOut_ - valueOut_;
5082    }
5083    return 0;
5084  } else {
5085    //theta_=0.0;
5086    return -1;
5087  }
5088}
5089// Restores bound to original bound
5090void 
5091ClpSimplexOther::originalBound(int iSequence, double theta, 
5092                               const double * lowerChange,
5093                               const double * upperChange)
5094{
5095     if (getFakeBound(iSequence) != noFake) {
5096          numberFake_--;
5097          setFakeBound(iSequence, noFake);
5098          if (iSequence >= numberColumns_) {
5099               // rows
5100               int iRow = iSequence - numberColumns_;
5101               rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
5102               rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
5103               if (rowScale_) {
5104                    if (rowLowerWork_[iRow] > -1.0e50)
5105                         rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5106                    if (rowUpperWork_[iRow] < 1.0e50)
5107                         rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5108               } else if (rhsScale_ != 1.0) {
5109                    if (rowLowerWork_[iRow] > -1.0e50)
5110                         rowLowerWork_[iRow] *= rhsScale_;
5111                    if (rowUpperWork_[iRow] < 1.0e50)
5112                         rowUpperWork_[iRow] *= rhsScale_;
5113               }
5114          } else {
5115               // columns
5116               columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
5117               columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
5118               if (rowScale_) {
5119                    double multiplier = 1.0 * inverseColumnScale_[iSequence];
5120                    if (columnLowerWork_[iSequence] > -1.0e50)
5121                         columnLowerWork_[iSequence] *= multiplier * rhsScale_;
5122                    if (columnUpperWork_[iSequence] < 1.0e50)
5123                         columnUpperWork_[iSequence] *= multiplier * rhsScale_;
5124               } else if (rhsScale_ != 1.0) {
5125                    if (columnLowerWork_[iSequence] > -1.0e50)
5126                         columnLowerWork_[iSequence] *= rhsScale_;
5127                    if (columnUpperWork_[iSequence] < 1.0e50)
5128                         columnUpperWork_[iSequence] *= rhsScale_;
5129               }
5130          }
5131     }
5132}
5133/* Expands out all possible combinations for a knapsack
5134   If buildObj NULL then just computes space needed - returns number elements
5135   On entry numberOutput is maximum allowed, on exit it is number needed or
5136   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
5137   least space to return values which reconstruct input.
5138   Rows returned will be original rows but no entries will be returned for
5139   any rows all of whose entries are in knapsack.  So up to user to allow for this.
5140   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
5141   in expanded knapsack.  Values in buildRow and buildElement;
5142*/
5143int
5144ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
5145                                double * buildObj, CoinBigIndex * buildStart,
5146                                int * buildRow, double * buildElement, int reConstruct) const
5147{
5148     int iRow;
5149     int iColumn;
5150     // Get column copy
5151     CoinPackedMatrix * columnCopy = matrix();
5152     // Get a row copy in standard format
5153     CoinPackedMatrix matrixByRow;
5154     matrixByRow.reverseOrderedCopyOf(*columnCopy);
5155     const double * elementByRow = matrixByRow.getElements();
5156     const int * column = matrixByRow.getIndices();
5157     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5158     const int * rowLength = matrixByRow.getVectorLengths();
5159     CoinBigIndex j;
5160     int * whichColumn = new int [numberColumns_];
5161     int * whichRow = new int [numberRows_];
5162     int numJ = 0;
5163     // Get what other columns can compensate for
5164     double * lo = new double [numberRows_];
5165     double * high = new double [numberRows_];
5166     {
5167          // Use to get tight column bounds
5168          ClpSimplex tempModel(*this);
5169          tempModel.tightenPrimalBounds(0.0, 0, true);
5170          // Now another model without knapsacks
5171          int nCol = 0;
5172          for (iRow = 0; iRow < numberRows_; iRow++) {
5173               whichRow[iRow] = iRow;
5174          }
5175          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
5176               whichColumn[iColumn] = -1;
5177          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
5178               int iColumn = column[j];
5179               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
5180                    whichColumn[iColumn] = 0;
5181               } else {
5182                    assert (!columnLower_[iColumn]); // fix later
5183               }
5184          }
5185          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5186               if (whichColumn[iColumn] < 0)
5187                    whichColumn[nCol++] = iColumn;
5188          }
5189          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
5190          // Row copy
5191          CoinPackedMatrix matrixByRow;
5192          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
5193          const double * elementByRow = matrixByRow.getElements();
5194          const int * column = matrixByRow.getIndices();
5195          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
5196          const int * rowLength = matrixByRow.getVectorLengths();
5197          const double * columnLower = tempModel2.getColLower();
5198          const double * columnUpper = tempModel2.getColUpper();
5199          for (iRow = 0; iRow < numberRows_; iRow++) {
5200               lo[iRow] = -COIN_DBL_MAX;
5201               high[iRow] = COIN_DBL_MAX;
5202               if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
5203
5204                    // possible row
5205                    int infiniteUpper = 0;
5206                    int infiniteLower = 0;
5207                    double maximumUp = 0.0;
5208                    double maximumDown = 0.0;
5209                    CoinBigIndex rStart = rowStart[iRow];
5210                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
5211                    CoinBigIndex j;
5212                    // Compute possible lower and upper ranges
5213
5214                    for (j = rStart; j < rEnd; ++j) {
5215                         double value = elementByRow[j];
5216                         iColumn = column[j];
5217                         if (value > 0.0) {
5218                              if (columnUpper[iColumn] >= 1.0e20) {
5219                                   ++infiniteUpper;
5220                              } else {
5221                                   maximumUp += columnUpper[iColumn] * value;
5222                              }
5223                              if (columnLower[iColumn] <= -1.0e20) {
5224                                   ++infiniteLower;
5225                              } else {
5226                                   maximumDown += columnLower[iColumn] * value;
5227                              }
5228                         } else if (valu