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

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

As Clp trunk now seems to use CoinUtils? stable, I
have had to duplicate code from there for the moment.

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