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

Last change on this file since 1784 was 1784, checked in by forrest, 9 years ago

minor change for null objective

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