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

Last change on this file since 1330 was 1330, checked in by forrest, 11 years ago

improve dualize

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 103.5 KB
Line 
1// Copyright (C) 2004, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5
6#include <math.h>
7
8#include "CoinHelperFunctions.hpp"
9#include "ClpSimplexOther.hpp"
10#include "ClpSimplexDual.hpp"
11#include "ClpSimplexPrimal.hpp"
12#include "ClpEventHandler.hpp"
13#include "ClpHelperFunctions.hpp"
14#include "ClpFactorization.hpp"
15#include "ClpDualRowDantzig.hpp"
16#include "CoinPackedMatrix.hpp"
17#include "CoinIndexedVector.hpp"
18#include "CoinBuild.hpp"
19#include "CoinMpsIO.hpp"
20#include "ClpMessage.hpp"
21#include <cfloat>
22#include <cassert>
23#include <string>
24#include <stdio.h>
25#include <iostream>
26/* Dual ranging.
27   This computes increase/decrease in cost for each given variable and corresponding
28   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
29   and numberColumns.. for artificials/slacks.
30   For non-basic variables the sequence number will be that of the non-basic variables.
31   
32   Up to user to provide correct length arrays.
33   
34*/
35void ClpSimplexOther::dualRanging(int numberCheck,const int * which,
36                            double * costIncreased, int * sequenceIncreased,
37                                  double * costDecreased, int * sequenceDecreased,
38                                  double * valueIncrease, double * valueDecrease)
39{
40  rowArray_[1]->clear();
41  columnArray_[1]->clear();
42  // long enough for rows+columns
43  assert(rowArray_[3]->capacity()>=numberRows_+numberColumns_);
44  rowArray_[3]->clear();
45  int * backPivot = rowArray_[3]->getIndices();
46  int i;
47  for ( i=0;i<numberRows_+numberColumns_;i++) {
48    backPivot[i]=-1;
49  }
50  for (i=0;i<numberRows_;i++) {
51    int iSequence = pivotVariable_[i];
52    backPivot[iSequence]=i;
53  }
54  // dualTolerance may be zero if from CBC.  In fact use that fact
55  bool inCBC = !dualTolerance_;
56  if (inCBC)
57    assert (integerType_);
58  dualTolerance_ = dblParam_[ClpDualTolerance];
59  for ( i=0;i<numberCheck;i++) {
60    rowArray_[0]->clear();
61    //rowArray_[0]->checkClear();
62    //rowArray_[1]->checkClear();
63    //columnArray_[1]->checkClear();
64    columnArray_[0]->clear();
65    //columnArray_[0]->checkClear();
66    int iSequence = which[i];
67    double costIncrease=COIN_DBL_MAX;
68    double costDecrease=COIN_DBL_MAX;
69    int sequenceIncrease=-1;
70    int sequenceDecrease=-1;
71    if (valueIncrease) {
72      assert (valueDecrease);
73      valueIncrease[i]=iSequence<numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
74      valueDecrease[i]=valueIncrease[i];
75    }
76   
77    switch(getStatus(iSequence)) {
78     
79    case basic:
80      {
81        // non-trvial
82        // Get pivot row
83        int iRow=backPivot[iSequence];
84        assert (iRow>=0);
85        double plusOne=1.0;
86        rowArray_[0]->createPacked(1,&iRow,&plusOne);
87        factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
88        // put row of tableau in rowArray[0] and columnArray[0]
89        matrix_->transposeTimes(this,-1.0,
90                                rowArray_[0],columnArray_[1],columnArray_[0]);
91        double alphaIncrease;
92        double alphaDecrease;
93        // do ratio test up and down
94        checkDualRatios(rowArray_[0],columnArray_[0],costIncrease,sequenceIncrease,alphaIncrease,
95                    costDecrease,sequenceDecrease,alphaDecrease);
96        if (valueIncrease) {
97          if (sequenceIncrease>=0)
98            valueIncrease[i] = primalRanging1(sequenceIncrease,iSequence);
99          if (sequenceDecrease>=0)
100            valueDecrease[i] = primalRanging1(sequenceDecrease,iSequence);
101        }
102        if (inCBC) { 
103          if (sequenceIncrease>=0) {
104            double djValue = dj_[sequenceIncrease];
105            if (fabs(djValue)>10.0*dualTolerance_) {
106              // we are going to use for cutoff so be exact
107              costIncrease = fabs(djValue/alphaIncrease); 
108              /* Not sure this is good idea as I don't think correct e.g.
109                 suppose a continuous variable has dj slightly greater. */
110              if(false&&sequenceIncrease<numberColumns_&&integerType_[sequenceIncrease]) {
111                // can improve
112                double movement = (columnScale_==NULL) ? 1.0 : 
113                  rhsScale_/columnScale_[sequenceIncrease];
114                costIncrease = CoinMax(fabs(djValue*movement),costIncrease);
115              }
116            } else {
117              costIncrease=0.0;
118            }
119          }
120          if (sequenceDecrease>=0) {
121            double djValue = dj_[sequenceDecrease];
122            if (fabs(djValue)>10.0*dualTolerance_) {
123              // we are going to use for cutoff so be exact
124              costDecrease = fabs(djValue/alphaDecrease); 
125              if(sequenceDecrease<numberColumns_&&integerType_[sequenceDecrease]) {
126                // can improve
127                double movement = (columnScale_==NULL) ? 1.0 : 
128                  rhsScale_/columnScale_[sequenceDecrease];
129                costDecrease = CoinMax(fabs(djValue*movement),costDecrease);
130              }
131            } else {
132              costDecrease=0.0;
133            }
134          }
135        }
136      }
137      break;
138    case isFixed:
139      break;
140    case isFree:
141    case superBasic:
142      costIncrease=0.0;
143      costDecrease=0.0;
144      sequenceIncrease=iSequence;
145      sequenceDecrease=iSequence;
146      break;
147    case atUpperBound:
148      costIncrease = CoinMax(0.0,-dj_[iSequence]);
149      sequenceIncrease = iSequence;
150      if (valueIncrease) 
151        valueIncrease[i] = primalRanging1(iSequence,iSequence);
152      break;
153    case atLowerBound:
154      costDecrease = CoinMax(0.0,dj_[iSequence]);
155      sequenceDecrease = iSequence;
156      if (valueIncrease) 
157        valueDecrease[i] = primalRanging1(iSequence,iSequence);
158      break;
159    }
160    double scaleFactor;
161#ifdef CLP_AUXILIARY_MODEL
162    if (!auxiliaryModel_) {
163#endif
164      if (rowScale_) {
165        if (iSequence<numberColumns_) 
166          scaleFactor = 1.0/(objectiveScale_*columnScale_[iSequence]);
167        else
168          scaleFactor = rowScale_[iSequence-numberColumns_]/objectiveScale_;
169      } else {
170        scaleFactor = 1.0/objectiveScale_;
171      }
172#ifdef CLP_AUXILIARY_MODEL
173    } else {
174      if (auxiliaryModel_->rowScale()) {
175        if (iSequence<numberColumns_) 
176          scaleFactor = 1.0/(objectiveScale_*auxiliaryModel_->columnScale()[iSequence]);
177        else
178          scaleFactor = auxiliaryModel_->rowScale()[iSequence-numberColumns_]/objectiveScale_;
179      } else {
180        scaleFactor = 1.0/objectiveScale_;
181      }
182    }
183#endif
184    if (costIncrease<1.0e30)
185      costIncrease *= scaleFactor;
186    if (costDecrease<1.0e30)
187      costDecrease *= scaleFactor;
188    if (optimizationDirection_==1.0) {
189      costIncreased[i] = costIncrease;
190      sequenceIncreased[i] = sequenceIncrease;
191      costDecreased[i] = costDecrease;
192      sequenceDecreased[i] = sequenceDecrease;
193    } else if (optimizationDirection_==-1.0) {
194      costIncreased[i] = costDecrease;
195      sequenceIncreased[i] = sequenceDecrease;
196      costDecreased[i] = costIncrease;
197      sequenceDecreased[i] = sequenceIncrease;
198      if (valueIncrease) {
199        double temp = valueIncrease[i];
200        valueIncrease[i]=valueDecrease[i];
201        valueDecrease[i]=temp;
202      }
203    } else if (optimizationDirection_==0.0) {
204      // !!!!!! ???
205      costIncreased[i] = COIN_DBL_MAX;
206      sequenceIncreased[i] = -1;
207      costDecreased[i] = COIN_DBL_MAX;
208      sequenceDecreased[i] = -1;
209    } else {
210      abort();
211    }
212  }
213  //rowArray_[0]->clear();
214  //rowArray_[1]->clear();
215  //columnArray_[1]->clear();
216  //columnArray_[0]->clear();
217  //rowArray_[3]->clear();
218  if (!optimizationDirection_)
219    printf("*** ????? Ranging with zero optimization costs\n");
220}
221/*
222   Row array has row part of pivot row
223   Column array has column part.
224   This is used in dual ranging
225*/
226void
227ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
228                                 CoinIndexedVector * columnArray,
229                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
230                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
231{
232  double acceptablePivot = 1.0e-9;
233  double * work;
234  int number;
235  int * which;
236  int iSection;
237
238  double thetaDown = 1.0e31;
239  double thetaUp = 1.0e31;
240  int sequenceDown =-1;
241  int sequenceUp = -1;
242  double alphaDown=0.0;
243  double alphaUp=0.0;
244
245  int addSequence;
246
247  for (iSection=0;iSection<2;iSection++) {
248
249    int i;
250    if (!iSection) {
251      work = rowArray->denseVector();
252      number = rowArray->getNumElements();
253      which = rowArray->getIndices();
254      addSequence = numberColumns_;
255    } else {
256      work = columnArray->denseVector();
257      number = columnArray->getNumElements();
258      which = columnArray->getIndices();
259      addSequence = 0;
260    }
261   
262    for (i=0;i<number;i++) {
263      int iSequence = which[i];
264      int iSequence2 = iSequence + addSequence;
265      double alpha=work[i];
266      if (fabs(alpha)<acceptablePivot)
267        continue;
268      double oldValue=dj_[iSequence2];
269
270      switch(getStatus(iSequence2)) {
271         
272      case basic: 
273        break;
274      case ClpSimplex::isFixed:
275        break;
276      case isFree:
277      case superBasic:
278        // treat dj as if zero
279        thetaDown=0.0;
280        thetaUp=0.0;
281        sequenceDown=iSequence2;
282        sequenceUp=iSequence2;
283        break;
284      case atUpperBound:
285        if (alpha>0.0) {
286          // test up
287          if (oldValue + thetaUp*alpha > dualTolerance_) {
288            thetaUp = (dualTolerance_-oldValue)/alpha;
289            sequenceUp = iSequence2;
290            alphaUp=alpha;
291          }
292        } else {
293          // test down
294          if (oldValue - thetaDown*alpha > dualTolerance_) {
295            thetaDown = -(dualTolerance_-oldValue)/alpha;
296            sequenceDown = iSequence2;
297            alphaDown=alpha;
298          }
299        }
300        break;
301      case atLowerBound:
302        if (alpha<0.0) {
303          // test up
304          if (oldValue + thetaUp*alpha <- dualTolerance_) {
305            thetaUp = -(dualTolerance_+oldValue)/alpha;
306            sequenceUp = iSequence2;
307            alphaUp=alpha;
308          }
309        } else {
310          // test down
311          if (oldValue - thetaDown*alpha < -dualTolerance_) {
312            thetaDown = (dualTolerance_+oldValue)/alpha;
313            sequenceDown = iSequence2;
314            alphaDown=alpha;
315          }
316        }
317        break;
318      }
319    }
320  }
321  if (sequenceUp>=0) {
322    costIncrease = thetaUp;
323    sequenceIncrease = sequenceUp;
324    alphaIncrease = alphaUp;
325  }
326  if (sequenceDown>=0) {
327    costDecrease = thetaDown;
328    sequenceDecrease = sequenceDown;
329    alphaDecrease = alphaDown;
330  }
331}
332/** Primal ranging.
333    This computes increase/decrease in value for each given variable and corresponding
334    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
335    and numberColumns.. for artificials/slacks.
336    For basic variables the sequence number will be that of the basic variables.
337   
338    Up to user to provide correct length arrays.
339   
340    When here - guaranteed optimal
341*/
342void 
343ClpSimplexOther::primalRanging(int numberCheck,const int * which,
344                  double * valueIncreased, int * sequenceIncreased,
345                  double * valueDecreased, int * sequenceDecreased)
346{
347  rowArray_[0]->clear();
348  rowArray_[1]->clear();
349  lowerIn_=-COIN_DBL_MAX;
350  upperIn_=COIN_DBL_MAX;
351  valueIn_ = 0.0;
352  for ( int i=0;i<numberCheck;i++) {
353    int iSequence = which[i];
354    double valueIncrease=COIN_DBL_MAX;
355    double valueDecrease=COIN_DBL_MAX;
356    int sequenceIncrease=-1;
357    int sequenceDecrease=-1;
358   
359    switch(getStatus(iSequence)) {
360     
361    case basic:
362    case isFree:
363    case superBasic:
364      // Easy
365      valueDecrease=CoinMax(0.0,upper_[iSequence]-solution_[iSequence]);
366      valueIncrease=CoinMax(0.0,solution_[iSequence]-lower_[iSequence]);
367      sequenceDecrease=iSequence;
368      sequenceIncrease=iSequence;
369      break;
370    case isFixed:
371    case atUpperBound:
372    case atLowerBound:
373      {
374        // Non trivial
375        // Other bound is ignored
376        unpackPacked(rowArray_[1],iSequence);
377        factorization_->updateColumn(rowArray_[2],rowArray_[1]);
378        // Get extra rows
379        matrix_->extendUpdated(this,rowArray_[1],0);
380        // do ratio test
381        checkPrimalRatios(rowArray_[1],1);
382        if (pivotRow_>=0) {
383          valueIncrease = theta_;
384          sequenceIncrease=pivotVariable_[pivotRow_];
385        }
386        checkPrimalRatios(rowArray_[1],-1);
387        if (pivotRow_>=0) {
388          valueDecrease = theta_;
389          sequenceDecrease=pivotVariable_[pivotRow_];
390        }
391        rowArray_[1]->clear();
392      }
393      break;
394    }
395    double scaleFactor;
396    if (rowScale_) {
397      if (iSequence<numberColumns_) 
398        scaleFactor = columnScale_[iSequence]/rhsScale_;
399      else
400        scaleFactor = 1.0/(rowScale_[iSequence-numberColumns_]*rhsScale_);
401    } else {
402      scaleFactor = 1.0/rhsScale_;
403    }
404    if (valueIncrease<1.0e30)
405      valueIncrease *= scaleFactor;
406    else
407      valueIncrease = COIN_DBL_MAX;
408    if (valueDecrease<1.0e30)
409      valueDecrease *= scaleFactor;
410    else
411      valueDecrease = COIN_DBL_MAX;
412    valueIncreased[i] = valueIncrease;
413    sequenceIncreased[i] = sequenceIncrease;
414    valueDecreased[i] = valueDecrease;
415    sequenceDecreased[i] = sequenceDecrease;
416  }
417}
418// Returns new value of whichOther when whichIn enters basis
419double 
420ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
421{
422  rowArray_[0]->clear();
423  rowArray_[1]->clear();
424  int iSequence = whichIn;
425  double newValue=solution_[whichOther];
426  double alphaOther=0.0;
427  Status status = getStatus(iSequence);
428  assert (status==atLowerBound||status==atUpperBound);
429  int wayIn = (status==atLowerBound) ? 1 : -1;
430 
431  switch(getStatus(iSequence)) {
432   
433  case basic:
434  case isFree:
435  case superBasic:
436    assert (whichIn==whichOther);
437    // Easy
438    newValue = wayIn>0 ? upper_[iSequence] : lower_[iSequence];
439    break;
440  case isFixed:
441  case atUpperBound:
442  case atLowerBound:
443    // Non trivial
444    {
445      // Other bound is ignored
446      unpackPacked(rowArray_[1],iSequence);
447      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
448      // Get extra rows
449      matrix_->extendUpdated(this,rowArray_[1],0);
450      // do ratio test
451      double acceptablePivot=1.0e-7;
452      double * work=rowArray_[1]->denseVector();
453      int number=rowArray_[1]->getNumElements();
454      int * which=rowArray_[1]->getIndices();
455     
456      // we may need to swap sign
457      double way = wayIn;
458      double theta = 1.0e30;
459      for (int iIndex=0;iIndex<number;iIndex++) {
460
461        int iRow = which[iIndex];
462        double alpha = work[iIndex]*way;
463        int iPivot=pivotVariable_[iRow];
464        if (iPivot==whichOther) {
465          alphaOther=alpha;
466          continue;
467        }
468        double oldValue = solution_[iPivot];
469        if (fabs(alpha)>acceptablePivot) {
470          if (alpha>0.0) {
471            // basic variable going towards lower bound
472            double bound = lower_[iPivot];
473            oldValue -= bound;
474            if (oldValue-theta*alpha<0.0) {
475              theta = CoinMax(0.0,oldValue/alpha);
476            }
477          } else {
478            // basic variable going towards upper bound
479            double bound = upper_[iPivot];
480            oldValue = oldValue-bound;
481            if (oldValue-theta*alpha>0.0) {
482              theta = CoinMax(0.0,oldValue/alpha);
483            }
484          }
485        }
486      }
487      if (whichIn!=whichOther) {
488        if (theta<1.0e30)
489          newValue -= theta*alphaOther;
490        else
491          newValue = alphaOther>0.0 ? -1.0e30 : 1.0e30;
492      } else {
493        newValue += theta*wayIn;
494      }
495    }
496    rowArray_[1]->clear();
497    break;
498  }
499  double scaleFactor;
500  if (rowScale_) {
501    if (whichOther<numberColumns_) 
502      scaleFactor = columnScale_[whichOther]/rhsScale_;
503    else
504      scaleFactor = 1.0/(rowScale_[whichOther-numberColumns_]*rhsScale_);
505  } else {
506    scaleFactor = 1.0/rhsScale_;
507  }
508  if (newValue<1.0e29)
509    if (newValue>-1.0e29)
510      newValue *= scaleFactor;
511    else
512      newValue = -COIN_DBL_MAX;
513  else
514    newValue = COIN_DBL_MAX;
515  return newValue;
516}
517/*
518   Row array has pivot column
519   This is used in primal ranging
520*/
521void 
522ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
523                         int direction)
524{
525  // sequence stays as row number until end
526  pivotRow_=-1;
527  double acceptablePivot=1.0e-7;
528  double * work=rowArray->denseVector();
529  int number=rowArray->getNumElements();
530  int * which=rowArray->getIndices();
531
532  // we need to swap sign if going down
533  double way = direction;
534  theta_ = 1.0e30;
535  for (int iIndex=0;iIndex<number;iIndex++) {
536   
537    int iRow = which[iIndex];
538    double alpha = work[iIndex]*way;
539    int iPivot=pivotVariable_[iRow];
540    double oldValue = solution_[iPivot];
541    if (fabs(alpha)>acceptablePivot) {
542      if (alpha>0.0) {
543        // basic variable going towards lower bound
544        double bound = lower_[iPivot];
545        oldValue -= bound;
546        if (oldValue-theta_*alpha<0.0) {
547          pivotRow_ = iRow;
548          theta_ = CoinMax(0.0,oldValue/alpha);
549        }
550      } else {
551        // basic variable going towards upper bound
552        double bound = upper_[iPivot];
553        oldValue = oldValue-bound;
554        if (oldValue-theta_*alpha>0.0) {
555          pivotRow_ = iRow;
556          theta_ = CoinMax(0.0,oldValue/alpha);
557        }
558      }
559    }
560  }
561}
562/* Write the basis in MPS format to the specified file.
563   If writeValues true writes values of structurals
564   (and adds VALUES to end of NAME card)
565   
566   Row and column names may be null.
567   formatType is
568   <ul>
569   <li> 0 - normal
570   <li> 1 - extra accuracy
571   <li> 2 - IEEE hex (later)
572   </ul>
573   
574   Returns non-zero on I/O error
575
576   This is based on code contributed by Thorsten Koch
577*/
578int 
579ClpSimplexOther::writeBasis(const char *filename,
580                            bool writeValues,
581                            int formatType) const
582{
583  formatType=CoinMax(0,formatType);
584  formatType=CoinMin(2,formatType);
585  if (!writeValues)
586    formatType=0;
587  // See if INTEL if IEEE
588  if (formatType==2) {
589    // test intel here and add 1 if not intel
590    double value=1.0;
591    char x[8];
592    memcpy(x,&value,8);
593    if (x[0]==63) {
594      formatType ++; // not intel
595    } else {
596      assert (x[0]==0);
597    }
598  }
599 
600  char number[20];
601  FILE * fp = fopen(filename,"w");
602  if (!fp)
603    return -1;
604   
605  // NAME card
606
607  if (strcmp(strParam_[ClpProbName].c_str(),"")==0) {
608    fprintf(fp, "NAME          BLANK      ");
609  } else {
610    fprintf(fp, "NAME          %s       ",strParam_[ClpProbName].c_str());
611  }
612  if (formatType>=2)
613    fprintf(fp,"IEEE");
614  else if (writeValues)
615    fprintf(fp,"VALUES");
616  // finish off name
617  fprintf(fp,"\n");
618  int iRow=0;
619  for(int iColumn =0; iColumn < numberColumns_; iColumn++) {
620    bool printit=false;
621    if( getColumnStatus(iColumn) == ClpSimplex::basic) {
622      printit=true;
623      // Find non basic row
624      for(; iRow < numberRows_; iRow++) {
625        if (getRowStatus(iRow) != ClpSimplex::basic) 
626          break;
627      }
628      if (lengthNames_) {
629        if (iRow!=numberRows_) {
630          fprintf(fp, " %s %-8s       %s",
631                  getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
632                  columnNames_[iColumn].c_str(),
633                  rowNames_[iRow].c_str());
634          iRow++;
635        } else {
636          // Allow for too many basics!
637          fprintf(fp, " BS %-8s       ",
638                  columnNames_[iColumn].c_str());
639          // Dummy row name if values
640          if (writeValues)
641            fprintf(fp,"      _dummy_");
642        }
643      } else {
644        // no names
645        if (iRow!=numberRows_) {
646          fprintf(fp, " %s C%7.7d     R%7.7d",
647                  getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
648                  iColumn,iRow);
649          iRow++;
650        } else {
651          // Allow for too many basics!
652          fprintf(fp, " BS C%7.7d",iColumn);
653          // Dummy row name if values
654          if (writeValues)
655            fprintf(fp,"      _dummy_");
656        }
657      }
658    } else  {
659      if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
660        printit=true;
661        if (lengthNames_) 
662          fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
663        else 
664          fprintf(fp, " UL C%7.7d", iColumn);
665        // Dummy row name if values
666        if (writeValues)
667          fprintf(fp,"      _dummy_");
668      }
669    }
670    if (printit&&writeValues) {
671      // add value
672      CoinConvertDouble(0,formatType,columnActivity_[iColumn],number);
673      fprintf(fp,"     %s",number);
674    }
675    if (printit)
676      fprintf(fp,"\n");
677  }
678  fprintf(fp, "ENDATA\n");
679  fclose(fp);
680  return 0;
681}
682// Read a basis from the given filename
683int 
684ClpSimplexOther::readBasis(const char *fileName)
685{
686  int status=0;
687  bool canOpen=false;
688  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
689    // stdin
690    canOpen=true;
691  } else {
692    FILE *fp=fopen(fileName,"r");
693    if (fp) {
694      // can open - lets go for it
695      fclose(fp);
696      canOpen=true;
697    } else {
698      handler_->message(CLP_UNABLE_OPEN,messages_)
699        <<fileName<<CoinMessageEol;
700      return -1;
701    }
702  }
703  CoinMpsIO m;
704  m.passInMessageHandler(handler_);
705  *m.messagesPointer()=coinMessages();
706  bool savePrefix =m.messageHandler()->prefix();
707  m.messageHandler()->setPrefix(handler_->prefix());
708  status=m.readBasis(fileName,"",columnActivity_,status_+numberColumns_,
709                     status_,
710                     columnNames_,numberColumns_,
711                     rowNames_,numberRows_);
712  m.messageHandler()->setPrefix(savePrefix);
713  if (status>=0) {
714    if (!status) {
715      // set values
716      int iColumn,iRow;
717      for (iRow=0;iRow<numberRows_;iRow++) {
718        if (getRowStatus(iRow)==atLowerBound)
719          rowActivity_[iRow]=rowLower_[iRow];
720        else if (getRowStatus(iRow)==atUpperBound)
721          rowActivity_[iRow]=rowUpper_[iRow];
722      }
723      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
724        if (getColumnStatus(iColumn)==atLowerBound)
725          columnActivity_[iColumn]=columnLower_[iColumn];
726        else if (getColumnStatus(iColumn)==atUpperBound)
727          columnActivity_[iColumn]=columnUpper_[iColumn];
728      }
729    } else {
730      memset(rowActivity_,0,numberRows_*sizeof(double));
731      matrix_->times(-1.0,columnActivity_,rowActivity_);
732    }
733  } else {
734    // errors
735    handler_->message(CLP_IMPORT_ERRORS,messages_)
736      <<status<<fileName<<CoinMessageEol;
737  }
738  return status;
739}
740/* Creates dual of a problem if looks plausible
741   (defaults will always create model)
742   fractionRowRanges is fraction of rows allowed to have ranges
743   fractionColumnRanges is fraction of columns allowed to have ranges
744*/
745ClpSimplex * 
746ClpSimplexOther::dualOfModel(double fractionRowRanges,double fractionColumnRanges) const
747{
748  const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
749  bool changed=false;
750  int numberChanged=0;
751  int iColumn;
752  // check if we need to change bounds to rows
753  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
754    if (columnUpper_[iColumn]<1.0e20&&
755        columnLower_[iColumn]>-1.0e20) {
756      changed=true;
757      numberChanged++;
758    }
759  }
760  int iRow;
761  int numberExtraRows=0;
762  if (numberChanged<=fractionColumnRanges*numberColumns_) {
763    for (iRow=0;iRow<numberRows_;iRow++) {
764      if (rowLower_[iRow]>-1.0e20&&
765          rowUpper_[iRow]<1.0e20) {
766        if (rowUpper_[iRow]!=rowLower_[iRow])
767          numberExtraRows++;
768      }
769    }
770    if (numberExtraRows>fractionRowRanges*numberRows_)
771      return NULL;
772  } else {
773    return NULL;
774  }
775  if (changed) {
776    ClpSimplex * model3 = new ClpSimplex(*model2);
777    CoinBuild build;
778    double one=1.0;
779    int numberColumns = model3->numberColumns();
780    const double * columnLower = model3->columnLower();
781    const double * columnUpper = model3->columnUpper();
782    for (iColumn=0;iColumn<numberColumns;iColumn++) {
783      if (columnUpper[iColumn]<1.0e20&&
784          columnLower[iColumn]>-1.0e20) {
785        if (fabs(columnLower[iColumn])<fabs(columnUpper[iColumn])) {
786          double value = columnUpper[iColumn];
787          model3->setColumnUpper(iColumn,COIN_DBL_MAX);
788          build.addRow(1,&iColumn,&one,-COIN_DBL_MAX,value);
789        } else {
790          double value = columnLower[iColumn];
791          model3->setColumnLower(iColumn,-COIN_DBL_MAX);
792          build.addRow(1,&iColumn,&one,value,COIN_DBL_MAX);
793        }
794      }
795    }
796    model3->addRows(build);
797    model2=model3;
798  }
799  int numberColumns = model2->numberColumns();
800  const double * columnLower = model2->columnLower();
801  const double * columnUpper = model2->columnUpper();
802  int numberRows = model2->numberRows();
803  double * rowLower = CoinCopyOfArray(model2->rowLower(),numberRows);
804  double * rowUpper = CoinCopyOfArray(model2->rowUpper(),numberRows);
805
806  const double * objective = model2->objective();
807  CoinPackedMatrix * matrix = model2->matrix();
808  // get transpose
809  CoinPackedMatrix rowCopy = *matrix;
810  const int * row = matrix->getIndices();
811  const int * columnLength = matrix->getVectorLengths();
812  const CoinBigIndex * columnStart = matrix->getVectorStarts();
813  const double * elementByColumn = matrix->getElements();
814  double objOffset=0.0;
815  for (iColumn=0;iColumn<numberColumns;iColumn++) {
816    double offset=0.0;
817    double objValue =optimizationDirection_*objective[iColumn];
818    if (columnUpper[iColumn]>1.0e20) {
819      if (columnLower[iColumn]>-1.0e20)
820        offset=columnLower[iColumn];
821    } else if (columnLower[iColumn]<-1.0e20) {
822      offset=columnUpper[iColumn];
823    } else {
824      // taken care of before
825      abort();
826    }
827    if (offset) {
828      objOffset += offset*objValue;
829      for (CoinBigIndex j=columnStart[iColumn];
830           j<columnStart[iColumn]+columnLength[iColumn];j++) {
831        int iRow = row[j];
832        if (rowLower[iRow]>-1.0e20)
833          rowLower[iRow] -= offset*elementByColumn[j];
834        if (rowUpper[iRow]<1.0e20)
835          rowUpper[iRow] -= offset*elementByColumn[j];
836      }
837    }
838  }
839  int * which = new int[numberRows+numberExtraRows];
840  rowCopy.reverseOrdering();
841  rowCopy.transpose();
842  double * fromRowsLower = new double[numberRows+numberExtraRows];
843  double * fromRowsUpper = new double[numberRows+numberExtraRows];
844  double * newObjective = new double[numberRows+numberExtraRows];
845  double * fromColumnsLower = new double[numberColumns];
846  double * fromColumnsUpper = new double[numberColumns];
847  for (iColumn=0;iColumn<numberColumns;iColumn++) {
848    double objValue =optimizationDirection_*objective[iColumn];
849    // Offset is already in
850    if (columnUpper[iColumn]>1.0e20) {
851      if (columnLower[iColumn]>-1.0e20) {
852        fromColumnsLower[iColumn]=-COIN_DBL_MAX;
853        fromColumnsUpper[iColumn]=objValue;
854      } else {
855        // free
856        fromColumnsLower[iColumn]=objValue;
857        fromColumnsUpper[iColumn]=objValue;
858      }
859    } else if (columnLower[iColumn]<-1.0e20) {
860      fromColumnsLower[iColumn]=objValue;
861      fromColumnsUpper[iColumn]=COIN_DBL_MAX;
862    } else {
863      abort();
864    }
865  }
866  int kRow=0;
867  int kExtraRow=numberRows;
868  for (iRow=0;iRow<numberRows;iRow++) {
869    if (rowLower[iRow]<-1.0e20) {
870      assert (rowUpper[iRow]<1.0e20);
871      newObjective[kRow]=-rowUpper[iRow];
872      fromRowsLower[kRow]=-COIN_DBL_MAX;
873      fromRowsUpper[kRow]=0.0;
874      which[kRow]=iRow;
875      kRow++;
876    } else if (rowUpper[iRow]>1.0e20) {
877      newObjective[kRow]=-rowLower[iRow];
878      fromRowsLower[kRow]=0.0;
879      fromRowsUpper[kRow]=COIN_DBL_MAX;
880      which[kRow]=iRow;
881      kRow++;
882    } else {
883      if (rowUpper[iRow]==rowLower[iRow]) {
884        newObjective[kRow]=-rowLower[iRow];
885        fromRowsLower[kRow]=-COIN_DBL_MAX;;
886        fromRowsUpper[kRow]=COIN_DBL_MAX;
887        which[kRow]=iRow;
888        kRow++;
889      } else {
890        // range
891        newObjective[kRow]=-rowUpper[iRow];
892        fromRowsLower[kRow]=-COIN_DBL_MAX;
893        fromRowsUpper[kRow]=0.0;
894        which[kRow]=iRow;
895        kRow++;
896        newObjective[kExtraRow]=-rowLower[iRow];
897        fromRowsLower[kExtraRow]=0.0;
898        fromRowsUpper[kExtraRow]=COIN_DBL_MAX;
899        which[kExtraRow]=iRow;
900        kExtraRow++;
901      }
902    }
903  }
904  if (numberExtraRows) {
905    CoinPackedMatrix newCopy;
906    newCopy.setExtraGap(0.0);
907    newCopy.setExtraMajor(0.0);
908    newCopy.submatrixOfWithDuplicates(rowCopy,kExtraRow,which);
909    rowCopy=newCopy;
910  }
911  ClpSimplex * modelDual = new ClpSimplex();
912  modelDual->loadProblem(rowCopy,fromRowsLower,fromRowsUpper,newObjective,
913                        fromColumnsLower,fromColumnsUpper);
914  modelDual->setObjectiveOffset(objOffset);
915  modelDual->setDualBound(model2->dualBound());
916  modelDual->setInfeasibilityCost(model2->infeasibilityCost());
917  modelDual->setDualTolerance(model2->dualTolerance());
918  modelDual->setPrimalTolerance(model2->primalTolerance());
919  modelDual->setPerturbation(model2->perturbation());
920  modelDual->setSpecialOptions(model2->specialOptions());
921  modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
922  delete [] fromRowsLower;
923  delete [] fromRowsUpper;
924  delete [] fromColumnsLower;
925  delete [] fromColumnsUpper;
926  delete [] newObjective;
927  delete [] which;
928  delete [] rowLower;
929  delete [] rowUpper;
930  if (changed)
931    delete model2;
932  modelDual->createStatus();
933  return modelDual;
934}
935// Restores solution from dualized problem
936int
937ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
938{
939  int returnCode=0;;
940  createStatus();
941  // Number of rows in dual problem was original number of columns
942  assert (numberColumns_==dualProblem->numberRows());
943  // If slack on d-row basic then column at bound otherwise column basic
944  // If d-column basic then rhs tight
945  int numberBasic=0;
946  int iRow,iColumn=0;
947  // Get number of extra rows from ranges
948  int numberExtraRows=0;
949  for (iRow=0;iRow<numberRows_;iRow++) {
950    if (rowLower_[iRow]>-1.0e20&&
951        rowUpper_[iRow]<1.0e20) {
952      if (rowUpper_[iRow]!=rowLower_[iRow])
953         numberExtraRows++;
954    }
955  }
956  const double * objective = this->objective();
957  const double * dualDual = dualProblem->dualRowSolution();
958  const double * dualDj = dualProblem->dualColumnSolution();
959  const double * dualSol = dualProblem->primalColumnSolution();
960  const double * dualActs = dualProblem->primalRowSolution();
961#if 0
962  ClpSimplex thisCopy=*this;
963  thisCopy.dual(); // for testing
964  const double * primalDual = thisCopy.dualRowSolution();
965  const double * primalDj = thisCopy.dualColumnSolution();
966  const double * primalSol = thisCopy.primalColumnSolution();
967  const double * primalActs = thisCopy.primalRowSolution();
968  char ss[]={'F','B','U','L','S','F'};
969  printf ("Dual problem row info %d rows\n",dualProblem->numberRows());
970  for (iRow=0;iRow<dualProblem->numberRows();iRow++)
971    printf("%d at %c primal %g dual %g\n",
972           iRow,ss[dualProblem->getRowStatus(iRow)],
973           dualActs[iRow],dualDual[iRow]);
974  printf ("Dual problem column info %d columns\n",dualProblem->numberColumns());
975  for (iColumn=0;iColumn<dualProblem->numberColumns();iColumn++)
976    printf("%d at %c primal %g dual %g\n",
977           iColumn,ss[dualProblem->getColumnStatus(iColumn)],
978           dualSol[iColumn],dualDj[iColumn]);
979  printf ("Primal problem row info %d rows\n",thisCopy.numberRows());
980  for (iRow=0;iRow<thisCopy.numberRows();iRow++)
981    printf("%d at %c primal %g dual %g\n",
982           iRow,ss[thisCopy.getRowStatus(iRow)],
983           primalActs[iRow],primalDual[iRow]);
984  printf ("Primal problem column info %d columns\n",thisCopy.numberColumns());
985  for (iColumn=0;iColumn<thisCopy.numberColumns();iColumn++)
986    printf("%d at %c primal %g dual %g\n",
987           iColumn,ss[thisCopy.getColumnStatus(iColumn)],
988           primalSol[iColumn],primalDj[iColumn]);
989#endif
990  // position at bound information
991  int jColumn=numberRows_;
992  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
993    double objValue =optimizationDirection_*objective[iColumn];
994    Status status = dualProblem->getRowStatus(iColumn);
995    double otherValue = COIN_DBL_MAX;
996    if (columnUpper_[iColumn]<1.0e20&&
997        columnLower_[iColumn]>-1.0e20) {
998        if (fabs(columnLower_[iColumn])<fabs(columnUpper_[iColumn])) {
999          otherValue = columnUpper_[iColumn]+dualDj[jColumn];
1000        } else {
1001          otherValue = columnLower_[iColumn]+dualDj[jColumn];
1002        }
1003        jColumn++;
1004    }
1005    if (status==basic) {
1006      // column is at bound
1007      if (otherValue==COIN_DBL_MAX) {
1008        reducedCost_[iColumn]=objValue-dualActs[iColumn];
1009        if (columnUpper_[iColumn]>1.0e20) {
1010          setColumnStatus(iColumn,atLowerBound);
1011          columnActivity_[iColumn]=columnLower_[iColumn];
1012        } else {
1013          setColumnStatus(iColumn,atUpperBound);
1014          columnActivity_[iColumn]=columnUpper_[iColumn];
1015        }
1016      } else {
1017        reducedCost_[iColumn]=objValue-dualActs[iColumn];
1018        //printf("other dual sol %g\n",otherValue);
1019        if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
1020          setColumnStatus(iColumn,atLowerBound);
1021          columnActivity_[iColumn]=columnLower_[iColumn];
1022        } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
1023          setColumnStatus(iColumn,atUpperBound);
1024          columnActivity_[iColumn]=columnUpper_[iColumn];
1025        } else {
1026          abort();
1027        }
1028      }
1029    } else {
1030      if (otherValue==COIN_DBL_MAX) {
1031        // column basic
1032        setColumnStatus(iColumn,basic);
1033        numberBasic++;
1034        if (columnLower_[iColumn]>-1.0e20) {
1035          columnActivity_[iColumn]=-dualDual[iColumn] + columnLower_[iColumn];
1036        } else if (columnUpper_[iColumn]<1.0e20) {
1037          columnActivity_[iColumn]=-dualDual[iColumn] + columnUpper_[iColumn];
1038        } else {
1039          columnActivity_[iColumn]=-dualDual[iColumn];
1040        }
1041        reducedCost_[iColumn]=0.0;
1042      } else {
1043        // may be at other bound
1044        //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1045        if (dualProblem->getColumnStatus(jColumn-1)!=basic) {
1046          // column basic
1047          setColumnStatus(iColumn,basic);
1048          numberBasic++;
1049          //printf("Col %d otherV %g dualDual %g\n",iColumn,
1050          // otherValue,dualDual[iColumn]);
1051          columnActivity_[iColumn]=-dualDual[iColumn];
1052          columnActivity_[iColumn]=otherValue;
1053          reducedCost_[iColumn]=0.0;
1054        } else {
1055          reducedCost_[iColumn]=objValue-dualActs[iColumn];
1056          if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
1057            setColumnStatus(iColumn,atLowerBound);
1058            columnActivity_[iColumn]=columnLower_[iColumn];
1059          } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
1060            setColumnStatus(iColumn,atUpperBound);
1061            columnActivity_[iColumn]=columnUpper_[iColumn];
1062          } else {
1063            abort();
1064          }
1065        }
1066      }
1067    }
1068  }
1069  // now rows
1070  int kRow=0;
1071  int kExtraRow=jColumn;
1072  int numberRanges=0;
1073  for (iRow=0;iRow<numberRows_;iRow++) {
1074    Status status = dualProblem->getColumnStatus(kRow);
1075    if (status==basic) {
1076      // row is at bound
1077      dual_[iRow]=dualSol[kRow];;
1078    } else {
1079      // row basic
1080      setRowStatus(iRow,basic);
1081      numberBasic++;
1082      dual_[iRow]=0.0;
1083    }
1084    if (rowLower_[iRow]<-1.0e20) {
1085      if (status==basic) {
1086        rowActivity_[iRow]=rowUpper_[iRow];
1087        setRowStatus(iRow,atUpperBound);
1088      } else {
1089        rowActivity_[iRow]=rowUpper_[iRow]+dualSol[kRow];
1090      }       
1091      kRow++;
1092    } else if (rowUpper_[iRow]>1.0e20) {
1093      if (status==basic) {
1094        rowActivity_[iRow]=rowLower_[iRow];
1095        setRowStatus(iRow,atLowerBound);
1096      } else {
1097        rowActivity_[iRow]=rowLower_[iRow]+dualSol[kRow];
1098      }       
1099      kRow++;
1100    } else {
1101      if (rowUpper_[iRow]==rowLower_[iRow]) {
1102        rowActivity_[iRow]=rowLower_[iRow];
1103        if (status==basic) {
1104          setRowStatus(iRow,atLowerBound);
1105        }       
1106        kRow++;
1107      } else {
1108        // range
1109        numberRanges++;
1110        Status statusL = dualProblem->getColumnStatus(kExtraRow);
1111        kExtraRow++;
1112        if (status==basic) {
1113          assert (statusL!=basic);
1114          rowActivity_[iRow]=rowUpper_[iRow];
1115          setRowStatus(iRow,atUpperBound);
1116        } else if (statusL==basic) {
1117          rowActivity_[iRow]=rowLower_[iRow];
1118          setRowStatus(iRow,atLowerBound);
1119        } else {
1120          rowActivity_[iRow]=rowLower_[iRow]+dualSol[kRow];
1121          // row basic
1122          setRowStatus(iRow,basic);
1123          numberBasic++;
1124          dual_[iRow]=0.0;
1125        }
1126        kRow++;
1127      }
1128    }
1129  }
1130  if (numberBasic!=numberRows_) {
1131    printf("Bad basis - ranges - coding needed\n");
1132    assert (numberRanges);
1133    returnCode=1;
1134  }
1135  if (optimizationDirection_<0.0) {
1136    for (iRow=0;iRow<numberRows_;iRow++) {
1137      dual_[iRow]=-dual_[iRow];
1138    }
1139  }
1140  // redo row activities
1141  memset(rowActivity_,0,numberRows_*sizeof(double));
1142  matrix_->times(1.0,columnActivity_,rowActivity_);
1143  // redo reduced costs
1144  memcpy(reducedCost_,this->objective(),numberColumns_*sizeof(double));
1145  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
1146  checkSolutionInternal();
1147  // Below will go to ..DEBUG later
1148  if (returnCode) 
1149    return 1; // Skip checks as will fail at present
1150#ifndef NDEBUG
1151  // Check if correct
1152  double * columnActivity = CoinCopyOfArray(columnActivity_,numberColumns_);
1153  double * rowActivity = CoinCopyOfArray(rowActivity_,numberRows_);
1154  double * reducedCost = CoinCopyOfArray(reducedCost_,numberColumns_);
1155  double * dual = CoinCopyOfArray(dual_,numberRows_);
1156  this->dual();
1157  CoinRelFltEq eq(1.0e-5);
1158  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1159    assert(eq(columnActivity[iColumn],columnActivity_[iColumn]));
1160  }
1161  for (iRow=0;iRow<numberRows_;iRow++) {
1162    assert(eq(dual[iRow],dual_[iRow]));
1163  }
1164  for (iRow=0;iRow<numberRows_;iRow++) {
1165    assert(eq(rowActivity[iRow],rowActivity_[iRow]));
1166  }
1167  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1168    assert(eq(reducedCost[iColumn],reducedCost_[iColumn]));
1169  }
1170  delete [] columnActivity;
1171  delete [] rowActivity;
1172  delete [] reducedCost;
1173  delete [] dual;
1174#endif
1175  return returnCode;
1176}
1177/* Does very cursory presolve.
1178   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1179*/
1180ClpSimplex * 
1181ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1182                        int & nBound, bool moreBounds, bool tightenBounds)
1183{
1184#if 0
1185  //#ifndef NDEBUG
1186  {
1187    int n=0;
1188    int i;
1189    for (i=0;i<numberColumns_;i++)
1190      if (getColumnStatus(i)==ClpSimplex::basic)
1191        n++;
1192    for (i=0;i<numberRows_;i++)
1193      if (getRowStatus(i)==ClpSimplex::basic)
1194        n++;
1195    assert (n==numberRows_);
1196  }
1197#endif
1198 
1199  const double * element = matrix_->getElements();
1200  const int * row = matrix_->getIndices();
1201  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1202  const int * columnLength = matrix_->getVectorLengths();
1203
1204  CoinZeroN(rhs,numberRows_);
1205  int iColumn;
1206  int iRow;
1207  CoinZeroN(whichRow,numberRows_);
1208  int * backColumn = whichColumn+numberColumns_;
1209  int numberRows2=0;
1210  int numberColumns2=0;
1211  double offset=0.0;
1212  const double * objective = this->objective();
1213  double * solution = columnActivity_;
1214  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1215    double lower = columnLower_[iColumn];
1216    double upper = columnUpper_[iColumn];
1217    if (upper>lower||getColumnStatus(iColumn)==ClpSimplex::basic) {
1218      backColumn[iColumn]=numberColumns2;
1219      whichColumn[numberColumns2++]=iColumn;
1220      for (CoinBigIndex j = columnStart[iColumn];
1221           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1222        int iRow = row[j];
1223        int n=whichRow[iRow];
1224        if (n==0&&element[j])
1225          whichRow[iRow]=-iColumn-1;
1226        else if (n<0) 
1227          whichRow[iRow]=2;
1228      }
1229    } else {
1230      // fixed
1231      backColumn[iColumn]=-1;
1232      solution[iColumn]=upper;
1233      if (upper) {
1234        offset += objective[iColumn]*upper;
1235        for (CoinBigIndex j = columnStart[iColumn];
1236             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1237          int iRow = row[j];
1238          double value = element[j];
1239          rhs[iRow] += upper*value;
1240        }
1241      }
1242    }
1243  }
1244  int returnCode=0;
1245  double tolerance = primalTolerance();
1246  nBound=2*numberRows_;
1247  for (iRow=0;iRow<numberRows_;iRow++) {
1248    int n=whichRow[iRow];
1249    if (n>0) {
1250      whichRow[numberRows2++]=iRow;
1251    } else if (n<0) {
1252      //whichRow[numberRows2++]=iRow;
1253      //continue;
1254      // Can only do in certain circumstances as we don't know current value
1255      if (rowLower_[iRow]==rowUpper_[iRow]||getRowStatus(iRow)==ClpSimplex::basic) {
1256        // save row and column for bound
1257        whichRow[--nBound]=iRow;
1258        whichRow[nBound+numberRows_]=-n-1;
1259      } else if (moreBounds) {
1260        // save row and column for bound
1261        whichRow[--nBound]=iRow;
1262        whichRow[nBound+numberRows_]=-n-1;
1263      } else {
1264        whichRow[numberRows2++]=iRow;
1265      }
1266    } else {
1267      // empty
1268      double rhsValue = rhs[iRow];
1269      if (rhsValue<rowLower_[iRow]-tolerance||rhsValue>rowUpper_[iRow]+tolerance) {
1270        returnCode=1; // infeasible
1271      }
1272    }
1273  }
1274  ClpSimplex * small=NULL;
1275  if (!returnCode) {
1276    small = new ClpSimplex(this,numberRows2,whichRow,
1277                     numberColumns2,whichColumn,true,false);
1278    // Set some stuff
1279    small->setDualBound(dualBound_);
1280    small->setInfeasibilityCost(infeasibilityCost_);
1281    small->setSpecialOptions(specialOptions_);
1282    small->setPerturbation(perturbation_);
1283    small->defaultFactorizationFrequency();
1284    small->setAlphaAccuracy(alphaAccuracy_);
1285    // If no rows left then no tightening!
1286    if (!numberRows2||!numberColumns2) 
1287      tightenBounds=false;
1288
1289    int numberElements=getNumElements();
1290    int numberElements2=small->getNumElements();
1291    small->setObjectiveOffset(objectiveOffset()-offset);
1292    handler_->message(CLP_CRUNCH_STATS,messages_)
1293      <<numberRows2<< -(numberRows_ - numberRows2)
1294      <<numberColumns2<< -(numberColumns_ - numberColumns2)
1295      <<numberElements2<< -(numberElements - numberElements2)
1296      <<CoinMessageEol;
1297    // And set objective value to match
1298    small->setObjectiveValue(this->objectiveValue());
1299    double * rowLower2 = small->rowLower();
1300    double * rowUpper2 = small->rowUpper();
1301    int jRow;
1302    for (jRow=0;jRow<numberRows2;jRow++) {
1303      iRow = whichRow[jRow];
1304      if (rowLower2[jRow]>-1.0e20)
1305        rowLower2[jRow] -= rhs[iRow];
1306      if (rowUpper2[jRow]<1.0e20)
1307        rowUpper2[jRow] -= rhs[iRow];
1308    }
1309    // and bounds
1310    double * columnLower2 = small->columnLower();
1311    double * columnUpper2 = small->columnUpper();
1312    const char * integerInformation = integerType_;
1313    for (jRow=nBound;jRow<2*numberRows_;jRow++) {
1314      iRow = whichRow[jRow];
1315      iColumn = whichRow[jRow+numberRows_];
1316      double lowerRow = rowLower_[iRow];
1317      if (lowerRow>-1.0e20)
1318        lowerRow -= rhs[iRow];
1319      double upperRow = rowUpper_[iRow];
1320      if (upperRow<1.0e20)
1321        upperRow -= rhs[iRow];
1322      int jColumn = backColumn[iColumn];
1323      double lower = columnLower2[jColumn];
1324      double upper = columnUpper2[jColumn];
1325      double value=0.0;
1326      for (CoinBigIndex j = columnStart[iColumn];
1327           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1328        if (iRow==row[j]) {
1329          value=element[j];
1330          break;
1331        }
1332      }
1333      assert (value);
1334      // convert rowLower and Upper to implied bounds on column
1335      double newLower=-COIN_DBL_MAX;
1336      double newUpper=COIN_DBL_MAX;
1337      if (value>0.0) {
1338        if (lowerRow>-1.0e20)
1339          newLower = lowerRow/value;
1340        if (upperRow<1.0e20)
1341          newUpper = upperRow/value;
1342      } else {
1343        if (upperRow<1.0e20)
1344          newLower = upperRow/value;
1345        if (lowerRow>-1.0e20)
1346          newUpper = lowerRow/value;
1347      }
1348      if (integerInformation&&integerInformation[iColumn]) {
1349        if (newLower-floor(newLower)<10.0*tolerance) 
1350          newLower=floor(newLower);
1351        else
1352          newLower=ceil(newLower);
1353        if (ceil(newUpper)-newUpper<10.0*tolerance) 
1354          newUpper=ceil(newUpper);
1355        else
1356          newUpper=floor(newUpper);
1357      }
1358      newLower = CoinMax(lower,newLower);
1359      newUpper = CoinMin(upper,newUpper);
1360      if (newLower>newUpper+tolerance) {
1361        //printf("XXYY inf on bound\n");
1362        returnCode=1;
1363      }
1364      columnLower2[jColumn]=newLower;
1365      columnUpper2[jColumn]=CoinMax(newLower,newUpper);
1366      if (getRowStatus(iRow)!=ClpSimplex::basic) {
1367        if (getColumnStatus(iColumn)==ClpSimplex::basic) {
1368          if (columnLower2[jColumn]==columnUpper2[jColumn]) {
1369            // can only get here if will be fixed
1370            small->setColumnStatus(jColumn,ClpSimplex::isFixed);
1371          } else {
1372            // solution is valid
1373            if (fabs(columnActivity_[iColumn]-columnLower2[jColumn])<
1374                fabs(columnActivity_[iColumn]-columnUpper2[jColumn]))
1375              small->setColumnStatus(jColumn,ClpSimplex::atLowerBound);
1376            else
1377              small->setColumnStatus(jColumn,ClpSimplex::atUpperBound);
1378          }
1379        } else {
1380          //printf("what now neither basic\n");
1381        }
1382      }
1383    }
1384    if (returnCode) {
1385      delete small;
1386      small=NULL;
1387    } else if (tightenBounds&&integerInformation) {
1388      // See if we can tighten any bounds
1389      // use rhs for upper and small duals for lower
1390      double * up = rhs;
1391      double * lo = small->dualRowSolution();
1392      const double * element = small->clpMatrix()->getElements();
1393      const int * row = small->clpMatrix()->getIndices();
1394      const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1395      //const int * columnLength = small->clpMatrix()->getVectorLengths();
1396      CoinZeroN(lo,numberRows2);
1397      CoinZeroN(up,numberRows2);
1398      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1399        double upper=columnUpper2[iColumn];
1400        double lower=columnLower2[iColumn];
1401        //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1402        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1403          int iRow=row[j];
1404          double value = element[j];
1405          if (value>0.0) {
1406            if (upper<1.0e20)
1407              up[iRow] += upper*value;
1408            else
1409              up[iRow] = COIN_DBL_MAX;
1410            if (lower>-1.0e20)
1411              lo[iRow] += lower*value;
1412            else
1413              lo[iRow] = -COIN_DBL_MAX;
1414          } else {
1415            if (upper<1.0e20)
1416              lo[iRow] += upper*value;
1417            else
1418              lo[iRow] = -COIN_DBL_MAX;
1419            if (lower>-1.0e20)
1420              up[iRow] += lower*value;
1421            else
1422              up[iRow] = COIN_DBL_MAX;
1423          }
1424        }
1425      }
1426      double * rowLower2 = small->rowLower();
1427      double * rowUpper2 = small->rowUpper();
1428      bool feasible=true;
1429      // make safer
1430      for (int iRow=0;iRow<numberRows2;iRow++) {
1431        double lower = lo[iRow];
1432        if (lower>rowUpper2[iRow]+tolerance) {
1433          feasible=false;
1434          break;
1435        } else {
1436          lo[iRow] = CoinMin(lower-rowUpper2[iRow],0.0)-tolerance;
1437        }
1438        double upper = up[iRow];
1439        if (upper<rowLower2[iRow]-tolerance) {
1440          feasible=false;
1441          break;
1442        } else {
1443          up[iRow] = CoinMax(upper-rowLower2[iRow],0.0)+tolerance;
1444        }
1445      }
1446      if (!feasible) {
1447        delete small;
1448        small=NULL;
1449      } else {
1450        // and tighten
1451        for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1452          if (integerInformation[whichColumn[iColumn]]) {
1453            double upper=columnUpper2[iColumn];
1454            double lower=columnLower2[iColumn];
1455            double newUpper = upper;
1456            double newLower = lower;
1457            double difference = upper-lower;
1458            if (lower>-1000.0&&upper<1000.0) {
1459              for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1460                int iRow=row[j];
1461                double value = element[j];
1462                if (value>0.0) {
1463                  double upWithOut = up[iRow] - value*difference;
1464                  if (upWithOut<0.0) {
1465                    newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
1466                  }
1467                  double lowWithOut = lo[iRow] + value*difference;
1468                  if (lowWithOut>0.0) {
1469                    newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
1470                  }
1471                } else {
1472                  double upWithOut = up[iRow] + value*difference;
1473                  if (upWithOut<0.0) {
1474                    newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
1475                  }
1476                  double lowWithOut = lo[iRow] - value*difference;
1477                  if (lowWithOut>0.0) {
1478                    newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
1479                  }
1480                }
1481              }
1482              if (newLower>lower||newUpper<upper) {
1483                if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
1484                  newUpper = floor(newUpper);
1485                else
1486                  newUpper = floor(newUpper+0.5);
1487                if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
1488                  newLower = ceil(newLower);
1489                else
1490                  newLower = ceil(newLower-0.5);
1491                // change may be too small - check
1492                if (newLower>lower||newUpper<upper) {
1493                  if (newUpper>=newLower) {
1494                    // Could also tighten in this
1495                    //printf("%d bounds %g %g tightened to %g %g\n",
1496                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1497                    //     newLower,newUpper);
1498#if 1
1499                    columnUpper2[iColumn]=newUpper;
1500                    columnLower2[iColumn]=newLower;
1501                    columnUpper_[whichColumn[iColumn]]=newUpper;
1502                    columnLower_[whichColumn[iColumn]]=newLower;
1503#endif
1504                    // and adjust bounds on rows
1505                    newUpper -= upper;
1506                    newLower -= lower;
1507                    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1508                      int iRow=row[j];
1509                      double value = element[j];
1510                      if (value>0.0) {
1511                        up[iRow] += newUpper*value;
1512                        lo[iRow] += newLower*value;
1513                      } else {
1514                        lo[iRow] += newUpper*value;
1515                        up[iRow] += newLower*value;
1516                      }
1517                    }
1518                  } else {
1519                    // infeasible
1520                    //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1521                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1522                    //     newLower,newUpper);
1523#if 1
1524                    delete small;
1525                    small=NULL;
1526                    break;
1527#endif
1528                  }
1529                }
1530              }
1531            }
1532          }
1533        }
1534      }
1535    }
1536  }
1537  return small;
1538}
1539/* After very cursory presolve.
1540   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1541*/
1542void 
1543ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1544                             const int * whichRow, 
1545                             const int * whichColumn, int nBound)
1546{
1547  getbackSolution(small,whichRow,whichColumn);
1548  // and deal with status for bounds
1549  const double * element = matrix_->getElements();
1550  const int * row = matrix_->getIndices();
1551  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1552  const int * columnLength = matrix_->getVectorLengths();
1553  double tolerance = primalTolerance();
1554  double djTolerance = dualTolerance();
1555  for (int jRow=nBound;jRow<2*numberRows_;jRow++) {
1556    int iRow = whichRow[jRow];
1557    int iColumn = whichRow[jRow+numberRows_];
1558    if (getColumnStatus(iColumn)!=ClpSimplex::basic) {
1559      double lower = columnLower_[iColumn];
1560      double upper = columnUpper_[iColumn];
1561      double value = columnActivity_[iColumn];
1562      double djValue = reducedCost_[iColumn];
1563      dual_[iRow]=0.0;
1564      if (upper>lower) {
1565        if (value<lower+tolerance&&djValue>-djTolerance) {
1566          setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1567          setRowStatus(iRow,ClpSimplex::basic);
1568        } else if (value>upper-tolerance&&djValue<djTolerance) {
1569          setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1570          setRowStatus(iRow,ClpSimplex::basic);
1571        } else {
1572          // has to be basic
1573          setColumnStatus(iColumn,ClpSimplex::basic);
1574          reducedCost_[iColumn] = 0.0;
1575          double value=0.0;
1576          for (CoinBigIndex j = columnStart[iColumn];
1577               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1578            if (iRow==row[j]) {
1579              value=element[j];
1580              break;
1581            }
1582          }
1583          dual_[iRow]=djValue/value;
1584          if (rowUpper_[iRow]>rowLower_[iRow]) {
1585            if (fabs(rowActivity_[iRow]-rowLower_[iRow])<
1586                fabs(rowActivity_[iRow]-rowUpper_[iRow]))
1587              setRowStatus(iRow,ClpSimplex::atLowerBound);
1588            else
1589              setRowStatus(iRow,ClpSimplex::atUpperBound);
1590          } else {
1591            setRowStatus(iRow,ClpSimplex::isFixed);
1592          }
1593        }
1594      } else {
1595        // row can always be basic
1596        setRowStatus(iRow,ClpSimplex::basic);
1597      }
1598    } else {
1599      // row can always be basic
1600      setRowStatus(iRow,ClpSimplex::basic);
1601    }
1602  }
1603  //#ifndef NDEBUG
1604#if 0
1605  if  (small.status()==0) {
1606    int n=0;
1607    int i;
1608    for (i=0;i<numberColumns;i++)
1609      if (getColumnStatus(i)==ClpSimplex::basic)
1610        n++;
1611    for (i=0;i<numberRows;i++)
1612      if (getRowStatus(i)==ClpSimplex::basic)
1613        n++;
1614    assert (n==numberRows);
1615  }
1616#endif
1617}
1618/* Tightens integer bounds - returns number tightened or -1 if infeasible
1619 */
1620int 
1621ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1622{
1623  // See if we can tighten any bounds
1624  // use rhs for upper and small duals for lower
1625  double * up = rhsSpace;
1626  double * lo = dual_;
1627  const double * element = matrix_->getElements();
1628  const int * row = matrix_->getIndices();
1629  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1630  const int * columnLength = matrix_->getVectorLengths();
1631  CoinZeroN(lo,numberRows_);
1632  CoinZeroN(up,numberRows_);
1633  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1634    double upper=columnUpper_[iColumn];
1635    double lower=columnLower_[iColumn];
1636    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1637    for (CoinBigIndex j=columnStart[iColumn];
1638         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1639      int iRow=row[j];
1640      double value = element[j];
1641      if (value>0.0) {
1642        if (upper<1.0e20)
1643          up[iRow] += upper*value;
1644        else
1645          up[iRow] = COIN_DBL_MAX;
1646        if (lower>-1.0e20)
1647          lo[iRow] += lower*value;
1648        else
1649          lo[iRow] = -COIN_DBL_MAX;
1650      } else {
1651        if (upper<1.0e20)
1652          lo[iRow] += upper*value;
1653        else
1654          lo[iRow] = -COIN_DBL_MAX;
1655        if (lower>-1.0e20)
1656          up[iRow] += lower*value;
1657        else
1658          up[iRow] = COIN_DBL_MAX;
1659      }
1660    }
1661  }
1662  bool feasible=true;
1663  // make safer
1664  double tolerance = primalTolerance();
1665  for (int iRow=0;iRow<numberRows_;iRow++) {
1666    double lower = lo[iRow];
1667    if (lower>rowUpper_[iRow]+tolerance) {
1668      feasible=false;
1669      break;
1670    } else {
1671      lo[iRow] = CoinMin(lower-rowUpper_[iRow],0.0)-tolerance;
1672    }
1673    double upper = up[iRow];
1674    if (upper<rowLower_[iRow]-tolerance) {
1675      feasible=false;
1676      break;
1677    } else {
1678      up[iRow] = CoinMax(upper-rowLower_[iRow],0.0)+tolerance;
1679    }
1680  }
1681  int numberTightened=0;
1682  if (!feasible) {
1683    return -1;
1684  } else if (integerType_) {
1685    // and tighten
1686    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1687      if (integerType_[iColumn]) {
1688        double upper=columnUpper_[iColumn];
1689        double lower=columnLower_[iColumn];
1690        double newUpper = upper;
1691        double newLower = lower;
1692        double difference = upper-lower;
1693        if (lower>-1000.0&&upper<1000.0) {
1694          for (CoinBigIndex j=columnStart[iColumn];
1695               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1696            int iRow=row[j];
1697            double value = element[j];
1698            if (value>0.0) {
1699              double upWithOut = up[iRow] - value*difference;
1700              if (upWithOut<0.0) {
1701                newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
1702              }
1703              double lowWithOut = lo[iRow] + value*difference;
1704              if (lowWithOut>0.0) {
1705                newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
1706              }
1707            } else {
1708              double upWithOut = up[iRow] + value*difference;
1709              if (upWithOut<0.0) {
1710                newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
1711              }
1712              double lowWithOut = lo[iRow] - value*difference;
1713              if (lowWithOut>0.0) {
1714                newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
1715              }
1716            }
1717          }
1718          if (newLower>lower||newUpper<upper) {
1719            if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
1720              newUpper = floor(newUpper);
1721            else
1722              newUpper = floor(newUpper+0.5);
1723            if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
1724              newLower = ceil(newLower);
1725            else
1726              newLower = ceil(newLower-0.5);
1727            // change may be too small - check
1728            if (newLower>lower||newUpper<upper) {
1729              if (newUpper>=newLower) {
1730                numberTightened++;
1731                //printf("%d bounds %g %g tightened to %g %g\n",
1732                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1733                //     newLower,newUpper);
1734                columnUpper_[iColumn]=newUpper;
1735                columnLower_[iColumn]=newLower;
1736                // and adjust bounds on rows
1737                newUpper -= upper;
1738                newLower -= lower;
1739                for (CoinBigIndex j=columnStart[iColumn];
1740                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
1741                  int iRow=row[j];
1742                  double value = element[j];
1743                  if (value>0.0) {
1744                    up[iRow] += newUpper*value;
1745                    lo[iRow] += newLower*value;
1746                  } else {
1747                    lo[iRow] += newUpper*value;
1748                    up[iRow] += newLower*value;
1749                  }
1750                }
1751              } else {
1752                // infeasible
1753                //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1754                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1755                //     newLower,newUpper);
1756                return -1;
1757              }
1758            }
1759          }
1760        }
1761      }
1762    }
1763  }
1764  return numberTightened;
1765}
1766/* Parametrics
1767   This is an initial slow version.
1768   The code uses current bounds + theta * change (if change array not NULL)
1769   and similarly for objective.
1770   It starts at startingTheta and returns ending theta in endingTheta.
1771   If reportIncrement 0.0 it will report on any movement
1772   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1773   If it can not reach input endingTheta return code will be 1 for infeasible,
1774   2 for unbounded, if error on ranges -1,  otherwise 0.
1775   Normal report is just theta and objective but
1776   if event handler exists it may do more
1777   On exit endingTheta is maximum reached (can be used for next startingTheta)
1778*/
1779int 
1780ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,double reportIncrement,
1781                             const double * changeLowerBound, const double * changeUpperBound,
1782                             const double * changeLowerRhs, const double * changeUpperRhs,
1783                             const double * changeObjective)
1784{
1785  bool needToDoSomething=true;
1786  bool canTryQuick = (reportIncrement) ? true : false;
1787  // Save copy of model
1788  ClpSimplex copyModel = *this;
1789  int savePerturbation = perturbation_;
1790  perturbation_=102; // switch off
1791  while (needToDoSomething) {
1792    needToDoSomething=false;
1793    algorithm_ = -1;
1794   
1795    // save data
1796    ClpDataSave data = saveData();
1797    int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0,NULL,0);
1798    int iRow,iColumn;
1799    double * chgUpper=NULL;
1800    double * chgLower=NULL;
1801    double * chgObjective=NULL;
1802   
1803    // Dantzig (as will not be used) (out later)
1804    ClpDualRowPivot * savePivot = dualRowPivot_;
1805    dualRowPivot_ = new ClpDualRowDantzig();
1806   
1807    if (!returnCode) {
1808      // Find theta when bounds will cross over and create arrays
1809      int numberTotal = numberRows_+numberColumns_;
1810      chgLower = new double[numberTotal];
1811      memset(chgLower,0,numberTotal*sizeof(double));
1812      chgUpper = new double[numberTotal];
1813      memset(chgUpper,0,numberTotal*sizeof(double));
1814      chgObjective = new double[numberTotal];
1815      memset(chgObjective,0,numberTotal*sizeof(double));
1816      assert (!rowScale_);
1817      double maxTheta=1.0e50;
1818      if (changeLowerRhs||changeUpperRhs) {
1819        for (iRow=0;iRow<numberRows_;iRow++) {
1820          double lower = rowLower_[iRow];
1821          double upper = rowUpper_[iRow];
1822          if (lower>upper) {
1823            maxTheta=-1.0;
1824            break;
1825          }
1826          double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
1827          double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
1828          if (lower>-1.0e20&&upper<1.0e20) {
1829            if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
1830              maxTheta = (upper-lower)/(changeLower-changeUpper);
1831            }
1832          }
1833          if (lower>-1.0e20) {
1834            lower_[numberColumns_+iRow] += startingTheta*changeLower;
1835            chgLower[numberColumns_+iRow]=changeLower;
1836          }
1837          if (upper<1.0e20) {
1838            upper_[numberColumns_+iRow] += startingTheta*changeUpper;
1839            chgUpper[numberColumns_+iRow]=changeUpper;
1840          }
1841        }
1842      }
1843      if (maxTheta>0.0) {
1844        if (changeLowerBound||changeUpperBound) {
1845          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1846            double lower = columnLower_[iColumn];
1847            double upper = columnUpper_[iColumn];
1848            if (lower>upper) {
1849              maxTheta=-1.0;
1850              break;
1851            }
1852            double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
1853            double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
1854            if (lower>-1.0e20&&upper<1.0e20) {
1855              if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
1856                maxTheta = (upper-lower)/(changeLower-changeUpper);
1857              }
1858            }
1859            if (lower>-1.0e20) {
1860              lower_[iColumn] += startingTheta*changeLower;
1861              chgLower[iColumn]=changeLower;
1862            }
1863            if (upper<1.0e20) {
1864              upper_[iColumn] += startingTheta*changeUpper;
1865              chgUpper[iColumn]=changeUpper;
1866            }
1867          }
1868        }
1869        if (maxTheta==1.0e50)
1870          maxTheta = COIN_DBL_MAX;
1871      }
1872      if (maxTheta<0.0) {
1873        // bad ranges or initial
1874        returnCode = -1;
1875      }
1876      endingTheta = CoinMin(endingTheta,maxTheta);
1877      if (endingTheta<startingTheta) {
1878        // bad initial
1879        returnCode = -2;
1880      }
1881    }
1882    double saveEndingTheta=endingTheta;
1883    if (!returnCode) {
1884      if (changeObjective) {
1885        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1886          chgObjective[iColumn] = changeObjective[iColumn];
1887          cost_[iColumn] += startingTheta*changeObjective[iColumn];
1888        }
1889      }
1890      double * saveDuals=NULL;
1891      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0,saveDuals,-1,data);
1892      assert (!problemStatus_);
1893      // Now do parametrics
1894      printf("at starting theta of %g objective value is %g\n",startingTheta,
1895             objectiveValue());
1896      while (!returnCode) {
1897        assert (reportIncrement);
1898        returnCode = parametricsLoop(startingTheta,endingTheta,reportIncrement,
1899                                     chgLower,chgUpper,chgObjective,data,
1900                                     canTryQuick);
1901        if (!returnCode) {
1902          //double change = endingTheta-startingTheta;
1903          startingTheta=endingTheta;
1904          endingTheta = saveEndingTheta;
1905          //for (int i=0;i<numberTotal;i++) {
1906          //lower_[i] += change*chgLower[i];
1907          //upper_[i] += change*chgUpper[i];
1908          //cost_[i] += change*chgObjective[i];
1909          //}
1910          printf("at theta of %g objective value is %g\n",startingTheta,
1911                 objectiveValue());
1912          if (startingTheta>=endingTheta)
1913            break;
1914        } else if (returnCode==-1) {
1915          // trouble - do external solve
1916          needToDoSomething=true;
1917        } else {
1918          abort();
1919        }
1920      }
1921    }
1922    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
1923   
1924    delete dualRowPivot_;
1925    dualRowPivot_ = savePivot;
1926    // Restore any saved stuff
1927    restoreData(data);
1928    if (needToDoSomething) {
1929      double saveStartingTheta=startingTheta; // known to be feasible
1930      int cleanedUp=1;
1931      while (cleanedUp) {
1932        // tweak
1933        if (cleanedUp==1) {
1934          if (!reportIncrement)
1935            startingTheta = CoinMin(startingTheta+1.0e-5,saveEndingTheta);
1936          else
1937            startingTheta = CoinMin(startingTheta+reportIncrement,saveEndingTheta);
1938        } else {
1939          // restoring to go slowly
1940          startingTheta=saveStartingTheta;
1941        }
1942        // only works if not scaled
1943        int i;
1944        const double * obj1 = objective();
1945        double * obj2 = copyModel.objective();
1946        const double * lower1 = columnLower_;
1947        double * lower2 = copyModel.columnLower();
1948        const double * upper1 = columnUpper_;
1949        double * upper2 = copyModel.columnUpper();
1950        for (i=0;i<numberColumns_;i++) {
1951          obj2[i] = obj1[i] + startingTheta*chgObjective[i];
1952          lower2[i] = lower1[i] + startingTheta*chgLower[i];
1953          upper2[i] = upper1[i] + startingTheta*chgUpper[i];
1954        }
1955        lower1 = rowLower_;
1956        lower2 = copyModel.rowLower();
1957        upper1 = rowUpper_;
1958        upper2 = copyModel.rowUpper();
1959        for (i=0;i<numberRows_;i++) {
1960          lower2[i] = lower1[i] + startingTheta*chgLower[i+numberColumns_];
1961          upper2[i] = upper1[i] + startingTheta*chgUpper[i+numberColumns_];
1962        }
1963        copyModel.dual();
1964        if (copyModel.problemStatus()) {
1965          printf("Can not get to theta of %g\n",startingTheta);
1966          canTryQuick=false; // do slowly to get exact amount
1967          // back to last known good
1968          if (cleanedUp==1)
1969            cleanedUp=2;
1970          else
1971            abort();
1972        } else {
1973          // and move stuff back
1974          int numberTotal = numberRows_+numberColumns_;
1975          CoinMemcpyN(copyModel.statusArray(),numberTotal,status_);
1976          CoinMemcpyN(copyModel.primalColumnSolution(),numberColumns_,columnActivity_);
1977          CoinMemcpyN(copyModel.primalRowSolution(),numberRows_,rowActivity_);
1978          cleanedUp=0;
1979        }
1980      }
1981    }
1982    delete [] chgLower;
1983    delete [] chgUpper;
1984    delete [] chgObjective;
1985  }
1986  perturbation_ = savePerturbation;
1987  return problemStatus_;
1988}
1989int 
1990ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta,double reportIncrement,
1991                                 const double * changeLower, const double * changeUpper,
1992                                 const double * changeObjective, ClpDataSave & data,
1993                                 bool canTryQuick)
1994{
1995  // stuff is already at starting
1996  // For this crude version just try and go to end
1997  double change=0.0;
1998  if (reportIncrement&&canTryQuick) { 
1999    endingTheta = CoinMin(endingTheta,startingTheta+reportIncrement);
2000    change = endingTheta-startingTheta;
2001  }
2002  int numberTotal = numberRows_+numberColumns_;
2003  int i;
2004  for ( i=0;i<numberTotal;i++) {
2005    lower_[i] += change*changeLower[i];
2006    upper_[i] += change*changeUpper[i];
2007    switch(getStatus(i)) {
2008     
2009    case basic:
2010    case isFree:
2011    case superBasic:
2012      break;
2013    case isFixed:
2014    case atUpperBound:
2015      solution_[i]=upper_[i];
2016      break;
2017    case atLowerBound:
2018      solution_[i]=lower_[i];
2019      break;
2020    }
2021    cost_[i] += change*changeObjective[i];
2022  }
2023  problemStatus_=-1;
2024 
2025  // This says whether to restore things etc
2026  // startup will have factorized so can skip
2027  int factorType=0;
2028  // Start check for cycles
2029  progress_.startCheck();
2030  // Say change made on first iteration
2031  changeMade_=1;
2032  /*
2033    Status of problem:
2034    0 - optimal
2035    1 - infeasible
2036    2 - unbounded
2037    -1 - iterating
2038    -2 - factorization wanted
2039    -3 - redo checking without factorization
2040    -4 - looks infeasible
2041  */
2042  while (problemStatus_<0) {
2043    int iRow, iColumn;
2044    // clear
2045    for (iRow=0;iRow<4;iRow++) {
2046      rowArray_[iRow]->clear();
2047    }   
2048   
2049    for (iColumn=0;iColumn<2;iColumn++) {
2050      columnArray_[iColumn]->clear();
2051    }   
2052   
2053    // give matrix (and model costs and bounds a chance to be
2054    // refreshed (normally null)
2055    matrix_->refresh(this);
2056    // may factorize, checks if problem finished
2057    statusOfProblemInParametrics(factorType,data);
2058    // Say good factorization
2059    factorType=1;
2060    if (data.sparseThreshold_) {
2061      // use default at present
2062      factorization_->sparseThreshold(0);
2063      factorization_->goSparse();
2064    }
2065   
2066    // exit if victory declared
2067    if (problemStatus_>=0)
2068      break;
2069   
2070    // test for maximum iterations
2071    if (hitMaximumIterations()) {
2072      problemStatus_=3;
2073      break;
2074    }
2075    // Check event
2076    {
2077      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2078      if (status>=0) {
2079        problemStatus_=5;
2080        secondaryStatus_=ClpEventHandler::endOfFactorization;
2081        break;
2082      }
2083    }
2084    // Do iterations
2085    if (canTryQuick) {
2086      double * saveDuals=NULL;
2087      reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals,0);
2088    } else {
2089      whileIterating(startingTheta,  endingTheta, reportIncrement,
2090                     changeLower, changeUpper,
2091                     changeObjective);
2092    }
2093  }
2094  if (!problemStatus_) {
2095    theta_=change+startingTheta;
2096    eventHandler_->event(ClpEventHandler::theta);
2097    return 0;
2098  } else if (problemStatus_==10) {
2099    return -1;
2100  } else {
2101    return problemStatus_;
2102  }
2103}
2104/* Checks if finished.  Updates status */
2105void 
2106ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
2107{
2108  if (type==2) {
2109    // trouble - go to recovery
2110    problemStatus_=10;
2111    return;
2112  }
2113  if (problemStatus_>-3||factorization_->pivots()) {
2114    // factorize
2115    // later on we will need to recover from singularities
2116    // also we could skip if first time
2117    if (type) {
2118      // is factorization okay?
2119      if (internalFactorize(1)) {
2120        // trouble - go to recovery
2121        problemStatus_=10;
2122        return;
2123      }
2124    }
2125    if (problemStatus_!=-4||factorization_->pivots()>10)
2126      problemStatus_=-3;
2127  }
2128  // at this stage status is -3 or -4 if looks infeasible
2129  // get primal and dual solutions
2130  gutsOfSolution(NULL,NULL);
2131  double realDualInfeasibilities=sumDualInfeasibilities_;
2132  // If bad accuracy treat as singular
2133  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
2134    // trouble - go to recovery
2135    problemStatus_=10;
2136    return;
2137  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
2138    // Can reduce tolerance
2139    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
2140    factorization_->pivotTolerance(newTolerance);
2141  } 
2142  // Check if looping
2143  int loop;
2144  if (type!=2) 
2145    loop = progress_.looping();
2146  else
2147    loop=-1;
2148  if (loop>=0) {
2149    problemStatus_ = loop; //exit if in loop
2150    if (!problemStatus_) {
2151      // declaring victory
2152      numberPrimalInfeasibilities_ = 0;
2153      sumPrimalInfeasibilities_ = 0.0;
2154    } else {
2155      problemStatus_ = 10; // instead - try other algorithm
2156    }
2157    return;
2158  } else if (loop<-1) {
2159    // something may have changed
2160    gutsOfSolution(NULL,NULL);
2161  }
2162  progressFlag_ = 0; //reset progress flag
2163  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
2164    handler_->message(CLP_SIMPLEX_STATUS,messages_)
2165      <<numberIterations_<<objectiveValue();
2166    handler_->printing(sumPrimalInfeasibilities_>0.0)
2167      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
2168    handler_->printing(sumDualInfeasibilities_>0.0)
2169      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
2170    handler_->printing(numberDualInfeasibilitiesWithoutFree_
2171                       <numberDualInfeasibilities_)
2172                         <<numberDualInfeasibilitiesWithoutFree_;
2173    handler_->message()<<CoinMessageEol;
2174  }
2175  /* If we are primal feasible and any dual infeasibilities are on
2176     free variables then it is better to go to primal */
2177  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
2178      numberDualInfeasibilities_) {
2179    problemStatus_=10;
2180    return;
2181  }
2182 
2183  // check optimal
2184  // give code benefit of doubt
2185  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
2186      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2187    // say optimal (with these bounds etc)
2188    numberDualInfeasibilities_ = 0;
2189    sumDualInfeasibilities_ = 0.0;
2190    numberPrimalInfeasibilities_ = 0;
2191    sumPrimalInfeasibilities_ = 0.0;
2192  }
2193  if (dualFeasible()||problemStatus_==-4) {
2194    progress_.modifyObjective(objectiveValue_
2195                               -sumDualInfeasibilities_*dualBound_);
2196  }
2197  if (numberPrimalInfeasibilities_) {
2198    if (problemStatus_==-4||problemStatus_==-5) {
2199      problemStatus_=1; // infeasible
2200    }
2201  } else if (numberDualInfeasibilities_) {
2202    // clean up
2203    problemStatus_=10;
2204  } else {
2205    problemStatus_=0;
2206  }
2207  lastGoodIteration_ = numberIterations_;
2208  if (problemStatus_<0) {
2209    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
2210    if (sumDualInfeasibilities_)
2211      numberDualInfeasibilities_=1;
2212  }
2213  // Allow matrices to be sorted etc
2214  int fake=-999; // signal sort
2215  matrix_->correctSequence(this,fake,fake);
2216}
2217/* This has the flow between re-factorizations
2218   Reasons to come out:
2219   -1 iterations etc
2220   -2 inaccuracy
2221   -3 slight inaccuracy (and done iterations)
2222   +0 looks optimal (might be unbounded - but we will investigate)
2223   +1 looks infeasible
2224   +3 max iterations
2225*/
2226int 
2227ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta,double reportIncrement,
2228                                const double * changeLower, const double * changeUpper,
2229                                const double * changeObjective)
2230{
2231  {
2232    int i;
2233    for (i=0;i<4;i++) {
2234      rowArray_[i]->clear();
2235    }   
2236    for (i=0;i<2;i++) {
2237      columnArray_[i]->clear();
2238    }   
2239  }     
2240  // if can't trust much and long way from optimal then relax
2241  if (largestPrimalError_>10.0)
2242    factorization_->relaxAccuracyCheck(CoinMin(1.0e2,largestPrimalError_/10.0));
2243  else
2244    factorization_->relaxAccuracyCheck(1.0);
2245  // status stays at -1 while iterating, >=0 finished, -2 to invert
2246  // status -3 to go to top without an invert
2247  int returnCode = -1;
2248  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
2249  double useTheta = startingTheta;
2250  double * primalChange = new double[numberRows_];
2251  double * dualChange = new double[numberColumns_];
2252  int numberTotal = numberColumns_+numberRows_;
2253  int iSequence;
2254  // See if bounds
2255  int type=0;
2256  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2257    if (changeLower[iSequence]||changeUpper[iSequence]) {
2258      type=1;
2259      break;
2260    }
2261  }
2262  // See if objective
2263  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2264    if (changeObjective[iSequence]) {
2265      type |= 2;
2266      break;
2267    }
2268  }
2269  assert (type);
2270  while (problemStatus_==-1) {
2271    double increaseTheta = CoinMin(endingTheta-useTheta,1.0e50);
2272   
2273    // Get theta for bounds - we know can't crossover
2274    int pivotType = nextTheta(type,increaseTheta,primalChange,dualChange,
2275                              changeLower,changeUpper,changeObjective);
2276    if (pivotType)
2277      abort();
2278    // choose row to go out
2279    // dualRow will go to virtual row pivot choice algorithm
2280    reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
2281    if (pivotRow_>=0) {
2282      // we found a pivot row
2283      if (handler_->detail(CLP_SIMPLEX_PIVOTROW,messages_)<100) {
2284        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
2285          <<pivotRow_
2286          <<CoinMessageEol;
2287      }
2288      // check accuracy of weights
2289      dualRowPivot_->checkAccuracy();
2290      // Get good size for pivot
2291      // Allow first few iterations to take tiny
2292      double acceptablePivot=1.0e-9;
2293      if (numberIterations_>100)
2294        acceptablePivot=1.0e-8;
2295      if (factorization_->pivots()>10||
2296          (factorization_->pivots()&&saveSumDual))
2297        acceptablePivot=1.0e-5; // if we have iterated be more strict
2298      else if (factorization_->pivots()>5)
2299        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
2300      else if (factorization_->pivots())
2301        acceptablePivot=1.0e-8; // relax
2302      double bestPossiblePivot=1.0;
2303      // get sign for finding row of tableau
2304      // normal iteration
2305      // create as packed
2306      double direction=directionOut_;
2307      rowArray_[0]->createPacked(1,&pivotRow_,&direction);
2308      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
2309      // put row of tableau in rowArray[0] and columnArray[0]
2310      matrix_->transposeTimes(this,-1.0,
2311                              rowArray_[0],rowArray_[3],columnArray_[0]);
2312      // do ratio test for normal iteration
2313      bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)->dualColumn(rowArray_[0],
2314                                                                columnArray_[0],columnArray_[1],
2315                                                                rowArray_[3],acceptablePivot,NULL);
2316      if (sequenceIn_>=0) {
2317        // normal iteration
2318        // update the incoming column
2319        double btranAlpha = -alpha_*directionOut_; // for check
2320        unpackPacked(rowArray_[1]);
2321        // moved into updateWeights factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2322        // and update dual weights (can do in parallel - with extra array)
2323        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
2324                                              rowArray_[2],
2325                                              rowArray_[3],
2326                                              rowArray_[1]);
2327        // see if update stable
2328#ifdef CLP_DEBUG
2329        if ((handler_->logLevel()&32))
2330          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
2331#endif
2332        double checkValue=1.0e-7;
2333        // if can't trust much and long way from optimal then relax
2334        if (largestPrimalError_>10.0)
2335          checkValue = CoinMin(1.0e-4,1.0e-8*largestPrimalError_);
2336        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
2337            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
2338          handler_->message(CLP_DUAL_CHECK,messages_)
2339            <<btranAlpha
2340            <<alpha_
2341            <<CoinMessageEol;
2342          if (factorization_->pivots()) {
2343            dualRowPivot_->unrollWeights();
2344            problemStatus_=-2; // factorize now
2345            rowArray_[0]->clear();
2346            rowArray_[1]->clear();
2347            columnArray_[0]->clear();
2348            returnCode=-2;
2349            break;
2350          } else {
2351            // take on more relaxed criterion
2352            double test;
2353            if (fabs(btranAlpha)<1.0e-8||fabs(alpha_)<1.0e-8)
2354              test = 1.0e-1*fabs(alpha_);
2355            else
2356              test = 1.0e-4*(1.0+fabs(alpha_));
2357            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
2358                fabs(btranAlpha-alpha_)>test) {
2359              dualRowPivot_->unrollWeights();
2360              // need to reject something
2361              char x = isColumn(sequenceOut_) ? 'C' :'R';
2362              handler_->message(CLP_SIMPLEX_FLAG,messages_)
2363                <<x<<sequenceWithin(sequenceOut_)
2364                <<CoinMessageEol;
2365              setFlagged(sequenceOut_);
2366              progress_.clearBadTimes();
2367              lastBadIteration_ = numberIterations_; // say be more cautious
2368              rowArray_[0]->clear();
2369              rowArray_[1]->clear();
2370              columnArray_[0]->clear();
2371              if (fabs(alpha_)<1.0e-10&&fabs(btranAlpha)<1.0e-8&&numberIterations_>100) {
2372                //printf("I think should declare infeasible\n");
2373                problemStatus_=1;
2374                returnCode=1;
2375                break;
2376              }
2377              continue;
2378            }
2379          }
2380        }
2381        // update duals BEFORE replaceColumn so can do updateColumn
2382        double objectiveChange=0.0;
2383        // do duals first as variables may flip bounds
2384        // rowArray_[0] and columnArray_[0] may have flips
2385        // so use rowArray_[3] for work array from here on
2386        int nswapped = 0;
2387        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
2388        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
2389        nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],
2390                                     rowArray_[2],theta_,
2391                                     objectiveChange,false);
2392
2393        // which will change basic solution
2394        if (nswapped) {
2395          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
2396          dualRowPivot_->updatePrimalSolution(rowArray_[2],
2397                                              1.0,objectiveChange);
2398          // recompute dualOut_
2399          valueOut_ = solution_[sequenceOut_];
2400          if (directionOut_<0) {
2401            dualOut_ = valueOut_ - upperOut_;
2402          } else {
2403            dualOut_ = lowerOut_ - valueOut_;
2404          }
2405        }
2406        // amount primal will move
2407        double movement = -dualOut_*directionOut_/alpha_;
2408        // so objective should increase by fabs(dj)*movement
2409        // but we already have objective change - so check will be good
2410        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
2411#ifdef CLP_DEBUG
2412          if (handler_->logLevel()&32)
2413            printf("movement %g, swap change %g, rest %g  * %g\n",
2414                   objectiveChange+fabs(movement*dualIn_),
2415                   objectiveChange,movement,dualIn_);
2416#endif
2417          if(factorization_->pivots()) {
2418            // going backwards - factorize
2419            dualRowPivot_->unrollWeights();
2420            problemStatus_=-2; // factorize now
2421            returnCode=-2;
2422            break;
2423          }
2424        }
2425        CoinAssert(fabs(dualOut_)<1.0e50);
2426        // if stable replace in basis
2427        int updateStatus = factorization_->replaceColumn(this,
2428                                                         rowArray_[2],
2429                                                         rowArray_[1],
2430                                                         pivotRow_,
2431                                                         alpha_);
2432        // if no pivots, bad update but reasonable alpha - take and invert
2433        if (updateStatus==2&&
2434                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
2435          updateStatus=4;
2436        if (updateStatus==1||updateStatus==4) {
2437          // slight error
2438          if (factorization_->pivots()>5||updateStatus==4) {
2439            problemStatus_=-2; // factorize now
2440            returnCode=-3;
2441          }
2442        } else if (updateStatus==2) {
2443          // major error
2444          dualRowPivot_->unrollWeights();
2445          // later we may need to unwind more e.g. fake bounds
2446          if (factorization_->pivots()) {
2447            problemStatus_=-2; // factorize now
2448            returnCode=-2;
2449            break;
2450          } else {
2451            // need to reject something
2452            char x = isColumn(sequenceOut_) ? 'C' :'R';
2453            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2454              <<x<<sequenceWithin(sequenceOut_)
2455              <<CoinMessageEol;
2456            setFlagged(sequenceOut_);
2457            progress_.clearBadTimes();
2458            lastBadIteration_ = numberIterations_; // say be more cautious
2459            rowArray_[0]->clear();
2460            rowArray_[1]->clear();
2461            columnArray_[0]->clear();
2462            // make sure dual feasible
2463            // look at all rows and columns
2464            double objectiveChange=0.0;
2465            reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
2466                              0.0,objectiveChange,true);
2467            continue;
2468          }
2469        } else if (updateStatus==3) {
2470          // out of memory
2471          // increase space if not many iterations
2472          if (factorization_->pivots()<
2473              0.5*factorization_->maximumPivots()&&
2474              factorization_->pivots()<200)
2475            factorization_->areaFactor(
2476                                       factorization_->areaFactor() * 1.1);
2477          problemStatus_=-2; // factorize now
2478        } else if (updateStatus==5) {
2479          problemStatus_=-2; // factorize now
2480        }
2481        // update primal solution
2482        if (theta_<0.0) {
2483#ifdef CLP_DEBUG
2484          if (handler_->logLevel()&32)
2485            printf("negative theta %g\n",theta_);
2486#endif
2487          theta_=0.0;
2488        }
2489        // do actual flips
2490        reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0],columnArray_[0],theta_);
2491        //rowArray_[1]->expand();
2492        dualRowPivot_->updatePrimalSolution(rowArray_[1],
2493                                            movement,
2494                                            objectiveChange);
2495        // modify dualout
2496        dualOut_ /= alpha_;
2497        dualOut_ *= -directionOut_;
2498        //setStatus(sequenceIn_,basic);
2499        dj_[sequenceIn_]=0.0;
2500        double oldValue=valueIn_;
2501        if (directionIn_==-1) {
2502          // as if from upper bound
2503          valueIn_ = upperIn_+dualOut_;
2504        } else {
2505          // as if from lower bound
2506          valueIn_ = lowerIn_+dualOut_;
2507        }
2508        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
2509        // outgoing
2510        // set dj to zero unless values pass
2511        if (directionOut_>0) {
2512          valueOut_ = lowerOut_;
2513          dj_[sequenceOut_] = theta_;
2514        } else {
2515          valueOut_ = upperOut_;
2516          dj_[sequenceOut_] = -theta_;
2517        }
2518        solution_[sequenceOut_]=valueOut_;
2519        int whatNext=housekeeping(objectiveChange);
2520        // and set bounds correctly
2521        reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_); 
2522        reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
2523        if (whatNext==1) {
2524          problemStatus_ =-2; // refactorize
2525        } else if (whatNext==2) {
2526          // maximum iterations or equivalent
2527          problemStatus_= 3;
2528          returnCode=3;
2529          break;
2530        }
2531        // Check event
2532        {
2533          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2534          if (status>=0) {
2535            problemStatus_=5;
2536            secondaryStatus_=ClpEventHandler::endOfIteration;
2537            returnCode=4;
2538            break;
2539          }
2540        }
2541      } else {
2542        // no incoming column is valid
2543        pivotRow_=-1;
2544#ifdef CLP_DEBUG
2545        if (handler_->logLevel()&32)
2546          printf("** no column pivot\n");
2547#endif
2548        if (factorization_->pivots()<5) {
2549          // If not in branch and bound etc save ray
2550          if ((specialOptions_&(1024|4096))==0) {
2551            // create ray anyway
2552            delete [] ray_;
2553            ray_ = new double [ numberRows_];
2554            rowArray_[0]->expand(); // in case packed
2555            ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
2556          }
2557          // If we have just factorized and infeasibility reasonable say infeas
2558          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
2559            if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
2560                || (specialOptions_&64)==0) {
2561              // say infeasible
2562              problemStatus_=1;
2563              // unless primal feasible!!!!
2564              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
2565              //     numberDualInfeasibilities_,sumDualInfeasibilities_);
2566              if (numberDualInfeasibilities_)
2567                problemStatus_=10;
2568              rowArray_[0]->clear();
2569              columnArray_[0]->clear();
2570              returnCode=1;
2571              break;
2572            }
2573          }
2574          // If special option set - put off as long as possible
2575          if ((specialOptions_&64)==0) {
2576            problemStatus_=-4; //say looks infeasible
2577          } else {
2578            // flag
2579            char x = isColumn(sequenceOut_) ? 'C' :'R';
2580            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2581              <<x<<sequenceWithin(sequenceOut_)
2582              <<CoinMessageEol;
2583            setFlagged(sequenceOut_);
2584            if (!factorization_->pivots()) {
2585              rowArray_[0]->clear();
2586              columnArray_[0]->clear();
2587              continue;
2588            }
2589          }
2590        }
2591        rowArray_[0]->clear();
2592        columnArray_[0]->clear();
2593        returnCode=1;
2594        break;
2595      }
2596    } else {
2597      // no pivot row
2598#ifdef CLP_DEBUG
2599      if (handler_->logLevel()&32)
2600        printf("** no row pivot\n");
2601#endif
2602      int numberPivots = factorization_->pivots();
2603      bool specialCase;
2604      int useNumberFake;
2605      returnCode=0;
2606      if (numberPivots<20&&
2607          (specialOptions_&2048)!=0&&!numberChanged_&&perturbation_>=100
2608          &&dualBound_>1.0e8) {
2609        specialCase=true;
2610        // as dual bound high - should be okay
2611        useNumberFake=0;
2612      } else {
2613        specialCase=false;
2614        useNumberFake=numberFake_;
2615      }
2616      if (!numberPivots||specialCase) {
2617        // may have crept through - so may be optimal
2618        // check any flagged variables
2619        int iRow;
2620        for (iRow=0;iRow<numberRows_;iRow++) {
2621          int iPivot=pivotVariable_[iRow];
2622          if (flagged(iPivot))
2623            break;
2624        }
2625        if (iRow<numberRows_&&numberPivots) {
2626          // try factorization
2627          returnCode=-2;
2628        }
2629       
2630        if (useNumberFake||numberDualInfeasibilities_) {
2631          // may be dual infeasible
2632          problemStatus_=-5;
2633        } else {
2634          if (iRow<numberRows_) {
2635            problemStatus_=-5;
2636          } else {
2637            if (numberPivots) {
2638              // objective may be wrong
2639              objectiveValue_ = innerProduct(cost_,
2640                                                                        numberColumns_+numberRows_,
2641                                                                        solution_);
2642              objectiveValue_ += objective_->nonlinearOffset();
2643              objectiveValue_ /= (objectiveScale_*rhsScale_);
2644              if ((specialOptions_&16384)==0) {
2645                // and dual_ may be wrong (i.e. for fixed or basic)
2646                CoinIndexedVector * arrayVector = rowArray_[1];
2647                arrayVector->clear();
2648                int iRow;
2649                double * array = arrayVector->denseVector();
2650                /* Use dual_ instead of array
2651                   Even though dual_ is only numberRows_ long this is
2652                   okay as gets permuted to longer rowArray_[2]
2653                */
2654                arrayVector->setDenseVector(dual_);
2655                int * index = arrayVector->getIndices();
2656                int number=0;
2657                for (iRow=0;iRow<numberRows_;iRow++) {
2658                  int iPivot=pivotVariable_[iRow];
2659                  double value = cost_[iPivot];
2660                  dual_[iRow]=value;
2661                  if (value) {
2662                    index[number++]=iRow;
2663                  }
2664                }
2665                arrayVector->setNumElements(number);
2666                // Extended duals before "updateTranspose"
2667                matrix_->dualExpanded(this,arrayVector,NULL,0);
2668                // Btran basic costs
2669                rowArray_[2]->clear();
2670                factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
2671                // and return vector
2672                arrayVector->setDenseVector(array);
2673              }
2674            }
2675            problemStatus_=0;
2676            sumPrimalInfeasibilities_=0.0;
2677            if ((specialOptions_&(1024+16384))!=0) {
2678              CoinIndexedVector * arrayVector = rowArray_[1];
2679              arrayVector->clear();
2680              double * rhs = arrayVector->denseVector();
2681              times(1.0,solution_,rhs);
2682              bool bad2=false;
2683              int i;
2684              for ( i=0;i<numberRows_;i++) {
2685                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
2686                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
2687                  bad2=true;
2688                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
2689                }
2690                rhs[i]=0.0;
2691              }
2692              for ( i=0;i<numberColumns_;i++) {
2693                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
2694                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
2695                  bad2=true;
2696                }
2697              }
2698              if (bad2) {
2699                problemStatus_=-3;
2700                returnCode=-2;
2701                // Force to re-factorize early next time
2702                int numberPivots = factorization_->pivots();
2703                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
2704              }
2705            }
2706          }
2707        }
2708      } else {
2709        problemStatus_=-3;
2710        returnCode=-2;
2711        // Force to re-factorize early next time
2712        int numberPivots = factorization_->pivots();
2713        forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
2714      }
2715      break;
2716    }
2717  }
2718  delete [] primalChange;
2719  delete [] dualChange;
2720  return returnCode;
2721}
2722// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
2723int 
2724ClpSimplexOther::nextTheta(int type, double maxTheta, double * primalChange, double * dualChange,
2725                           const double * changeLower, const double * changeUpper,
2726                           const double * changeObjective)
2727{
2728  int numberTotal = numberColumns_+numberRows_;
2729  int iSequence;
2730  int iRow;
2731  theta_=maxTheta;
2732  bool toLower=false;
2733  if ((type&1)!=0) {
2734    // get change
2735    for (iSequence=0;iSequence<numberTotal;iSequence++) {
2736      primalChange[iSequence]=0.0;
2737      switch(getStatus(iSequence)) {
2738       
2739      case basic:
2740      case isFree:
2741      case superBasic:
2742        break;
2743      case isFixed:
2744      case atUpperBound:
2745        primalChange[iSequence]=changeUpper[iSequence];
2746        break;
2747      case atLowerBound:
2748        primalChange[iSequence]=changeLower[iSequence];
2749        break;
2750      }
2751    }
2752    // use array
2753    double * array = rowArray_[1]->denseVector();
2754    times(1.0,primalChange,array);
2755    int * index = rowArray_[1]->getIndices();
2756    int number=0;
2757    for (iRow=0;iRow<numberRows_;iRow++) {
2758      double value = array[iRow];
2759      if (value) {
2760        array[iRow]=value;
2761        index[number++]=iRow;
2762      }
2763    }
2764    // ftran it
2765    rowArray_[1]->setNumElements(number);
2766    factorization_->updateColumn(rowArray_[0],rowArray_[1]);
2767    number=rowArray_[1]->getNumElements();
2768    pivotRow_=-1;
2769    for (iRow=0;iRow<number;iRow++) {
2770      int iPivot = index[iRow];
2771      iSequence = pivotVariable_[iPivot];
2772      // solution value will be sol - theta*alpha
2773      // bounds will be bounds + change *theta
2774      double currentSolution = solution_[iSequence];
2775      double currentLower = lower_[iSequence];
2776      double currentUpper = upper_[iSequence];
2777      double alpha = array[iPivot];
2778      assert (currentSolution>=currentLower-primalTolerance_);
2779      assert (currentSolution<=currentUpper+primalTolerance_);
2780      double thetaCoefficient;
2781      double hitsLower = COIN_DBL_MAX;
2782      thetaCoefficient = changeLower[iSequence]+alpha;
2783      if (fabs(thetaCoefficient)>1.0e-8)
2784        hitsLower = (currentSolution-currentLower)/thetaCoefficient;
2785      if (hitsLower<0.0) {
2786        // does not hit - but should we check further
2787        hitsLower=COIN_DBL_MAX;
2788      }
2789      double hitsUpper = COIN_DBL_MAX;
2790      thetaCoefficient = changeUpper[iSequence]+alpha;
2791      if (fabs(thetaCoefficient)>1.0e-8)
2792        hitsUpper = (currentSolution-currentUpper)/thetaCoefficient;
2793      if (hitsUpper<0.0) {
2794        // does not hit - but should we check further
2795        hitsUpper=COIN_DBL_MAX;
2796      }
2797      if (CoinMin(hitsLower,hitsUpper)<theta_) {
2798        theta_ = CoinMin(hitsLower,hitsUpper);
2799        toLower = hitsLower<hitsUpper;
2800        pivotRow_=iPivot;
2801      }
2802    }
2803  }
2804  if ((type&2)!=0) {
2805    abort();
2806  }
2807  if (pivotRow_>=0) {
2808    sequenceOut_ = pivotVariable_[pivotRow_];
2809    valueOut_ = solution_[sequenceOut_];
2810    lowerOut_ = lower_[sequenceOut_];
2811    upperOut_ = upper_[sequenceOut_];
2812    if (!toLower) {
2813      directionOut_ = -1;
2814      dualOut_ = valueOut_ - upperOut_;
2815    } else if (valueOut_<lowerOut_) {
2816      directionOut_ = 1;
2817      dualOut_ = lowerOut_ - valueOut_;
2818    }
2819    return 0;
2820  } else {
2821    return -1;
2822  }
2823}
2824/* Expands out all possible combinations for a knapsack
2825   If buildObj NULL then just computes space needed - returns number elements
2826   On entry numberOutput is maximum allowed, on exit it is number needed or
2827   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
2828   least space to return values which reconstruct input.
2829   Rows returned will be original rows but no entries will be returned for
2830   any rows all of whose entries are in knapsack.  So up to user to allow for this.
2831   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
2832   in expanded knapsack.  Values in buildRow and buildElement;
2833*/
2834int 
2835ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
2836                                double * buildObj, CoinBigIndex * buildStart,
2837                                int * buildRow, double * buildElement,int reConstruct) const
2838{
2839  int iRow;
2840  int iColumn;
2841  // Get column copy
2842  CoinPackedMatrix * columnCopy = matrix();
2843  // Get a row copy in standard format
2844  CoinPackedMatrix matrixByRow;
2845  matrixByRow.reverseOrderedCopyOf(*columnCopy);
2846  const double * elementByRow = matrixByRow.getElements();
2847  const int * column = matrixByRow.getIndices();
2848  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2849  const int * rowLength = matrixByRow.getVectorLengths();
2850  CoinBigIndex j;
2851  int * whichColumn = new int [numberColumns_];
2852  int * whichRow = new int [numberRows_];
2853  int numJ=0;
2854  // Get what other columns can compensate for
2855  double * lo = new double [numberRows_];
2856  double * high = new double [numberRows_];
2857  {
2858    // Use to get tight column bounds
2859    ClpSimplex tempModel(*this);
2860    tempModel.tightenPrimalBounds(0.0,0,true);
2861    // Now another model without knapsacks
2862    int nCol=0;
2863    for (iRow=0;iRow<numberRows_;iRow++) {
2864      whichRow[iRow]=iRow;
2865    }
2866    for (iColumn=0;iColumn<numberColumns_;iColumn++)
2867      whichColumn[iColumn]=-1;
2868    for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
2869      int iColumn = column[j];
2870      if (columnUpper_[iColumn]>columnLower_[iColumn]) {
2871        whichColumn[iColumn]=0;
2872      } else {
2873        assert (!columnLower_[iColumn]); // fix later
2874      }
2875    }
2876    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2877      if (whichColumn[iColumn]<0)
2878        whichColumn[nCol++]=iColumn;
2879    }
2880    ClpSimplex tempModel2(&tempModel,numberRows_,whichRow,nCol,whichColumn,false,false,false);
2881    // Row copy
2882    CoinPackedMatrix matrixByRow;
2883    matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
2884    const double * elementByRow = matrixByRow.getElements();
2885    const int * column = matrixByRow.getIndices();
2886    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2887    const int * rowLength = matrixByRow.getVectorLengths();
2888    const double * columnLower = tempModel2.getColLower();
2889    const double * columnUpper = tempModel2.getColUpper();
2890    for (iRow = 0; iRow < numberRows_; iRow++) {
2891      lo[iRow]=-COIN_DBL_MAX;
2892      high[iRow]=COIN_DBL_MAX;
2893      if (rowLower_[iRow]>-1.0e20||rowUpper_[iRow]<1.0e20) {
2894
2895        // possible row
2896        int infiniteUpper = 0;
2897        int infiniteLower = 0;
2898        double maximumUp = 0.0;
2899        double maximumDown = 0.0;
2900        CoinBigIndex rStart = rowStart[iRow];
2901        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2902        CoinBigIndex j;
2903        // Compute possible lower and upper ranges
2904
2905        for (j = rStart; j < rEnd; ++j) {
2906          double value=elementByRow[j];
2907          iColumn = column[j];
2908          if (value > 0.0) {
2909            if (columnUpper[iColumn] >= 1.0e20) {
2910              ++infiniteUpper;
2911            } else {
2912              maximumUp += columnUpper[iColumn] * value;
2913            }
2914            if (columnLower[iColumn] <= -1.0e20) {
2915              ++infiniteLower;
2916            } else {
2917              maximumDown += columnLower[iColumn] * value;
2918            }
2919          } else if (value<0.0) {
2920            if (columnUpper[iColumn] >= 1.0e20) {
2921              ++infiniteLower;
2922            } else {
2923              maximumDown += columnUpper[iColumn] * value;
2924            }
2925            if (columnLower[iColumn] <= -1.0e20) {
2926              ++infiniteUpper;
2927            } else {
2928              maximumUp += columnLower[iColumn] * value;
2929            }
2930          }
2931        }
2932        // Build in a margin of error
2933        maximumUp += 1.0e-8*fabs(maximumUp)+1.0e-7;
2934        maximumDown -= 1.0e-8*fabs(maximumDown)+1.0e-7;
2935        // we want to save effective rhs
2936        double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
2937        double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
2938        if (up==COIN_DBL_MAX||rowLower_[iRow]==-COIN_DBL_MAX) {
2939          // However low we go it doesn't matter
2940          lo[iRow]=-COIN_DBL_MAX;
2941        } else {
2942          // If we go below this then can not be feasible
2943          lo[iRow]=rowLower_[iRow]-up;
2944        }
2945        if (down==-COIN_DBL_MAX||rowUpper_[iRow]==COIN_DBL_MAX) {
2946          // However high we go it doesn't matter
2947          high[iRow]=COIN_DBL_MAX;
2948        } else {
2949          // If we go above this then can not be feasible
2950          high[iRow]=rowUpper_[iRow]-down;
2951        }
2952      }
2953    }
2954  }
2955  numJ =0;
2956  for (iColumn=0;iColumn<numberColumns_;iColumn++)
2957    whichColumn[iColumn]=-1;
2958  int * markRow = new int [numberRows_];
2959  for (iRow=0;iRow<numberRows_;iRow++) 
2960    markRow[iRow]=1;
2961  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
2962    int iColumn = column[j];
2963    if (columnUpper_[iColumn]>columnLower_[iColumn]) {
2964      whichColumn[iColumn]=numJ;
2965      numJ++;
2966    }
2967  }
2968  /* mark rows
2969     -n in knapsack and n other variables
2970     1 no entries
2971     n+1000 not involved in knapsack but n entries
2972     0 only in knapsack
2973  */
2974  for (iRow=0;iRow<numberRows_;iRow++) {
2975    int type=1;
2976    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
2977      int iColumn = column[j];
2978      if (whichColumn[iColumn]>=0) {
2979        if (type==1) {
2980          type=0;
2981        } else if (type>0) {
2982          assert (type>1000); 
2983          type = -(type-1000);
2984        }
2985      } else if (type==1) {
2986        type = 1001;
2987      } else if (type<0) {
2988        type --;
2989      } else if (type==0) {
2990        type = -1;
2991      } else {
2992        assert (type>1000);
2993        type++;
2994      }
2995    }
2996    markRow[iRow]=type;
2997    if (type<0&&type>-30&&false)
2998      printf("markrow on row %d is %d\n",iRow,markRow[iRow]);
2999  }
3000  int * bound = new int [numberColumns_+1];
3001  int * stack = new int [numberColumns_+1];
3002  int * flip = new int [numberColumns_+1];
3003  double * offset = new double[numberColumns_+1];
3004  double * size = new double [numberColumns_+1];
3005  double * rhsOffset = new double[numberRows_];
3006  int * build = new int[numberColumns_];
3007  int maxNumber=numberOutput;
3008  numJ=0;
3009  double minSize = rowLower_[knapsackRow];
3010  double maxSize = rowUpper_[knapsackRow];
3011  double knapsackOffset=0.0;
3012  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
3013    int iColumn = column[j];
3014    double lowerColumn=columnLower_[iColumn];
3015    double upperColumn=columnUpper_[iColumn];
3016    if (lowerColumn==upperColumn)
3017      continue;
3018    double gap = upperColumn-lowerColumn;
3019    if (gap>1.0e8)
3020      gap=1.0e8;
3021    assert (fabs(floor(gap+0.5)-gap)<1.0e-5);
3022    whichColumn[numJ]=iColumn;
3023    bound[numJ]=static_cast<int> (gap);
3024    if (elementByRow[j]>0.0) {
3025      flip[numJ]=1;
3026      offset[numJ]=lowerColumn;
3027      size[numJ++]=elementByRow[j];
3028    } else {
3029      flip[numJ]=-1;
3030      offset[numJ]=upperColumn;
3031      size[numJ++]=-elementByRow[j];
3032      lowerColumn = upperColumn;
3033    }
3034    knapsackOffset += elementByRow[j]*lowerColumn;
3035  } 
3036  int jRow;
3037  for (iRow=0;iRow<numberRows_;iRow++)
3038    whichRow[iRow]=iRow;
3039  ClpSimplex smallModel(this,numberRows_,whichRow,numJ,whichColumn,true,true,true);
3040  // modify rhs to allow for nonzero lower bounds
3041  //double * rowLower = smallModel.rowLower();
3042  //double * rowUpper = smallModel.rowUpper();
3043  //const double * columnLower = smallModel.columnLower();
3044  //const double * columnUpper = smallModel.columnUpper();
3045  const CoinPackedMatrix * matrix = smallModel.matrix();
3046  const double * element = matrix->getElements();
3047  const int * row = matrix->getIndices();
3048  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3049  const int * columnLength = matrix->getVectorLengths();
3050  const double * objective = smallModel.objective();
3051  //double objectiveOffset=0.0;
3052  // would use for fixed?
3053  CoinZeroN(rhsOffset,numberRows_);
3054  double * rowActivity = smallModel.primalRowSolution();
3055  CoinZeroN(rowActivity,numberRows_);
3056  maxSize -= knapsackOffset;
3057  minSize -= knapsackOffset;
3058  // now generate
3059  int i;
3060  int iStack=numJ;
3061  for (i=0;i<numJ;i++) {
3062    stack[i]=0;
3063  }
3064  double tooMuch = 10.0*maxSize+10000;
3065  stack[numJ]=1;
3066  size[numJ]=tooMuch;
3067  bound[numJ]=0;
3068  double sum=tooMuch;
3069  // allow for all zero being OK
3070  stack[numJ-1]=-1;
3071  sum -= size[numJ-1];
3072  numberOutput=0;
3073  int nelCreate=0;
3074  /* typeRun is - 0 for initial sizes
3075                  1 for build
3076                  2 for reconstruct
3077  */
3078  int typeRun = buildObj ? 1 : 0;
3079  if (reConstruct>=0) {
3080    assert (buildRow&&buildElement);
3081    typeRun=2;
3082  }
3083  if (typeRun==1)
3084    buildStart[0]=0;
3085  while (iStack>=0) {
3086    if (sum>=minSize&&sum<=maxSize) {
3087      double checkSize =0.0;
3088      bool good=true;
3089      int nRow=0;
3090      double obj=0.0;
3091      CoinZeroN(rowActivity,nRow);
3092      for (iColumn=0;iColumn<numJ;iColumn++) {
3093        int iValue = stack[iColumn];
3094        if (iValue>bound[iColumn]) {
3095          good=false;
3096          break;
3097        } else {
3098          double realValue = offset[iColumn] + flip[iColumn]*iValue;
3099          if (realValue) {
3100            obj += objective[iColumn]*realValue;
3101            for (CoinBigIndex j=columnStart[iColumn];
3102                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
3103              double value = element[j]*realValue;
3104              int kRow = row[j];
3105              if (rowActivity[kRow]) {
3106                rowActivity[kRow] += value;
3107                if (!rowActivity[kRow])
3108                  rowActivity[kRow]=1.0e-100;
3109              } else {
3110                build[nRow++]=kRow;
3111                rowActivity[kRow]=value;
3112              }
3113            }
3114          }
3115        }
3116      }
3117      if (good) {
3118        for (jRow=0;jRow<nRow;jRow++) {
3119          int kRow=build[jRow];
3120          double value = rowActivity[kRow];
3121          if (value>high[kRow]||value<lo[kRow]) {
3122            good=false;
3123            break;
3124          }
3125        }
3126      }
3127      if (good) {
3128        if (typeRun==1) {
3129          buildObj[numberOutput]=obj;
3130          for (jRow=0;jRow<nRow;jRow++) {
3131            int kRow=build[jRow];
3132            double value = rowActivity[kRow];
3133            if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
3134              buildElement[nelCreate]=value;
3135              buildRow[nelCreate++]=kRow;
3136            }
3137          }
3138          buildStart[numberOutput+1]=nelCreate;
3139        } else if (!typeRun) {
3140          for (jRow=0;jRow<nRow;jRow++) {
3141            int kRow=build[jRow];
3142            double value = rowActivity[kRow];
3143            if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
3144              nelCreate++;
3145            }
3146          }
3147        }
3148        if (typeRun==2&&reConstruct==numberOutput) {
3149          // build and exit
3150          nelCreate=0;
3151          for (iColumn=0;iColumn<numJ;iColumn++) {
3152            int iValue = stack[iColumn];
3153            double realValue = offset[iColumn] + flip[iColumn]*iValue;
3154            if (realValue) {
3155              buildRow[nelCreate]=whichColumn[iColumn];
3156              buildElement[nelCreate++]=realValue;
3157            }
3158          }
3159          numberOutput=1;
3160          for (i=0;i<numJ;i++) {
3161            bound[i]=0;
3162          }
3163          break;
3164        }
3165        numberOutput++;
3166        if (numberOutput>maxNumber) {
3167          nelCreate=-numberOutput;
3168          numberOutput=-1;
3169          for (i=0;i<numJ;i++) {
3170            bound[i]=0;
3171          }
3172          break;
3173        } else if (typeRun==1&&numberOutput==maxNumber) {
3174          // On second run
3175          for (i=0;i<numJ;i++) {
3176            bound[i]=0;
3177          }
3178          break;
3179        } 
3180        for (int j=0;j<numJ;j++) {
3181          checkSize += stack[j]*size[j];
3182        }
3183        assert (fabs(sum-checkSize)<1.0e-3);
3184      }
3185      for (jRow=0;jRow<nRow;jRow++) {
3186        int kRow=build[jRow];
3187        rowActivity[kRow]=0.0;
3188      }
3189    }
3190    if (sum>maxSize||stack[iStack]>bound[iStack]) {
3191      sum -= size[iStack]*stack[iStack]; 
3192      stack[iStack--]=0;
3193      if (iStack>=0) {
3194        stack[iStack] ++;
3195        sum += size[iStack];
3196      }
3197    } else {
3198      // must be less
3199      // add to last possible
3200      iStack = numJ-1;
3201      sum += size[iStack];
3202      stack[iStack]++;
3203    }
3204  }
3205  //printf("%d will be created\n",numberOutput);
3206  delete [] whichColumn;
3207  delete [] whichRow;
3208  delete [] bound;
3209  delete [] stack;
3210  delete [] flip;
3211  delete [] size;
3212  delete [] offset;
3213  delete [] rhsOffset;
3214  delete [] build;
3215  delete [] markRow;
3216  delete [] lo;
3217  delete [] high;
3218  return nelCreate;
3219}
Note: See TracBrowser for help on using the repository browser.