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

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

add ClpTraceDebug?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 260.8 KB
Line 
1/* $Id: ClpSimplexOther.cpp 1817 2011-11-03 09:26:23Z 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      double saveDualBound=dualBound_;
2973      dualBound_=CoinMax(dualBound_,1.0e15);
2974      swapped=true;
2975      double * temp;
2976      memcpy(saveLower,lower_,numberTotal*sizeof(double));
2977      temp=saveLower;
2978      saveLower=lower_;
2979      lower_=temp;
2980      //columnLowerWork_ = lower_;
2981      //rowLowerWork_ = lower_ + numberColumns_;
2982      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
2983      temp=saveUpper;
2984      saveUpper=upper_;
2985      upper_=temp;
2986      //columnUpperWork_ = upper_;
2987      //rowUpperWork_ = upper_ + numberColumns_;
2988      if (rowScale_) {
2989        // scale saved and change arrays
2990        double * lowerChange = lower_+numberTotal;
2991        double * upperChange = upper_+numberTotal;
2992        double * lowerSave = lowerChange+numberTotal;
2993        double * upperSave = upperChange+numberTotal;
2994        for (int i=0;i<numberColumns_;i++) {
2995          double multiplier = inverseColumnScale_[i];
2996          if (lowerSave[i]>-1.0e20)
2997            lowerSave[i] *= multiplier;
2998          if (upperSave[i]<1.0e20)
2999            upperSave[i] *= multiplier;
3000          lowerChange[i] *= multiplier;
3001          upperChange[i] *= multiplier;
3002        }
3003        lowerChange += numberColumns_;
3004        upperChange += numberColumns_;
3005        lowerSave += numberColumns_;
3006        upperSave += numberColumns_;
3007        for (int i=0;i<numberRows_;i++) {
3008          double multiplier = rowScale_[i];
3009          if (lowerSave[i]>-1.0e20)
3010            lowerSave[i] *= multiplier;
3011          if (upperSave[i]<1.0e20)
3012            upperSave[i] *= multiplier;
3013          lowerChange[i] *= multiplier;
3014          upperChange[i] *= multiplier;
3015        }
3016      }
3017      //double saveEndingTheta = endingTheta;
3018      double * saveDuals = NULL;
3019      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3020      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3021        // probably can get rid of this if we adjust every change in theta
3022        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3023        //   sumPrimalInfeasibilities_);
3024        int pass=100;
3025        while(sumPrimalInfeasibilities_) {
3026          pass--;
3027          if (!pass)
3028            break;
3029          problemStatus_=-1;
3030          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3031            double value=solution_[iSequence];
3032            // remember scaling
3033            if (value<lower_[iSequence]-1.0e-9) {
3034              lower_[iSequence]=value;
3035              lowerCopy[iSequence]=value;
3036            } else if (value>upper_[iSequence]+1.0e-9) {
3037              upper_[iSequence]=value;
3038              upperCopy[iSequence]=value;
3039            }
3040          }
3041          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3042        }
3043      }
3044      assert (!problemStatus_);
3045      if (nLowerChange||nUpperChange) {
3046#ifndef ALL_DANTZIG
3047        // Dantzig
3048        savePivot = dualRowPivot_;
3049        dualRowPivot_ = new ClpDualRowDantzig();
3050        dualRowPivot_->setModel(this);
3051#endif
3052        for (int i=0;i<numberRows_+numberColumns_;i++)
3053          setFakeBound(i, noFake);
3054        // Now do parametrics
3055        handler_->message(CLP_PARAMETRICS_STATS, messages_)
3056          << startingTheta << objectiveValue() << CoinMessageEol;
3057        bool canSkipFactorization=true;
3058        while (!returnCode) {
3059                 paramData.startingTheta=startingTheta;
3060                 paramData.endingTheta=endingTheta;
3061                 returnCode = parametricsLoop(paramData,
3062                                       data,canSkipFactorization);
3063                 startingTheta=paramData.startingTheta;
3064                 endingTheta=paramData.endingTheta;
3065          canSkipFactorization=false;
3066          if (!returnCode) {
3067            //startingTheta = endingTheta;
3068            //endingTheta = saveEndingTheta;
3069            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3070              << startingTheta << objectiveValue() << CoinMessageEol;
3071            if (startingTheta >= endingTheta-primalTolerance_
3072                ||problemStatus_==2)
3073              break;
3074          } else if (returnCode == -1) {
3075            // trouble - do external solve
3076            abort(); //needToDoSomething = true;
3077          } else if (problemStatus_==1) {
3078            // can't move any further
3079            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3080              << endingTheta << objectiveValue() << CoinMessageEol;
3081            problemStatus_=0;
3082          }
3083        }
3084      }
3085      dualBound_ = saveDualBound;
3086      //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3087    }
3088    if (problemStatus_==2) {
3089      delete [] ray_;
3090      ray_ = new double [numberColumns_];
3091    }
3092    if (swapped&&lower_) {
3093      double * temp=saveLower;
3094      saveLower=lower_;
3095      lower_=temp;
3096      temp=saveUpper;
3097      saveUpper=upper_;
3098      upper_=temp;
3099    }
3100    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
3101  }   
3102  if (!scalingFlag_) {
3103    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3104    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3105    memcpy(rowLower_,lowerCopy+numberColumns_,
3106           numberRows_*sizeof(double));
3107    memcpy(rowUpper_,upperCopy+numberColumns_,
3108           numberRows_*sizeof(double));
3109  } else {
3110    //  extra for unscaled
3111    double * unscaledCopy;
3112    unscaledCopy = lowerCopy + numberTotal; 
3113    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3114    memcpy(rowLower_,unscaledCopy+numberColumns_,
3115           numberRows_*sizeof(double));
3116    unscaledCopy = upperCopy + numberTotal; 
3117    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3118    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3119           numberRows_*sizeof(double));
3120  }
3121  delete [] saveLower;
3122  delete [] saveUpper;
3123#ifdef ALL_DANTZIG
3124  if (savePivot) {
3125#endif
3126    delete dualRowPivot_;
3127    dualRowPivot_ = savePivot;
3128#ifdef ALL_DANTZIG
3129  }
3130#endif
3131  // Restore any saved stuff
3132  restoreData(data);
3133  perturbation_ = savePerturbation;
3134  delete rowArray_[4];
3135  rowArray_[4]=NULL;
3136  delete rowArray_[5];
3137  rowArray_[5]=NULL;
3138  char line[100];
3139  sprintf(line,"Ending theta %g\n", endingTheta);
3140  handler_->message(CLP_GENERAL,messages_)
3141    << line << CoinMessageEol;
3142  return problemStatus_;
3143}
3144int
3145ClpSimplexOther::parametricsLoop(parametricsData & paramData,
3146                                 ClpDataSave & data,bool canSkipFactorization)
3147{
3148  double & startingTheta = paramData.startingTheta;
3149  double & endingTheta = paramData.endingTheta;
3150  int numberTotal = numberRows_+numberColumns_;
3151  // stuff is already at starting
3152  const int * lowerList = paramData.lowerList;
3153  const int * upperList = paramData.upperList;
3154  problemStatus_ = -1;
3155  //double saveEndingTheta=endingTheta;
3156
3157  // This says whether to restore things etc
3158  // startup will have factorized so can skip
3159  int factorType = 0;
3160  // Start check for cycles
3161  progress_.startCheck();
3162  // Say change made on first iteration
3163  changeMade_ = 1;
3164  /*
3165    Status of problem:
3166    0 - optimal
3167    1 - infeasible
3168    2 - unbounded
3169    -1 - iterating
3170    -2 - factorization wanted
3171    -3 - redo checking without factorization
3172    -4 - looks infeasible
3173  */
3174  while (problemStatus_ < 0) {
3175    int iRow, iColumn;
3176    // clear
3177    for (iRow = 0; iRow < 6; iRow++) {
3178      rowArray_[iRow]->clear();
3179    }
3180   
3181    for (iColumn = 0; iColumn < 2; iColumn++) {
3182      columnArray_[iColumn]->clear();
3183    }
3184   
3185    // give matrix (and model costs and bounds a chance to be
3186    // refreshed (normally null)
3187    matrix_->refresh(this);
3188    // may factorize, checks if problem finished
3189    if (!canSkipFactorization)
3190      statusOfProblemInParametrics(factorType, data);
3191    canSkipFactorization=false;
3192    if (numberPrimalInfeasibilities_) {
3193      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3194        // treat as success
3195        problemStatus_=0;
3196        endingTheta=startingTheta;
3197        break;
3198      }
3199      // probably can get rid of this if we adjust every change in theta
3200      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3201      //     sumPrimalInfeasibilities_);
3202      const double * lowerChange = lower_+numberTotal;
3203      const double * upperChange = upper_+numberTotal;
3204      const double * startLower = lowerChange+numberTotal;
3205      const double * startUpper = upperChange+numberTotal;
3206      //startingTheta -= 1.0e-7;
3207      int nLowerChange = lowerList[-1];
3208      for (int i = 0; i < nLowerChange; i++) {
3209        int iSequence = lowerList[i];
3210        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3211      }
3212      int nUpperChange = upperList[-1];
3213      for (int i = 0; i < nUpperChange; i++) {
3214        int iSequence = upperList[i];
3215        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3216      }
3217      // adjust rhs in case dual uses
3218      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3219      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3220      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3221      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3222      if (rowScale_) {
3223        for (int i=0;i<numberColumns_;i++) {
3224          double multiplier = columnScale_[i];
3225          if (columnLower_[i]>-1.0e20)
3226            columnLower_[i] *= multiplier;
3227          if (columnUpper_[i]<1.0e20)
3228            columnUpper_[i] *= multiplier;
3229        }
3230        for (int i=0;i<numberRows_;i++) {
3231          double multiplier = inverseRowScale_[i];
3232          if (rowLower_[i]>-1.0e20)
3233            rowLower_[i] *= multiplier;
3234          if (rowUpper_[i]<1.0e20)
3235            rowUpper_[i] *= multiplier;
3236        }
3237      }
3238      double * saveDuals = NULL;
3239      problemStatus_=-1;
3240      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3241      int pass=100;
3242      double moved=0.0;
3243      while(sumPrimalInfeasibilities_) {
3244        pass--;
3245        if (!pass)
3246          break;
3247        problemStatus_=-1;
3248        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3249          double value=solution_[iSequence];
3250          if (value<lower_[iSequence]-1.0e-9) {
3251            moved += lower_[iSequence]-value;
3252            lower_[iSequence]=value;
3253          } else if (value>upper_[iSequence]+1.0e-9) {
3254            moved += upper_[iSequence]-value;
3255            upper_[iSequence]=value;
3256          }
3257        }
3258        if (!moved) {
3259          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3260            double value=solution_[iSequence];
3261            if (value<lower_[iSequence]-1.0e-9) {
3262              moved += lower_[iSequence]-value;
3263              lower_[iSequence]=value;
3264            } else if (value>upper_[iSequence]+1.0e-9) {
3265              moved += upper_[iSequence]-value;
3266              upper_[iSequence]=value;
3267            }
3268          }
3269        }
3270        assert (moved);
3271        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3272      }
3273      // adjust
3274      //printf("Should adjust - moved %g\n",moved);
3275    }
3276    // Say good factorization
3277    factorType = 1;
3278    if (data.sparseThreshold_) {
3279      // use default at present
3280      factorization_->sparseThreshold(0);
3281      factorization_->goSparse();
3282    }
3283   
3284    // exit if victory declared
3285    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3286      break;
3287   
3288    // test for maximum iterations
3289    if (hitMaximumIterations()) {
3290      problemStatus_ = 3;
3291      break;
3292    }
3293#ifdef CLP_USER_DRIVEN
3294    // Check event
3295    {
3296      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3297      if (status >= 0) {
3298        problemStatus_ = 5;
3299        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3300        break;
3301      }
3302    }
3303#endif
3304    // Do iterations
3305    problemStatus_=-1;
3306    whileIterating(paramData, 0.0,
3307                   NULL);
3308    //startingTheta = endingTheta;
3309    //endingTheta = saveEndingTheta;
3310  }
3311  if (!problemStatus_/*||problemStatus_==2*/) {
3312    theta_ = endingTheta;
3313#ifdef CLP_USER_DRIVEN
3314    {
3315      double saveTheta=theta_;
3316      theta_ = endingTheta;
3317      int status=eventHandler_->event(ClpEventHandler::theta);
3318      if (status>=0&&status<10) {
3319        endingTheta=theta_;
3320        theta_=saveTheta;
3321        problemStatus_=-1;
3322      } else {
3323        if (status>=10) {
3324          problemStatus_=status-10;
3325          startingTheta=endingTheta;
3326        }
3327        theta_=saveTheta;
3328      }
3329    }
3330#endif
3331    return 0;
3332  } else if (problemStatus_ == 10) {
3333    return -1;
3334  } else {
3335    return problemStatus_;
3336  }
3337}
3338/* Checks if finished.  Updates status */
3339void
3340ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3341{
3342     if (type == 2) {
3343          // trouble - go to recovery
3344          problemStatus_ = 10;
3345          return;
3346     }
3347     if (problemStatus_ > -3 || factorization_->pivots()) {
3348          // factorize
3349          // later on we will need to recover from singularities
3350          // also we could skip if first time
3351          if (type) {
3352               // is factorization okay?
3353               if (internalFactorize(1)) {
3354                    // trouble - go to recovery
3355                    problemStatus_ = 10;
3356                    return;
3357               }
3358          }
3359          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3360               problemStatus_ = -3;
3361     }
3362     // at this stage status is -3 or -4 if looks infeasible
3363     // get primal and dual solutions
3364     gutsOfSolution(NULL, NULL);
3365     double realDualInfeasibilities = sumDualInfeasibilities_;
3366     // If bad accuracy treat as singular
3367     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3368          // trouble - go to recovery
3369          problemStatus_ = 10;
3370          return;
3371     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3372          // Can reduce tolerance
3373          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3374          factorization_->pivotTolerance(newTolerance);
3375     }
3376     // Check if looping
3377     int loop;
3378     if (type != 2)
3379          loop = progress_.looping();
3380     else
3381          loop = -1;
3382     if (loop >= 0) {
3383          problemStatus_ = loop; //exit if in loop
3384          if (!problemStatus_) {
3385               // declaring victory
3386               numberPrimalInfeasibilities_ = 0;
3387               sumPrimalInfeasibilities_ = 0.0;
3388          } else {
3389               problemStatus_ = 10; // instead - try other algorithm
3390          }
3391          return;
3392     } else if (loop < -1) {
3393          // something may have changed
3394          gutsOfSolution(NULL, NULL);
3395     }
3396     progressFlag_ = 0; //reset progress flag
3397     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3398          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3399                    << numberIterations_ << objectiveValue();
3400          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3401                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3402          handler_->printing(sumDualInfeasibilities_ > 0.0)
3403                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3404          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3405                             < numberDualInfeasibilities_)
3406                    << numberDualInfeasibilitiesWithoutFree_;
3407          handler_->message() << CoinMessageEol;
3408     }
3409#ifdef CLP_USER_DRIVEN
3410     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3411       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3412       if (status>=0) {
3413         // fix up
3414         for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
3415           double value=solution_[iSequence];
3416           if (value<=lower_[iSequence]-primalTolerance_) {
3417             lower_[iSequence]=value;
3418           } else if (value>=upper_[iSequence]+primalTolerance_) {
3419             upper_[iSequence]=value;
3420           }
3421         }
3422         numberPrimalInfeasibilities_ = 0;
3423         sumPrimalInfeasibilities_ = 0.0;
3424       }
3425     }
3426#endif
3427     /* If we are primal feasible and any dual infeasibilities are on
3428        free variables then it is better to go to primal */
3429     if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3430               numberDualInfeasibilities_) {
3431          problemStatus_ = 10;
3432          return;
3433     }
3434
3435     // check optimal
3436     // give code benefit of doubt
3437     if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3438               sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3439          // say optimal (with these bounds etc)
3440          numberDualInfeasibilities_ = 0;
3441          sumDualInfeasibilities_ = 0.0;
3442          numberPrimalInfeasibilities_ = 0;
3443          sumPrimalInfeasibilities_ = 0.0;
3444     }
3445     if (dualFeasible() || problemStatus_ == -4) {
3446          progress_.modifyObjective(objectiveValue_
3447                                    - sumDualInfeasibilities_ * dualBound_);
3448     }
3449     if (numberPrimalInfeasibilities_) {
3450          if (problemStatus_ == -4 || problemStatus_ == -5) {
3451               problemStatus_ = 1; // infeasible
3452          }
3453     } else if (numberDualInfeasibilities_) {
3454          // clean up
3455          problemStatus_ = 10;
3456     } else {
3457          problemStatus_ = 0;
3458     }
3459     lastGoodIteration_ = numberIterations_;
3460     if (problemStatus_ < 0) {
3461          sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3462          if (sumDualInfeasibilities_)
3463               numberDualInfeasibilities_ = 1;
3464     }
3465     // Allow matrices to be sorted etc
3466     int fake = -999; // signal sort
3467     matrix_->correctSequence(this, fake, fake);
3468}
3469//static double lastThetaX=0.0;
3470/* This has the flow between re-factorizations
3471   Reasons to come out:
3472   -1 iterations etc
3473   -2 inaccuracy
3474   -3 slight inaccuracy (and done iterations)
3475   +0 looks optimal (might be unbounded - but we will investigate)
3476   +1 looks infeasible
3477   +3 max iterations
3478   +4 accuracy problems
3479*/
3480int
3481ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
3482                                const double * /*changeObjective*/)
3483{
3484  double & startingTheta = paramData.startingTheta;
3485  double & endingTheta = paramData.endingTheta;
3486  const double * lowerChange = paramData.lowerChange;
3487  const double * upperChange = paramData.upperChange;
3488  int numberTotal = numberColumns_ + numberRows_;
3489  const int * lowerList = paramData.lowerList;
3490  const int * upperList = paramData.upperList;
3491  // do basic pointers
3492  int * backwardBasic = paramData.backwardBasic;
3493  for (int i=0;i<numberTotal;i++)
3494    backwardBasic[i]=-1;
3495  for (int i=0;i<numberRows_;i++) {
3496    int iPivot=pivotVariable_[i];
3497    backwardBasic[iPivot]=i;
3498  }
3499     {
3500          int i;
3501          for (i = 0; i < 4; i++) {
3502               rowArray_[i]->clear();
3503          }
3504          for (i = 0; i < 2; i++) {
3505               columnArray_[i]->clear();
3506          }
3507     }
3508     // if can't trust much and long way from optimal then relax
3509     if (largestPrimalError_ > 10.0)
3510          factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3511     else
3512          factorization_->relaxAccuracyCheck(1.0);
3513     // status stays at -1 while iterating, >=0 finished, -2 to invert
3514     // status -3 to go to top without an invert
3515     int returnCode = -1;
3516     double lastTheta = startingTheta;
3517     double useTheta = startingTheta;
3518     while (problemStatus_ == -1) {
3519          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3520          // Get theta for bounds - we know can't crossover
3521          int pivotType = nextTheta(1, increaseTheta, paramData,
3522                                     NULL);
3523          useTheta += theta_;
3524          double change = useTheta - lastTheta;
3525          if (change>1.0e-14) {
3526            int n;
3527            n=lowerList[-1];
3528            for (int i=0;i<n;i++) {
3529              int iSequence = lowerList[i];
3530              lower_[iSequence] += change * lowerChange[iSequence];
3531              if(getStatus(iSequence)==atLowerBound) {
3532                solution_[iSequence] = lower_[iSequence];
3533              }
3534            }
3535            n=upperList[-1];
3536            for (int i=0;i<n;i++) {
3537              int iSequence = upperList[i];
3538              upper_[iSequence] += change * upperChange[iSequence];
3539              if(getStatus(iSequence)==atUpperBound||
3540                 getStatus(iSequence)==isFixed) {
3541                solution_[iSequence] = upper_[iSequence];
3542              }
3543            }
3544          }
3545          sequenceIn_=-1;
3546          if (pivotType) {
3547            if (useTheta>lastTheta+1.0e-9) {
3548              handler_->message(CLP_PARAMETRICS_STATS, messages_)
3549                << useTheta << objectiveValue() << CoinMessageEol;
3550              lastTheta = useTheta;
3551            }
3552            problemStatus_ = -2;
3553            if (!factorization_->pivots()&&pivotRow_<0)
3554              problemStatus_=2;
3555#ifdef CLP_USER_DRIVEN
3556            {
3557              double saveTheta=theta_;
3558              theta_ = endingTheta;
3559              int status=eventHandler_->event(ClpEventHandler::theta);
3560              if (status>=0&&status<10) {
3561                endingTheta=theta_;
3562                theta_=saveTheta;
3563                problemStatus_=-1;
3564                continue;
3565              } else {
3566                if (status>=10)
3567                  problemStatus_=status-10;
3568                if (status<0)
3569                  startingTheta = useTheta;
3570                theta_=saveTheta;
3571              }
3572            }
3573#else
3574            startingTheta = useTheta;
3575#endif
3576            return 4;
3577          }
3578          // choose row to go out
3579          //reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3580          if (pivotRow_ >= 0) {
3581               // we found a pivot row
3582               if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3583                    handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3584                              << pivotRow_
3585                              << CoinMessageEol;
3586               }
3587               // check accuracy of weights
3588               dualRowPivot_->checkAccuracy();
3589               // do ratio test for normal iteration
3590               double bestPossiblePivot = bestPivot();
3591               if (sequenceIn_ >= 0) {
3592                    // normal iteration
3593                    // update the incoming column
3594                    double btranAlpha = -alpha_ * directionOut_; // for check
3595                    unpackPacked(rowArray_[1]);
3596                    // and update dual weights (can do in parallel - with extra array)
3597                    rowArray_[2]->clear();
3598                    alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3599                                                          rowArray_[2],
3600                                                          rowArray_[3],
3601                                                          rowArray_[1]);
3602                    // see if update stable
3603#ifdef CLP_DEBUG
3604                    if ((handler_->logLevel() & 32))
3605                         printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3606#endif
3607                    double checkValue = 1.0e-7;
3608                    // if can't trust much and long way from optimal then relax
3609                    if (largestPrimalError_ > 10.0)
3610                         checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3611                    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3612                              fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3613                         handler_->message(CLP_DUAL_CHECK, messages_)
3614                                   << btranAlpha
3615                                   << alpha_
3616                                   << CoinMessageEol;
3617                         // clear arrays
3618                         rowArray_[4]->clear();
3619                         if (factorization_->pivots()) {
3620                              dualRowPivot_->unrollWeights();
3621                              problemStatus_ = -2; // factorize now
3622                              rowArray_[0]->clear();
3623                              rowArray_[1]->clear();
3624                              columnArray_[0]->clear();
3625                              returnCode = -2;
3626                              break;
3627                         } else {
3628                              // take on more relaxed criterion
3629                              double test;
3630                              if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
3631                                   test = 1.0e-1 * fabs(alpha_);
3632                              else
3633                                   test = 1.0e-4 * (1.0 + fabs(alpha_));
3634                              if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
3635                                        fabs(btranAlpha - alpha_) > test) {
3636                                   dualRowPivot_->unrollWeights();
3637                                   // need to reject something
3638                                   char x = isColumn(sequenceOut_) ? 'C' : 'R';
3639                                   handler_->message(CLP_SIMPLEX_FLAG, messages_)
3640                                             << x << sequenceWithin(sequenceOut_)
3641                                             << CoinMessageEol;
3642                                   setFlagged(sequenceOut_);
3643                                   progress_.clearBadTimes();
3644                                   lastBadIteration_ = numberIterations_; // say be more cautious
3645                                   rowArray_[0]->clear();
3646                                   rowArray_[1]->clear();
3647                                   columnArray_[0]->clear();
3648                                   if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3649                                        //printf("I think should declare infeasible\n");
3650                                        problemStatus_ = 1;
3651                                        returnCode = 1;
3652                                        break;
3653                                   }
3654                                   continue;
3655                              }
3656                         }
3657                    }
3658                    // update duals BEFORE replaceColumn so can do updateColumn
3659                    double objectiveChange = 0.0;
3660                    // do duals first as variables may flip bounds
3661                    // rowArray_[0] and columnArray_[0] may have flips
3662                    // so use rowArray_[3] for work array from here on
3663                    int nswapped = 0;
3664                    //rowArray_[0]->cleanAndPackSafe(1.0e-60);
3665                    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
3666#if CLP_CAN_HAVE_ZERO_OBJ
3667                    if ((specialOptions_&2097152)==0) {
3668#endif
3669                      nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3670                                                                                               rowArray_[2], theta_,
3671                                                                                               objectiveChange, false);
3672                      assert (!nswapped);
3673#if CLP_CAN_HAVE_ZERO_OBJ
3674                    } else {
3675                      rowArray_[0]->clear();
3676                      rowArray_[2]->clear();
3677                      columnArray_[0]->clear();
3678                    }
3679#endif
3680                    // which will change basic solution
3681                    if (nswapped) {
3682                         factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3683                         dualRowPivot_->updatePrimalSolution(rowArray_[2],
3684                                                             1.0, objectiveChange);
3685                         // recompute dualOut_
3686                         valueOut_ = solution_[sequenceOut_];
3687                         if (directionOut_ < 0) {
3688                              dualOut_ = valueOut_ - upperOut_;
3689                         } else {
3690                              dualOut_ = lowerOut_ - valueOut_;
3691                         }
3692                    }
3693                    // amount primal will move
3694                    double movement = -dualOut_ * directionOut_ / alpha_;
3695                    // so objective should increase by fabs(dj)*movement
3696                    // but we already have objective change - so check will be good
3697                    if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3698#ifdef CLP_DEBUG
3699                         if (handler_->logLevel() & 32)
3700                              printf("movement %g, swap change %g, rest %g  * %g\n",
3701                                     objectiveChange + fabs(movement * dualIn_),
3702                                     objectiveChange, movement, dualIn_);
3703#endif
3704                         assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3705                         if(factorization_->pivots()) {
3706                              // going backwards - factorize
3707                              dualRowPivot_->unrollWeights();
3708                              problemStatus_ = -2; // factorize now
3709                              returnCode = -2;
3710                              break;
3711                         }
3712                    }
3713                    CoinAssert(fabs(dualOut_) < 1.0e50);
3714                    // if stable replace in basis
3715                    int updateStatus = factorization_->replaceColumn(this,
3716                                       rowArray_[2],
3717                                       rowArray_[1],
3718                                       pivotRow_,
3719                                       alpha_);
3720                    // if no pivots, bad update but reasonable alpha - take and invert
3721                    if (updateStatus == 2 &&
3722                              !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3723                         updateStatus = 4;
3724                    if (updateStatus == 1 || updateStatus == 4) {
3725                         // slight error
3726                         if (factorization_->pivots() > 5 || updateStatus == 4) {
3727                              problemStatus_ = -2; // factorize now
3728                              returnCode = -3;
3729                         }
3730                    } else if (updateStatus == 2) {
3731                         // major error
3732                         dualRowPivot_->unrollWeights();
3733                         // later we may need to unwind more e.g. fake bounds
3734                         if (factorization_->pivots()) {
3735                              problemStatus_ = -2; // factorize now
3736                              returnCode = -2;
3737                              break;
3738                         } else {
3739                              // need to reject something
3740                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3741                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3742                                        << x << sequenceWithin(sequenceOut_)
3743                                        << CoinMessageEol;
3744                              setFlagged(sequenceOut_);
3745                              progress_.clearBadTimes();
3746                              lastBadIteration_ = numberIterations_; // say be more cautious
3747                              rowArray_[0]->clear();
3748                              rowArray_[1]->clear();
3749                              columnArray_[0]->clear();
3750                              // make sure dual feasible
3751                              // look at all rows and columns
3752                              double objectiveChange = 0.0;
3753                              reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3754                                        0.0, objectiveChange, true);
3755                              continue;
3756                         }
3757                    } else if (updateStatus == 3) {
3758                         // out of memory
3759                         // increase space if not many iterations
3760                         if (factorization_->pivots() <
3761                                   0.5 * factorization_->maximumPivots() &&
3762                                   factorization_->pivots() < 200)
3763                              factorization_->areaFactor(
3764                                   factorization_->areaFactor() * 1.1);
3765                         problemStatus_ = -2; // factorize now
3766                    } else if (updateStatus == 5) {
3767                         problemStatus_ = -2; // factorize now
3768                    }
3769                    // update change vector
3770                    {
3771                      double * work = rowArray_[1]->denseVector();
3772                      int number = rowArray_[1]->getNumElements();
3773                      int * which = rowArray_[1]->getIndices();
3774                      assert (!rowArray_[4]->packedMode());
3775                      assert (rowArray_[1]->packedMode());
3776                      double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3777                      double multiplier = -pivotValue/alpha_;
3778                      if (multiplier) {
3779                        for (int i = 0; i < number; i++) {
3780                          int iRow = which[i];
3781                          rowArray_[4]->quickAddNonZero(iRow,multiplier*work[i]);
3782                        }
3783                      }
3784                      pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3785                      // we want pivot to be -multiplier
3786                      rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3787                    }
3788                    // update primal solution
3789#if CLP_CAN_HAVE_ZERO_OBJ
3790                    if ((specialOptions_&2097152)!=0) 
3791                      theta_=0.0;
3792#endif
3793                    if (theta_ < 0.0) {
3794#ifdef CLP_DEBUG
3795                         if (handler_->logLevel() & 32)
3796                              printf("negative theta %g\n", theta_);
3797#endif
3798                         theta_ = 0.0;
3799                    }
3800                    // do actual flips
3801                    reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
3802                    //rowArray_[1]->expand();
3803                    dualRowPivot_->updatePrimalSolution(rowArray_[1],
3804                                                        movement,
3805                                                        objectiveChange);
3806                    // modify dualout
3807                    dualOut_ /= alpha_;
3808                    dualOut_ *= -directionOut_;
3809                    //setStatus(sequenceIn_,basic);
3810                    dj_[sequenceIn_] = 0.0;
3811                    //double oldValue = valueIn_;
3812                    if (directionIn_ == -1) {
3813                         // as if from upper bound
3814                         valueIn_ = upperIn_ + dualOut_;
3815                    } else {
3816                         // as if from lower bound
3817                         valueIn_ = lowerIn_ + dualOut_;
3818                    }
3819                    objectiveChange = 0.0;
3820#if CLP_CAN_HAVE_ZERO_OBJ
3821                    if ((specialOptions_&2097152)==0) {
3822#endif
3823                      for (int i=0;i<numberTotal;i++)
3824                        objectiveChange += solution_[i]*cost_[i];
3825                      objectiveChange -= objectiveValue_;
3826#if CLP_CAN_HAVE_ZERO_OBJ
3827                    }
3828#endif
3829                    // outgoing
3830                    originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
3831                    lowerOut_=lower_[sequenceOut_];
3832                    upperOut_=upper_[sequenceOut_];
3833                    // set dj to zero unless values pass
3834                    if (directionOut_ > 0) {
3835                         valueOut_ = lowerOut_;
3836                         dj_[sequenceOut_] = theta_;
3837#if CLP_CAN_HAVE_ZERO_OBJ>1
3838#ifdef COIN_REUSE_RANDOM
3839                         if ((specialOptions_&2097152)!=0) {
3840                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
3841                         }
3842#endif
3843#endif
3844                    } else {
3845                         valueOut_ = upperOut_;
3846                         dj_[sequenceOut_] = -theta_;
3847#if CLP_CAN_HAVE_ZERO_OBJ>1
3848#ifdef COIN_REUSE_RANDOM
3849                         if ((specialOptions_&2097152)!=0) {
3850                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
3851                         }
3852#endif
3853#endif
3854                    }
3855                    solution_[sequenceOut_] = valueOut_;
3856                    int whatNext = housekeeping(objectiveChange);
3857                    assert (backwardBasic[sequenceOut_]==pivotRow_);
3858                    backwardBasic[sequenceOut_]=-1;
3859                    backwardBasic[sequenceIn_]=pivotRow_;
3860                    {
3861                      char in[200],out[200];
3862                      int iSequence=sequenceIn_;
3863                      if (iSequence<numberColumns_) {
3864                        if (lengthNames_) 
3865                          strcpy(in,columnNames_[iSequence].c_str());
3866                         else 
3867                          sprintf(in,"C%7.7d",iSequence);
3868                      } else {
3869                        iSequence -= numberColumns_;
3870                        if (lengthNames_) 
3871                          strcpy(in,rowNames_[iSequence].c_str());
3872                         else 
3873                          sprintf(in,"R%7.7d",iSequence);
3874                      }
3875                      iSequence=sequenceOut_;
3876                      if (iSequence<numberColumns_) {
3877                        if (lengthNames_) 
3878                          strcpy(out,columnNames_[iSequence].c_str());
3879                         else 
3880                          sprintf(out,"C%7.7d",iSequence);
3881                      } else {
3882                        iSequence -= numberColumns_;
3883                        if (lengthNames_) 
3884                          strcpy(out,rowNames_[iSequence].c_str());
3885                         else 
3886                          sprintf(out,"R%7.7d",iSequence);
3887                      }
3888                      handler_->message(CLP_PARAMETRICS_STATS2, messages_)
3889                        << useTheta << objectiveValue() 
3890                        << in << out << CoinMessageEol;
3891                    }
3892                    if (useTheta>lastTheta+1.0e-9) {
3893                      handler_->message(CLP_PARAMETRICS_STATS, messages_)
3894                        << useTheta << objectiveValue() << CoinMessageEol;
3895                      lastTheta = useTheta;
3896                    }
3897                    // and set bounds correctly
3898                    originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
3899                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
3900                    if (whatNext == 1) {
3901                         problemStatus_ = -2; // refactorize
3902                    } else if (whatNext == 2) {
3903                         // maximum iterations or equivalent
3904                         problemStatus_ = 3;
3905                         returnCode = 3;
3906                         break;
3907                    }
3908#ifdef CLP_USER_DRIVEN
3909                    // Check event
3910                    {
3911                         int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3912                         if (status >= 0) {
3913                              problemStatus_ = 5;
3914                              secondaryStatus_ = ClpEventHandler::endOfIteration;
3915                              returnCode = 4;
3916                              break;
3917                         }
3918                    }
3919#endif
3920               } else {
3921                    // no incoming column is valid
3922#ifdef CLP_USER_DRIVEN
3923                 rowArray_[0]->clear();
3924                 columnArray_[0]->clear();
3925                 theta_ = useTheta;
3926                 lastTheta = useTheta;
3927                 int action = eventHandler_->event(ClpEventHandler::noTheta);
3928                 if (action>=0) {
3929                   endingTheta=theta_;
3930                   theta_ = 0.0;
3931                   //adjust [4] from handler - but
3932                   //rowArray_[4]->clear(); // temp
3933                   if (action>=0&&action<10)
3934                     problemStatus_=-1; // carry on
3935                   else if (action==15)
3936                     problemStatus_ =5; // say stopped
3937                   returnCode = 1;
3938                   if (action==0||action>=10) 
3939                     break;
3940                   else
3941                     continue;
3942                 } else {
3943                 theta_ = 0.0;
3944                 }
3945#endif
3946                    pivotRow_ = -1;
3947#ifdef CLP_DEBUG
3948                    if (handler_->logLevel() & 32)
3949                         printf("** no column pivot\n");
3950#endif
3951                    if (factorization_->pivots() < 10) { 
3952                         // If we have just factorized and infeasibility reasonable say infeas
3953                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
3954                              if (valueOut_ > upperOut_ + 1.0e-3 || valueOut_ < lowerOut_ - 1.0e-3
3955                                        || (specialOptions_ & 64) == 0) {
3956                                   // say infeasible
3957                                   problemStatus_ = 1;
3958                                   // unless primal feasible!!!!
3959                                   //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
3960                                   //     numberDualInfeasibilities_,sumDualInfeasibilities_);
3961                                   if (numberDualInfeasibilities_)
3962                                        problemStatus_ = 10;
3963                                   rowArray_[0]->clear();
3964                                   columnArray_[0]->clear();
3965                              }
3966                         }
3967                         // If special option set - put off as long as possible
3968                         if ((specialOptions_ & 64) == 0) {
3969                              problemStatus_ = -4; //say looks infeasible
3970                         } else {
3971                              // flag
3972                              char x = isColumn(sequenceOut_) ? 'C' : 'R';
3973                              handler_->message(CLP_SIMPLEX_FLAG, messages_)
3974                                        << x << sequenceWithin(sequenceOut_)
3975                                        << CoinMessageEol;
3976                              setFlagged(sequenceOut_);
3977                              if (!factorization_->pivots()) {
3978                                   rowArray_[0]->clear();
3979                                   columnArray_[0]->clear();
3980                                   continue;
3981                              }
3982                         }
3983                    }
3984                    rowArray_[0]->clear();
3985                    columnArray_[0]->clear();
3986                    returnCode = 1;
3987                    break;
3988               }
3989          } else {
3990               // no pivot row
3991#ifdef CLP_USER_DRIVEN
3992            {
3993              double saveTheta=theta_;
3994              theta_ = endingTheta;
3995              int status=eventHandler_->event(ClpEventHandler::theta);
3996              if (status>=0&&status<10) {
3997                endingTheta=theta_;
3998                theta_=saveTheta;
3999                continue;
4000              } else {
4001                theta_=saveTheta;
4002              }
4003            }
4004#endif
4005#ifdef CLP_DEBUG
4006               if (handler_->logLevel() & 32)
4007                    printf("** no row pivot\n");
4008#endif
4009               int numberPivots = factorization_->pivots();
4010               bool specialCase;
4011               int useNumberFake;
4012               returnCode = 0;
4013               if (numberPivots < 20 &&
4014                         (specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
4015                         && dualBound_ > 1.0e8) {
4016                    specialCase = true;
4017                    // as dual bound high - should be okay
4018                    useNumberFake = 0;
4019               } else {
4020                    specialCase = false;
4021                    useNumberFake = numberFake_;
4022               }
4023               if (!numberPivots || specialCase) {
4024                    // may have crept through - so may be optimal
4025                    // check any flagged variables
4026                    int iRow;
4027                    for (iRow = 0; iRow < numberRows_; iRow++) {
4028                         int iPivot = pivotVariable_[iRow];
4029                         if (flagged(iPivot))
4030                              break;
4031                    }
4032                    if (iRow < numberRows_ && numberPivots) {
4033                         // try factorization
4034                         returnCode = -2;
4035                    }
4036
4037                    if (useNumberFake || numberDualInfeasibilities_) {
4038                         // may be dual infeasible
4039                         problemStatus_ = -5;
4040                    } else {
4041                         if (iRow < numberRows_) {
4042                              problemStatus_ = -5;
4043                         } else {
4044                              if (numberPivots) {
4045                                   // objective may be wrong
4046                                   objectiveValue_ = innerProduct(cost_,
4047                                                                  numberColumns_ + numberRows_,
4048                                                                  solution_);
4049                                   objectiveValue_ += objective_->nonlinearOffset();
4050                                   objectiveValue_ /= (objectiveScale_ * rhsScale_);
4051                                   if ((specialOptions_ & 16384) == 0) {
4052                                        // and dual_ may be wrong (i.e. for fixed or basic)
4053                                        CoinIndexedVector * arrayVector = rowArray_[1];
4054                                        arrayVector->clear();
4055                                        int iRow;
4056                                        double * array = arrayVector->denseVector();
4057                                        /* Use dual_ instead of array
4058                                           Even though dual_ is only numberRows_ long this is
4059                                           okay as gets permuted to longer rowArray_[2]
4060                                        */
4061                                        arrayVector->setDenseVector(dual_);
4062                                        int * index = arrayVector->getIndices();
4063                                        int number = 0;
4064                                        for (iRow = 0; iRow < numberRows_; iRow++) {
4065                                             int iPivot = pivotVariable_[iRow];
4066                                             double value = cost_[iPivot];
4067                                             dual_[iRow] = value;
4068                                             if (value) {
4069                                                  index[number++] = iRow;
4070                                             }
4071                                        }
4072                                        arrayVector->setNumElements(number);
4073                                        // Extended duals before "updateTranspose"
4074                                        matrix_->dualExpanded(this, arrayVector, NULL, 0);
4075                                        // Btran basic costs
4076                                        rowArray_[2]->clear();
4077                                        factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
4078                                        // and return vector
4079                                        arrayVector->setDenseVector(array);
4080                                   }
4081                              }
4082                              problemStatus_ = 0;
4083                              sumPrimalInfeasibilities_ = 0.0;
4084                              if ((specialOptions_&(1024 + 16384)) != 0) {
4085                                   CoinIndexedVector * arrayVector = rowArray_[1];
4086                                   arrayVector->clear();
4087                                   double * rhs = arrayVector->denseVector();
4088                                   times(1.0, solution_, rhs);
4089                                   bool bad2 = false;
4090                                   int i;
4091                                   for ( i = 0; i < numberRows_; i++) {
4092                                        if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
4093                                                  rhs[i] > rowUpperWork_[i] + primalTolerance_) {
4094                                             bad2 = true;
4095                                        } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
4096                                        }
4097                                        rhs[i] = 0.0;
4098                                   }
4099                                   for ( i = 0; i < numberColumns_; i++) {
4100                                        if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
4101                                                  solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4102                                             bad2 = true;
4103                                        }
4104                                   }
4105                                   if (bad2) {
4106                                        problemStatus_ = -3;
4107                                        returnCode = -2;
4108                                        // Force to re-factorize early next time
4109                                        int numberPivots = factorization_->pivots();
4110                                        forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4111                                   }
4112                              }
4113                         }
4114                    }
4115               } else {
4116                    problemStatus_ = -3;
4117                    returnCode = -2;
4118                    // Force to re-factorize early next time
4119                    int numberPivots = factorization_->pivots();
4120                    forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4121               }
4122               break;
4123          }
4124     }
4125     startingTheta = lastTheta+theta_;
4126     return returnCode;
4127}
4128// Finds best possible pivot
4129double 
4130ClpSimplexOther::bestPivot(bool justColumns)
4131{
4132  // Get good size for pivot
4133  // Allow first few iterations to take tiny
4134  double acceptablePivot = 1.0e-9;
4135  if (numberIterations_ > 100)
4136    acceptablePivot = 1.0e-8;
4137  if (factorization_->pivots() > 10 ||
4138      (factorization_->pivots() && sumDualInfeasibilities_))
4139    acceptablePivot = 1.0e-5; // if we have iterated be more strict
4140  else if (factorization_->pivots() > 5)
4141    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4142  else if (factorization_->pivots())
4143    acceptablePivot = 1.0e-8; // relax
4144  double bestPossiblePivot = 1.0;
4145  // get sign for finding row of tableau
4146  // normal iteration
4147  // create as packed
4148  double direction = directionOut_;
4149  rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4150  factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4151  // put row of tableau in rowArray[0] and columnArray[0]
4152  matrix_->transposeTimes(this, -1.0,
4153                          rowArray_[0], rowArray_[3], columnArray_[0]);
4154  sequenceIn_=-1;
4155  if (justColumns)
4156    rowArray_[0]->clear();
4157  // do ratio test for normal iteration
4158  bestPossiblePivot = 
4159    reinterpret_cast<ClpSimplexDual *> 
4160    ( this)->dualColumn(rowArray_[0],
4161                        columnArray_[0], columnArray_[1],
4162                        rowArray_[3], acceptablePivot, NULL);
4163  return bestPossiblePivot;
4164}
4165// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4166int
4167ClpSimplexOther::nextTheta(int type, double maxTheta, parametricsData & paramData,
4168                           const double * /*changeObjective*/)
4169{
4170  const double * lowerChange = paramData.lowerChange;
4171  const double * upperChange = paramData.upperChange;
4172  const int * lowerList = paramData.lowerList;
4173  const int * upperList = paramData.upperList;
4174  int iSequence;
4175  theta_ = maxTheta;
4176  bool toLower = false;
4177  assert (type==1);
4178  // may need to decide based on model?
4179  bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4180  double * array = rowArray_[4]->denseVector();
4181  const int * row = matrix_->getIndices();
4182  const int * columnLength = matrix_->getVectorLengths();
4183  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4184  const double * elementByColumn = matrix_->getElements();
4185#if 0
4186  double tempArray[5000];
4187  bool checkIt=false;
4188  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
4189    memcpy(tempArray,array,numberRows_*sizeof(double));
4190    checkIt=true;
4191    needFullUpdate=true;
4192  }
4193#endif
4194  if (!factorization_->pivots()||needFullUpdate) {
4195    rowArray_[4]->clear();
4196    // get change
4197    if (!rowScale_) {
4198      int n;
4199      n=lowerList[-2];
4200      int i;
4201      for (i=0;i<n;i++) {
4202        int iSequence = lowerList[i];
4203        assert (iSequence<numberColumns_);
4204        if (getColumnStatus(iSequence)==atLowerBound) {
4205          double value=lowerChange[iSequence];
4206          for (CoinBigIndex j = columnStart[iSequence];
4207               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4208            rowArray_[4]->quickAddNonZero(row[j], elementByColumn[j]*value);
4209          }
4210        }
4211      }
4212      n=lowerList[-1];
4213      const double * change = lowerChange+numberColumns_;
4214      for (;i<n;i++) {
4215        int iSequence = lowerList[i]-numberColumns_;
4216        assert (iSequence>=0);
4217        if (getRowStatus(iSequence)==atLowerBound) {
4218          double value=change[iSequence];
4219          rowArray_[4]->quickAddNonZero(iSequence, -value);
4220        }
4221      }
4222      n=upperList[-2];
4223      for (i=0;i<n;i++) {
4224        int iSequence = upperList[i];
4225        assert (iSequence<numberColumns_);
4226        if (getColumnStatus(iSequence)==atUpperBound) {
4227          double value=upperChange[iSequence];
4228          for (CoinBigIndex j = columnStart[iSequence];
4229               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4230            rowArray_[4]->quickAddNonZero(row[j], elementByColumn[j]*value);
4231          }
4232        }
4233      }
4234      n=upperList[-1];
4235      change = upperChange+numberColumns_;
4236      for (;i<n;i++) {
4237        int iSequence = upperList[i]-numberColumns_;
4238        assert (iSequence>=0);
4239        if (getRowStatus(iSequence)==atUpperBound) {
4240          double value=change[iSequence];
4241          rowArray_[4]->quickAddNonZero(iSequence, -value);
4242        }
4243      }
4244    } else {
4245      int n;
4246      n=lowerList[-2];
4247      int i;
4248      for (i=0;i<n;i++) {
4249        int iSequence = lowerList[i];
4250        assert (iSequence<numberColumns_);
4251        if (getColumnStatus(iSequence)==atLowerBound) {
4252          double value=lowerChange[iSequence];
4253          // apply scaling
4254          double scale = columnScale_[iSequence];
4255          for (CoinBigIndex j = columnStart[iSequence];
4256               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4257            int iRow = row[j];
4258            rowArray_[4]->quickAddNonZero(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4259          }
4260        }
4261      }
4262      n=lowerList[-1];
4263      const double * change = lowerChange+numberColumns_;
4264      for (;i<n;i++) {
4265        int iSequence = lowerList[i]-numberColumns_;
4266        assert (iSequence>=0);
4267        if (getRowStatus(iSequence)==atLowerBound) {
4268          double value=change[iSequence];
4269          rowArray_[4]->quickAddNonZero(iSequence, -value);
4270        }
4271      }
4272      n=upperList[-2];
4273      for (i=0;i<n;i++) {
4274        int iSequence = upperList[i];
4275        assert (iSequence<numberColumns_);
4276        if (getColumnStatus(iSequence)==atUpperBound) {
4277          double value=upperChange[iSequence];
4278          // apply scaling
4279          double scale = columnScale_[iSequence];
4280          for (CoinBigIndex j = columnStart[iSequence];
4281               j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4282            int iRow = row[j];
4283            rowArray_[4]->quickAddNonZero(iRow, elementByColumn[j]*scale * rowScale_[iRow]*value);
4284          }
4285        }
4286      }
4287      n=upperList[-1];
4288      change = upperChange+numberColumns_;
4289      for (;i<n;i++) {
4290        int iSequence = upperList[i]-numberColumns_;
4291        assert (iSequence>=0);
4292        if (getRowStatus(iSequence)==atUpperBound) {
4293          double value=change[iSequence];
4294          rowArray_[4]->quickAddNonZero(iSequence, -value);
4295        }
4296      }
4297    }
4298    // ftran it
4299    factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4300#if 0
4301    if (checkIt) {
4302      for (int i=0;i<numberRows_;i++) {
4303        assert (fabs(tempArray[i]-array[i])<1.0e-8);
4304      }
4305    }
4306#endif
4307  } else if (sequenceIn_>=0) {
4308    //assert (sequenceIn_>=0);
4309    assert (sequenceOut_>=0);
4310    assert (sequenceIn_!=sequenceOut_);
4311    double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4312    int needed=0;
4313    assert (!rowArray_[5]->getNumElements());
4314    if (change) {
4315      if (sequenceIn_<numberColumns_) {
4316        if (!rowScale_) {
4317          for (CoinBigIndex i = columnStart[sequenceIn_];
4318               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4319            rowArray_[5]->quickAddNonZero(row[i], elementByColumn[i]*change);
4320          }
4321        } else {
4322          // apply scaling
4323          double scale = columnScale_[sequenceIn_];
4324          for (CoinBigIndex i = columnStart[sequenceIn_];
4325               i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4326            int iRow = row[i];
4327            rowArray_[5]->quickAddNonZero(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4328          }
4329        }
4330      } else {
4331        rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4332      }
4333      needed++;
4334    }
4335    if (getStatus(sequenceOut_)==atLowerBound)
4336      change=lowerChange[sequenceOut_];
4337    else
4338      change=upperChange[sequenceOut_];
4339    if (change) {
4340      if (sequenceOut_<numberColumns_) {
4341        if (!rowScale_) {
4342          for (CoinBigIndex i = columnStart[sequenceOut_];
4343               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4344            rowArray_[5]->quickAddNonZero(row[i], elementByColumn[i]*change);
4345          }
4346        } else {
4347          // apply scaling
4348          double scale = columnScale_[sequenceOut_];
4349          for (CoinBigIndex i = columnStart[sequenceOut_];
4350               i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4351            int iRow = row[i];
4352            rowArray_[5]->quickAddNonZero(iRow, elementByColumn[i]*scale * rowScale_[iRow]*change);
4353          }
4354        }
4355      } else {
4356        rowArray_[5]->quickAddNonZero(sequenceOut_-numberColumns_,-change);
4357      }
4358      needed++;
4359    }
4360    //printf("seqin %d seqout %d needed %d\n",
4361    //     sequenceIn_,sequenceOut_,needed);
4362    if (needed) {
4363      // ftran it
4364      factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4365      // add
4366      double * array5 = rowArray_[5]->denseVector();
4367      int * index5 = rowArray_[5]->getIndices();
4368      int number5 = rowArray_[5]->getNumElements();
4369      for (int i = 0; i < number5; i++) {
4370        int iPivot = index5[i];
4371        rowArray_[4]->quickAddNonZero(iPivot,array5[iPivot]);
4372        array5[iPivot]=0.0;
4373      }
4374      rowArray_[5]->setNumElements(0);
4375    }
4376  }
4377  pivotRow_ = -1;
4378  const int * index = rowArray_[4]->getIndices();
4379  int number = rowArray_[4]->getNumElements();
4380  char * markDone = paramData.markDone;
4381  const int * backwardBasic = paramData.backwardBasic;
4382  // first ones with alpha
4383  for (int i=0;i<number;i++) {
4384    int iPivot=index[i];
4385    iSequence = pivotVariable_[iPivot];
4386    assert(!markDone[iSequence]);
4387    markDone[iSequence]=1;
4388    // solution value will be sol - theta*alpha
4389    // bounds will be bounds + change *theta
4390    double currentSolution = solution_[iSequence];
4391    double alpha = array[iPivot];
4392    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
4393    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
4394    if (thetaCoefficientLower > 1.0e-8) {
4395      double currentLower = lower_[iSequence];
4396      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4397      double gap=currentSolution-currentLower;
4398      if (thetaCoefficientLower*theta_>gap) {
4399        theta_ = gap/thetaCoefficientLower;
4400        toLower=true;
4401        pivotRow_=iPivot;
4402      }
4403    }
4404    if (thetaCoefficientUpper < -1.0e-8) {
4405      double currentUpper = upper_[iSequence];
4406      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4407      double gap=currentSolution-currentUpper; //negative
4408      if (thetaCoefficientUpper*theta_<gap) {
4409        theta_ = gap/thetaCoefficientUpper;
4410        toLower=false;
4411        pivotRow_=iPivot;
4412      }
4413    }
4414  }
4415  // now others
4416  int nLook=lowerList[-1];
4417  for (int i=0;i<nLook;i++) {
4418    int iSequence = lowerList[i];
4419    if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) {
4420      double currentSolution = solution_[iSequence];
4421      double currentLower = lower_[iSequence];
4422      ClpTraceDebug (currentSolution >= currentLower - 100.0*primalTolerance_);
4423      double thetaCoefficient = lowerChange[iSequence];
4424      if (thetaCoefficient > 0.0) {
4425        double gap=currentSolution-currentLower;
4426        if (thetaCoefficient*theta_>gap) {
4427          theta_ = gap/thetaCoefficient;
4428          toLower=true;
4429          pivotRow_ = backwardBasic[iSequence];
4430        }
4431      }
4432    }
4433  }
4434  nLook=upperList[-1];
4435  for (int i=0;i<nLook;i++) {
4436    int iSequence = upperList[i];
4437    if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) {
4438      double currentSolution = solution_[iSequence];
4439      double currentUpper = upper_[iSequence];
4440      ClpTraceDebug (currentSolution <= currentUpper + 100.0*primalTolerance_);
4441      double thetaCoefficient = upperChange[iSequence];
4442      if (thetaCoefficient < 0) {
4443        double gap=currentSolution-currentUpper; //negative
4444        if (thetaCoefficient*theta_<gap) {
4445          theta_ = gap/thetaCoefficient;
4446          toLower=false;
4447          pivotRow_ = backwardBasic[iSequence];
4448        }
4449      }
4450    }
4451  }
4452  theta_ = CoinMax(theta_,0.0);
4453  if (theta_>1.0e-15) {
4454    // update solution
4455    for (int iRow = 0; iRow < number; iRow++) {
4456      int iPivot = index[iRow];
4457      iSequence = pivotVariable_[iPivot];
4458      markDone[iSequence]=0;
4459      // solution value will be sol - theta*alpha
4460      double alpha = array[iPivot];
4461      solution_[iSequence] -= theta_ * alpha;
4462    }
4463  } else {
4464    for (int iRow = 0; iRow < number; iRow++) {
4465      int iPivot = index[iRow];
4466      iSequence = pivotVariable_[iPivot];
4467      markDone[iSequence]=0;
4468    }
4469  }
4470#if 0
4471  for (int i=0;i<numberTotal;i++)
4472    assert(!markDone[i]);
4473#endif
4474  if (pivotRow_ >= 0) {
4475    sequenceOut_ = pivotVariable_[pivotRow_];
4476    valueOut_ = solution_[sequenceOut_];
4477    lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
4478    upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
4479    if (!toLower) {
4480      directionOut_ = -1;
4481      dualOut_ = valueOut_ - upperOut_;
4482    } else {
4483      directionOut_ = 1;
4484      dualOut_ = lowerOut_ - valueOut_;
4485    }
4486    return 0;
4487  } else {
4488    //theta_=0.0;
4489    return -1;
4490  }
4491}
4492// Restores bound to original bound
4493void 
4494ClpSimplexOther::originalBound(int iSequence, double theta, 
4495                               const double * lowerChange,
4496                               const double * upperChange)
4497{
4498     if (getFakeBound(iSequence) != noFake) {
4499          numberFake_--;
4500          setFakeBound(iSequence, noFake);
4501          if (iSequence >= numberColumns_) {
4502               // rows
4503               int iRow = iSequence - numberColumns_;
4504               rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
4505               rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
4506               if (rowScale_) {
4507                    if (rowLowerWork_[iRow] > -1.0e50)
4508                         rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
4509                    if (rowUpperWork_[iRow] < 1.0e50)
4510                         rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
4511               } else if (rhsScale_ != 1.0) {
4512                    if (rowLowerWork_[iRow] > -1.0e50)
4513                         rowLowerWork_[iRow] *= rhsScale_;
4514                    if (rowUpperWork_[iRow] < 1.0e50)
4515                         rowUpperWork_[iRow] *= rhsScale_;
4516               }
4517          } else {
4518               // columns
4519               columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
4520               columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
4521               if (rowScale_) {
4522                    double multiplier = 1.0 * inverseColumnScale_[iSequence];
4523                    if (columnLowerWork_[iSequence] > -1.0e50)
4524                         columnLowerWork_[iSequence] *= multiplier * rhsScale_;
4525                    if (columnUpperWork_[iSequence] < 1.0e50)
4526                         columnUpperWork_[iSequence] *= multiplier * rhsScale_;
4527               } else if (rhsScale_ != 1.0) {
4528                    if (columnLowerWork_[iSequence] > -1.0e50)
4529                         columnLowerWork_[iSequence] *= rhsScale_;
4530                    if (columnUpperWork_[iSequence] < 1.0e50)
4531                         columnUpperWork_[iSequence] *= rhsScale_;
4532               }
4533          }
4534     }
4535}
4536/* Expands out all possible combinations for a knapsack
4537   If buildObj NULL then just computes space needed - returns number elements
4538   On entry numberOutput is maximum allowed, on exit it is number needed or
4539   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
4540   least space to return values which reconstruct input.
4541   Rows returned will be original rows but no entries will be returned for
4542   any rows all of whose entries are in knapsack.  So up to user to allow for this.
4543   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
4544   in expanded knapsack.  Values in buildRow and buildElement;
4545*/
4546int
4547ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
4548                                double * buildObj, CoinBigIndex * buildStart,
4549                                int * buildRow, double * buildElement, int reConstruct) const
4550{
4551     int iRow;
4552     int iColumn;
4553     // Get column copy
4554     CoinPackedMatrix * columnCopy = matrix();
4555     // Get a row copy in standard format
4556     CoinPackedMatrix matrixByRow;
4557     matrixByRow.reverseOrderedCopyOf(*columnCopy);
4558     const double * elementByRow = matrixByRow.getElements();
4559     const int * column = matrixByRow.getIndices();
4560     const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
4561     const int * rowLength = matrixByRow.getVectorLengths();
4562     CoinBigIndex j;
4563     int * whichColumn = new int [numberColumns_];
4564     int * whichRow = new int [numberRows_];
4565     int numJ = 0;
4566     // Get what other columns can compensate for
4567     double * lo = new double [numberRows_];
4568     double * high = new double [numberRows_];
4569     {
4570          // Use to get tight column bounds
4571          ClpSimplex tempModel(*this);
4572          tempModel.tightenPrimalBounds(0.0, 0, true);
4573          // Now another model without knapsacks
4574          int nCol = 0;
4575          for (iRow = 0; iRow < numberRows_; iRow++) {
4576               whichRow[iRow] = iRow;
4577          }
4578          for (iColumn = 0; iColumn < numberColumns_; iColumn++)
4579               whichColumn[iColumn] = -1;
4580          for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4581               int iColumn = column[j];
4582               if (columnUpper_[iColumn] > columnLower_[iColumn]) {
4583                    whichColumn[iColumn] = 0;
4584               } else {
4585                    assert (!columnLower_[iColumn]); // fix later
4586               }
4587          }
4588          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4589               if (whichColumn[iColumn] < 0)
4590                    whichColumn[nCol++] = iColumn;
4591          }
4592          ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
4593          // Row copy
4594          CoinPackedMatrix matrixByRow;
4595          matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
4596          const double * elementByRow = matrixByRow.getElements();
4597          const int * column = matrixByRow.getIndices();
4598          const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
4599          const int * rowLength = matrixByRow.getVectorLengths();
4600          const double * columnLower = tempModel2.getColLower();
4601          const double * columnUpper = tempModel2.getColUpper();
4602          for (iRow = 0; iRow < numberRows_; iRow++) {
4603               lo[iRow] = -COIN_DBL_MAX;
4604               high[iRow] = COIN_DBL_MAX;
4605               if (rowLower_[iRow] > -1.0e20 || rowUpper_[iRow] < 1.0e20) {
4606
4607                    // possible row
4608                    int infiniteUpper = 0;
4609                    int infiniteLower = 0;
4610                    double maximumUp = 0.0;
4611                    double maximumDown = 0.0;
4612                    CoinBigIndex rStart = rowStart[iRow];
4613                    CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4614                    CoinBigIndex j;
4615                    // Compute possible lower and upper ranges
4616
4617                    for (j = rStart; j < rEnd; ++j) {
4618                         double value = elementByRow[j];
4619                         iColumn = column[j];
4620                         if (value > 0.0) {
4621                              if (columnUpper[iColumn] >= 1.0e20) {
4622                                   ++infiniteUpper;
4623                              } else {
4624                                   maximumUp += columnUpper[iColumn] * value;
4625                              }
4626                              if (columnLower[iColumn] <= -1.0e20) {
4627                                   ++infiniteLower;
4628                              } else {
4629                                   maximumDown += columnLower[iColumn] * value;
4630                              }
4631                         } else if (value < 0.0) {
4632                              if (columnUpper[iColumn] >= 1.0e20) {
4633                                   ++infiniteLower;
4634                              } else {
4635                                   maximumDown += columnUpper[iColumn] * value;
4636                              }
4637                              if (columnLower[iColumn] <= -1.0e20) {
4638                                   ++infiniteUpper;
4639                              } else {
4640                                   maximumUp += columnLower[iColumn] * value;
4641                              }
4642                         }
4643                    }
4644                    // Build in a margin of error
4645                    maximumUp += 1.0e-8 * fabs(maximumUp) + 1.0e-7;
4646                    maximumDown -= 1.0e-8 * fabs(maximumDown) + 1.0e-7;
4647                    // we want to save effective rhs
4648                    double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
4649                    double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
4650                    if (up == COIN_DBL_MAX || rowLower_[iRow] == -COIN_DBL_MAX) {
4651                         // However low we go it doesn't matter
4652                         lo[iRow] = -COIN_DBL_MAX;
4653                    } else {
4654                         // If we go below this then can not be feasible
4655                         lo[iRow] = rowLower_[iRow] - up;
4656                    }
4657                    if (down == -COIN_DBL_MAX || rowUpper_[iRow] == COIN_DBL_MAX) {
4658                         // However high we go it doesn't matter
4659                         high[iRow] = COIN_DBL_MAX;
4660                    } else {
4661                         // If we go above this then can not be feasible
4662                         high[iRow] = rowUpper_[iRow] - down;
4663                    }
4664               }
4665          }
4666     }
4667     numJ = 0;
4668     for (iColumn = 0; iColumn < numberColumns_; iColumn++)
4669          whichColumn[iColumn] = -1;
4670     int * markRow = new int [numberRows_];
4671     for (iRow = 0; iRow < numberRows_; iRow++)
4672          markRow[iRow] = 1;
4673     for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4674          int iColumn = column[j];
4675          if (columnUpper_[iColumn] > columnLower_[iColumn]) {
4676               whichColumn[iColumn] = numJ;
4677               numJ++;
4678          }
4679     }
4680     /* mark rows
4681        -n in knapsack and n other variables
4682        1 no entries
4683        n+1000 not involved in knapsack but n entries
4684        0 only in knapsack
4685     */
4686     for (iRow = 0; iRow < numberRows_; iRow++) {
4687          int type = 1;
4688          for (j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4689               int iColumn = column[j];
4690               if (whichColumn[iColumn] >= 0) {
4691                    if (type == 1) {
4692                         type = 0;
4693                    } else if (type > 0) {
4694                         assert (type > 1000);
4695                         type = -(type - 1000);
4696                    }
4697               } else if (type == 1) {
4698                    type = 1001;
4699               } else if (type < 0) {
4700                    type --;
4701               } else if (type == 0) {
4702                    type = -1;
4703               } else {
4704                    assert (type > 1000);
4705                    type++;
4706               }
4707          }
4708          markRow[iRow] = type;
4709          if (type < 0 && type > -30 && false)
4710               printf("markrow on row %d is %d\n", iRow, markRow[iRow]);
4711     }
4712     int * bound = new int [numberColumns_+1];
4713     int * stack = new int [numberColumns_+1];
4714     int * flip = new int [numberColumns_+1];
4715     double * offset = new double[numberColumns_+1];
4716     double * size = new double [numberColumns_+1];
4717     double * rhsOffset = new double[numberRows_];
4718     int * build = new int[numberColumns_];
4719     int maxNumber = numberOutput;
4720     numJ = 0;
4721     double minSize = rowLower_[knapsackRow];
4722     double maxSize = rowUpper_[knapsackRow];
4723     double knapsackOffset = 0.0;
4724     for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4725          int iColumn = column[j];
4726          double lowerColumn = columnLower_[iColumn];
4727          double upperColumn = columnUpper_[iColumn];
4728          if (lowerColumn == upperColumn)
4729               continue;
4730          double gap = upperColumn - lowerColumn;
4731          if (gap > 1.0e8)
4732               gap = 1.0e8;
4733          assert (fabs(floor(gap + 0.5) - gap) < 1.0e-5);
4734          whichColumn[numJ] = iColumn;
4735          bound[numJ] = static_cast<int> (gap);
4736          if (elementByRow[j] > 0.0) {
4737               flip[numJ] = 1;
4738               offset[numJ] = lowerColumn;
4739               size[numJ++] = elementByRow[j];
4740          } else {
4741               flip[numJ] = -1;
4742               offset[numJ] = upperColumn;
4743               size[numJ++] = -elementByRow[j];
4744               lowerColumn = upperColumn;
4745          }
4746          knapsackOffset += elementByRow[j] * lowerColumn;
4747     }
4748     int jRow;
4749     for (iRow = 0; iRow < numberRows_; iRow++)
4750          whichRow[iRow] = iRow;
4751     ClpSimplex smallModel(this, numberRows_, whichRow, numJ, whichColumn, true, true, true);
4752     // modify rhs to allow for nonzero lower bounds
4753     //double * rowLower = smallModel.rowLower();
4754     //double * rowUpper = smallModel.rowUpper();
4755     //const double * columnLower = smallModel.columnLower();
4756     //const double * columnUpper = smallModel.columnUpper();
4757     const CoinPackedMatrix * matrix = smallModel.matrix();
4758     const double * element = matrix->getElements();
4759     const int * row = matrix->getIndices();
4760     const CoinBigIndex * columnStart = matrix->getVectorStarts();
4761     const int * columnLength = matrix->getVectorLengths();
4762     const double * objective = smallModel.objective();
4763     //double objectiveOffset=0.0;
4764     // would use for fixed?
4765     CoinZeroN(rhsOffset, numberRows_);
4766     double * rowActivity = smallModel.primalRowSolution();
4767     CoinZeroN(rowActivity, numberRows_);
4768     maxSize -= knapsackOffset;
4769     minSize -= knapsackOffset;
4770     // now generate
4771     int i;
4772     int iStack = numJ;
4773     for (i = 0; i < numJ; i++) {
4774          stack[i] = 0;
4775     }
4776     double tooMuch = 10.0 * maxSize + 10000;
4777     stack[numJ] = 1;
4778     size[numJ] = tooMuch;
4779     bound[numJ] = 0;
4780     double sum = tooMuch;
4781     // allow for all zero being OK
4782     stack[numJ-1] = -1;
4783     sum -= size[numJ-1];
4784     numberOutput = 0;
4785     int nelCreate = 0;
4786     /* typeRun is - 0 for initial sizes
4787                     1 for build
4788          2 for reconstruct
4789     */
4790     int typeRun = buildObj ? 1 : 0;
4791     if (reConstruct >= 0) {
4792          assert (buildRow && buildElement);
4793          typeRun = 2;
4794     }
4795     if (typeRun == 1)
4796          buildStart[0] = 0;
4797     while (iStack >= 0) {
4798          if (sum >= minSize && sum <= maxSize) {
4799               double checkSize = 0.0;
4800               bool good = true;
4801               int nRow = 0;
4802               double obj = 0.0;
4803               CoinZeroN(rowActivity, numberRows_);
4804               for (iColumn = 0; iColumn < numJ; iColumn++) {
4805                    int iValue = stack[iColumn];
4806                    if (iValue > bound[iColumn]) {
4807                         good = false;
4808                         break;
4809                    } else {
4810                         double realValue = offset[iColumn] + flip[iColumn] * iValue;
4811                         if (realValue) {
4812                              obj += objective[iColumn] * realValue;
4813                              for (CoinBigIndex j = columnStart[iColumn];
4814                                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4815                                   double value = element[j] * realValue;
4816                                   int kRow = row[j];
4817                                   if (rowActivity[kRow]) {
4818                                        rowActivity[kRow] += value;
4819                                        if (!rowActivity[kRow])
4820                                             rowActivity[kRow] = 1.0e-100;
4821                                   } else {
4822                                        build[nRow++] = kRow;
4823                                        rowActivity[kRow] = value;
4824                                   }
4825                              }
4826                         }
4827                    }
4828               }
4829               if (good) {
4830                    for (jRow = 0; jRow < nRow; jRow++) {
4831                         int kRow = build[jRow];
4832                         double value = rowActivity[kRow];
4833                         if (value > high[kRow] || value < lo[kRow]) {
4834                              good = false;
4835                              break;
4836                         }
4837                    }
4838               }
4839               if (good) {
4840                    if (typeRun == 1) {
4841                         buildObj[numberOutput] = obj;
4842                         for (jRow = 0; jRow < nRow; jRow++) {
4843                              int kRow = build[jRow];
4844                              double value = rowActivity[kRow];
4845                              if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
4846                                   buildElement[nelCreate] = value;
4847                                   buildRow[nelCreate++] = kRow;
4848                              }
4849                         }
4850                         buildStart[numberOutput+1] = nelCreate;
4851                    } else if (!typeRun) {
4852                         for (jRow = 0; jRow < nRow; jRow++) {
4853                              int kRow = build[jRow];
4854                              double value = rowActivity[kRow];
4855                              if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
4856                                   nelCreate++;
4857                              }
4858                         }
4859                    }
4860                    if (typeRun == 2 && reConstruct == numberOutput) {
4861                         // build and exit
4862                         nelCreate = 0;
4863                         for (iColumn = 0; iColumn < numJ; iColumn++) {
4864                              int iValue = stack[iColumn];
4865                              double realValue = offset[iColumn] + flip[iColumn] * iValue;
4866                              if (realValue) {
4867                                   buildRow[nelCreate] = whichColumn[iColumn];
4868                                   buildElement[nelCreate++] = realValue;
4869                              }
4870                         }
4871                         numberOutput = 1;
4872                         for (i = 0; i < numJ; i++) {
4873                              bound[i] = 0;
4874                         }
4875                         break;
4876                    }
4877                    numberOutput++;
4878                    if (numberOutput > maxNumber) {
4879                         nelCreate = -numberOutput;
4880                         numberOutput = -1;
4881                         for (i = 0; i < numJ; i++) {
4882                              bound[i] = 0;
4883                         }
4884                         break;
4885                    } else if (typeRun == 1 && numberOutput == maxNumber) {
4886                         // On second run
4887                         for (i = 0; i < numJ; i++) {
4888                              bound[i] = 0;
4889                         }
4890                         break;
4891                    }
4892                    for (int j = 0; j < numJ; j++) {
4893                         checkSize += stack[j] * size[j];
4894                    }
4895                    assert (fabs(sum - checkSize) < 1.0e-3);
4896               }
4897               for (jRow = 0; jRow < nRow; jRow++) {
4898                    int kRow = build[jRow];
4899                    rowActivity[kRow] = 0.0;
4900               }
4901          }
4902          if (sum > maxSize || stack[iStack] > bound[iStack]) {
4903               sum -= size[iStack] * stack[iStack];
4904               stack[iStack--] = 0;
4905               if (iStack >= 0) {
4906                    stack[iStack] ++;
4907                    sum += size[iStack];
4908               }
4909          } else {
4910               // must be less
4911               // add to last possible
4912               iStack = numJ - 1;
4913               sum += size[iStack];
4914               stack[iStack]++;
4915          }
4916     }
4917     //printf("%d will be created\n",numberOutput);
4918     delete [] whichColumn;
4919     delete [] whichRow;
4920     delete [] bound;
4921     delete [] stack;
4922     delete [] flip;
4923     delete [] size;
4924     delete [] offset;
4925     delete [] rhsOffset;
4926     delete [] build;
4927     delete [] markRow;
4928     delete [] lo;
4929     delete [] high;
4930     return nelCreate;
4931}
4932// Quick try at cleaning up duals if postsolve gets wrong
4933void
4934ClpSimplexOther::cleanupAfterPostsolve()
4935{
4936     // First mark singleton equality rows
4937     char * mark = new char [ numberRows_];
4938     memset(mark, 0, numberRows_);
4939     const int * row = matrix_->getIndices();
4940     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4941     const int * columnLength = matrix_->getVectorLengths();
4942     const double * element = matrix_->getElements();
4943     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4944          for (CoinBigIndex j = columnStart[iColumn];
4945                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4946               int iRow = row[j];
4947               if (mark[iRow])
4948                    mark[iRow] = 2;
4949               else
4950                    mark[iRow] = 1;
4951          }
4952     }
4953     // for now just == rows
4954     for (int iRow = 0; iRow < numberRows_; iRow++) {
4955          if (rowUpper_[iRow] > rowLower_[iRow])
4956               mark[iRow] = 3;
4957     }
4958     double dualTolerance = dblParam_[ClpDualTolerance];
4959     double primalTolerance = dblParam_[ClpPrimalTolerance];
4960     int numberCleaned = 0;
4961     double maxmin = optimizationDirection_;
4962     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4963          double dualValue = reducedCost_[iColumn] * maxmin;
4964          double primalValue = columnActivity_[iColumn];
4965          double lower = columnLower_[iColumn];
4966          double upper = columnUpper_[iColumn];
4967          int way = 0;
4968          switch(getColumnStatus(iColumn)) {
4969
4970          case basic:
4971               // dual should be zero
4972               if (dualValue > dualTolerance) {
4973                    way = -1;
4974               } else if (dualValue < -dualTolerance) {
4975                    way = 1;
4976               }
4977               break;
4978          case ClpSimplex::isFixed:
4979               break;
4980          case atUpperBound:
4981               // dual should not be positive
4982               if (dualValue > dualTolerance) {
4983                    way = -1;
4984               }
4985               break;
4986          case atLowerBound:
4987               // dual should not be negative
4988               if (dualValue < -dualTolerance) {
4989                    way = 1;
4990               }
4991               break;
4992          case superBasic:
4993          case isFree:
4994               if (primalValue < upper - primalTolerance) {
4995                    // dual should not be negative
4996                    if (dualValue < -dualTolerance) {
4997                         way = 1;
4998                    }
4999               }
5000               if (primalValue > lower + primalTolerance) {
5001                    // dual should not be positive
5002                    if (dualValue > dualTolerance) {
5003                         way = -1;
5004                    }
5005               }
5006               break;
5007          }
5008          if (way) {
5009               // see if can find singleton row
5010               for (CoinBigIndex j = columnStart[iColumn];
5011                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
5012                    int iRow = row[j];
5013                    if (mark[iRow] == 1) {
5014                         double value = element[j];
5015                         // dj - addDual*value == 0.0
5016                         double addDual = dualValue / value;
5017                         dual_[iRow] += addDual;
5018                         reducedCost_[iColumn] = 0.0;
5019                         numberCleaned++;
5020                         break;
5021                    }
5022               }
5023          }
5024     }
5025     delete [] mark;
5026#ifdef CLP_INVESTIGATE
5027     printf("cleanupAfterPostsolve cleaned up %d columns\n", numberCleaned);
5028#endif
5029     // Redo
5030     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
5031     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
5032     checkSolutionInternal();
5033}
5034// Returns gub version of model or NULL
5035ClpSimplex * 
5036ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns,
5037                            int neededGub,
5038                            int factorizationFrequency)
5039{
5040  // find gub
5041  int numberRows = this->numberRows();
5042  int numberColumns = this->numberColumns();
5043  int iRow, iColumn;
5044  int * columnIsGub = new int [numberColumns];
5045  const double * columnLower = this->columnLower();
5046  const double * columnUpper = this->columnUpper();
5047  int numberFixed=0;
5048  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5049    if (columnUpper[iColumn] == columnLower[iColumn]) {
5050      columnIsGub[iColumn]=-2;
5051      numberFixed++;
5052    } else if (columnLower[iColumn]>=0) {
5053      columnIsGub[iColumn]=-1;
5054    } else {
5055      columnIsGub[iColumn]=-3;
5056    }
5057  }
5058  CoinPackedMatrix * matrix = this->matrix();
5059  // get row copy
5060  CoinPackedMatrix rowCopy = *matrix;
5061  rowCopy.reverseOrdering();
5062  const int * column = rowCopy.getIndices();
5063  const int * rowLength = rowCopy.getVectorLengths();
5064  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
5065  const double * element = rowCopy.getElements();
5066  int numberNonGub = 0;
5067  int numberEmpty = numberRows;
5068  int * rowIsGub = new int [numberRows];
5069  int smallestGubRow=-1;
5070  int count=numberColumns+1;
5071  double * rowLower = this->rowLower();
5072  double * rowUpper = this->rowUpper();
5073  // make sure we can get rid of upper bounds
5074  double * fixedRow = new double [numberRows];
5075  for (iRow = 0 ; iRow < numberRows ; iRow++) {
5076    double sumFixed=0.0;
5077    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5078      int iColumn = column[j];
5079      double value = columnLower[iColumn];
5080      if (value) 
5081        sumFixed += element[j] * value;
5082    }
5083    fixedRow[iRow]=rowUpper[iRow]-sumFixed;
5084  }
5085  for (iRow = numberRows - 1; iRow >= 0; iRow--) {
5086    bool gubRow = true;
5087    int numberInRow=0;
5088    double sumFixed=0.0;
5089    double gap = fixedRow[iRow]-1.0e-12;
5090    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5091      int iColumn = column[j];
5092      if (columnIsGub[iColumn]!=-2) {
5093        if (element[j] != 1.0||columnIsGub[iColumn]==-3||
5094            columnUpper[iColumn]-columnLower[iColumn]<gap) {
5095          gubRow = false;
5096          break;
5097        } else {
5098          numberInRow++;
5099          if (columnIsGub[iColumn] >= 0) {
5100            gubRow = false;
5101            break;
5102          }
5103        }
5104      } else {
5105        sumFixed += columnLower[iColumn]*element[j];
5106      }
5107    }
5108    if (!gubRow) {
5109      whichRows[numberNonGub++] = iRow;
5110      rowIsGub[iRow] = -1;
5111    } else if (numberInRow) {
5112      if (numberInRow<count) {
5113        count = numberInRow;
5114        smallestGubRow=iRow;
5115      }
5116      for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5117        int iColumn = column[j];
5118        if (columnIsGub[iColumn]!=-2) 
5119          columnIsGub[iColumn] = iRow;
5120      }
5121      rowIsGub[iRow] = 0;
5122    } else {
5123      // empty row!
5124      whichRows[--numberEmpty] = iRow;
5125      rowIsGub[iRow] = -2;
5126      if (sumFixed>rowUpper[iRow]+1.0e-4||
5127          sumFixed<rowLower[iRow]-1.0e-4) {
5128        fprintf(stderr,"******** No infeasible empty rows - please!\n");
5129        abort();
5130      }
5131    }
5132  }
5133  delete [] fixedRow;
5134  char message[100];
5135  int numberGub = numberEmpty - numberNonGub;