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

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

patches to try and make parametrics faster

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