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

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

stop seg fault on infinite theta

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 260.7 KB
RevLine 
[1370]1/* $Id: ClpSimplexOther.cpp 1795 2011-09-15 08:28:21Z forrest $ */
[343]2// Copyright (C) 2004, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1665]4// This code is licensed under the terms of the Eclipse Public License (EPL).
[343]5
6#include "CoinPragma.hpp"
7
8#include <math.h>
9
10#include "CoinHelperFunctions.hpp"
11#include "ClpSimplexOther.hpp"
[636]12#include "ClpSimplexDual.hpp"
13#include "ClpSimplexPrimal.hpp"
14#include "ClpEventHandler.hpp"
15#include "ClpHelperFunctions.hpp"
[343]16#include "ClpFactorization.hpp"
[636]17#include "ClpDualRowDantzig.hpp"
[1761]18#include "ClpNonLinearCost.hpp"
[1676]19#include "ClpDynamicMatrix.hpp"
[343]20#include "CoinPackedMatrix.hpp"
21#include "CoinIndexedVector.hpp"
[618]22#include "CoinBuild.hpp"
[480]23#include "CoinMpsIO.hpp"
[1722]24#include "CoinFloatEqual.hpp"
[343]25#include "ClpMessage.hpp"
26#include <cfloat>
27#include <cassert>
28#include <string>
29#include <stdio.h>
[1643]30#include <iostream>
[343]31/* Dual ranging.
32   This computes increase/decrease in cost for each given variable and corresponding
[1525]33   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
[343]34   and numberColumns.. for artificials/slacks.
35   For non-basic variables the sequence number will be that of the non-basic variables.
[1525]36
[343]37   Up to user to provide correct length arrays.
[1525]38
[343]39*/
[1525]40void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
41                                  double * costIncreased, int * sequenceIncreased,
42                                  double * costDecreased, int * sequenceDecreased,
43                                  double * valueIncrease, double * valueDecrease)
[343]44{
[1525]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;
[631]79          }
[1525]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];
[631]88          }
[1525]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);
[1738]123                    //valueIncrease[i] = scale2;
[1525]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");
[343]228}
[1525]229/*
[343]230   Row array has row part of pivot row
231   Column array has column part.
232   This is used in dual ranging
233*/
234void
[389]235ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
[631]236                                 CoinIndexedVector * columnArray,
237                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
238                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
[343]239{
[1525]240     double acceptablePivot = 1.0e-9;
241     double * work;
242     int number;
243     int * which;
244     int iSection;
[343]245
[1525]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;
[343]252
[1525]253     int addSequence;
[343]254
[1525]255     for (iSection = 0; iSection < 2; iSection++) {
[343]256
[1525]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          }
[343]269
[1525]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     }
[343]339}
[389]340/** Primal ranging.
341    This computes increase/decrease in value for each given variable and corresponding
[1525]342    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
[389]343    and numberColumns.. for artificials/slacks.
344    For basic variables the sequence number will be that of the basic variables.
[1525]345
[389]346    Up to user to provide correct length arrays.
[1525]347
[389]348    When here - guaranteed optimal
349*/
[1525]350void
351ClpSimplexOther::primalRanging(int numberCheck, const int * which,
352                               double * valueIncreased, int * sequenceIncreased,
353                               double * valueDecreased, int * sequenceDecreased)
[389]354{
[1525]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     }
[389]424}
[912]425// Returns new value of whichOther when whichIn enters basis
[1525]426double
[912]427ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
428{
[1525]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;
[1197]437
[1525]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;
[912]523}
[1525]524/*
[389]525   Row array has pivot column
526   This is used in primal ranging
527*/
[1525]528void
[389]529ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
[1525]530                                   int direction)
[389]531{
[1525]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();
[389]538
[1525]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     }
[389]568}
[480]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)
[1525]572
[480]573   Row and column names may be null.
574   formatType is
575   <ul>
576   <li> 0 - normal
[1525]577   <li> 1 - extra accuracy
[480]578   <li> 2 - IEEE hex (later)
579   </ul>
[1525]580
[480]581   Returns non-zero on I/O error
582
583   This is based on code contributed by Thorsten Koch
584*/
[1525]585int
[480]586ClpSimplexOther::writeBasis(const char *filename,
[1525]587                            bool writeValues,
588                            int formatType) const
[480]589{
[1525]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     }
[480]606
[1525]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;
[480]688}
689// Read a basis from the given filename
[1525]690int
[480]691ClpSimplexOther::readBasis(const char *fileName)
692{
[1525]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;
[480]746}
[1111]747/* Creates dual of a problem if looks plausible
748   (defaults will always create model)
[1525]749   fractionRowRanges is fraction of rows allowed to have ranges
750   fractionColumnRanges is fraction of columns allowed to have ranges
[1111]751*/
[1525]752ClpSimplex *
753ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
[618]754{
[1525]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);
[618]812
[1525]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;
[618]941}
942// Restores solution from dualized problem
[1034]943int
[1785]944ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
945                                 bool checkAccuracy)
[618]946{
[1525]947     int returnCode = 0;;
948     createStatus();
949     // Number of rows in dual problem was original number of columns
950     assert (numberColumns_ == dualProblem->numberRows());
951     // If slack on d-row basic then column at bound otherwise column basic
952     // If d-column basic then rhs tight
953     int numberBasic = 0;
954     int iRow, iColumn = 0;
955     // Get number of extra rows from ranges
956     int numberExtraRows = 0;
957     for (iRow = 0; iRow < numberRows_; iRow++) {
958          if (rowLower_[iRow] > -1.0e20 &&
959                    rowUpper_[iRow] < 1.0e20) {
960               if (rowUpper_[iRow] != rowLower_[iRow])
961                    numberExtraRows++;
962          }
963     }
964     const double * objective = this->objective();
965     const double * dualDual = dualProblem->dualRowSolution();
966     const double * dualDj = dualProblem->dualColumnSolution();
967     const double * dualSol = dualProblem->primalColumnSolution();
968     const double * dualActs = dualProblem->primalRowSolution();
[618]969#if 0
[1525]970     ClpSimplex thisCopy = *this;
971     thisCopy.dual(); // for testing
972     const double * primalDual = thisCopy.dualRowSolution();
973     const double * primalDj = thisCopy.dualColumnSolution();
974     const double * primalSol = thisCopy.primalColumnSolution();
975     const double * primalActs = thisCopy.primalRowSolution();
976     char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
977     printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
978     for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
979          printf("%d at %c primal %g dual %g\n",
980                 iRow, ss[dualProblem->getRowStatus(iRow)],
981                 dualActs[iRow], dualDual[iRow]);
982     printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
983     for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
984          printf("%d at %c primal %g dual %g\n",
985                 iColumn, ss[dualProblem->getColumnStatus(iColumn)],
986                 dualSol[iColumn], dualDj[iColumn]);
987     printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
988     for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
989          printf("%d at %c primal %g dual %g\n",
990                 iRow, ss[thisCopy.getRowStatus(iRow)],
991                 primalActs[iRow], primalDual[iRow]);
992     printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
993     for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
994          printf("%d at %c primal %g dual %g\n",
995                 iColumn, ss[thisCopy.getColumnStatus(iColumn)],
996                 primalSol[iColumn], primalDj[iColumn]);
[618]997#endif
[1525]998     // position at bound information
999     int jColumn = numberRows_;
1000     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1001          double objValue = optimizationDirection_ * objective[iColumn];
1002          Status status = dualProblem->getRowStatus(iColumn);
1003          double otherValue = COIN_DBL_MAX;
1004          if (columnUpper_[iColumn] < 1.0e20 &&
1005                    columnLower_[iColumn] > -1.0e20) {
1006               if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1007                    otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1008               } else {
1009                    otherValue = columnLower_[iColumn] + dualDj[jColumn];
1010               }
1011               jColumn++;
1012          }
1013          if (status == basic) {
1014               // column is at bound
1015               if (otherValue == COIN_DBL_MAX) {
1016                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1017                    if (columnUpper_[iColumn] > 1.0e20) {
1018                         if (columnLower_[iColumn] > -1.0e20) {
1019                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1020                                   setColumnStatus(iColumn, atLowerBound);
1021                              else
1022                                   setColumnStatus(iColumn, isFixed);
1023                              columnActivity_[iColumn] = columnLower_[iColumn];
1024                         } else {
1025                              // free
1026                              setColumnStatus(iColumn, isFree);
1027                              columnActivity_[iColumn] = 0.0;
1028                         }
1029                    } else {
1030                         setColumnStatus(iColumn, atUpperBound);
1031                         columnActivity_[iColumn] = columnUpper_[iColumn];
1032                    }
1033               } else {
1034                    reducedCost_[iColumn] = objValue - dualActs[iColumn];
1035                    //printf("other dual sol %g\n",otherValue);
1036                    if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1037                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1038                              setColumnStatus(iColumn, atLowerBound);
1039                         else
1040                              setColumnStatus(iColumn, isFixed);
1041                         columnActivity_[iColumn] = columnLower_[iColumn];
1042                    } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1043                         if (columnUpper_[iColumn] > columnLower_[iColumn])
1044                              setColumnStatus(iColumn, atUpperBound);
1045                         else
1046                              setColumnStatus(iColumn, isFixed);
1047                         columnActivity_[iColumn] = columnUpper_[iColumn];
1048                    } else {
1049                         abort();
1050                    }
1051               }
[618]1052          } else {
[1525]1053               if (otherValue == COIN_DBL_MAX) {
1054                    // column basic
1055                    setColumnStatus(iColumn, basic);
1056                    numberBasic++;
1057                    if (columnLower_[iColumn] > -1.0e20) {
1058                         columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1059                    } else if (columnUpper_[iColumn] < 1.0e20) {
1060                         columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1061                    } else {
1062                         columnActivity_[iColumn] = -dualDual[iColumn];
1063                    }
1064                    reducedCost_[iColumn] = 0.0;
1065               } else {
1066                    // may be at other bound
1067                    //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1068                    if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1069                         // column basic
1070                         setColumnStatus(iColumn, basic);
1071                         numberBasic++;
1072                         //printf("Col %d otherV %g dualDual %g\n",iColumn,
1073                         // otherValue,dualDual[iColumn]);
1074                         columnActivity_[iColumn] = -dualDual[iColumn];
1075                         columnActivity_[iColumn] = otherValue;
1076                         reducedCost_[iColumn] = 0.0;
1077                    } else {
1078                         reducedCost_[iColumn] = objValue - dualActs[iColumn];
1079                         if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1080                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1081                                   setColumnStatus(iColumn, atLowerBound);
1082                              else
1083                                   setColumnStatus(iColumn, isFixed);
1084                              columnActivity_[iColumn] = columnLower_[iColumn];
1085                         } else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1086                              if (columnUpper_[iColumn] > columnLower_[iColumn])
1087                                   setColumnStatus(iColumn, atUpperBound);
1088                              else
1089                                   setColumnStatus(iColumn, isFixed);
1090                              columnActivity_[iColumn] = columnUpper_[iColumn];
1091                         } else {
1092                              abort();
1093                         }
1094                    }
1095               }
[618]1096          }
[1525]1097     }
1098     // now rows
1099     int kExtraRow = jColumn;
1100     int numberRanges = 0;
1101     for (iRow = 0; iRow < numberRows_; iRow++) {
1102          Status status = dualProblem->getColumnStatus(iRow);
1103          if (status == basic) {
1104               // row is at bound
1105               dual_[iRow] = dualSol[iRow];;
1106          } else {
1107               // row basic
1108               setRowStatus(iRow, basic);
1109               numberBasic++;
1110               dual_[iRow] = 0.0;
1111          }
1112          if (rowLower_[iRow] < -1.0e20) {
1113               if (status == basic) {
1114                    rowActivity_[iRow] = rowUpper_[iRow];
1115                    setRowStatus(iRow, atUpperBound);
1116               } else {
1117                    assert (dualDj[iRow] < 1.0e-5);
1118                    rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1119               }
1120          } else if (rowUpper_[iRow] > 1.0e20) {
1121               if (status == basic) {
1122                    rowActivity_[iRow] = rowLower_[iRow];
1123                    setRowStatus(iRow, atLowerBound);
1124               } else {
1125                    rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1126                    assert (dualDj[iRow] > -1.0e-5);
1127               }
1128          } else {
1129               if (rowUpper_[iRow] == rowLower_[iRow]) {
1130                    rowActivity_[iRow] = rowLower_[iRow];
1131                    if (status == basic) {
1132                         setRowStatus(iRow, isFixed);
1133                    }
1134               } else {
1135                    // range
1136                    numberRanges++;
1137                    Status statusL = dualProblem->getColumnStatus(kExtraRow);
1138                    //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1139                    //     iRow,status,kExtraRow,statusL, dualSol[iRow],
1140                    //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1141                    if (status == basic) {
1142                         assert (statusL != basic);
1143                         rowActivity_[iRow] = rowUpper_[iRow];
1144                         setRowStatus(iRow, atUpperBound);
1145                    } else if (statusL == basic) {
1146                         numberBasic--; // already counted
1147                         rowActivity_[iRow] = rowLower_[iRow];
1148                         setRowStatus(iRow, atLowerBound);
1149                         dual_[iRow] = dualSol[kExtraRow];;
1150                    } else {
1151                         rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1152                         assert (dualDj[iRow] < 1.0e-5);
1153                         // row basic
1154                         //setRowStatus(iRow,basic);
1155                         //numberBasic++;
1156                         dual_[iRow] = 0.0;
1157                    }
1158                    kExtraRow++;
1159               }
1160          }
1161     }
1162     if (numberBasic != numberRows_) {
1163          printf("Bad basis - ranges - coding needed\n");
1164          assert (numberRanges);
1165          abort();
1166     }
1167     if (optimizationDirection_ < 0.0) {
1168          for (iRow = 0; iRow < numberRows_; iRow++) {
1169               dual_[iRow] = -dual_[iRow];
1170          }
1171     }
1172     // redo row activities
1173     memset(rowActivity_, 0, numberRows_ * sizeof(double));
1174     matrix_->times(1.0, columnActivity_, rowActivity_);
1175     // redo reduced costs
1176     memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1177     matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1178     checkSolutionInternal();
1179     if (sumDualInfeasibilities_ > 1.0e-5 || sumPrimalInfeasibilities_ > 1.0e-5) {
1180          returnCode = 1;
[1344]1181#ifdef CLP_INVESTIGATE
[1525]1182          printf("There are %d dual infeasibilities summing to %g ",
1183                 numberDualInfeasibilities_, sumDualInfeasibilities_);
1184          printf("and %d primal infeasibilities summing to %g\n",
1185                 numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
[1344]1186#endif
[1525]1187     }
1188     // Below will go to ..DEBUG later
[1344]1189#if 1 //ndef NDEBUG
[1785]1190     if (checkAccuracy) {
1191       // Check if correct
1192       double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1193       double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1194       double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1195       double * dual = CoinCopyOfArray(dual_, numberRows_);
1196       this->dual(); //primal();
1197       CoinRelFltEq eq(1.0e-5);
1198       for (iRow = 0; iRow < numberRows_; iRow++) {
1199         assert(eq(dual[iRow], dual_[iRow]));
1200       }
1201       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1202         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1203       }
1204       for (iRow = 0; iRow < numberRows_; iRow++) {
1205         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1206       }
1207       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1208         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1209       }
1210       delete [] columnActivity;
1211       delete [] rowActivity;
1212       delete [] reducedCost;
1213       delete [] dual;
[1525]1214     }
[1330]1215#endif
[1525]1216     return returnCode;
[618]1217}
1218/* Does very cursory presolve.
1219   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1220*/
[1525]1221ClpSimplex *
[618]1222ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
[636]1223                        int & nBound, bool moreBounds, bool tightenBounds)
[618]1224{
[1525]1225     //#define CHECK_STATUS
[1429]1226#ifdef CHECK_STATUS
[1525]1227     {
1228          int n = 0;
1229          int i;
1230          for (i = 0; i < numberColumns_; i++)
1231               if (getColumnStatus(i) == ClpSimplex::basic)
1232                    n++;
1233          for (i = 0; i < numberRows_; i++)
1234               if (getRowStatus(i) == ClpSimplex::basic)
1235                    n++;
1236          assert (n == numberRows_);
1237     }
[618]1238#endif
1239
[1525]1240     const double * element = matrix_->getElements();
1241     const int * row = matrix_->getIndices();
1242     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1243     const int * columnLength = matrix_->getVectorLengths();
[689]1244
[1525]1245     CoinZeroN(rhs, numberRows_);
1246     int iColumn;
1247     int iRow;
1248     CoinZeroN(whichRow, numberRows_);
1249     int * backColumn = whichColumn + numberColumns_;
1250     int numberRows2 = 0;
1251     int numberColumns2 = 0;
1252     double offset = 0.0;
1253     const double * objective = this->objective();
1254     double * solution = columnActivity_;
1255     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1256          double lower = columnLower_[iColumn];
1257          double upper = columnUpper_[iColumn];
1258          if (upper > lower || getColumnStatus(iColumn) == ClpSimplex::basic) {
1259               backColumn[iColumn] = numberColumns2;
1260               whichColumn[numberColumns2++] = iColumn;
1261               for (CoinBigIndex j = columnStart[iColumn];
1262                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1263                    int iRow = row[j];
1264                    int n = whichRow[iRow];
1265                    if (n == 0 && element[j])
1266                         whichRow[iRow] = -iColumn - 1;
1267                    else if (n < 0)
1268                         whichRow[iRow] = 2;
1269               }
[618]1270          } else {
[1525]1271               // fixed
1272               backColumn[iColumn] = -1;
1273               solution[iColumn] = upper;
1274               if (upper) {
1275                    offset += objective[iColumn] * upper;
1276                    for (CoinBigIndex j = columnStart[iColumn];
1277                              j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1278                         int iRow = row[j];
1279                         double value = element[j];
1280                         rhs[iRow] += upper * value;
1281                    }
1282               }
[618]1283          }
[1525]1284     }
1285     int returnCode = 0;
1286     double tolerance = primalTolerance();
1287     nBound = 2 * numberRows_;
1288     for (iRow = 0; iRow < numberRows_; iRow++) {
1289          int n = whichRow[iRow];
1290          if (n > 0) {
1291               whichRow[numberRows2++] = iRow;
1292          } else if (n < 0) {
1293               //whichRow[numberRows2++]=iRow;
1294               //continue;
1295               // Can only do in certain circumstances as we don't know current value
1296               if (rowLower_[iRow] == rowUpper_[iRow] || getRowStatus(iRow) == ClpSimplex::basic) {
1297                    // save row and column for bound
1298                    whichRow[--nBound] = iRow;
1299                    whichRow[nBound+numberRows_] = -n - 1;
1300               } else if (moreBounds) {
1301                    // save row and column for bound
1302                    whichRow[--nBound] = iRow;
1303                    whichRow[nBound+numberRows_] = -n - 1;
1304               } else {
1305                    whichRow[numberRows2++] = iRow;
1306               }
[636]1307          } else {
[1525]1308               // empty
1309               double rhsValue = rhs[iRow];
1310               if (rhsValue < rowLower_[iRow] - tolerance || rhsValue > rowUpper_[iRow] + tolerance) {
1311                    returnCode = 1; // infeasible
1312               }
[636]1313          }
[1525]1314     }
1315     ClpSimplex * small = NULL;
1316     if (!returnCode) {
[1660]1317       //printf("CRUNCH from (%d,%d) to (%d,%d)\n",
1318       //     numberRows_,numberColumns_,numberRows2,numberColumns2);
[1525]1319          small = new ClpSimplex(this, numberRows2, whichRow,
1320                                 numberColumns2, whichColumn, true, false);
1321#if 0
1322          ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1323          if (rowCopy) {
1324               assert(!small->rowCopy());
1325               small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1326                                    numberColumns2, whichColumn));
1327          }
[636]1328#endif
[1525]1329          // Set some stuff
1330          small->setDualBound(dualBound_);
1331          small->setInfeasibilityCost(infeasibilityCost_);
1332          small->setSpecialOptions(specialOptions_);
1333          small->setPerturbation(perturbation_);
1334          small->defaultFactorizationFrequency();
1335          small->setAlphaAccuracy(alphaAccuracy_);
1336          // If no rows left then no tightening!
1337          if (!numberRows2 || !numberColumns2)
1338               tightenBounds = false;
1339
1340          int numberElements = getNumElements();
1341          int numberElements2 = small->getNumElements();
1342          small->setObjectiveOffset(objectiveOffset() - offset);
1343          handler_->message(CLP_CRUNCH_STATS, messages_)
1344                    << numberRows2 << -(numberRows_ - numberRows2)
1345                    << numberColumns2 << -(numberColumns_ - numberColumns2)
1346                    << numberElements2 << -(numberElements - numberElements2)
1347                    << CoinMessageEol;
1348          // And set objective value to match
1349          small->setObjectiveValue(this->objectiveValue());
1350          double * rowLower2 = small->rowLower();
1351          double * rowUpper2 = small->rowUpper();
1352          int jRow;
1353          for (jRow = 0; jRow < numberRows2; jRow++) {
1354               iRow = whichRow[jRow];
1355               if (rowLower2[jRow] > -1.0e20)
1356                    rowLower2[jRow] -= rhs[iRow];
1357               if (rowUpper2[jRow] < 1.0e20)
1358                    rowUpper2[jRow] -= rhs[iRow];
1359          }
1360          // and bounds
1361          double * columnLower2 = small->columnLower();
1362          double * columnUpper2 = small->columnUpper();
1363          const char * integerInformation = integerType_;
1364          for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1365               iRow = whichRow[jRow];
1366               iColumn = whichRow[jRow+numberRows_];
1367               double lowerRow = rowLower_[iRow];
1368               if (lowerRow > -1.0e20)
1369                    lowerRow -= rhs[iRow];
1370               double upperRow = rowUpper_[iRow];
1371               if (upperRow < 1.0e20)
1372                    upperRow -= rhs[iRow];
1373               int jColumn = backColumn[iColumn];
1374               double lower = columnLower2[jColumn];
1375               double upper = columnUpper2[jColumn];
1376               double value = 0.0;
1377               for (CoinBigIndex j = columnStart[iColumn];
1378                         j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1379                    if (iRow == row[j]) {
1380                         value = element[j];
1381                         break;
[636]1382                    }
[1525]1383               }
1384               assert (value);
1385               // convert rowLower and Upper to implied bounds on column
1386               double newLower = -COIN_DBL_MAX;
1387               double newUpper = COIN_DBL_MAX;
1388               if (value > 0.0) {
1389                    if (lowerRow > -1.0e20)
1390                         newLower = lowerRow / value;
1391                    if (upperRow < 1.0e20)
1392                         newUpper = upperRow / value;
1393               } else {
1394                    if (upperRow < 1.0e20)
1395                         newLower = upperRow / value;
1396                    if (lowerRow > -1.0e20)
1397                         newUpper = lowerRow / value;
1398               }
1399               if (integerInformation && integerInformation[iColumn]) {
1400                    if (newLower - floor(newLower) < 10.0 * tolerance)
1401                         newLower = floor(newLower);
1402                    else
1403                         newLower = ceil(newLower);
1404                    if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1405                         newUpper = ceil(newUpper);
1406                    else
1407                         newUpper = floor(newUpper);
1408               }
1409               newLower = CoinMax(lower, newLower);
1410               newUpper = CoinMin(upper, newUpper);
1411               if (newLower > newUpper + tolerance) {
1412                    //printf("XXYY inf on bound\n");
1413                    returnCode = 1;
1414               }
1415               columnLower2[jColumn] = newLower;
1416               columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1417               if (getRowStatus(iRow) != ClpSimplex::basic) {
1418                    if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1419                         if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1420                              // can only get here if will be fixed
1421                              small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1422                         } else {
1423                              // solution is valid
1424                              if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1425                                        fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1426                                   small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1427                              else
1428                                   small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1429                         }
1430                    } else {
1431                         //printf("what now neither basic\n");
1432                    }
1433               }
1434          }
1435          if (returnCode) {
1436               delete small;
1437               small = NULL;
1438          } else if (tightenBounds && integerInformation) {
1439               // See if we can tighten any bounds
1440               // use rhs for upper and small duals for lower
1441               double * up = rhs;
1442               double * lo = small->dualRowSolution();
1443               const double * element = small->clpMatrix()->getElements();
1444               const int * row = small->clpMatrix()->getIndices();
1445               const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1446               //const int * columnLength = small->clpMatrix()->getVectorLengths();
1447               CoinZeroN(lo, numberRows2);
1448               CoinZeroN(up, numberRows2);
1449               for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1450                    double upper = columnUpper2[iColumn];
1451                    double lower = columnLower2[iColumn];
1452                    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1453                    for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1454                         int iRow = row[j];
1455                         double value = element[j];
1456                         if (value > 0.0) {
1457                              if (upper < 1.0e20)
1458                                   up[iRow] += upper * value;
1459                              else
1460                                   up[iRow] = COIN_DBL_MAX;
1461                              if (lower > -1.0e20)
1462                                   lo[iRow] += lower * value;
1463                              else
1464                                   lo[iRow] = -COIN_DBL_MAX;
1465                         } else {
1466                              if (upper < 1.0e20)
1467                                   lo[iRow] += upper * value;
1468                              else
1469                                   lo[iRow] = -COIN_DBL_MAX;
1470                              if (lower > -1.0e20)
1471                                   up[iRow] += lower * value;
1472                              else
1473                                   up[iRow] = COIN_DBL_MAX;
1474                         }
1475                    }
1476               }
1477               double * rowLower2 = small->rowLower();
1478               double * rowUpper2 = small->rowUpper();
1479               bool feasible = true;
1480               // make safer
1481               for (int iRow = 0; iRow < numberRows2; iRow++) {
1482                    double lower = lo[iRow];
1483                    if (lower > rowUpper2[iRow] + tolerance) {
1484                         feasible = false;
1485                         break;
1486                    } else {
1487                         lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1488                    }
1489                    double upper = up[iRow];
1490                    if (upper < rowLower2[iRow] - tolerance) {
1491                         feasible = false;
1492                         break;
1493                    } else {
1494                         up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1495                    }
1496               }
1497               if (!feasible) {
1498                    delete small;
1499                    small = NULL;
1500               } else {
1501                    // and tighten
1502                    for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1503                         if (integerInformation[whichColumn[iColumn]]) {
1504                              double upper = columnUpper2[iColumn];
1505                              double lower = columnLower2[iColumn];
1506                              double newUpper = upper;
1507                              double newLower = lower;
1508                              double difference = upper - lower;
1509                              if (lower > -1000.0 && upper < 1000.0) {
1510                                   for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1511                                        int iRow = row[j];
1512                                        double value = element[j];
1513                                        if (value > 0.0) {
1514                                             double upWithOut = up[iRow] - value * difference;
1515                                             if (upWithOut < 0.0) {
1516                                                  newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1517                                             }
1518                                             double lowWithOut = lo[iRow] + value * difference;
1519                                             if (lowWithOut > 0.0) {
1520                                                  newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1521                                             }
1522                                        } else {
1523                                             double upWithOut = up[iRow] + value * difference;
1524                                             if (upWithOut < 0.0) {
1525                                                  newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1526                                             }
1527                                             double lowWithOut = lo[iRow] - value * difference;
1528                                             if (lowWithOut > 0.0) {
1529                                                  newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1530                                             }
1531                                        }
1532                                   }
1533                                   if (newLower > lower || newUpper < upper) {
1534                                        if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1535                                             newUpper = floor(newUpper);
1536                                        else
1537                                             newUpper = floor(newUpper + 0.5);
1538                                        if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1539                                             newLower = ceil(newLower);
1540                                        else
1541                                             newLower = ceil(newLower - 0.5);
1542                                        // change may be too small - check
1543                                        if (newLower > lower || newUpper < upper) {
1544                                             if (newUpper >= newLower) {
1545                                                  // Could also tighten in this
1546                                                  //printf("%d bounds %g %g tightened to %g %g\n",
1547                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1548                                                  //     newLower,newUpper);
[636]1549#if 1
[1525]1550                                                  columnUpper2[iColumn] = newUpper;
1551                                                  columnLower2[iColumn] = newLower;
1552                                                  columnUpper_[whichColumn[iColumn]] = newUpper;
1553                                                  columnLower_[whichColumn[iColumn]] = newLower;
[636]1554#endif
[1525]1555                                                  // and adjust bounds on rows
1556                                                  newUpper -= upper;
1557                                                  newLower -= lower;
1558                                                  for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1559                                                       int iRow = row[j];
1560                                                       double value = element[j];
1561                                                       if (value > 0.0) {
1562                                                            up[iRow] += newUpper * value;
1563                                                            lo[iRow] += newLower * value;
1564                                                       } else {
1565                                                            lo[iRow] += newUpper * value;
1566                                                            up[iRow] += newLower * value;
1567                                                       }
1568                                                  }
1569                                             } else {
1570                                                  // infeasible
1571                                                  //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1572                                                  //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1573                                                  //     newLower,newUpper);
1574#if 1
1575                                                  delete small;
1576                                                  small = NULL;
1577                                                  break;
1578#endif
1579                                             }
1580                                        }
1581                                   }
1582                              }
1583                         }
1584                    }
1585               }
[636]1586          }
[1525]1587     }
[1344]1588#if 0
[1525]1589     if (small) {
1590          static int which = 0;
1591          which++;
1592          char xxxx[20];
1593          sprintf(xxxx, "bad%d.mps", which);
1594          small->writeMps(xxxx, 0, 1);
1595          sprintf(xxxx, "largebad%d.mps", which);
1596          writeMps(xxxx, 0, 1);
1597          printf("bad%d %x old size %d %d new %d %d\n", which, small,
1598                 numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
[1344]1599#if 0
[1525]1600          for (int i = 0; i < numberColumns_; i++)
1601               printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1602          for (int i = 0; i < numberRows_; i++)
1603               printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
[1344]1604#endif
[1525]1605     }
[1344]1606#endif
[1429]1607#ifdef CHECK_STATUS
[1525]1608     {
1609          int n = 0;
1610          int i;
1611          for (i = 0; i < small->numberColumns(); i++)
1612               if (small->getColumnStatus(i) == ClpSimplex::basic)
1613                    n++;
1614          for (i = 0; i < small->numberRows(); i++)
1615               if (small->getRowStatus(i) == ClpSimplex::basic)
1616                    n++;
1617          assert (n == small->numberRows());
1618     }
[1429]1619#endif
[1525]1620     return small;
[618]1621}
1622/* After very cursory presolve.
1623   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1624*/
[1525]1625void
[618]1626ClpSimplexOther::afterCrunch(const ClpSimplex & small,
[1525]1627                             const int * whichRow,
[618]1628                             const int * whichColumn, int nBound)
1629{
[1344]1630#ifndef NDEBUG
[1525]1631     for (int i = 0; i < small.numberRows(); i++)
1632          assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1633     for (int i = 0; i < small.numberColumns(); i++)
1634          assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
[1344]1635#endif
[1525]1636     getbackSolution(small, whichRow, whichColumn);
1637     // and deal with status for bounds
1638     const double * element = matrix_->getElements();
1639     const int * row = matrix_->getIndices();
1640     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1641     const int * columnLength = matrix_->getVectorLengths();
1642     double tolerance = primalTolerance();
1643     double djTolerance = dualTolerance();
1644     for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1645          int iRow = whichRow[jRow];
1646          int iColumn = whichRow[jRow+numberRows_];
1647          if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1648               double lower = columnLower_[iColumn];
1649               double upper = columnUpper_[iColumn];
1650               double value = columnActivity_[iColumn];
1651               double djValue = reducedCost_[iColumn];
1652               dual_[iRow] = 0.0;
1653               if (upper > lower) {
1654                    if (value < lower + tolerance && djValue > -djTolerance) {
1655                         setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1656                         setRowStatus(iRow, ClpSimplex::basic);
1657                    } else if (value > upper - tolerance && djValue < djTolerance) {
1658                         setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1659                         setRowStatus(iRow, ClpSimplex::basic);
1660                    } else {
1661                         // has to be basic
1662                         setColumnStatus(iColumn, ClpSimplex::basic);
1663                         reducedCost_[iColumn] = 0.0;
1664                         double value = 0.0;
1665                         for (CoinBigIndex j = columnStart[iColumn];
1666                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1667                              if (iRow == row[j]) {
1668                                   value = element[j];
1669                                   break;
1670                              }
1671                         }
1672                         dual_[iRow] = djValue / value;
1673                         if (rowUpper_[iRow] > rowLower_[iRow]) {
1674                              if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1675                                        fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1676                                   setRowStatus(iRow, ClpSimplex::atLowerBound);
1677                              else
1678                                   setRowStatus(iRow, ClpSimplex::atUpperBound);
1679                         } else {
1680                              setRowStatus(iRow, ClpSimplex::isFixed);
1681                         }
1682                    }
1683               } else {
1684                    // row can always be basic
1685                    setRowStatus(iRow, ClpSimplex::basic);
1686               }
[618]1687          } else {
[1525]1688               // row can always be basic
1689               setRowStatus(iRow, ClpSimplex::basic);
[618]1690          }
[1525]1691     }
1692     //#ifndef NDEBUG
[618]1693#if 0
[1525]1694     if  (small.status() == 0) {
1695          int n = 0;
1696          int i;
1697          for (i = 0; i < numberColumns; i++)
1698               if (getColumnStatus(i) == ClpSimplex::basic)
1699                    n++;
1700          for (i = 0; i < numberRows; i++)
1701               if (getRowStatus(i) == ClpSimplex::basic)
1702                    n++;
1703          assert (n == numberRows);
1704     }
[618]1705#endif
1706}
[636]1707/* Tightens integer bounds - returns number tightened or -1 if infeasible
1708 */
[1525]1709int
[636]1710ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1711{
[1525]1712     // See if we can tighten any bounds
1713     // use rhs for upper and small duals for lower
1714     double * up = rhsSpace;
1715     double * lo = dual_;
1716     const double * element = matrix_->getElements();
1717     const int * row = matrix_->getIndices();
1718     const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1719     const int * columnLength = matrix_->getVectorLengths();
1720     CoinZeroN(lo, numberRows_);
1721     CoinZeroN(up, numberRows_);
1722     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1723          double upper = columnUpper_[iColumn];
1724          double lower = columnLower_[iColumn];
1725          //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1726          for (CoinBigIndex j = columnStart[iColumn];
1727                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1728               int iRow = row[j];
1729               double value = element[j];
1730               if (value > 0.0) {
1731                    if (upper < 1.0e20)
1732                         up[iRow] += upper * value;
1733                    else
1734                         up[iRow] = COIN_DBL_MAX;
1735                    if (lower > -1.0e20)
1736                         lo[iRow] += lower * value;
1737                    else
1738                         lo[iRow] = -COIN_DBL_MAX;
1739               } else {
1740                    if (upper < 1.0e20)
1741                         lo[iRow] += upper * value;
1742                    else
1743                         lo[iRow] = -COIN_DBL_MAX;
1744                    if (lower > -1.0e20)
1745                         up[iRow] += lower * value;
1746                    else
1747                         up[iRow] = COIN_DBL_MAX;
1748               }
[636]1749          }
[1525]1750     }
1751     bool feasible = true;
1752     // make safer
1753     double tolerance = primalTolerance();
1754     for (int iRow = 0; iRow < numberRows_; iRow++) {
1755          double lower = lo[iRow];
1756          if (lower > rowUpper_[iRow] + tolerance) {
1757               feasible = false;
1758               break;
1759          } else {
1760               lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
[636]1761          }
[1525]1762          double upper = up[iRow];
1763          if (upper < rowLower_[iRow] - tolerance) {
1764               feasible = false;
1765               break;
1766          } else {
1767               up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1768          }
1769     }
1770     int numberTightened = 0;
1771     if (!feasible) {
1772          return -1;
1773     } else if (integerType_) {
1774          // and tighten
1775          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1776               if (integerType_[iColumn]) {
1777                    double upper = columnUpper_[iColumn];
1778                    double lower = columnLower_[iColumn];
1779                    double newUpper = upper;
1780                    double newLower = lower;
1781                    double difference = upper - lower;
1782                    if (lower > -1000.0 && upper < 1000.0) {
1783                         for (CoinBigIndex j = columnStart[iColumn];
1784                                   j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1785                              int iRow = row[j];
1786                              double value = element[j];
1787                              if (value > 0.0) {
1788                                   double upWithOut = up[iRow] - value * difference;
1789                                   if (upWithOut < 0.0) {
1790                                        newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1791                                   }
1792                                   double lowWithOut = lo[iRow] + value * difference;
1793                                   if (lowWithOut > 0.0) {
1794                                        newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1795                                   }
1796                              } else {
1797                                   double upWithOut = up[iRow] + value * difference;
1798                                   if (upWithOut < 0.0) {
1799                                        newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1800                                   }
1801                                   double lowWithOut = lo[iRow] - value * difference;
1802                                   if (lowWithOut > 0.0) {
1803                                        newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1804                                   }
1805                              }
1806                         }
1807                         if (newLower > lower || newUpper < upper) {
1808                              if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1809                                   newUpper = floor(newUpper);
1810                              else
1811                                   newUpper = floor(newUpper + 0.5);
1812                              if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1813                                   newLower = ceil(newLower);
1814                              else
1815                                   newLower = ceil(newLower - 0.5);
1816                              // change may be too small - check
1817                              if (newLower > lower || newUpper < upper) {
1818                                   if (newUpper >= newLower) {
1819                                        numberTightened++;
1820                                        //printf("%d bounds %g %g tightened to %g %g\n",
1821                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1822                                        //     newLower,newUpper);
1823                                        columnUpper_[iColumn] = newUpper;
1824                                        columnLower_[iColumn] = newLower;
1825                                        // and adjust bounds on rows
1826                                        newUpper -= upper;
1827                                        newLower -= lower;
1828                                        for (CoinBigIndex j = columnStart[iColumn];
1829                                                  j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1830                                             int iRow = row[j];
1831                                             double value = element[j];
1832                                             if (value > 0.0) {
1833                                                  up[iRow] += newUpper * value;
1834                                                  lo[iRow] += newLower * value;
1835                                             } else {
1836                                                  lo[iRow] += newUpper * value;
1837                                                  up[iRow] += newLower * value;
1838                                             }
1839                                        }
1840                                   } else {
1841                                        // infeasible
1842                                        //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1843                                        //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1844                                        //     newLower,newUpper);
1845                                        return -1;
1846                                   }
1847                              }
1848                         }
1849                    }
1850               }
1851          }
1852     }
1853     return numberTightened;
[636]1854}
1855/* Parametrics
1856   This is an initial slow version.
1857   The code uses current bounds + theta * change (if change array not NULL)
1858   and similarly for objective.
1859   It starts at startingTheta and returns ending theta in endingTheta.
1860   If reportIncrement 0.0 it will report on any movement
1861   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1862   If it can not reach input endingTheta return code will be 1 for infeasible,
1863   2 for unbounded, if error on ranges -1,  otherwise 0.
1864   Normal report is just theta and objective but
1865   if event handler exists it may do more
1866   On exit endingTheta is maximum reached (can be used for next startingTheta)
1867*/
[1525]1868int
1869ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
[1780]1870                             const double * lowerChangeBound, const double * upperChangeBound,
1871                             const double * lowerChangeRhs, const double * upperChangeRhs,
[636]1872                             const double * changeObjective)
1873{
[1525]1874     bool needToDoSomething = true;
1875     bool canTryQuick = (reportIncrement) ? true : false;
1876     // Save copy of model
1877     ClpSimplex copyModel = *this;
1878     int savePerturbation = perturbation_;
1879     perturbation_ = 102; // switch off
1880     while (needToDoSomething) {
1881          needToDoSomething = false;
1882          algorithm_ = -1;
1883
1884          // save data
1885          ClpDataSave data = saveData();
[1780]1886          // Dantzig
1887          ClpDualRowPivot * savePivot = dualRowPivot_;
1888          dualRowPivot_ = new ClpDualRowDantzig();
1889          dualRowPivot_->setModel(this);
[1525]1890          int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1891          int iRow, iColumn;
1892          double * chgUpper = NULL;
1893          double * chgLower = NULL;
1894          double * chgObjective = NULL;
1895
1896
1897          if (!returnCode) {
1898               // Find theta when bounds will cross over and create arrays
1899               int numberTotal = numberRows_ + numberColumns_;
1900               chgLower = new double[numberTotal];
1901               memset(chgLower, 0, numberTotal * sizeof(double));
1902               chgUpper = new double[numberTotal];
1903               memset(chgUpper, 0, numberTotal * sizeof(double));
1904               chgObjective = new double[numberTotal];
1905               memset(chgObjective, 0, numberTotal * sizeof(double));
1906               assert (!rowScale_);
1907               double maxTheta = 1.0e50;
[1780]1908               if (lowerChangeRhs || upperChangeRhs) {
[1525]1909                    for (iRow = 0; iRow < numberRows_; iRow++) {
1910                         double lower = rowLower_[iRow];
1911                         double upper = rowUpper_[iRow];
1912                         if (lower > upper) {
1913                              maxTheta = -1.0;
1914                              break;
1915                         }
[1780]1916                         double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0;
1917                         double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0;
[1525]1918                         if (lower > -1.0e20 && upper < 1.0e20) {
[1780]1919                              if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1920                                   maxTheta = (upper - lower) / (lowerChange - upperChange);
[1525]1921                              }
1922                         }
1923                         if (lower > -1.0e20) {
[1780]1924                              lower_[numberColumns_+iRow] += startingTheta * lowerChange;
1925                              chgLower[numberColumns_+iRow] = lowerChange;
[1525]1926                         }
1927                         if (upper < 1.0e20) {
[1780]1928                              upper_[numberColumns_+iRow] += startingTheta * upperChange;
1929                              chgUpper[numberColumns_+iRow] = upperChange;
[1525]1930                         }
1931                    }
1932               }
1933               if (maxTheta > 0.0) {
[1780]1934                    if (lowerChangeBound || upperChangeBound) {
[1525]1935                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1936                              double lower = columnLower_[iColumn];
1937                              double upper = columnUpper_[iColumn];
1938                              if (lower > upper) {
1939                                   maxTheta = -1.0;
1940                                   break;
1941                              }
[1780]1942                              double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0;
1943                              double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0;
[1525]1944                              if (lower > -1.0e20 && upper < 1.0e20) {
[1780]1945                                   if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1946                                        maxTheta = (upper - lower) / (lowerChange - upperChange);
[1525]1947                                   }
1948                              }
1949                              if (lower > -1.0e20) {
[1780]1950                                   lower_[iColumn] += startingTheta * lowerChange;
1951                                   chgLower[iColumn] = lowerChange;
[1525]1952                              }
1953                              if (upper < 1.0e20) {
[1780]1954                                   upper_[iColumn] += startingTheta * upperChange;
1955                                   chgUpper[iColumn] = upperChange;
[1525]1956                              }
1957                         }
1958                    }
1959                    if (maxTheta == 1.0e50)
1960                         maxTheta = COIN_DBL_MAX;
1961               }
1962               if (maxTheta < 0.0) {
1963                    // bad ranges or initial
1964                    returnCode = -1;
1965               }
[1643]1966               if (maxTheta < endingTheta) {
1967                 char line[100];
1968                 sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
1969                         endingTheta,maxTheta);
1970                 handler_->message(CLP_GENERAL,messages_)
1971                   << line << CoinMessageEol;
1972                 endingTheta = maxTheta;
1973               }
[1525]1974               if (endingTheta < startingTheta) {
1975                    // bad initial
1976                    returnCode = -2;
1977               }
[636]1978          }
[1525]1979          double saveEndingTheta = endingTheta;
1980          if (!returnCode) {
1981               if (changeObjective) {
1982                    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1983                         chgObjective[iColumn] = changeObjective[iColumn];
1984                         cost_[iColumn] += startingTheta * changeObjective[iColumn];
1985                    }
1986               }
1987               double * saveDuals = NULL;
1988               reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
1989               assert (!problemStatus_);
[1761]1990               for (int i=0;i<numberRows_+numberColumns_;i++)
1991                 setFakeBound(i, noFake);
[1525]1992               // Now do parametrics
[1643]1993               handler_->message(CLP_PARAMETRICS_STATS, messages_)
1994                 << startingTheta << objectiveValue() << CoinMessageEol;
[1525]1995               while (!returnCode) {
[1643]1996                    //assert (reportIncrement);
[1785]1997                 parametricsData paramData;
1998                 paramData.startingTheta=startingTheta;
1999                 paramData.endingTheta=endingTheta;
2000                 paramData.lowerChange = chgLower;
2001                 paramData.upperChange = chgUpper;
2002                    returnCode = parametricsLoop(paramData, reportIncrement,
[1525]2003                                                 chgLower, chgUpper, chgObjective, data,
2004                                                 canTryQuick);
[1785]2005                 startingTheta=paramData.startingTheta;
2006                 endingTheta=paramData.endingTheta;
[1525]2007                    if (!returnCode) {
2008                         //double change = endingTheta-startingTheta;
2009                         startingTheta = endingTheta;
2010                         endingTheta = saveEndingTheta;
2011                         //for (int i=0;i<numberTotal;i++) {
2012                         //lower_[i] += change*chgLower[i];
2013                         //upper_[i] += change*chgUpper[i];
2014                         //cost_[i] += change*chgObjective[i];
2015                         //}
[1643]2016                         handler_->message(CLP_PARAMETRICS_STATS, messages_)
2017                           << startingTheta << objectiveValue() << CoinMessageEol;
[1525]2018                         if (startingTheta >= endingTheta)
2019                              break;
2020                    } else if (returnCode == -1) {
2021                         // trouble - do external solve
2022                         needToDoSomething = true;
[1643]2023                    } else if (problemStatus_==1) {
2024                      // can't move any further
2025                      if (!canTryQuick) {
2026                        handler_->message(CLP_PARAMETRICS_STATS, messages_)
2027                          << endingTheta << objectiveValue() << CoinMessageEol;
2028                        problemStatus_=0;
2029                      }
[1525]2030                    } else {
2031                         abort();
2032                    }
2033               }
[636]2034          }
[1525]2035          reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2036
2037          delete dualRowPivot_;
2038          dualRowPivot_ = savePivot;
2039          // Restore any saved stuff
2040          restoreData(data);
2041          if (needToDoSomething) {
2042               double saveStartingTheta = startingTheta; // known to be feasible
2043               int cleanedUp = 1;
2044               while (cleanedUp) {
2045                    // tweak
2046                    if (cleanedUp == 1) {
2047                         if (!reportIncrement)
2048                              startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2049                         else
2050                              startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2051                    } else {
2052                         // restoring to go slowly
2053                         startingTheta = saveStartingTheta;
2054                    }
2055                    // only works if not scaled
2056                    int i;
2057                    const double * obj1 = objective();
2058                    double * obj2 = copyModel.objective();
2059                    const double * lower1 = columnLower_;
2060                    double * lower2 = copyModel.columnLower();
2061                    const double * upper1 = columnUpper_;
2062                    double * upper2 = copyModel.columnUpper();
2063                    for (i = 0; i < numberColumns_; i++) {
2064                         obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2065                         lower2[i] = lower1[i] + startingTheta * chgLower[i];
2066                         upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2067                    }
2068                    lower1 = rowLower_;
2069                    lower2 = copyModel.rowLower();
2070                    upper1 = rowUpper_;
2071                    upper2 = copyModel.rowUpper();
2072                    for (i = 0; i < numberRows_; i++) {
2073                         lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2074                         upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2075                    }
2076                    copyModel.dual();
2077                    if (copyModel.problemStatus()) {
[1643]2078                      char line[100];
2079                      sprintf(line,"Can not get to theta of %g\n", startingTheta);
2080                      handler_->message(CLP_GENERAL,messages_)
2081                        << line << CoinMessageEol;
[1525]2082                         canTryQuick = false; // do slowly to get exact amount
2083                         // back to last known good
2084                         if (cleanedUp == 1)
2085                              cleanedUp = 2;
2086                         else
2087                              abort();
2088                    } else {
2089                         // and move stuff back
2090                         int numberTotal = numberRows_ + numberColumns_;
2091                         CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2092                         CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2093                         CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2094                         cleanedUp = 0;
2095                    }
2096               }
[636]2097          }
[1525]2098          delete [] chgLower;
2099          delete [] chgUpper;
2100          delete [] chgObjective;
2101     }
2102     perturbation_ = savePerturbation;
[1643]2103     char line[100];
2104     sprintf(line,"Ending theta %g\n", endingTheta);
2105     handler_->message(CLP_GENERAL,messages_)
2106       << line << CoinMessageEol;
[1525]2107     return problemStatus_;
[636]2108}
[1643]2109/* Version of parametrics which reads from file
2110   See CbcClpParam.cpp for details of format
2111   Returns -2 if unable to open file */
2112int 
2113ClpSimplexOther::parametrics(const char * dataFile)
2114{
2115  int returnCode=-2;
2116  FILE *fp = fopen(dataFile, "r");
2117  char line[200];
2118  if (!fp) {
2119    handler_->message(CLP_UNABLE_OPEN, messages_)
2120      << dataFile << CoinMessageEol;
2121    return -2;
2122  }
2123
2124  if (!fgets(line, 200, fp)) {
2125    sprintf(line,"Empty parametrics file %s?",dataFile);
2126    handler_->message(CLP_GENERAL,messages_)
2127      << line << CoinMessageEol;
2128    fclose(fp);
2129    return -2;
2130  }
2131  char * pos = line;
2132  char * put = line;
2133  while (*pos >= ' ' && *pos != '\n') {
2134    if (*pos != ' ' && *pos != '\t') {
2135      *put = static_cast<char>(tolower(*pos));
2136      put++;
2137    }
2138    pos++;
2139  }
2140  *put = '\0';
2141  pos = line;
2142  double startTheta=0.0;
2143  double endTheta=0.0;
2144  double intervalTheta=COIN_DBL_MAX;
2145  int detail=0;
2146  bool good = true;
2147  while (good) {
2148    good=false;
2149    // check ROWS
2150    char * comma = strchr(pos, ',');
2151    if (!comma)
2152      break;
2153    *comma = '\0';
2154    if (strcmp(pos,"rows"))
2155      break;
2156    *comma = ',';
2157    pos = comma+1;
2158    // check lower theta
2159    comma = strchr(pos, ',');
2160    if (!comma)
2161      break;
2162    *comma = '\0';
2163    startTheta = atof(pos);
2164    *comma = ',';
2165    pos = comma+1;
2166    // check upper theta
2167    comma = strchr(pos, ',');
2168    good=true;
2169    if (comma)
2170      *comma = '\0';
2171    endTheta = atof(pos);
2172    if (comma) {
2173      *comma = ',';
2174      pos = comma+1;
2175      comma = strchr(pos, ',');
2176      if (comma)
2177        *comma = '\0';
2178      intervalTheta = atof(pos);
2179      if (comma) {
2180        *comma = ',';
2181        pos = comma+1;
2182        comma = strchr(pos, ',');
2183        if (comma)
2184          *comma = '\0';
2185        detail = atoi(pos);
2186        if (comma) 
2187        *comma = ',';
2188      }
2189    }
2190    break;
2191  }
2192  if (good) {
2193    if (startTheta<0.0||
2194        startTheta>endTheta||
2195        intervalTheta<0.0)
2196      good=false;
2197    if (detail<0||detail>1)
2198      good=false;
2199  }
2200  if (intervalTheta>=endTheta)
2201    intervalTheta=0.0;
2202  if (!good) {
2203    sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2204    handler_->message(CLP_GENERAL,messages_)
2205      << line << CoinMessageEol;
2206    fclose(fp);
2207    return -2;
2208  }
2209  if (!fgets(line, 200, fp)) {
2210    sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2211    handler_->message(CLP_GENERAL,messages_)
2212      << line << CoinMessageEol;
2213    fclose(fp);
2214    return -2;
2215  }
2216  double * lowerRowMove = NULL;
2217  double * upperRowMove = NULL;
2218  double * lowerColumnMove = NULL;
2219  double * upperColumnMove = NULL;
2220  double * objectiveMove = NULL;
2221  char saveLine[200];
2222  saveLine[0]='\0';
2223  std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2224  int gotRow[] = { -1, -1, -1, -1, -1};
2225  int orderRow[5];
2226  assert(sizeof(gotRow) == sizeof(orderRow));
2227  int nAcross = 0;
2228  pos = line;
2229  put = line;
2230  while (*pos >= ' ' && *pos != '\n') {
2231    if (*pos != ' ' && *pos != '\t') {
2232      *put = static_cast<char>(tolower(*pos));
2233      put++;
2234    }
2235    pos++;
2236  }
2237  *put = '\0';
2238  pos = line;
2239  int i;
2240  good = true;
2241  if (strncmp(line,"column",6)) {
2242    while (pos) {
2243      char * comma = strchr(pos, ',');
2244      if (comma)
2245        *comma = '\0';
2246      for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2247        if (headingsRow[i] == pos) {
2248          if (gotRow[i] < 0) {
2249            orderRow[nAcross] = i;
2250            gotRow[i] = nAcross++;
2251          } else {
2252            // duplicate
2253            good = false;
2254          }
2255          break;
2256        }
2257      }
2258      if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2259        good = false;
2260      if (comma) {
2261        *comma = ',';
2262        pos = comma + 1;
2263      } else {
2264        break;
2265      }
2266    }
2267    if (gotRow[0] < 0 && gotRow[1] < 0)
2268      good = false;
2269    if (gotRow[0] >= 0 && gotRow[1] >= 0)
2270      good = false;
2271    if (gotRow[0] >= 0 && !lengthNames())
2272      good = false;
2273    if (gotRow[4]<0) {
2274      if (gotRow[2] < 0 && gotRow[3] >= 0)
2275        good = false;
2276      else if (gotRow[3] < 0 && gotRow[2] >= 0)
2277        good = false;
2278    } else if (gotRow[2]>=0||gotRow[3]>=0) {
2279      good = false;
2280    }
2281    if (good) {
2282      char ** rowNames = new char * [numberRows_];
2283      int iRow;
2284      for (iRow = 0; iRow < numberRows_; iRow++) {
2285        rowNames[iRow] =
2286          CoinStrdup(rowName(iRow).c_str());
2287      }
2288      lowerRowMove = new double [numberRows_];
2289      memset(lowerRowMove,0,numberRows_*sizeof(double));
2290      upperRowMove = new double [numberRows_];
2291      memset(upperRowMove,0,numberRows_*sizeof(double));
2292      int nLine = 0;
2293      int nBadLine = 0;
2294      int nBadName = 0;
2295      bool goodLine=false;
2296      while (fgets(line, 200, fp)) {
2297        goodLine=true;
2298        if (!strncmp(line, "ENDATA", 6)||
2299            !strncmp(line, "COLUMN",6))
2300          break;
2301        goodLine=false;
2302        nLine++;
2303        iRow = -1;
2304        double upper = 0.0;
2305        double lower = 0.0;
2306        char * pos = line;
2307        char * put = line;
2308        while (*pos >= ' ' && *pos != '\n') {
2309          if (*pos != ' ' && *pos != '\t') {
2310            *put = *pos;
2311            put++;
2312          }
2313          pos++;
2314        }
2315        *put = '\0';
2316        pos = line;
2317        for (int i = 0; i < nAcross; i++) {
2318          char * comma = strchr(pos, ',');
2319          if (comma) {
2320            *comma = '\0';
2321          } else if (i < nAcross - 1) {
2322            nBadLine++;
2323            break;
2324          }
2325          switch (orderRow[i]) {
2326            // name
2327          case 0:
2328            // For large problems this could be slow
2329            for (iRow = 0; iRow < numberRows_; iRow++) {
2330              if (!strcmp(rowNames[iRow], pos))
2331                break;
2332            }
2333            if (iRow == numberRows_)
2334              iRow = -1;
2335            break;
2336            // number
2337          case 1:
2338            iRow = atoi(pos);
2339            if (iRow < 0 || iRow >= numberRows_)
2340              iRow = -1;
2341            break;
2342            // lower
2343          case 2:
2344            upper = atof(pos);
2345            break;
2346            // upper
2347          case 3:
2348            lower = atof(pos);
2349            break;
2350            // rhs
2351          case 4:
2352            lower = atof(pos);
2353            upper = lower;
2354            break;
2355          }
2356          if (comma) {
2357            *comma = ',';
2358            pos = comma + 1;
2359          }
2360        }
2361        if (iRow >= 0) {
2362          if (rowLower_[iRow]>-1.0e20)
2363            lowerRowMove[iRow] = lower;
2364          else
2365            lowerRowMove[iRow]=0.0;
2366          if (rowUpper_[iRow]<1.0e20)
2367            upperRowMove[iRow] = upper;
2368          else
2369            upperRowMove[iRow] = lower;
2370        } else {
2371          nBadName++;
2372          if(saveLine[0]=='\0')
2373            strcpy(saveLine,line);
2374        }
2375      }
2376      sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2377      handler_->message(CLP_GENERAL,messages_)
2378        << line << CoinMessageEol;
2379      if (nBadName) {
2380        sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2381        handler_->message(CLP_GENERAL,messages_)
2382          << line << CoinMessageEol;
2383        returnCode=-1;
2384        good=false;
2385      }
2386      for (iRow = 0; iRow < numberRows_; iRow++) {
2387        free(rowNames[iRow]);
2388      }
2389      delete [] rowNames;
2390    } else {
2391      sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2392      handler_->message(CLP_GENERAL,messages_)
2393        << line << CoinMessageEol;
2394      returnCode=-1;
2395      good=false;
2396    }
2397  }
2398  if (good&&(!strncmp(line, "COLUMN",6)||!strncmp(line, "column",6))) {
2399    if (!fgets(line, 200, fp)) {
2400      sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2401      handler_->message(CLP_GENERAL,messages_)
2402        << line << CoinMessageEol;
2403      fclose(fp);
2404      return -2;
2405    }
2406    std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2407    saveLine[0]='\0';
2408    int gotColumn[] = { -1, -1, -1, -1, -1};
2409    int orderColumn[5];
2410    assert(sizeof(gotColumn) == sizeof(orderColumn));
2411    nAcross = 0;
2412    pos = line;
2413    put = line;
2414    while (*pos >= ' ' && *pos != '\n') {
2415      if (*pos != ' ' && *pos != '\t') {
2416        *put = static_cast<char>(tolower(*pos));
2417        put++;
2418      }
2419      pos++;
2420    }
2421    *put = '\0';
2422    pos = line;
2423    int i;
2424    if (strncmp(line,"endata",6)&&good) {
2425      while (pos) {
2426        char * comma = strchr(pos, ',');
2427        if (comma)
2428          *comma = '\0';
2429        for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2430          if (headingsColumn[i] == pos) {
2431            if (gotColumn[i] < 0) {
2432              orderColumn[nAcross] = i;
2433              gotColumn[i] = nAcross++;
2434            } else {
2435              // duplicate
2436              good = false;
2437            }
2438            break;
2439          }
2440        }
2441        if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2442          good = false;
2443        if (comma) {
2444          *comma = ',';
2445          pos = comma + 1;
2446        } else {
2447          break;
2448        }
2449      }
2450      if (gotColumn[0] < 0 && gotColumn[1] < 0)
2451        good = false;
2452      if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2453        good = false;
2454      if (gotColumn[0] >= 0 && !lengthNames())
2455        good = false;
2456      if (good) {
2457        char ** columnNames = new char * [numberColumns_];
2458        int iColumn;
2459        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2460          columnNames[iColumn] =
2461            CoinStrdup(columnName(iColumn).c_str());
2462        }
2463        lowerColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2464        memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2465        upperColumnMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2466        memset(upperColumnMove,0,numberColumns_*sizeof(double));
2467        objectiveMove = reinterpret_cast<double *> (malloc(numberColumns_ * sizeof(double)));
2468        memset(objectiveMove,0,numberColumns_*sizeof(double));
2469        int nLine = 0;
2470        int nBadLine = 0;
2471        int nBadName = 0;
2472        bool goodLine=false;
2473        while (fgets(line, 200, fp)) {
2474          goodLine=true;
2475          if (!strncmp(line, "ENDATA", 6))
2476            break;
2477          goodLine=false;
2478          nLine++;
2479          iColumn = -1;
2480          double upper = 0.0;
2481          double lower = 0.0;
2482          double obj =0.0;
2483          char * pos = line;
2484          char * put = line;
2485          while (*pos >= ' ' && *pos != '\n') {
2486            if (*pos != ' ' && *pos != '\t') {
2487              *put = *pos;
2488              put++;
2489            }
2490            pos++;
2491          }
2492          *put = '\0';
2493          pos = line;
2494          for (int i = 0; i < nAcross; i++) {
2495            char * comma = strchr(pos, ',');
2496            if (comma) {
2497              *comma = '\0';
2498            } else if (i < nAcross - 1) {
2499              nBadLine++;
2500              break;
2501            }
2502            switch (orderColumn[i]) {
2503              // name
2504            case 0:
2505              // For large problems this could be slow
2506              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2507                if (!strcmp(columnNames[iColumn], pos))
2508                  break;
2509              }
2510              if (iColumn == numberColumns_)
2511                iColumn = -1;
2512              break;
2513              // number
2514            case 1:
2515              iColumn = atoi(pos);
2516              if (iColumn < 0 || iColumn >= numberColumns_)
2517                iColumn = -1;
2518              break;
2519              // lower
2520            case 2:
2521              upper = atof(pos);
2522              break;
2523              // upper
2524            case 3:
2525              lower = atof(pos);
2526              break;
2527              // objective
2528            case 4:
2529              obj = atof(pos);
2530              upper = lower;
2531              break;
2532            }
2533            if (comma) {
2534              *comma = ',';
2535              pos = comma + 1;
2536            }
2537          }
2538          if (iColumn >= 0) {
2539            if (columnLower_[iColumn]>-1.0e20)
2540              lowerColumnMove[iColumn] = lower;
2541            else
2542              lowerColumnMove[iColumn]=0.0;
2543            if (columnUpper_[iColumn]<1.0e20)
2544              upperColumnMove[iColumn] = upper;
2545            else
2546              upperColumnMove[iColumn] = lower;
2547            objectiveMove[iColumn] = obj;
2548          } else {
2549            nBadName++;
2550            if(saveLine[0]=='\0')
2551              strcpy(saveLine,line);
2552          }
2553        }
2554        sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2555        handler_->message(CLP_GENERAL,messages_)
2556          << line << CoinMessageEol;
2557        if (nBadName) {
2558          sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2559          handler_->message(CLP_GENERAL,messages_)
2560            << line << CoinMessageEol;
2561          returnCode=-1;
2562          good=false;
2563        }
2564        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2565          free(columnNames[iColumn]);
2566        }
2567        delete [] columnNames;
2568      } else {
2569        sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2570        handler_->message(CLP_GENERAL,messages_)
2571          << line << CoinMessageEol;
2572        returnCode=-1;
2573        good=false;
2574      }
2575    }
2576  }
2577  returnCode=-1;
2578  if (good) {
2579    // clean arrays
2580    if (lowerRowMove) {
2581      bool empty=true;
2582      for (int i=0;i<numberRows_;i++) {
2583        if (lowerRowMove[i]) {
2584          empty=false;
2585        break;
2586        }
2587      }
2588      if (empty) {
2589        delete [] lowerRowMove;
2590        lowerRowMove=NULL;
2591      }
2592    }
2593    if (upperRowMove) {
2594      bool empty=true;
2595      for (int i=0;i<numberRows_;i++) {
2596        if (upperRowMove[i]) {
2597          empty=false;
2598        break;
2599        }
2600      }
2601      if (empty) {
2602        delete [] upperRowMove;
2603        upperRowMove=NULL;
2604      }
2605    }
2606    if (lowerColumnMove) {
2607      bool empty=true;
2608      for (int i=0;i<numberColumns_;i++) {
2609        if (lowerColumnMove[i]) {
2610          empty=false;
2611        break;
2612        }
2613      }
2614      if (empty) {
2615        delete [] lowerColumnMove;
2616        lowerColumnMove=NULL;
2617      }
2618    }
2619    if (upperColumnMove) {
2620      bool empty=true;
2621      for (int i=0;i<numberColumns_;i++) {
2622        if (upperColumnMove[i]) {
2623          empty=false;
2624        break;
2625        }
2626      }
2627      if (empty) {
2628        delete [] upperColumnMove;
2629        upperColumnMove=NULL;
2630      }
2631    }
2632    if (objectiveMove) {
2633      bool empty=true;
2634      for (int i=0;i<numberColumns_;i++) {
2635        if (objectiveMove[i]) {
2636          empty=false;
2637        break;
2638        }
2639      }
2640      if (empty) {
2641        delete [] objectiveMove;
2642        objectiveMove=NULL;
2643      }
2644    }
2645    int saveScaling = scalingFlag_;
2646    scalingFlag_ = 0;
2647    int saveLogLevel = handler_->logLevel();
2648    if (detail>0&&!intervalTheta)
2649      handler_->setLogLevel(3);
2650    else
2651      handler_->setLogLevel(1);
2652    returnCode = parametrics(startTheta,endTheta,intervalTheta,
2653                             lowerColumnMove,upperColumnMove,
2654                             lowerRowMove,upperRowMove,
2655                             objectiveMove);
2656    scalingFlag_ = saveScaling;
2657    handler_->setLogLevel(saveLogLevel);
2658  }
2659  delete [] lowerRowMove;
2660  delete [] upperRowMove;
2661  delete [] lowerColumnMove;
2662  delete [] upperColumnMove;
2663  delete [] objectiveMove;
2664  fclose(fp);
2665  return returnCode;
2666}
[1525]2667int
[1785]2668ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
[1780]2669                                 const double * lowerChange, const double * upperChange,
[636]2670                                 const double * changeObjective, ClpDataSave & data,
2671                                 bool canTryQuick)
2672{
[1785]2673  double startingTheta = paramData.startingTheta;
2674  double & endingTheta = paramData.endingTheta;
[1525]2675     // stuff is already at starting
2676     // For this crude version just try and go to end
2677     double change = 0.0;
2678     if (reportIncrement && canTryQuick) {
2679          endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2680          change = endingTheta - startingTheta;
2681     }
2682     int numberTotal = numberRows_ + numberColumns_;
2683     int i;
2684     for ( i = 0; i < numberTotal; i++) {
[1780]2685          lower_[i] += change * lowerChange[i];
2686          upper_[i] += change * upperChange[i];
[1525]2687          switch(getStatus(i)) {
2688
2689          case basic:
2690          case isFree:
2691          case superBasic:
2692               break;
2693          case isFixed:
2694          case atUpperBound:
2695               solution_[i] = upper_[i];
2696               break;
2697          case atLowerBound:
2698               solution_[i] = lower_[i];
2699               break;
2700          }
2701          cost_[i] += change * changeObjective[i];
2702     }
2703     problemStatus_ = -1;
2704
2705     // This says whether to restore things etc
2706     // startup will have factorized so can skip
2707     int factorType = 0;
2708     // Start check for cycles
2709     progress_.startCheck();
2710     // Say change made on first iteration
2711     changeMade_ = 1;
2712     /*
2713       Status of problem:
2714       0 - optimal
2715       1 - infeasible
2716       2 - unbounded
2717       -1 - iterating
2718       -2 - factorization wanted
2719       -3 - redo checking without factorization
2720       -4 - looks infeasible
2721     */
2722     while (problemStatus_ < 0) {
2723          int iRow, iColumn;
2724          // clear
2725          for (iRow = 0; iRow < 4; iRow++) {
2726               rowArray_[iRow]->clear();
2727          }
2728
2729          for (iColumn = 0; iColumn < 2; iColumn++) {
2730               columnArray_[iColumn]->clear();
2731          }
2732
2733          // give matrix (and model costs and bounds a chance to be
2734          // refreshed (normally null)
2735          matrix_->refresh(this);
2736          // may factorize, checks if problem finished
2737          statusOfProblemInParametrics(factorType, data);
2738          // Say good factorization
2739          factorType = 1;
2740          if (data.sparseThreshold_) {
2741               // use default at present
2742               factorization_->sparseThreshold(0);
2743               factorization_->goSparse();
2744          }
2745
2746          // exit if victory declared
[1643]2747          if (problemStatus_ >= 0 && 
2748              (canTryQuick || startingTheta>=endingTheta-1.0e-7) )
[1525]2749               break;
2750
2751          // test for maximum iterations
2752          if (hitMaximumIterations()) {
2753               problemStatus_ = 3;
2754               break;
2755          }
2756          // Check event
2757          {
2758               int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2759               if (status >= 0) {
2760                    problemStatus_ = 5;
2761                    secondaryStatus_ = ClpEventHandler::endOfFactorization;
2762                    break;
2763               }
2764          }
2765          // Do iterations
[1643]2766          problemStatus_=-1;
[1525]2767          if (canTryQuick) {
2768               double * saveDuals = NULL;
2769               reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2770          } else {
[1785]2771               whileIterating(paramData, reportIncrement,
[1525]2772                              changeObjective);
[1643]2773               startingTheta = endingTheta;
[1525]2774          }
2775     }
2776     if (!problemStatus_) {
2777          theta_ = change + startingTheta;
2778          eventHandler_->event(ClpEventHandler::theta);
2779          return 0;
2780     } else if (problemStatus_ == 10) {
2781          return -1;
2782     } else {
2783          return problemStatus_;
2784     }
[636]2785}
[1780]2786/* Parametrics
2787   The code uses current bounds + theta * change (if change array not NULL)
2788   It starts at startingTheta and returns ending theta in endingTheta.
2789   If it can not reach input endingTheta return code will be 1 for infeasible,
2790   2 for unbounded, if error on ranges -1,  otherwise 0.
2791   Event handler may do more
2792   On exit endingTheta is maximum reached (can be used for next startingTheta)
2793*/
[1761]2794int
[1780]2795ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
2796                             const double * lowerChangeBound, const double * upperChangeBound,
2797                             const double * lowerChangeRhs, const double * upperChangeRhs)
[1761]2798{
[1780]2799  int savePerturbation = perturbation_;
2800  perturbation_ = 102; // switch off
2801  algorithm_ = -1;
2802  // extra region
2803  int maximumPivots = factorization_->maximumPivots();
2804  int numberDense = factorization_->numberDense();
2805  int length = numberRows_ + numberDense + maximumPivots;
2806  assert (!rowArray_[4]);
2807  rowArray_[4]=new CoinIndexedVector(length);
2808  assert (!rowArray_[5]);
2809  rowArray_[5]=new CoinIndexedVector(length);
2810
2811  // save data
2812  ClpDataSave data = saveData();
2813  int numberTotal = numberRows_ + numberColumns_;
2814  int ratio = (2*sizeof(int))/sizeof(double);
2815  assert (ratio==1||ratio==2);
2816  // allow for unscaled - even if not needed
[1785]2817  int lengthArrays = 4*numberTotal+(2*numberTotal+2)*ratio;
[1780]2818  /*
2819    Save information and modify
2820  */
2821  double * saveLower = new double [lengthArrays];
2822  double * saveUpper = new double [lengthArrays];
2823  double * lowerCopy = saveLower+2*numberTotal;
2824  double * upperCopy = saveUpper+2*numberTotal;
2825  double * lowerChange = saveLower+numberTotal;
2826  double * upperChange = saveUpper+numberTotal;
2827  int * lowerList = (reinterpret_cast<int *>(saveLower+4*numberTotal))+2;
2828  int * upperList = (reinterpret_cast<int *>(saveUpper+4*numberTotal))+2;
2829  // To mark as odd
2830  char * markDone = reinterpret_cast<char *>(lowerList+numberTotal);
2831  memset(markDone,0,numberTotal);
[1785]2832  int * backwardBasic = lowerList+2*numberTotal;
2833  parametricsData paramData;
2834  paramData.lowerChange = lowerChange;
2835  paramData.lowerList=lowerList;
2836  paramData.upperChange = upperChange;
2837  paramData.upperList=upperList;
2838  paramData.markDone=markDone;
2839  paramData.backwardBasic=backwardBasic;
[1780]2840  // Find theta when bounds will cross over and create arrays
2841  memset(lowerChange, 0, numberTotal * sizeof(double));
2842  memset(upperChange, 0, numberTotal * sizeof(double));
2843  if (lowerChangeBound)
2844    memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
2845  if (upperChangeBound)
2846    memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
2847  if (lowerChangeRhs)
2848    memcpy(lowerChange+numberColumns_,
2849           lowerChangeRhs,numberRows_*sizeof(double));
2850  if (upperChangeRhs)
2851    memcpy(upperChange+numberColumns_,
2852           upperChangeRhs,numberRows_*sizeof(double));
2853  int nLowerChange=0;
2854  int nUpperChange=0;
2855  for (int i=0;i<numberColumns_;i++) {
2856    if (lowerChange[i]) { 
2857      lowerList[nLowerChange++]=i;
2858    }
2859    if (upperChange[i]) { 
2860      upperList[nUpperChange++]=i;
2861    }
2862  }
2863  lowerList[-2]=nLowerChange;
2864  upperList[-2]=nUpperChange;
2865  for (int i=numberColumns_;i<numberTotal;i++) {
2866    if (lowerChange[i]) { 
2867      lowerList[nLowerChange++]=i;
2868    }
2869    if (upperChange[i]) { 
2870      upperList[nUpperChange++]=i;
2871    }
2872  }
2873  lowerList[-1]=nLowerChange;
2874  upperList[-1]=nUpperChange;
2875  memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
2876  memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
2877  memcpy(lowerCopy+numberColumns_,
2878         rowLower_,numberRows_*sizeof(double));
2879  memcpy(upperCopy+numberColumns_,
2880         rowUpper_,numberRows_*sizeof(double));
2881  {
2882    //  extra for unscaled
2883    double * unscaledCopy;
2884    unscaledCopy = lowerCopy + numberTotal; 
2885    memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
2886    memcpy(unscaledCopy+numberColumns_,
2887           rowLower_,numberRows_*sizeof(double));
2888    unscaledCopy = upperCopy + numberTotal; 
2889    memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
2890    memcpy(unscaledCopy+numberColumns_,
2891           rowUpper_,numberRows_*sizeof(double));
2892  }
2893  double maxTheta = 1.0e50;
2894  for (int iRow = 0; iRow < numberRows_; iRow++) {
2895    double lower = rowLower_[iRow];
2896    double upper = rowUpper_[iRow];
2897    if (lower<-1.0e30)
2898      lowerChange[numberColumns_+iRow]=0.0;
2899    double chgLower = lowerChange[numberColumns_+iRow];
2900    if (upper>1.0e30)
2901      upperChange[numberColumns_+iRow]=0.0;
2902    double chgUpper = upperChange[numberColumns_+iRow];
2903    if (lower > -1.0e30 && upper < 1.0e30) {
2904      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2905        maxTheta = (upper - lower) / (chgLower - chgUpper);
2906      }
2907    }
2908    lower+=startingTheta*chgLower;
2909    upper+=startingTheta*chgUpper;
2910    if (lower > upper) {
2911      maxTheta = -1.0;
[1761]2912      break;
[1780]2913    }
2914    rowLower_[iRow]=lower;
2915    rowUpper_[iRow]=upper;
2916  }
2917  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2918    double lower = columnLower_[iColumn];
2919    double upper = columnUpper_[iColumn];
2920    if (lower<-1.0e30)
2921      lowerChange[iColumn]=0.0;
2922    double chgLower = lowerChange[iColumn];
2923    if (upper>1.0e30)
2924      upperChange[iColumn]=0.0;
2925    double chgUpper = upperChange[iColumn];
2926    if (lower > -1.0e30 && upper < 1.0e30) {
2927      if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2928        maxTheta = (upper - lower) / (chgLower - chgUpper);
2929      }
2930    }
2931    lower+=startingTheta*chgLower;
2932    upper+=startingTheta*chgUpper;
2933    if (lower > upper) {
2934      maxTheta = -1.0;
[1761]2935      break;
2936    }
[1780]2937    columnLower_[iColumn]=lower;
2938    columnUpper_[iColumn]=upper;
[1761]2939  }
[1780]2940  if (maxTheta == 1.0e50)
2941    maxTheta = COIN_DBL_MAX;
2942  int returnCode=0;
2943  if (maxTheta < 0.0) {
2944    // bad ranges or initial
2945    returnCode = -1;
2946  }
2947  if (maxTheta < endingTheta) {
2948    char line[100];
2949    sprintf(line,"Crossover considerations reduce ending  theta from %g to %g\n", 
2950            endingTheta,maxTheta);
2951    handler_->message(CLP_GENERAL,messages_)
2952      << line << CoinMessageEol;
2953    endingTheta = maxTheta;
2954  }
2955  if (endingTheta < startingTheta) {
2956    // bad initial
2957    returnCode = -2;
2958  }
2959  bool swapped=false;
2960  // Dantzig
2961#define ALL_DANTZIG
2962#ifdef ALL_DANTZIG
2963  ClpDualRowPivot * savePivot = dualRowPivot_;
2964  dualRowPivot_ = new ClpDualRowDantzig();
2965  dualRowPivot_->setModel(this);
2966#else
2967  ClpDualRowPivot * savePivot = NULL;
2968#endif
2969  if (!returnCode) {
2970    returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
2971    if (!returnCode) {
2972      swapped=true;
2973      double * temp;
2974      memcpy(saveLower,lower_,numberTotal*sizeof(double));
2975      temp=saveLower;
2976      saveLower=lower_;
2977      lower_=temp;
2978      //columnLowerWork_ = lower_;
2979      //rowLowerWork_ = lower_ + numberColumns_;
2980      memcpy(saveUpper,upper_,numberTotal*sizeof(double));
2981      temp=saveUpper;
2982      saveUpper=upper_;
2983      upper_=temp;
2984      //columnUpperWork_ = upper_;
2985      //rowUpperWork_ = upper_ + numberColumns_;
2986      if (rowScale_) {
2987        // scale saved and change arrays
2988        double * lowerChange = lower_+numberTotal;
2989        double * upperChange = upper_+numberTotal;
2990        double * lowerSave = lowerChange+numberTotal;
2991        double * upperSave = upperChange+numberTotal;
2992        for (int i=0;i<numberColumns_;i++) {
2993          double multiplier = inverseColumnScale_[i];
2994          if (lowerSave[i]>-1.0e20)
2995            lowerSave[i] *= multiplier;
2996          if (upperSave[i]<1.0e20)
2997            upperSave[i] *= multiplier;
2998          lowerChange[i] *= multiplier;
2999          upperChange[i] *= multiplier;
3000        }
3001        lowerChange += numberColumns_;
3002        upperChange += numberColumns_;
3003        lowerSave += numberColumns_;
3004        upperSave += numberColumns_;
3005        for (int i=0;i<numberRows_;i++) {
3006          double multiplier = rowScale_[i];
3007          if (lowerSave[i]>-1.0e20)
3008            lowerSave[i] *= multiplier;
3009          if (upperSave[i]<1.0e20)
3010            upperSave[i] *= multiplier;
3011          lowerChange[i] *= multiplier;
3012          upperChange[i] *= multiplier;
3013        }
3014      }
3015      //double saveEndingTheta = endingTheta;
3016      double * saveDuals = NULL;
3017      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3018      if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3019        // probably can get rid of this if we adjust every change in theta
3020        //printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3021        //   sumPrimalInfeasibilities_);
3022        int pass=100;
3023        while(sumPrimalInfeasibilities_) {
3024          pass--;
3025          if (!pass)
3026            break;
3027          problemStatus_=-1;
3028          for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3029            double value=solution_[iSequence];
3030            // remember scaling
3031            if (value<lower_[iSequence]-1.0e-9) {
3032              lower_[iSequence]=value;
3033              lowerCopy[iSequence]=value;
3034            } else if (value>upper_[iSequence]+1.0e-9) {
3035              upper_[iSequence]=value;
3036              upperCopy[iSequence]=value;
3037            }
3038          }
3039          reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3040        }
3041      }
3042      assert (!problemStatus_);
3043      if (nLowerChange||nUpperChange) {
3044#ifndef ALL_DANTZIG
3045        // Dantzig
3046        savePivot = dualRowPivot_;
3047        dualRowPivot_ = new ClpDualRowDantzig();
3048        dualRowPivot_->setModel(this);
3049#endif
3050        for (int i=0;i<numberRows_+numberColumns_;i++)
3051          setFakeBound(i, noFake);
3052        // Now do parametrics
3053        handler_->message(CLP_PARAMETRICS_STATS, messages_)
3054          << startingTheta << objectiveValue() << CoinMessageEol;
3055        bool canSkipFactorization=true;
3056        while (!returnCode) {
[1785]3057                 paramData.startingTheta=startingTheta;
3058                 paramData.endingTheta=endingTheta;
3059                 returnCode = parametricsLoop(paramData,
[1780]3060                                       data,canSkipFactorization);
[1785]3061                 startingTheta=paramData.startingTheta;
3062                 endingTheta=paramData.endingTheta;
[1780]3063          canSkipFactorization=false;
3064          if (!returnCode) {
3065            //startingTheta = endingTheta;
3066            //endingTheta = saveEndingTheta;
3067            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3068              << startingTheta << objectiveValue() << CoinMessageEol;
3069            if (startingTheta >= endingTheta-primalTolerance_
3070                ||problemStatus_==2)
3071              break;
3072          } else if (returnCode == -1) {
3073            // trouble - do external solve
3074            abort(); //needToDoSomething = true;
3075          } else if (problemStatus_==1) {
3076            // can't move any further
3077            handler_->message(CLP_PARAMETRICS_STATS, messages_)
3078              << endingTheta << objectiveValue() << CoinMessageEol;
3079            problemStatus_=0;
3080          }
3081        }
3082      }
3083      //reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3084    }
3085    if (problemStatus_==2) {
3086      delete [] ray_;
3087      ray_ = new double [numberColumns_];
3088    }
3089    if (swapped&&lower_) {
3090      double * temp=saveLower;
3091      saveLower=lower_;
3092      lower_=temp;
3093      temp=saveUpper;
3094      saveUpper=upper_;
3095      upper_=temp;
3096    }
[1785]3097    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
[1780]3098  }   
3099  if (!scalingFlag_) {
3100    memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3101    memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3102    memcpy(rowLower_,lowerCopy+numberColumns_,
3103           numberRows_*sizeof(double));
3104    memcpy(rowUpper_,upperCopy+numberColumns_,
3105           numberRows_*sizeof(double));
3106  } else {
3107    //  extra for unscaled
3108    double * unscaledCopy;
3109    unscaledCopy = lowerCopy + numberTotal; 
3110    memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3111    memcpy(rowLower_,unscaledCopy+numberColumns_,
3112           numberRows_*sizeof(double));
3113    unscaledCopy = upperCopy + numberTotal; 
3114    memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3115    memcpy(rowUpper_,unscaledCopy+numberColumns_,
3116           numberRows_*sizeof(double));
3117  }
3118  delete [] saveLower;
3119  delete [] saveUpper;
3120#ifdef ALL_DANTZIG
3121  if (savePivot) {
3122#endif
3123    delete dualRowPivot_;
3124    dualRowPivot_ = savePivot;
3125#ifdef ALL_DANTZIG
3126  }
3127#endif
3128  // Restore any saved stuff
3129  restoreData(data);
3130  perturbation_ = savePerturbation;
3131  delete rowArray_[4];
3132  rowArray_[4]=NULL;
3133  delete rowArray_[5];
3134  rowArray_[5]=NULL;
3135  char line[100];
3136  sprintf(line,"Ending theta %g\n", endingTheta);
3137  handler_->message(CLP_GENERAL,messages_)
3138    << line << CoinMessageEol;
3139  return problemStatus_;
3140}
3141int
[1785]3142ClpSimplexOther::parametricsLoop(parametricsData & paramData,
[1780]3143                                 ClpDataSave & data,bool canSkipFactorization)
3144{
[1785]3145  double & startingTheta = paramData.startingTheta;
3146  double & endingTheta = paramData.endingTheta;
[1780]3147  int numberTotal = numberRows_+numberColumns_;
3148  // stuff is already at starting
[1785]3149  const int * lowerList = paramData.lowerList;
3150  const int * upperList = paramData.upperList;
[1761]3151  problemStatus_ = -1;
[1780]3152  //double saveEndingTheta=endingTheta;
[1761]3153
3154  // This says whether to restore things etc
3155  // startup will have factorized so can skip
3156  int factorType = 0;
3157  // Start check for cycles
3158  progress_.startCheck();
3159  // Say change made on first iteration
3160  changeMade_ = 1;
3161  /*
3162    Status of problem:
3163    0 - optimal
3164    1 - infeasible
3165    2 - unbounded
3166    -1 - iterating
3167    -2 - factorization wanted
3168    -3 - redo checking without factorization
3169    -4 - looks infeasible
3170  */
3171  while (problemStatus_ < 0) {
3172    int iRow, iColumn;
3173    // clear
[1780]3174    for (iRow = 0; iRow < 6; iRow++) {
[1761]3175      rowArray_[iRow]->clear();
3176    }
3177   
3178    for (iColumn = 0; iColumn < 2; iColumn++) {
3179      columnArray_[iColumn]->clear();
3180    }
3181   
3182    // give matrix (and model costs and bounds a chance to be
3183    // refreshed (normally null)
3184    matrix_->refresh(this);
3185    // may factorize, checks if problem finished
[1780]3186    if (!canSkipFactorization)
3187      statusOfProblemInParametrics(factorType, data);
3188    canSkipFactorization=false;
3189    if (numberPrimalInfeasibilities_) {
[1795]3190      if (largestPrimalError_>1.0e3&&startingTheta>1.0e10) {
3191        // treat as success
3192        problemStatus_=0;
3193        endingTheta=startingTheta;
3194        break;
3195      }
[1780]3196      // probably can get rid of this if we adjust every change in theta
3197      //printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3198      //     sumPrimalInfeasibilities_);
3199      const double * lowerChange = lower_+numberTotal;
3200      const double * upperChange = upper_+numberTotal;
3201      const double * startLower = lowerChange+numberTotal;
3202      const double * startUpper = upperChange+numberTotal;
3203      //startingTheta -= 1.0e-7;
3204      int nLowerChange = lowerList[-1];
3205      for (int i = 0; i < nLowerChange; i++) {
3206        int iSequence = lowerList[i];
3207        lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3208      }
3209      int nUpperChange = upperList[-1];
3210      for (int i = 0; i < nUpperChange; i++) {
3211        int iSequence = upperList[i];
3212        upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3213      }
3214      // adjust rhs in case dual uses
3215      memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3216      memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3217      memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3218      memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3219      if (rowScale_) {
3220        for (int i=0;i<numberColumns_;i++) {
3221          double multiplier = columnScale_[i];
3222          if (columnLower_[i]>-1.0e20)
3223            columnLower_[i] *= multiplier;
3224          if (columnUpper_[i]<1.0e20)
3225            columnUpper_[i] *= multiplier;
3226        }
3227        for (int i=0;i<numberRows_;i++) {
3228          double multiplier = inverseRowScale_[i];
3229          if (rowLower_[i]>-1.0e20)
3230            rowLower_[i] *= multiplier;
3231          if (rowUpper_[i]<1.0e20)
3232            rowUpper_[i] *= multiplier;
3233        }
3234      }
3235      double * saveDuals = NULL;
3236      problemStatus_=-1;
3237      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3238      int pass=100;
3239      double moved=0.0;
3240      while(sumPrimalInfeasibilities_) {
3241        pass--;
3242        if (!pass)
3243          break;
3244        problemStatus_=-1;
3245        for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3246          double value=solution_[iSequence];
3247          if (value<lower_[iSequence]-1.0e-9) {
3248            moved += lower_[iSequence]-value;
3249            lower_[iSequence]=value;
3250          } else if (value>upper_[iSequence]+1.0e-9) {
[1790]3251            moved += upper_[iSequence]-value;
[1780]3252            upper_[iSequence]=value;
3253          }
3254        }
[1790]3255        if (!moved) {
3256          for (int iSequence=0;iSequence<numberColumns_;iSequence++) {
3257            double value=solution_[iSequence];
3258            if (value<lower_[iSequence]-1.0e-9) {
3259              moved += lower_[iSequence]-value;
3260              lower_[iSequence]=value;
3261            } else if (value>upper_[iSequence]+1.0e-9) {
3262              moved += upper_[iSequence]-value;
3263              upper_[iSequence]=value;
3264            }
3265          }
3266        }
3267        assert (moved);
[1780]3268        reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3269      }
3270      // adjust
3271      //printf("Should adjust - moved %g\n",moved);
3272    }
[1761]3273    // Say good factorization
3274    factorType = 1;
3275    if (data.sparseThreshold_) {
3276      // use default at present
3277      factorization_->sparseThreshold(0);
3278      factorization_->goSparse();
3279    }
3280   
3281    // exit if victory declared
3282    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3283      break;
3284   
3285    // test for maximum iterations
3286    if (hitMaximumIterations()) {
3287      problemStatus_ = 3;
3288      break;
3289    }
[1780]3290#ifdef CLP_USER_DRIVEN
[1761]3291    // Check event
3292    {
3293      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3294      if (status >= 0) {
3295        problemStatus_ = 5;
3296        secondaryStatus_ = ClpEventHandler::endOfFactorization;
3297        break;
3298      }
3299    }
[1780]3300#endif
[1761]3301    // Do iterations
3302    problemStatus_=-1;
[1785]3303    whileIterating(paramData, 0.0,
[1761]3304                   NULL);
[1780]3305    //startingTheta = endingTheta;
3306    //endingTheta = saveEndingTheta;
[1761]3307  }
[1780]3308  if (!problemStatus_/*||problemStatus_==2*/) {
3309    theta_ = endingTheta;
3310#ifdef CLP_USER_DRIVEN
3311    {
3312      double saveTheta=theta_;
3313      theta_ = endingTheta;
3314      int status=eventHandler_->event(ClpEventHandler::theta);
3315      if (status>=0&&status<10) {
3316        endingTheta=theta_;
3317        theta_=saveTheta;
3318        problemStatus_=-1;
3319      } else {
3320        if (status>=10) {
3321          problemStatus_=status-10;
3322          startingTheta=endingTheta;
3323        }
3324        theta_=saveTheta;
3325      }
3326    }
3327#endif
[1761]3328    return 0;
3329  } else if (problemStatus_ == 10) {
3330    return -1;
3331  } else {
3332    return problemStatus_;
3333  }
3334}
[636]3335/* Checks if finished.  Updates status */
[1525]3336void
[636]3337ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3338{
[1525]3339     if (type == 2) {
3340          // trouble - go to recovery
3341          problemStatus_ = 10;
3342          return;
3343     }
3344     if (problemStatus_ > -3 || factorization_->pivots()) {
3345          // factorize
3346          // later on we will need to recover from singularities
3347          // also we could skip if first time
3348          if (type) {
3349               // is factorization okay?
3350               if (internalFactorize(1)) {
3351                    // trouble - go to recovery
3352                    problemStatus_ = 10;
3353                    return;
3354               }
3355          }
3356          if (problemStatus_ != -4 || factorization_->pivots() > 10)
3357               problemStatus_ = -3;
3358     }
3359     // at this stage status is -3 or -4 if looks infeasible
3360     // get primal and dual solutions
3361     gutsOfSolution(NULL, NULL);
3362     double realDualInfeasibilities = sumDualInfeasibilities_;
3363     // If bad accuracy treat as singular
3364     if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
3365          // trouble - go to recovery
3366          problemStatus_ = 10;
3367          return;
3368     } else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3369          // Can reduce tolerance
3370          double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3371          factorization_->pivotTolerance(newTolerance);
3372     }
3373     // Check if looping
3374     int loop;
3375     if (type != 2)
3376          loop = progress_.looping();
3377     else
3378          loop = -1;
3379     if (loop >= 0) {
3380          problemStatus_ = loop; //exit if in loop
3381          if (!problemStatus_) {
3382               // declaring victory
3383               numberPrimalInfeasibilities_ = 0;
3384               sumPrimalInfeasibilities_ = 0.0;
3385          } else {
3386               problemStatus_ = 10; // instead - try other algorithm
3387          }
3388          return;
3389     } else if (loop < -1) {
3390          // something may have changed
3391          gutsOfSolution(NULL, NULL);
3392     }
3393     progressFlag_ = 0; //reset progress flag
3394     if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3395          handler_->message(CLP_SIMPLEX_STATUS, messages_)
3396                    << numberIterations_ << objectiveValue();
3397          handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3398                    << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3399          handler_->printing(sumDualInfeasibilities_ > 0.0)
3400                    << sumDualInfeasibilities_ << numberDualInfeasibilities_;
3401          handler_->printing(numberDualInfeasibilitiesWithoutFree_
3402                             < numberDualInfeasibilities_)
3403                    << numberDualInfeasibilitiesWithoutFree_;
3404          handler_->message() << CoinMessageEol;
3405     }
[1790]3406#ifdef CLP_USER_DRIVEN
3407     if (sumPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-7) {
3408       int status=eventHandler_->event(ClpEventHandler::slightlyInfeasible);
3409       if (status>=0) {
3410         // fix up