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

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

For GUB

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