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

Last change on this file since 1780 was 1780, checked in by forrest, 10 years ago

stuff for parametrics

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