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

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

improve robustness of parametrics

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