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

Last change on this file since 1722 was 1722, checked in by stefan, 10 years ago

adjust to changes in CoinUtils? header files

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