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

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

stop seg fault on infinite theta

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