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

Last change on this file since 1331 was 1331, 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: 104.3 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          if (columnLower_[iColumn]>-1.0e20) {
1011            if (columnUpper_[iColumn]>columnLower_[iColumn])
1012              setColumnStatus(iColumn,atLowerBound);
1013            else
1014              setColumnStatus(iColumn,isFixed);
1015            columnActivity_[iColumn]=columnLower_[iColumn];
1016          } else {
1017            // free
1018            setColumnStatus(iColumn,isFree);
1019            columnActivity_[iColumn]=0.0;
1020          }
1021        } else {
1022          setColumnStatus(iColumn,atUpperBound);
1023          columnActivity_[iColumn]=columnUpper_[iColumn];
1024        }
1025      } else {
1026        reducedCost_[iColumn]=objValue-dualActs[iColumn];
1027        //printf("other dual sol %g\n",otherValue);
1028        if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
1029          if (columnUpper_[iColumn]>columnLower_[iColumn])
1030            setColumnStatus(iColumn,atLowerBound);
1031          else
1032            setColumnStatus(iColumn,isFixed);
1033          columnActivity_[iColumn]=columnLower_[iColumn];
1034        } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
1035          if (columnUpper_[iColumn]>columnLower_[iColumn])
1036            setColumnStatus(iColumn,atUpperBound);
1037          else
1038            setColumnStatus(iColumn,isFixed);
1039          columnActivity_[iColumn]=columnUpper_[iColumn];
1040        } else {
1041          abort();
1042        }
1043      }
1044    } else {
1045      if (otherValue==COIN_DBL_MAX) {
1046        // column basic
1047        setColumnStatus(iColumn,basic);
1048        numberBasic++;
1049        if (columnLower_[iColumn]>-1.0e20) {
1050          columnActivity_[iColumn]=-dualDual[iColumn] + columnLower_[iColumn];
1051        } else if (columnUpper_[iColumn]<1.0e20) {
1052          columnActivity_[iColumn]=-dualDual[iColumn] + columnUpper_[iColumn];
1053        } else {
1054          columnActivity_[iColumn]=-dualDual[iColumn];
1055        }
1056        reducedCost_[iColumn]=0.0;
1057      } else {
1058        // may be at other bound
1059        //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1060        if (dualProblem->getColumnStatus(jColumn-1)!=basic) {
1061          // column basic
1062          setColumnStatus(iColumn,basic);
1063          numberBasic++;
1064          //printf("Col %d otherV %g dualDual %g\n",iColumn,
1065          // otherValue,dualDual[iColumn]);
1066          columnActivity_[iColumn]=-dualDual[iColumn];
1067          columnActivity_[iColumn]=otherValue;
1068          reducedCost_[iColumn]=0.0;
1069        } else {
1070          reducedCost_[iColumn]=objValue-dualActs[iColumn];
1071          if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
1072            if (columnUpper_[iColumn]>columnLower_[iColumn])
1073              setColumnStatus(iColumn,atLowerBound);
1074            else
1075              setColumnStatus(iColumn,isFixed);
1076            columnActivity_[iColumn]=columnLower_[iColumn];
1077          } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
1078            if (columnUpper_[iColumn]>columnLower_[iColumn])
1079              setColumnStatus(iColumn,atUpperBound);
1080            else
1081              setColumnStatus(iColumn,isFixed);
1082            columnActivity_[iColumn]=columnUpper_[iColumn];
1083          } else {
1084            abort();
1085          }
1086        }
1087      }
1088    }
1089  }
1090  // now rows
1091  int kExtraRow=jColumn;
1092  int numberRanges=0;
1093  for (iRow=0;iRow<numberRows_;iRow++) {
1094    Status status = dualProblem->getColumnStatus(iRow);
1095    if (status==basic) {
1096      // row is at bound
1097      dual_[iRow]=dualSol[iRow];;
1098    } else {
1099      // row basic
1100      setRowStatus(iRow,basic);
1101      numberBasic++;
1102      dual_[iRow]=0.0;
1103    }
1104    if (rowLower_[iRow]<-1.0e20) {
1105      if (status==basic) {
1106        rowActivity_[iRow]=rowUpper_[iRow];
1107        setRowStatus(iRow,atUpperBound);
1108      } else {
1109        assert (dualDj[iRow]<1.0e-5);
1110        rowActivity_[iRow]=rowUpper_[iRow]+dualDj[iRow];
1111      }       
1112    } else if (rowUpper_[iRow]>1.0e20) {
1113      if (status==basic) {
1114        rowActivity_[iRow]=rowLower_[iRow];
1115        setRowStatus(iRow,atLowerBound);
1116      } else {
1117        rowActivity_[iRow]=rowLower_[iRow]+dualDj[iRow];
1118        assert (dualDj[iRow]>-1.0e-5);
1119      }       
1120    } else {
1121      if (rowUpper_[iRow]==rowLower_[iRow]) {
1122        rowActivity_[iRow]=rowLower_[iRow];
1123        if (status==basic) {
1124          setRowStatus(iRow,isFixed);
1125        }       
1126      } else {
1127        // range
1128        numberRanges++;
1129        Status statusL = dualProblem->getColumnStatus(kExtraRow);
1130        //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1131        //     iRow,status,kExtraRow,statusL, dualSol[iRow],
1132        //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1133        if (status==basic) {
1134          assert (statusL!=basic);
1135          rowActivity_[iRow]=rowUpper_[iRow];
1136          setRowStatus(iRow,atUpperBound);
1137        } else if (statusL==basic) {
1138          numberBasic--; // already counted
1139          rowActivity_[iRow]=rowLower_[iRow];
1140          setRowStatus(iRow,atLowerBound);
1141          dual_[iRow]=dualSol[kExtraRow];;
1142        } else {
1143          rowActivity_[iRow]=rowLower_[iRow]-dualDj[iRow];
1144          assert (dualDj[iRow]<1.0e-5);
1145          // row basic
1146          //setRowStatus(iRow,basic);
1147          //numberBasic++;
1148          dual_[iRow]=0.0;
1149        }
1150        kExtraRow++;
1151      }
1152    }
1153  }
1154  if (numberBasic!=numberRows_) {
1155    printf("Bad basis - ranges - coding needed\n");
1156    assert (numberRanges);
1157    abort();
1158  }
1159  if (optimizationDirection_<0.0) {
1160    for (iRow=0;iRow<numberRows_;iRow++) {
1161      dual_[iRow]=-dual_[iRow];
1162    }
1163  }
1164  // redo row activities
1165  memset(rowActivity_,0,numberRows_*sizeof(double));
1166  matrix_->times(1.0,columnActivity_,rowActivity_);
1167  // redo reduced costs
1168  memcpy(reducedCost_,this->objective(),numberColumns_*sizeof(double));
1169  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
1170  checkSolutionInternal();
1171  // Below will go to ..DEBUG later
1172#if 0 //ndef NDEBUG
1173  // Check if correct
1174  double * columnActivity = CoinCopyOfArray(columnActivity_,numberColumns_);
1175  double * rowActivity = CoinCopyOfArray(rowActivity_,numberRows_);
1176  double * reducedCost = CoinCopyOfArray(reducedCost_,numberColumns_);
1177  double * dual = CoinCopyOfArray(dual_,numberRows_);
1178  this->primal();
1179  CoinRelFltEq eq(1.0e-5);
1180  for (iRow=0;iRow<numberRows_;iRow++) {
1181    assert(eq(dual[iRow],dual_[iRow]));
1182  }
1183  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1184    assert(eq(columnActivity[iColumn],columnActivity_[iColumn]));
1185  }
1186  for (iRow=0;iRow<numberRows_;iRow++) {
1187    assert(eq(rowActivity[iRow],rowActivity_[iRow]));
1188  }
1189  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1190    assert(eq(reducedCost[iColumn],reducedCost_[iColumn]));
1191  }
1192  delete [] columnActivity;
1193  delete [] rowActivity;
1194  delete [] reducedCost;
1195  delete [] dual;
1196#endif
1197  return returnCode;
1198}
1199/* Does very cursory presolve.
1200   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
1201*/
1202ClpSimplex * 
1203ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1204                        int & nBound, bool moreBounds, bool tightenBounds)
1205{
1206#if 0
1207  //#ifndef NDEBUG
1208  {
1209    int n=0;
1210    int i;
1211    for (i=0;i<numberColumns_;i++)
1212      if (getColumnStatus(i)==ClpSimplex::basic)
1213        n++;
1214    for (i=0;i<numberRows_;i++)
1215      if (getRowStatus(i)==ClpSimplex::basic)
1216        n++;
1217    assert (n==numberRows_);
1218  }
1219#endif
1220 
1221  const double * element = matrix_->getElements();
1222  const int * row = matrix_->getIndices();
1223  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1224  const int * columnLength = matrix_->getVectorLengths();
1225
1226  CoinZeroN(rhs,numberRows_);
1227  int iColumn;
1228  int iRow;
1229  CoinZeroN(whichRow,numberRows_);
1230  int * backColumn = whichColumn+numberColumns_;
1231  int numberRows2=0;
1232  int numberColumns2=0;
1233  double offset=0.0;
1234  const double * objective = this->objective();
1235  double * solution = columnActivity_;
1236  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1237    double lower = columnLower_[iColumn];
1238    double upper = columnUpper_[iColumn];
1239    if (upper>lower||getColumnStatus(iColumn)==ClpSimplex::basic) {
1240      backColumn[iColumn]=numberColumns2;
1241      whichColumn[numberColumns2++]=iColumn;
1242      for (CoinBigIndex j = columnStart[iColumn];
1243           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1244        int iRow = row[j];
1245        int n=whichRow[iRow];
1246        if (n==0&&element[j])
1247          whichRow[iRow]=-iColumn-1;
1248        else if (n<0) 
1249          whichRow[iRow]=2;
1250      }
1251    } else {
1252      // fixed
1253      backColumn[iColumn]=-1;
1254      solution[iColumn]=upper;
1255      if (upper) {
1256        offset += objective[iColumn]*upper;
1257        for (CoinBigIndex j = columnStart[iColumn];
1258             j<columnStart[iColumn]+columnLength[iColumn];j++) {
1259          int iRow = row[j];
1260          double value = element[j];
1261          rhs[iRow] += upper*value;
1262        }
1263      }
1264    }
1265  }
1266  int returnCode=0;
1267  double tolerance = primalTolerance();
1268  nBound=2*numberRows_;
1269  for (iRow=0;iRow<numberRows_;iRow++) {
1270    int n=whichRow[iRow];
1271    if (n>0) {
1272      whichRow[numberRows2++]=iRow;
1273    } else if (n<0) {
1274      //whichRow[numberRows2++]=iRow;
1275      //continue;
1276      // Can only do in certain circumstances as we don't know current value
1277      if (rowLower_[iRow]==rowUpper_[iRow]||getRowStatus(iRow)==ClpSimplex::basic) {
1278        // save row and column for bound
1279        whichRow[--nBound]=iRow;
1280        whichRow[nBound+numberRows_]=-n-1;
1281      } else if (moreBounds) {
1282        // save row and column for bound
1283        whichRow[--nBound]=iRow;
1284        whichRow[nBound+numberRows_]=-n-1;
1285      } else {
1286        whichRow[numberRows2++]=iRow;
1287      }
1288    } else {
1289      // empty
1290      double rhsValue = rhs[iRow];
1291      if (rhsValue<rowLower_[iRow]-tolerance||rhsValue>rowUpper_[iRow]+tolerance) {
1292        returnCode=1; // infeasible
1293      }
1294    }
1295  }
1296  ClpSimplex * small=NULL;
1297  if (!returnCode) {
1298    small = new ClpSimplex(this,numberRows2,whichRow,
1299                     numberColumns2,whichColumn,true,false);
1300    // Set some stuff
1301    small->setDualBound(dualBound_);
1302    small->setInfeasibilityCost(infeasibilityCost_);
1303    small->setSpecialOptions(specialOptions_);
1304    small->setPerturbation(perturbation_);
1305    small->defaultFactorizationFrequency();
1306    small->setAlphaAccuracy(alphaAccuracy_);
1307    // If no rows left then no tightening!
1308    if (!numberRows2||!numberColumns2) 
1309      tightenBounds=false;
1310
1311    int numberElements=getNumElements();
1312    int numberElements2=small->getNumElements();
1313    small->setObjectiveOffset(objectiveOffset()-offset);
1314    handler_->message(CLP_CRUNCH_STATS,messages_)
1315      <<numberRows2<< -(numberRows_ - numberRows2)
1316      <<numberColumns2<< -(numberColumns_ - numberColumns2)
1317      <<numberElements2<< -(numberElements - numberElements2)
1318      <<CoinMessageEol;
1319    // And set objective value to match
1320    small->setObjectiveValue(this->objectiveValue());
1321    double * rowLower2 = small->rowLower();
1322    double * rowUpper2 = small->rowUpper();
1323    int jRow;
1324    for (jRow=0;jRow<numberRows2;jRow++) {
1325      iRow = whichRow[jRow];
1326      if (rowLower2[jRow]>-1.0e20)
1327        rowLower2[jRow] -= rhs[iRow];
1328      if (rowUpper2[jRow]<1.0e20)
1329        rowUpper2[jRow] -= rhs[iRow];
1330    }
1331    // and bounds
1332    double * columnLower2 = small->columnLower();
1333    double * columnUpper2 = small->columnUpper();
1334    const char * integerInformation = integerType_;
1335    for (jRow=nBound;jRow<2*numberRows_;jRow++) {
1336      iRow = whichRow[jRow];
1337      iColumn = whichRow[jRow+numberRows_];
1338      double lowerRow = rowLower_[iRow];
1339      if (lowerRow>-1.0e20)
1340        lowerRow -= rhs[iRow];
1341      double upperRow = rowUpper_[iRow];
1342      if (upperRow<1.0e20)
1343        upperRow -= rhs[iRow];
1344      int jColumn = backColumn[iColumn];
1345      double lower = columnLower2[jColumn];
1346      double upper = columnUpper2[jColumn];
1347      double value=0.0;
1348      for (CoinBigIndex j = columnStart[iColumn];
1349           j<columnStart[iColumn]+columnLength[iColumn];j++) {
1350        if (iRow==row[j]) {
1351          value=element[j];
1352          break;
1353        }
1354      }
1355      assert (value);
1356      // convert rowLower and Upper to implied bounds on column
1357      double newLower=-COIN_DBL_MAX;
1358      double newUpper=COIN_DBL_MAX;
1359      if (value>0.0) {
1360        if (lowerRow>-1.0e20)
1361          newLower = lowerRow/value;
1362        if (upperRow<1.0e20)
1363          newUpper = upperRow/value;
1364      } else {
1365        if (upperRow<1.0e20)
1366          newLower = upperRow/value;
1367        if (lowerRow>-1.0e20)
1368          newUpper = lowerRow/value;
1369      }
1370      if (integerInformation&&integerInformation[iColumn]) {
1371        if (newLower-floor(newLower)<10.0*tolerance) 
1372          newLower=floor(newLower);
1373        else
1374          newLower=ceil(newLower);
1375        if (ceil(newUpper)-newUpper<10.0*tolerance) 
1376          newUpper=ceil(newUpper);
1377        else
1378          newUpper=floor(newUpper);
1379      }
1380      newLower = CoinMax(lower,newLower);
1381      newUpper = CoinMin(upper,newUpper);
1382      if (newLower>newUpper+tolerance) {
1383        //printf("XXYY inf on bound\n");
1384        returnCode=1;
1385      }
1386      columnLower2[jColumn]=newLower;
1387      columnUpper2[jColumn]=CoinMax(newLower,newUpper);
1388      if (getRowStatus(iRow)!=ClpSimplex::basic) {
1389        if (getColumnStatus(iColumn)==ClpSimplex::basic) {
1390          if (columnLower2[jColumn]==columnUpper2[jColumn]) {
1391            // can only get here if will be fixed
1392            small->setColumnStatus(jColumn,ClpSimplex::isFixed);
1393          } else {
1394            // solution is valid
1395            if (fabs(columnActivity_[iColumn]-columnLower2[jColumn])<
1396                fabs(columnActivity_[iColumn]-columnUpper2[jColumn]))
1397              small->setColumnStatus(jColumn,ClpSimplex::atLowerBound);
1398            else
1399              small->setColumnStatus(jColumn,ClpSimplex::atUpperBound);
1400          }
1401        } else {
1402          //printf("what now neither basic\n");
1403        }
1404      }
1405    }
1406    if (returnCode) {
1407      delete small;
1408      small=NULL;
1409    } else if (tightenBounds&&integerInformation) {
1410      // See if we can tighten any bounds
1411      // use rhs for upper and small duals for lower
1412      double * up = rhs;
1413      double * lo = small->dualRowSolution();
1414      const double * element = small->clpMatrix()->getElements();
1415      const int * row = small->clpMatrix()->getIndices();
1416      const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1417      //const int * columnLength = small->clpMatrix()->getVectorLengths();
1418      CoinZeroN(lo,numberRows2);
1419      CoinZeroN(up,numberRows2);
1420      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1421        double upper=columnUpper2[iColumn];
1422        double lower=columnLower2[iColumn];
1423        //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1424        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1425          int iRow=row[j];
1426          double value = element[j];
1427          if (value>0.0) {
1428            if (upper<1.0e20)
1429              up[iRow] += upper*value;
1430            else
1431              up[iRow] = COIN_DBL_MAX;
1432            if (lower>-1.0e20)
1433              lo[iRow] += lower*value;
1434            else
1435              lo[iRow] = -COIN_DBL_MAX;
1436          } else {
1437            if (upper<1.0e20)
1438              lo[iRow] += upper*value;
1439            else
1440              lo[iRow] = -COIN_DBL_MAX;
1441            if (lower>-1.0e20)
1442              up[iRow] += lower*value;
1443            else
1444              up[iRow] = COIN_DBL_MAX;
1445          }
1446        }
1447      }
1448      double * rowLower2 = small->rowLower();
1449      double * rowUpper2 = small->rowUpper();
1450      bool feasible=true;
1451      // make safer
1452      for (int iRow=0;iRow<numberRows2;iRow++) {
1453        double lower = lo[iRow];
1454        if (lower>rowUpper2[iRow]+tolerance) {
1455          feasible=false;
1456          break;
1457        } else {
1458          lo[iRow] = CoinMin(lower-rowUpper2[iRow],0.0)-tolerance;
1459        }
1460        double upper = up[iRow];
1461        if (upper<rowLower2[iRow]-tolerance) {
1462          feasible=false;
1463          break;
1464        } else {
1465          up[iRow] = CoinMax(upper-rowLower2[iRow],0.0)+tolerance;
1466        }
1467      }
1468      if (!feasible) {
1469        delete small;
1470        small=NULL;
1471      } else {
1472        // and tighten
1473        for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
1474          if (integerInformation[whichColumn[iColumn]]) {
1475            double upper=columnUpper2[iColumn];
1476            double lower=columnLower2[iColumn];
1477            double newUpper = upper;
1478            double newLower = lower;
1479            double difference = upper-lower;
1480            if (lower>-1000.0&&upper<1000.0) {
1481              for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1482                int iRow=row[j];
1483                double value = element[j];
1484                if (value>0.0) {
1485                  double upWithOut = up[iRow] - value*difference;
1486                  if (upWithOut<0.0) {
1487                    newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
1488                  }
1489                  double lowWithOut = lo[iRow] + value*difference;
1490                  if (lowWithOut>0.0) {
1491                    newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
1492                  }
1493                } else {
1494                  double upWithOut = up[iRow] + value*difference;
1495                  if (upWithOut<0.0) {
1496                    newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
1497                  }
1498                  double lowWithOut = lo[iRow] - value*difference;
1499                  if (lowWithOut>0.0) {
1500                    newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
1501                  }
1502                }
1503              }
1504              if (newLower>lower||newUpper<upper) {
1505                if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
1506                  newUpper = floor(newUpper);
1507                else
1508                  newUpper = floor(newUpper+0.5);
1509                if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
1510                  newLower = ceil(newLower);
1511                else
1512                  newLower = ceil(newLower-0.5);
1513                // change may be too small - check
1514                if (newLower>lower||newUpper<upper) {
1515                  if (newUpper>=newLower) {
1516                    // Could also tighten in this
1517                    //printf("%d bounds %g %g tightened to %g %g\n",
1518                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1519                    //     newLower,newUpper);
1520#if 1
1521                    columnUpper2[iColumn]=newUpper;
1522                    columnLower2[iColumn]=newLower;
1523                    columnUpper_[whichColumn[iColumn]]=newUpper;
1524                    columnLower_[whichColumn[iColumn]]=newLower;
1525#endif
1526                    // and adjust bounds on rows
1527                    newUpper -= upper;
1528                    newLower -= lower;
1529                    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1530                      int iRow=row[j];
1531                      double value = element[j];
1532                      if (value>0.0) {
1533                        up[iRow] += newUpper*value;
1534                        lo[iRow] += newLower*value;
1535                      } else {
1536                        lo[iRow] += newUpper*value;
1537                        up[iRow] += newLower*value;
1538                      }
1539                    }
1540                  } else {
1541                    // infeasible
1542                    //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1543                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1544                    //     newLower,newUpper);
1545#if 1
1546                    delete small;
1547                    small=NULL;
1548                    break;
1549#endif
1550                  }
1551                }
1552              }
1553            }
1554          }
1555        }
1556      }
1557    }
1558  }
1559  return small;
1560}
1561/* After very cursory presolve.
1562   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
1563*/
1564void 
1565ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1566                             const int * whichRow, 
1567                             const int * whichColumn, int nBound)
1568{
1569  getbackSolution(small,whichRow,whichColumn);
1570  // and deal with status for bounds
1571  const double * element = matrix_->getElements();
1572  const int * row = matrix_->getIndices();
1573  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1574  const int * columnLength = matrix_->getVectorLengths();
1575  double tolerance = primalTolerance();
1576  double djTolerance = dualTolerance();
1577  for (int jRow=nBound;jRow<2*numberRows_;jRow++) {
1578    int iRow = whichRow[jRow];
1579    int iColumn = whichRow[jRow+numberRows_];
1580    if (getColumnStatus(iColumn)!=ClpSimplex::basic) {
1581      double lower = columnLower_[iColumn];
1582      double upper = columnUpper_[iColumn];
1583      double value = columnActivity_[iColumn];
1584      double djValue = reducedCost_[iColumn];
1585      dual_[iRow]=0.0;
1586      if (upper>lower) {
1587        if (value<lower+tolerance&&djValue>-djTolerance) {
1588          setColumnStatus(iColumn,ClpSimplex::atLowerBound);
1589          setRowStatus(iRow,ClpSimplex::basic);
1590        } else if (value>upper-tolerance&&djValue<djTolerance) {
1591          setColumnStatus(iColumn,ClpSimplex::atUpperBound);
1592          setRowStatus(iRow,ClpSimplex::basic);
1593        } else {
1594          // has to be basic
1595          setColumnStatus(iColumn,ClpSimplex::basic);
1596          reducedCost_[iColumn] = 0.0;
1597          double value=0.0;
1598          for (CoinBigIndex j = columnStart[iColumn];
1599               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1600            if (iRow==row[j]) {
1601              value=element[j];
1602              break;
1603            }
1604          }
1605          dual_[iRow]=djValue/value;
1606          if (rowUpper_[iRow]>rowLower_[iRow]) {
1607            if (fabs(rowActivity_[iRow]-rowLower_[iRow])<
1608                fabs(rowActivity_[iRow]-rowUpper_[iRow]))
1609              setRowStatus(iRow,ClpSimplex::atLowerBound);
1610            else
1611              setRowStatus(iRow,ClpSimplex::atUpperBound);
1612          } else {
1613            setRowStatus(iRow,ClpSimplex::isFixed);
1614          }
1615        }
1616      } else {
1617        // row can always be basic
1618        setRowStatus(iRow,ClpSimplex::basic);
1619      }
1620    } else {
1621      // row can always be basic
1622      setRowStatus(iRow,ClpSimplex::basic);
1623    }
1624  }
1625  //#ifndef NDEBUG
1626#if 0
1627  if  (small.status()==0) {
1628    int n=0;
1629    int i;
1630    for (i=0;i<numberColumns;i++)
1631      if (getColumnStatus(i)==ClpSimplex::basic)
1632        n++;
1633    for (i=0;i<numberRows;i++)
1634      if (getRowStatus(i)==ClpSimplex::basic)
1635        n++;
1636    assert (n==numberRows);
1637  }
1638#endif
1639}
1640/* Tightens integer bounds - returns number tightened or -1 if infeasible
1641 */
1642int 
1643ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1644{
1645  // See if we can tighten any bounds
1646  // use rhs for upper and small duals for lower
1647  double * up = rhsSpace;
1648  double * lo = dual_;
1649  const double * element = matrix_->getElements();
1650  const int * row = matrix_->getIndices();
1651  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1652  const int * columnLength = matrix_->getVectorLengths();
1653  CoinZeroN(lo,numberRows_);
1654  CoinZeroN(up,numberRows_);
1655  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1656    double upper=columnUpper_[iColumn];
1657    double lower=columnLower_[iColumn];
1658    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1659    for (CoinBigIndex j=columnStart[iColumn];
1660         j<columnStart[iColumn]+columnLength[iColumn];j++) {
1661      int iRow=row[j];
1662      double value = element[j];
1663      if (value>0.0) {
1664        if (upper<1.0e20)
1665          up[iRow] += upper*value;
1666        else
1667          up[iRow] = COIN_DBL_MAX;
1668        if (lower>-1.0e20)
1669          lo[iRow] += lower*value;
1670        else
1671          lo[iRow] = -COIN_DBL_MAX;
1672      } else {
1673        if (upper<1.0e20)
1674          lo[iRow] += upper*value;
1675        else
1676          lo[iRow] = -COIN_DBL_MAX;
1677        if (lower>-1.0e20)
1678          up[iRow] += lower*value;
1679        else
1680          up[iRow] = COIN_DBL_MAX;
1681      }
1682    }
1683  }
1684  bool feasible=true;
1685  // make safer
1686  double tolerance = primalTolerance();
1687  for (int iRow=0;iRow<numberRows_;iRow++) {
1688    double lower = lo[iRow];
1689    if (lower>rowUpper_[iRow]+tolerance) {
1690      feasible=false;
1691      break;
1692    } else {
1693      lo[iRow] = CoinMin(lower-rowUpper_[iRow],0.0)-tolerance;
1694    }
1695    double upper = up[iRow];
1696    if (upper<rowLower_[iRow]-tolerance) {
1697      feasible=false;
1698      break;
1699    } else {
1700      up[iRow] = CoinMax(upper-rowLower_[iRow],0.0)+tolerance;
1701    }
1702  }
1703  int numberTightened=0;
1704  if (!feasible) {
1705    return -1;
1706  } else if (integerType_) {
1707    // and tighten
1708    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1709      if (integerType_[iColumn]) {
1710        double upper=columnUpper_[iColumn];
1711        double lower=columnLower_[iColumn];
1712        double newUpper = upper;
1713        double newLower = lower;
1714        double difference = upper-lower;
1715        if (lower>-1000.0&&upper<1000.0) {
1716          for (CoinBigIndex j=columnStart[iColumn];
1717               j<columnStart[iColumn]+columnLength[iColumn];j++) {
1718            int iRow=row[j];
1719            double value = element[j];
1720            if (value>0.0) {
1721              double upWithOut = up[iRow] - value*difference;
1722              if (upWithOut<0.0) {
1723                newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
1724              }
1725              double lowWithOut = lo[iRow] + value*difference;
1726              if (lowWithOut>0.0) {
1727                newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
1728              }
1729            } else {
1730              double upWithOut = up[iRow] + value*difference;
1731              if (upWithOut<0.0) {
1732                newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
1733              }
1734              double lowWithOut = lo[iRow] - value*difference;
1735              if (lowWithOut>0.0) {
1736                newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
1737              }
1738            }
1739          }
1740          if (newLower>lower||newUpper<upper) {
1741            if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
1742              newUpper = floor(newUpper);
1743            else
1744              newUpper = floor(newUpper+0.5);
1745            if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
1746              newLower = ceil(newLower);
1747            else
1748              newLower = ceil(newLower-0.5);
1749            // change may be too small - check
1750            if (newLower>lower||newUpper<upper) {
1751              if (newUpper>=newLower) {
1752                numberTightened++;
1753                //printf("%d bounds %g %g tightened to %g %g\n",
1754                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1755                //     newLower,newUpper);
1756                columnUpper_[iColumn]=newUpper;
1757                columnLower_[iColumn]=newLower;
1758                // and adjust bounds on rows
1759                newUpper -= upper;
1760                newLower -= lower;
1761                for (CoinBigIndex j=columnStart[iColumn];
1762                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
1763                  int iRow=row[j];
1764                  double value = element[j];
1765                  if (value>0.0) {
1766                    up[iRow] += newUpper*value;
1767                    lo[iRow] += newLower*value;
1768                  } else {
1769                    lo[iRow] += newUpper*value;
1770                    up[iRow] += newLower*value;
1771                  }
1772                }
1773              } else {
1774                // infeasible
1775                //printf("%d bounds infeasible %g %g tightened to %g %g\n",
1776                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1777                //     newLower,newUpper);
1778                return -1;
1779              }
1780            }
1781          }
1782        }
1783      }
1784    }
1785  }
1786  return numberTightened;
1787}
1788/* Parametrics
1789   This is an initial slow version.
1790   The code uses current bounds + theta * change (if change array not NULL)
1791   and similarly for objective.
1792   It starts at startingTheta and returns ending theta in endingTheta.
1793   If reportIncrement 0.0 it will report on any movement
1794   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1795   If it can not reach input endingTheta return code will be 1 for infeasible,
1796   2 for unbounded, if error on ranges -1,  otherwise 0.
1797   Normal report is just theta and objective but
1798   if event handler exists it may do more
1799   On exit endingTheta is maximum reached (can be used for next startingTheta)
1800*/
1801int 
1802ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,double reportIncrement,
1803                             const double * changeLowerBound, const double * changeUpperBound,
1804                             const double * changeLowerRhs, const double * changeUpperRhs,
1805                             const double * changeObjective)
1806{
1807  bool needToDoSomething=true;
1808  bool canTryQuick = (reportIncrement) ? true : false;
1809  // Save copy of model
1810  ClpSimplex copyModel = *this;
1811  int savePerturbation = perturbation_;
1812  perturbation_=102; // switch off
1813  while (needToDoSomething) {
1814    needToDoSomething=false;
1815    algorithm_ = -1;
1816   
1817    // save data
1818    ClpDataSave data = saveData();
1819    int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0,NULL,0);
1820    int iRow,iColumn;
1821    double * chgUpper=NULL;
1822    double * chgLower=NULL;
1823    double * chgObjective=NULL;
1824   
1825    // Dantzig (as will not be used) (out later)
1826    ClpDualRowPivot * savePivot = dualRowPivot_;
1827    dualRowPivot_ = new ClpDualRowDantzig();
1828   
1829    if (!returnCode) {
1830      // Find theta when bounds will cross over and create arrays
1831      int numberTotal = numberRows_+numberColumns_;
1832      chgLower = new double[numberTotal];
1833      memset(chgLower,0,numberTotal*sizeof(double));
1834      chgUpper = new double[numberTotal];
1835      memset(chgUpper,0,numberTotal*sizeof(double));
1836      chgObjective = new double[numberTotal];
1837      memset(chgObjective,0,numberTotal*sizeof(double));
1838      assert (!rowScale_);
1839      double maxTheta=1.0e50;
1840      if (changeLowerRhs||changeUpperRhs) {
1841        for (iRow=0;iRow<numberRows_;iRow++) {
1842          double lower = rowLower_[iRow];
1843          double upper = rowUpper_[iRow];
1844          if (lower>upper) {
1845            maxTheta=-1.0;
1846            break;
1847          }
1848          double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
1849          double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
1850          if (lower>-1.0e20&&upper<1.0e20) {
1851            if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
1852              maxTheta = (upper-lower)/(changeLower-changeUpper);
1853            }
1854          }
1855          if (lower>-1.0e20) {
1856            lower_[numberColumns_+iRow] += startingTheta*changeLower;
1857            chgLower[numberColumns_+iRow]=changeLower;
1858          }
1859          if (upper<1.0e20) {
1860            upper_[numberColumns_+iRow] += startingTheta*changeUpper;
1861            chgUpper[numberColumns_+iRow]=changeUpper;
1862          }
1863        }
1864      }
1865      if (maxTheta>0.0) {
1866        if (changeLowerBound||changeUpperBound) {
1867          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1868            double lower = columnLower_[iColumn];
1869            double upper = columnUpper_[iColumn];
1870            if (lower>upper) {
1871              maxTheta=-1.0;
1872              break;
1873            }
1874            double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
1875            double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
1876            if (lower>-1.0e20&&upper<1.0e20) {
1877              if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
1878                maxTheta = (upper-lower)/(changeLower-changeUpper);
1879              }
1880            }
1881            if (lower>-1.0e20) {
1882              lower_[iColumn] += startingTheta*changeLower;
1883              chgLower[iColumn]=changeLower;
1884            }
1885            if (upper<1.0e20) {
1886              upper_[iColumn] += startingTheta*changeUpper;
1887              chgUpper[iColumn]=changeUpper;
1888            }
1889          }
1890        }
1891        if (maxTheta==1.0e50)
1892          maxTheta = COIN_DBL_MAX;
1893      }
1894      if (maxTheta<0.0) {
1895        // bad ranges or initial
1896        returnCode = -1;
1897      }
1898      endingTheta = CoinMin(endingTheta,maxTheta);
1899      if (endingTheta<startingTheta) {
1900        // bad initial
1901        returnCode = -2;
1902      }
1903    }
1904    double saveEndingTheta=endingTheta;
1905    if (!returnCode) {
1906      if (changeObjective) {
1907        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1908          chgObjective[iColumn] = changeObjective[iColumn];
1909          cost_[iColumn] += startingTheta*changeObjective[iColumn];
1910        }
1911      }
1912      double * saveDuals=NULL;
1913      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0,saveDuals,-1,data);
1914      assert (!problemStatus_);
1915      // Now do parametrics
1916      printf("at starting theta of %g objective value is %g\n",startingTheta,
1917             objectiveValue());
1918      while (!returnCode) {
1919        assert (reportIncrement);
1920        returnCode = parametricsLoop(startingTheta,endingTheta,reportIncrement,
1921                                     chgLower,chgUpper,chgObjective,data,
1922                                     canTryQuick);
1923        if (!returnCode) {
1924          //double change = endingTheta-startingTheta;
1925          startingTheta=endingTheta;
1926          endingTheta = saveEndingTheta;
1927          //for (int i=0;i<numberTotal;i++) {
1928          //lower_[i] += change*chgLower[i];
1929          //upper_[i] += change*chgUpper[i];
1930          //cost_[i] += change*chgObjective[i];
1931          //}
1932          printf("at theta of %g objective value is %g\n",startingTheta,
1933                 objectiveValue());
1934          if (startingTheta>=endingTheta)
1935            break;
1936        } else if (returnCode==-1) {
1937          // trouble - do external solve
1938          needToDoSomething=true;
1939        } else {
1940          abort();
1941        }
1942      }
1943    }
1944    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
1945   
1946    delete dualRowPivot_;
1947    dualRowPivot_ = savePivot;
1948    // Restore any saved stuff
1949    restoreData(data);
1950    if (needToDoSomething) {
1951      double saveStartingTheta=startingTheta; // known to be feasible
1952      int cleanedUp=1;
1953      while (cleanedUp) {
1954        // tweak
1955        if (cleanedUp==1) {
1956          if (!reportIncrement)
1957            startingTheta = CoinMin(startingTheta+1.0e-5,saveEndingTheta);
1958          else
1959            startingTheta = CoinMin(startingTheta+reportIncrement,saveEndingTheta);
1960        } else {
1961          // restoring to go slowly
1962          startingTheta=saveStartingTheta;
1963        }
1964        // only works if not scaled
1965        int i;
1966        const double * obj1 = objective();
1967        double * obj2 = copyModel.objective();
1968        const double * lower1 = columnLower_;
1969        double * lower2 = copyModel.columnLower();
1970        const double * upper1 = columnUpper_;
1971        double * upper2 = copyModel.columnUpper();
1972        for (i=0;i<numberColumns_;i++) {
1973          obj2[i] = obj1[i] + startingTheta*chgObjective[i];
1974          lower2[i] = lower1[i] + startingTheta*chgLower[i];
1975          upper2[i] = upper1[i] + startingTheta*chgUpper[i];
1976        }
1977        lower1 = rowLower_;
1978        lower2 = copyModel.rowLower();
1979        upper1 = rowUpper_;
1980        upper2 = copyModel.rowUpper();
1981        for (i=0;i<numberRows_;i++) {
1982          lower2[i] = lower1[i] + startingTheta*chgLower[i+numberColumns_];
1983          upper2[i] = upper1[i] + startingTheta*chgUpper[i+numberColumns_];
1984        }
1985        copyModel.dual();
1986        if (copyModel.problemStatus()) {
1987          printf("Can not get to theta of %g\n",startingTheta);
1988          canTryQuick=false; // do slowly to get exact amount
1989          // back to last known good
1990          if (cleanedUp==1)
1991            cleanedUp=2;
1992          else
1993            abort();
1994        } else {
1995          // and move stuff back
1996          int numberTotal = numberRows_+numberColumns_;
1997          CoinMemcpyN(copyModel.statusArray(),numberTotal,status_);
1998          CoinMemcpyN(copyModel.primalColumnSolution(),numberColumns_,columnActivity_);
1999          CoinMemcpyN(copyModel.primalRowSolution(),numberRows_,rowActivity_);
2000          cleanedUp=0;
2001        }
2002      }
2003    }
2004    delete [] chgLower;
2005    delete [] chgUpper;
2006    delete [] chgObjective;
2007  }
2008  perturbation_ = savePerturbation;
2009  return problemStatus_;
2010}
2011int 
2012ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta,double reportIncrement,
2013                                 const double * changeLower, const double * changeUpper,
2014                                 const double * changeObjective, ClpDataSave & data,
2015                                 bool canTryQuick)
2016{
2017  // stuff is already at starting
2018  // For this crude version just try and go to end
2019  double change=0.0;
2020  if (reportIncrement&&canTryQuick) { 
2021    endingTheta = CoinMin(endingTheta,startingTheta+reportIncrement);
2022    change = endingTheta-startingTheta;
2023  }
2024  int numberTotal = numberRows_+numberColumns_;
2025  int i;
2026  for ( i=0;i<numberTotal;i++) {
2027    lower_[i] += change*changeLower[i];
2028    upper_[i] += change*changeUpper[i];
2029    switch(getStatus(i)) {
2030     
2031    case basic:
2032    case isFree:
2033    case superBasic:
2034      break;
2035    case isFixed:
2036    case atUpperBound:
2037      solution_[i]=upper_[i];
2038      break;
2039    case atLowerBound:
2040      solution_[i]=lower_[i];
2041      break;
2042    }
2043    cost_[i] += change*changeObjective[i];
2044  }
2045  problemStatus_=-1;
2046 
2047  // This says whether to restore things etc
2048  // startup will have factorized so can skip
2049  int factorType=0;
2050  // Start check for cycles
2051  progress_.startCheck();
2052  // Say change made on first iteration
2053  changeMade_=1;
2054  /*
2055    Status of problem:
2056    0 - optimal
2057    1 - infeasible
2058    2 - unbounded
2059    -1 - iterating
2060    -2 - factorization wanted
2061    -3 - redo checking without factorization
2062    -4 - looks infeasible
2063  */
2064  while (problemStatus_<0) {
2065    int iRow, iColumn;
2066    // clear
2067    for (iRow=0;iRow<4;iRow++) {
2068      rowArray_[iRow]->clear();
2069    }   
2070   
2071    for (iColumn=0;iColumn<2;iColumn++) {
2072      columnArray_[iColumn]->clear();
2073    }   
2074   
2075    // give matrix (and model costs and bounds a chance to be
2076    // refreshed (normally null)
2077    matrix_->refresh(this);
2078    // may factorize, checks if problem finished
2079    statusOfProblemInParametrics(factorType,data);
2080    // Say good factorization
2081    factorType=1;
2082    if (data.sparseThreshold_) {
2083      // use default at present
2084      factorization_->sparseThreshold(0);
2085      factorization_->goSparse();
2086    }
2087   
2088    // exit if victory declared
2089    if (problemStatus_>=0)
2090      break;
2091   
2092    // test for maximum iterations
2093    if (hitMaximumIterations()) {
2094      problemStatus_=3;
2095      break;
2096    }
2097    // Check event
2098    {
2099      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2100      if (status>=0) {
2101        problemStatus_=5;
2102        secondaryStatus_=ClpEventHandler::endOfFactorization;
2103        break;
2104      }
2105    }
2106    // Do iterations
2107    if (canTryQuick) {
2108      double * saveDuals=NULL;
2109      reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals,0);
2110    } else {
2111      whileIterating(startingTheta,  endingTheta, reportIncrement,
2112                     changeLower, changeUpper,
2113                     changeObjective);
2114    }
2115  }
2116  if (!problemStatus_) {
2117    theta_=change+startingTheta;
2118    eventHandler_->event(ClpEventHandler::theta);
2119    return 0;
2120  } else if (problemStatus_==10) {
2121    return -1;
2122  } else {
2123    return problemStatus_;
2124  }
2125}
2126/* Checks if finished.  Updates status */
2127void 
2128ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
2129{
2130  if (type==2) {
2131    // trouble - go to recovery
2132    problemStatus_=10;
2133    return;
2134  }
2135  if (problemStatus_>-3||factorization_->pivots()) {
2136    // factorize
2137    // later on we will need to recover from singularities
2138    // also we could skip if first time
2139    if (type) {
2140      // is factorization okay?
2141      if (internalFactorize(1)) {
2142        // trouble - go to recovery
2143        problemStatus_=10;
2144        return;
2145      }
2146    }
2147    if (problemStatus_!=-4||factorization_->pivots()>10)
2148      problemStatus_=-3;
2149  }
2150  // at this stage status is -3 or -4 if looks infeasible
2151  // get primal and dual solutions
2152  gutsOfSolution(NULL,NULL);
2153  double realDualInfeasibilities=sumDualInfeasibilities_;
2154  // If bad accuracy treat as singular
2155  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
2156    // trouble - go to recovery
2157    problemStatus_=10;
2158    return;
2159  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
2160    // Can reduce tolerance
2161    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
2162    factorization_->pivotTolerance(newTolerance);
2163  } 
2164  // Check if looping
2165  int loop;
2166  if (type!=2) 
2167    loop = progress_.looping();
2168  else
2169    loop=-1;
2170  if (loop>=0) {
2171    problemStatus_ = loop; //exit if in loop
2172    if (!problemStatus_) {
2173      // declaring victory
2174      numberPrimalInfeasibilities_ = 0;
2175      sumPrimalInfeasibilities_ = 0.0;
2176    } else {
2177      problemStatus_ = 10; // instead - try other algorithm
2178    }
2179    return;
2180  } else if (loop<-1) {
2181    // something may have changed
2182    gutsOfSolution(NULL,NULL);
2183  }
2184  progressFlag_ = 0; //reset progress flag
2185  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
2186    handler_->message(CLP_SIMPLEX_STATUS,messages_)
2187      <<numberIterations_<<objectiveValue();
2188    handler_->printing(sumPrimalInfeasibilities_>0.0)
2189      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
2190    handler_->printing(sumDualInfeasibilities_>0.0)
2191      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
2192    handler_->printing(numberDualInfeasibilitiesWithoutFree_
2193                       <numberDualInfeasibilities_)
2194                         <<numberDualInfeasibilitiesWithoutFree_;
2195    handler_->message()<<CoinMessageEol;
2196  }
2197  /* If we are primal feasible and any dual infeasibilities are on
2198     free variables then it is better to go to primal */
2199  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
2200      numberDualInfeasibilities_) {
2201    problemStatus_=10;
2202    return;
2203  }
2204 
2205  // check optimal
2206  // give code benefit of doubt
2207  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
2208      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
2209    // say optimal (with these bounds etc)
2210    numberDualInfeasibilities_ = 0;
2211    sumDualInfeasibilities_ = 0.0;
2212    numberPrimalInfeasibilities_ = 0;
2213    sumPrimalInfeasibilities_ = 0.0;
2214  }
2215  if (dualFeasible()||problemStatus_==-4) {
2216    progress_.modifyObjective(objectiveValue_
2217                               -sumDualInfeasibilities_*dualBound_);
2218  }
2219  if (numberPrimalInfeasibilities_) {
2220    if (problemStatus_==-4||problemStatus_==-5) {
2221      problemStatus_=1; // infeasible
2222    }
2223  } else if (numberDualInfeasibilities_) {
2224    // clean up
2225    problemStatus_=10;
2226  } else {
2227    problemStatus_=0;
2228  }
2229  lastGoodIteration_ = numberIterations_;
2230  if (problemStatus_<0) {
2231    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
2232    if (sumDualInfeasibilities_)
2233      numberDualInfeasibilities_=1;
2234  }
2235  // Allow matrices to be sorted etc
2236  int fake=-999; // signal sort
2237  matrix_->correctSequence(this,fake,fake);
2238}
2239/* This has the flow between re-factorizations
2240   Reasons to come out:
2241   -1 iterations etc
2242   -2 inaccuracy
2243   -3 slight inaccuracy (and done iterations)
2244   +0 looks optimal (might be unbounded - but we will investigate)
2245   +1 looks infeasible
2246   +3 max iterations
2247*/
2248int 
2249ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta,double reportIncrement,
2250                                const double * changeLower, const double * changeUpper,
2251                                const double * changeObjective)
2252{
2253  {
2254    int i;
2255    for (i=0;i<4;i++) {
2256      rowArray_[i]->clear();
2257    }   
2258    for (i=0;i<2;i++) {
2259      columnArray_[i]->clear();
2260    }   
2261  }     
2262  // if can't trust much and long way from optimal then relax
2263  if (largestPrimalError_>10.0)
2264    factorization_->relaxAccuracyCheck(CoinMin(1.0e2,largestPrimalError_/10.0));
2265  else
2266    factorization_->relaxAccuracyCheck(1.0);
2267  // status stays at -1 while iterating, >=0 finished, -2 to invert
2268  // status -3 to go to top without an invert
2269  int returnCode = -1;
2270  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
2271  double useTheta = startingTheta;
2272  double * primalChange = new double[numberRows_];
2273  double * dualChange = new double[numberColumns_];
2274  int numberTotal = numberColumns_+numberRows_;
2275  int iSequence;
2276  // See if bounds
2277  int type=0;
2278  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2279    if (changeLower[iSequence]||changeUpper[iSequence]) {
2280      type=1;
2281      break;
2282    }
2283  }
2284  // See if objective
2285  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2286    if (changeObjective[iSequence]) {
2287      type |= 2;
2288      break;
2289    }
2290  }
2291  assert (type);
2292  while (problemStatus_==-1) {
2293    double increaseTheta = CoinMin(endingTheta-useTheta,1.0e50);
2294   
2295    // Get theta for bounds - we know can't crossover
2296    int pivotType = nextTheta(type,increaseTheta,primalChange,dualChange,
2297                              changeLower,changeUpper,changeObjective);
2298    if (pivotType)
2299      abort();
2300    // choose row to go out
2301    // dualRow will go to virtual row pivot choice algorithm
2302    reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
2303    if (pivotRow_>=0) {
2304      // we found a pivot row
2305      if (handler_->detail(CLP_SIMPLEX_PIVOTROW,messages_)<100) {
2306        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
2307          <<pivotRow_
2308          <<CoinMessageEol;
2309      }
2310      // check accuracy of weights
2311      dualRowPivot_->checkAccuracy();
2312      // Get good size for pivot
2313      // Allow first few iterations to take tiny
2314      double acceptablePivot=1.0e-9;
2315      if (numberIterations_>100)
2316        acceptablePivot=1.0e-8;
2317      if (factorization_->pivots()>10||
2318          (factorization_->pivots()&&saveSumDual))
2319        acceptablePivot=1.0e-5; // if we have iterated be more strict
2320      else if (factorization_->pivots()>5)
2321        acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
2322      else if (factorization_->pivots())
2323        acceptablePivot=1.0e-8; // relax
2324      double bestPossiblePivot=1.0;
2325      // get sign for finding row of tableau
2326      // normal iteration
2327      // create as packed
2328      double direction=directionOut_;
2329      rowArray_[0]->createPacked(1,&pivotRow_,&direction);
2330      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
2331      // put row of tableau in rowArray[0] and columnArray[0]
2332      matrix_->transposeTimes(this,-1.0,
2333                              rowArray_[0],rowArray_[3],columnArray_[0]);
2334      // do ratio test for normal iteration
2335      bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)->dualColumn(rowArray_[0],
2336                                                                columnArray_[0],columnArray_[1],
2337                                                                rowArray_[3],acceptablePivot,NULL);
2338      if (sequenceIn_>=0) {
2339        // normal iteration
2340        // update the incoming column
2341        double btranAlpha = -alpha_*directionOut_; // for check
2342        unpackPacked(rowArray_[1]);
2343        // moved into updateWeights factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2344        // and update dual weights (can do in parallel - with extra array)
2345        alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
2346                                              rowArray_[2],
2347                                              rowArray_[3],
2348                                              rowArray_[1]);
2349        // see if update stable
2350#ifdef CLP_DEBUG
2351        if ((handler_->logLevel()&32))
2352          printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
2353#endif
2354        double checkValue=1.0e-7;
2355        // if can't trust much and long way from optimal then relax
2356        if (largestPrimalError_>10.0)
2357          checkValue = CoinMin(1.0e-4,1.0e-8*largestPrimalError_);
2358        if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
2359            fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
2360          handler_->message(CLP_DUAL_CHECK,messages_)
2361            <<btranAlpha
2362            <<alpha_
2363            <<CoinMessageEol;
2364          if (factorization_->pivots()) {
2365            dualRowPivot_->unrollWeights();
2366            problemStatus_=-2; // factorize now
2367            rowArray_[0]->clear();
2368            rowArray_[1]->clear();
2369            columnArray_[0]->clear();
2370            returnCode=-2;
2371            break;
2372          } else {
2373            // take on more relaxed criterion
2374            double test;
2375            if (fabs(btranAlpha)<1.0e-8||fabs(alpha_)<1.0e-8)
2376              test = 1.0e-1*fabs(alpha_);
2377            else
2378              test = 1.0e-4*(1.0+fabs(alpha_));
2379            if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
2380                fabs(btranAlpha-alpha_)>test) {
2381              dualRowPivot_->unrollWeights();
2382              // need to reject something
2383              char x = isColumn(sequenceOut_) ? 'C' :'R';
2384              handler_->message(CLP_SIMPLEX_FLAG,messages_)
2385                <<x<<sequenceWithin(sequenceOut_)
2386                <<CoinMessageEol;
2387              setFlagged(sequenceOut_);
2388              progress_.clearBadTimes();
2389              lastBadIteration_ = numberIterations_; // say be more cautious
2390              rowArray_[0]->clear();
2391              rowArray_[1]->clear();
2392              columnArray_[0]->clear();
2393              if (fabs(alpha_)<1.0e-10&&fabs(btranAlpha)<1.0e-8&&numberIterations_>100) {
2394                //printf("I think should declare infeasible\n");
2395                problemStatus_=1;
2396                returnCode=1;
2397                break;
2398              }
2399              continue;
2400            }
2401          }
2402        }
2403        // update duals BEFORE replaceColumn so can do updateColumn
2404        double objectiveChange=0.0;
2405        // do duals first as variables may flip bounds
2406        // rowArray_[0] and columnArray_[0] may have flips
2407        // so use rowArray_[3] for work array from here on
2408        int nswapped = 0;
2409        //rowArray_[0]->cleanAndPackSafe(1.0e-60);
2410        //columnArray_[0]->cleanAndPackSafe(1.0e-60);
2411        nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],
2412                                     rowArray_[2],theta_,
2413                                     objectiveChange,false);
2414
2415        // which will change basic solution
2416        if (nswapped) {
2417          factorization_->updateColumn(rowArray_[3],rowArray_[2]);
2418          dualRowPivot_->updatePrimalSolution(rowArray_[2],
2419                                              1.0,objectiveChange);
2420          // recompute dualOut_
2421          valueOut_ = solution_[sequenceOut_];
2422          if (directionOut_<0) {
2423            dualOut_ = valueOut_ - upperOut_;
2424          } else {
2425            dualOut_ = lowerOut_ - valueOut_;
2426          }
2427        }
2428        // amount primal will move
2429        double movement = -dualOut_*directionOut_/alpha_;
2430        // so objective should increase by fabs(dj)*movement
2431        // but we already have objective change - so check will be good
2432        if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
2433#ifdef CLP_DEBUG
2434          if (handler_->logLevel()&32)
2435            printf("movement %g, swap change %g, rest %g  * %g\n",
2436                   objectiveChange+fabs(movement*dualIn_),
2437                   objectiveChange,movement,dualIn_);
2438#endif
2439          if(factorization_->pivots()) {
2440            // going backwards - factorize
2441            dualRowPivot_->unrollWeights();
2442            problemStatus_=-2; // factorize now
2443            returnCode=-2;
2444            break;
2445          }
2446        }
2447        CoinAssert(fabs(dualOut_)<1.0e50);
2448        // if stable replace in basis
2449        int updateStatus = factorization_->replaceColumn(this,
2450                                                         rowArray_[2],
2451                                                         rowArray_[1],
2452                                                         pivotRow_,
2453                                                         alpha_);
2454        // if no pivots, bad update but reasonable alpha - take and invert
2455        if (updateStatus==2&&
2456                   !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
2457          updateStatus=4;
2458        if (updateStatus==1||updateStatus==4) {
2459          // slight error
2460          if (factorization_->pivots()>5||updateStatus==4) {
2461            problemStatus_=-2; // factorize now
2462            returnCode=-3;
2463          }
2464        } else if (updateStatus==2) {
2465          // major error
2466          dualRowPivot_->unrollWeights();
2467          // later we may need to unwind more e.g. fake bounds
2468          if (factorization_->pivots()) {
2469            problemStatus_=-2; // factorize now
2470            returnCode=-2;
2471            break;
2472          } else {
2473            // need to reject something
2474            char x = isColumn(sequenceOut_) ? 'C' :'R';
2475            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2476              <<x<<sequenceWithin(sequenceOut_)
2477              <<CoinMessageEol;
2478            setFlagged(sequenceOut_);
2479            progress_.clearBadTimes();
2480            lastBadIteration_ = numberIterations_; // say be more cautious
2481            rowArray_[0]->clear();
2482            rowArray_[1]->clear();
2483            columnArray_[0]->clear();
2484            // make sure dual feasible
2485            // look at all rows and columns
2486            double objectiveChange=0.0;
2487            reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
2488                              0.0,objectiveChange,true);
2489            continue;
2490          }
2491        } else if (updateStatus==3) {
2492          // out of memory
2493          // increase space if not many iterations
2494          if (factorization_->pivots()<
2495              0.5*factorization_->maximumPivots()&&
2496              factorization_->pivots()<200)
2497            factorization_->areaFactor(
2498                                       factorization_->areaFactor() * 1.1);
2499          problemStatus_=-2; // factorize now
2500        } else if (updateStatus==5) {
2501          problemStatus_=-2; // factorize now
2502        }
2503        // update primal solution
2504        if (theta_<0.0) {
2505#ifdef CLP_DEBUG
2506          if (handler_->logLevel()&32)
2507            printf("negative theta %g\n",theta_);
2508#endif
2509          theta_=0.0;
2510        }
2511        // do actual flips
2512        reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0],columnArray_[0],theta_);
2513        //rowArray_[1]->expand();
2514        dualRowPivot_->updatePrimalSolution(rowArray_[1],
2515                                            movement,
2516                                            objectiveChange);
2517        // modify dualout
2518        dualOut_ /= alpha_;
2519        dualOut_ *= -directionOut_;
2520        //setStatus(sequenceIn_,basic);
2521        dj_[sequenceIn_]=0.0;
2522        double oldValue=valueIn_;
2523        if (directionIn_==-1) {
2524          // as if from upper bound
2525          valueIn_ = upperIn_+dualOut_;
2526        } else {
2527          // as if from lower bound
2528          valueIn_ = lowerIn_+dualOut_;
2529        }
2530        objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
2531        // outgoing
2532        // set dj to zero unless values pass
2533        if (directionOut_>0) {
2534          valueOut_ = lowerOut_;
2535          dj_[sequenceOut_] = theta_;
2536        } else {
2537          valueOut_ = upperOut_;
2538          dj_[sequenceOut_] = -theta_;
2539        }
2540        solution_[sequenceOut_]=valueOut_;
2541        int whatNext=housekeeping(objectiveChange);
2542        // and set bounds correctly
2543        reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_); 
2544        reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
2545        if (whatNext==1) {
2546          problemStatus_ =-2; // refactorize
2547        } else if (whatNext==2) {
2548          // maximum iterations or equivalent
2549          problemStatus_= 3;
2550          returnCode=3;
2551          break;
2552        }
2553        // Check event
2554        {
2555          int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2556          if (status>=0) {
2557            problemStatus_=5;
2558            secondaryStatus_=ClpEventHandler::endOfIteration;
2559            returnCode=4;
2560            break;
2561          }
2562        }
2563      } else {
2564        // no incoming column is valid
2565        pivotRow_=-1;
2566#ifdef CLP_DEBUG
2567        if (handler_->logLevel()&32)
2568          printf("** no column pivot\n");
2569#endif
2570        if (factorization_->pivots()<5) {
2571          // If not in branch and bound etc save ray
2572          if ((specialOptions_&(1024|4096))==0) {
2573            // create ray anyway
2574            delete [] ray_;
2575            ray_ = new double [ numberRows_];
2576            rowArray_[0]->expand(); // in case packed
2577            ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
2578          }
2579          // If we have just factorized and infeasibility reasonable say infeas
2580          if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
2581            if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
2582                || (specialOptions_&64)==0) {
2583              // say infeasible
2584              problemStatus_=1;
2585              // unless primal feasible!!!!
2586              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
2587              //     numberDualInfeasibilities_,sumDualInfeasibilities_);
2588              if (numberDualInfeasibilities_)
2589                problemStatus_=10;
2590              rowArray_[0]->clear();
2591              columnArray_[0]->clear();
2592              returnCode=1;
2593              break;
2594            }
2595          }
2596          // If special option set - put off as long as possible
2597          if ((specialOptions_&64)==0) {
2598            problemStatus_=-4; //say looks infeasible
2599          } else {
2600            // flag
2601            char x = isColumn(sequenceOut_) ? 'C' :'R';
2602            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2603              <<x<<sequenceWithin(sequenceOut_)
2604              <<CoinMessageEol;
2605            setFlagged(sequenceOut_);
2606            if (!factorization_->pivots()) {
2607              rowArray_[0]->clear();
2608              columnArray_[0]->clear();
2609              continue;
2610            }
2611          }
2612        }
2613        rowArray_[0]->clear();
2614        columnArray_[0]->clear();
2615        returnCode=1;
2616        break;
2617      }
2618    } else {
2619      // no pivot row
2620#ifdef CLP_DEBUG
2621      if (handler_->logLevel()&32)
2622        printf("** no row pivot\n");
2623#endif
2624      int numberPivots = factorization_->pivots();
2625      bool specialCase;
2626      int useNumberFake;
2627      returnCode=0;
2628      if (numberPivots<20&&
2629          (specialOptions_&2048)!=0&&!numberChanged_&&perturbation_>=100
2630          &&dualBound_>1.0e8) {
2631        specialCase=true;
2632        // as dual bound high - should be okay
2633        useNumberFake=0;
2634      } else {
2635        specialCase=false;
2636        useNumberFake=numberFake_;
2637      }
2638      if (!numberPivots||specialCase) {
2639        // may have crept through - so may be optimal
2640        // check any flagged variables
2641        int iRow;
2642        for (iRow=0;iRow<numberRows_;iRow++) {
2643          int iPivot=pivotVariable_[iRow];
2644          if (flagged(iPivot))
2645            break;
2646        }
2647        if (iRow<numberRows_&&numberPivots) {
2648          // try factorization
2649          returnCode=-2;
2650        }
2651       
2652        if (useNumberFake||numberDualInfeasibilities_) {
2653          // may be dual infeasible
2654          problemStatus_=-5;
2655        } else {
2656          if (iRow<numberRows_) {
2657            problemStatus_=-5;
2658          } else {
2659            if (numberPivots) {
2660              // objective may be wrong
2661              objectiveValue_ = innerProduct(cost_,
2662                                                                        numberColumns_+numberRows_,
2663                                                                        solution_);
2664              objectiveValue_ += objective_->nonlinearOffset();
2665              objectiveValue_ /= (objectiveScale_*rhsScale_);
2666              if ((specialOptions_&16384)==0) {
2667                // and dual_ may be wrong (i.e. for fixed or basic)
2668                CoinIndexedVector * arrayVector = rowArray_[1];
2669                arrayVector->clear();
2670                int iRow;
2671                double * array = arrayVector->denseVector();
2672                /* Use dual_ instead of array
2673                   Even though dual_ is only numberRows_ long this is
2674                   okay as gets permuted to longer rowArray_[2]
2675                */
2676                arrayVector->setDenseVector(dual_);
2677                int * index = arrayVector->getIndices();
2678                int number=0;
2679                for (iRow=0;iRow<numberRows_;iRow++) {
2680                  int iPivot=pivotVariable_[iRow];
2681                  double value = cost_[iPivot];
2682                  dual_[iRow]=value;
2683                  if (value) {
2684                    index[number++]=iRow;
2685                  }
2686                }
2687                arrayVector->setNumElements(number);
2688                // Extended duals before "updateTranspose"
2689                matrix_->dualExpanded(this,arrayVector,NULL,0);
2690                // Btran basic costs
2691                rowArray_[2]->clear();
2692                factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
2693                // and return vector
2694                arrayVector->setDenseVector(array);
2695              }
2696            }
2697            problemStatus_=0;
2698            sumPrimalInfeasibilities_=0.0;
2699            if ((specialOptions_&(1024+16384))!=0) {
2700              CoinIndexedVector * arrayVector = rowArray_[1];
2701              arrayVector->clear();
2702              double * rhs = arrayVector->denseVector();
2703              times(1.0,solution_,rhs);
2704              bool bad2=false;
2705              int i;
2706              for ( i=0;i<numberRows_;i++) {
2707                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
2708                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
2709                  bad2=true;
2710                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
2711                }
2712                rhs[i]=0.0;
2713              }
2714              for ( i=0;i<numberColumns_;i++) {
2715                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
2716                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
2717                  bad2=true;
2718                }
2719              }
2720              if (bad2) {
2721                problemStatus_=-3;
2722                returnCode=-2;
2723                // Force to re-factorize early next time
2724                int numberPivots = factorization_->pivots();
2725                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
2726              }
2727            }
2728          }
2729        }
2730      } else {
2731        problemStatus_=-3;
2732        returnCode=-2;
2733        // Force to re-factorize early next time
2734        int numberPivots = factorization_->pivots();
2735        forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
2736      }
2737      break;
2738    }
2739  }
2740  delete [] primalChange;
2741  delete [] dualChange;
2742  return returnCode;
2743}
2744// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
2745int 
2746ClpSimplexOther::nextTheta(int type, double maxTheta, double * primalChange, double * dualChange,
2747                           const double * changeLower, const double * changeUpper,
2748                           const double * changeObjective)
2749{
2750  int numberTotal = numberColumns_+numberRows_;
2751  int iSequence;
2752  int iRow;
2753  theta_=maxTheta;
2754  bool toLower=false;
2755  if ((type&1)!=0) {
2756    // get change
2757    for (iSequence=0;iSequence<numberTotal;iSequence++) {
2758      primalChange[iSequence]=0.0;
2759      switch(getStatus(iSequence)) {
2760       
2761      case basic:
2762      case isFree:
2763      case superBasic:
2764        break;
2765      case isFixed:
2766      case atUpperBound:
2767        primalChange[iSequence]=changeUpper[iSequence];
2768        break;
2769      case atLowerBound:
2770        primalChange[iSequence]=changeLower[iSequence];
2771        break;
2772      }
2773    }
2774    // use array
2775    double * array = rowArray_[1]->denseVector();
2776    times(1.0,primalChange,array);
2777    int * index = rowArray_[1]->getIndices();
2778    int number=0;
2779    for (iRow=0;iRow<numberRows_;iRow++) {
2780      double value = array[iRow];
2781      if (value) {
2782        array[iRow]=value;
2783        index[number++]=iRow;
2784      }
2785    }
2786    // ftran it
2787    rowArray_[1]->setNumElements(number);
2788    factorization_->updateColumn(rowArray_[0],rowArray_[1]);
2789    number=rowArray_[1]->getNumElements();
2790    pivotRow_=-1;
2791    for (iRow=0;iRow<number;iRow++) {
2792      int iPivot = index[iRow];
2793      iSequence = pivotVariable_[iPivot];
2794      // solution value will be sol - theta*alpha
2795      // bounds will be bounds + change *theta
2796      double currentSolution = solution_[iSequence];
2797      double currentLower = lower_[iSequence];
2798      double currentUpper = upper_[iSequence];
2799      double alpha = array[iPivot];
2800      assert (currentSolution>=currentLower-primalTolerance_);
2801      assert (currentSolution<=currentUpper+primalTolerance_);
2802      double thetaCoefficient;
2803      double hitsLower = COIN_DBL_MAX;
2804      thetaCoefficient = changeLower[iSequence]+alpha;
2805      if (fabs(thetaCoefficient)>1.0e-8)
2806        hitsLower = (currentSolution-currentLower)/thetaCoefficient;
2807      if (hitsLower<0.0) {
2808        // does not hit - but should we check further
2809        hitsLower=COIN_DBL_MAX;
2810      }
2811      double hitsUpper = COIN_DBL_MAX;
2812      thetaCoefficient = changeUpper[iSequence]+alpha;
2813      if (fabs(thetaCoefficient)>1.0e-8)
2814        hitsUpper = (currentSolution-currentUpper)/thetaCoefficient;
2815      if (hitsUpper<0.0) {
2816        // does not hit - but should we check further
2817        hitsUpper=COIN_DBL_MAX;
2818      }
2819      if (CoinMin(hitsLower,hitsUpper)<theta_) {
2820        theta_ = CoinMin(hitsLower,hitsUpper);
2821        toLower = hitsLower<hitsUpper;
2822        pivotRow_=iPivot;
2823      }
2824    }
2825  }
2826  if ((type&2)!=0) {
2827    abort();
2828  }
2829  if (pivotRow_>=0) {
2830    sequenceOut_ = pivotVariable_[pivotRow_];
2831    valueOut_ = solution_[sequenceOut_];
2832    lowerOut_ = lower_[sequenceOut_];
2833    upperOut_ = upper_[sequenceOut_];
2834    if (!toLower) {
2835      directionOut_ = -1;
2836      dualOut_ = valueOut_ - upperOut_;
2837    } else if (valueOut_<lowerOut_) {
2838      directionOut_ = 1;
2839      dualOut_ = lowerOut_ - valueOut_;
2840    }
2841    return 0;
2842  } else {
2843    return -1;
2844  }
2845}
2846/* Expands out all possible combinations for a knapsack
2847   If buildObj NULL then just computes space needed - returns number elements
2848   On entry numberOutput is maximum allowed, on exit it is number needed or
2849   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
2850   least space to return values which reconstruct input.
2851   Rows returned will be original rows but no entries will be returned for
2852   any rows all of whose entries are in knapsack.  So up to user to allow for this.
2853   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
2854   in expanded knapsack.  Values in buildRow and buildElement;
2855*/
2856int 
2857ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
2858                                double * buildObj, CoinBigIndex * buildStart,
2859                                int * buildRow, double * buildElement,int reConstruct) const
2860{
2861  int iRow;
2862  int iColumn;
2863  // Get column copy
2864  CoinPackedMatrix * columnCopy = matrix();
2865  // Get a row copy in standard format
2866  CoinPackedMatrix matrixByRow;
2867  matrixByRow.reverseOrderedCopyOf(*columnCopy);
2868  const double * elementByRow = matrixByRow.getElements();
2869  const int * column = matrixByRow.getIndices();
2870  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2871  const int * rowLength = matrixByRow.getVectorLengths();
2872  CoinBigIndex j;
2873  int * whichColumn = new int [numberColumns_];
2874  int * whichRow = new int [numberRows_];
2875  int numJ=0;
2876  // Get what other columns can compensate for
2877  double * lo = new double [numberRows_];
2878  double * high = new double [numberRows_];
2879  {
2880    // Use to get tight column bounds
2881    ClpSimplex tempModel(*this);
2882    tempModel.tightenPrimalBounds(0.0,0,true);
2883    // Now another model without knapsacks
2884    int nCol=0;
2885    for (iRow=0;iRow<numberRows_;iRow++) {
2886      whichRow[iRow]=iRow;
2887    }
2888    for (iColumn=0;iColumn<numberColumns_;iColumn++)
2889      whichColumn[iColumn]=-1;
2890    for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
2891      int iColumn = column[j];
2892      if (columnUpper_[iColumn]>columnLower_[iColumn]) {
2893        whichColumn[iColumn]=0;
2894      } else {
2895        assert (!columnLower_[iColumn]); // fix later
2896      }
2897    }
2898    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2899      if (whichColumn[iColumn]<0)
2900        whichColumn[nCol++]=iColumn;
2901    }
2902    ClpSimplex tempModel2(&tempModel,numberRows_,whichRow,nCol,whichColumn,false,false,false);
2903    // Row copy
2904    CoinPackedMatrix matrixByRow;
2905    matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
2906    const double * elementByRow = matrixByRow.getElements();
2907    const int * column = matrixByRow.getIndices();
2908    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
2909    const int * rowLength = matrixByRow.getVectorLengths();
2910    const double * columnLower = tempModel2.getColLower();
2911    const double * columnUpper = tempModel2.getColUpper();
2912    for (iRow = 0; iRow < numberRows_; iRow++) {
2913      lo[iRow]=-COIN_DBL_MAX;
2914      high[iRow]=COIN_DBL_MAX;
2915      if (rowLower_[iRow]>-1.0e20||rowUpper_[iRow]<1.0e20) {
2916
2917        // possible row
2918        int infiniteUpper = 0;
2919        int infiniteLower = 0;
2920        double maximumUp = 0.0;
2921        double maximumDown = 0.0;
2922        CoinBigIndex rStart = rowStart[iRow];
2923        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
2924        CoinBigIndex j;
2925        // Compute possible lower and upper ranges
2926
2927        for (j = rStart; j < rEnd; ++j) {
2928          double value=elementByRow[j];
2929          iColumn = column[j];
2930          if (value > 0.0) {
2931            if (columnUpper[iColumn] >= 1.0e20) {
2932              ++infiniteUpper;
2933            } else {
2934              maximumUp += columnUpper[iColumn] * value;
2935            }
2936            if (columnLower[iColumn] <= -1.0e20) {
2937              ++infiniteLower;
2938            } else {
2939              maximumDown += columnLower[iColumn] * value;
2940            }
2941          } else if (value<0.0) {
2942            if (columnUpper[iColumn] >= 1.0e20) {
2943              ++infiniteLower;
2944            } else {
2945              maximumDown += columnUpper[iColumn] * value;
2946            }
2947            if (columnLower[iColumn] <= -1.0e20) {
2948              ++infiniteUpper;
2949            } else {
2950              maximumUp += columnLower[iColumn] * value;
2951            }
2952          }
2953        }
2954        // Build in a margin of error
2955        maximumUp += 1.0e-8*fabs(maximumUp)+1.0e-7;
2956        maximumDown -= 1.0e-8*fabs(maximumDown)+1.0e-7;
2957        // we want to save effective rhs
2958        double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
2959        double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
2960        if (up==COIN_DBL_MAX||rowLower_[iRow]==-COIN_DBL_MAX) {
2961          // However low we go it doesn't matter
2962          lo[iRow]=-COIN_DBL_MAX;
2963        } else {
2964          // If we go below this then can not be feasible
2965          lo[iRow]=rowLower_[iRow]-up;
2966        }
2967        if (down==-COIN_DBL_MAX||rowUpper_[iRow]==COIN_DBL_MAX) {
2968          // However high we go it doesn't matter
2969          high[iRow]=COIN_DBL_MAX;
2970        } else {
2971          // If we go above this then can not be feasible
2972          high[iRow]=rowUpper_[iRow]-down;
2973        }
2974      }
2975    }
2976  }
2977  numJ =0;
2978  for (iColumn=0;iColumn<numberColumns_;iColumn++)
2979    whichColumn[iColumn]=-1;
2980  int * markRow = new int [numberRows_];
2981  for (iRow=0;iRow<numberRows_;iRow++) 
2982    markRow[iRow]=1;
2983  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
2984    int iColumn = column[j];
2985    if (columnUpper_[iColumn]>columnLower_[iColumn]) {
2986      whichColumn[iColumn]=numJ;
2987      numJ++;
2988    }
2989  }
2990  /* mark rows
2991     -n in knapsack and n other variables
2992     1 no entries
2993     n+1000 not involved in knapsack but n entries
2994     0 only in knapsack
2995  */
2996  for (iRow=0;iRow<numberRows_;iRow++) {
2997    int type=1;
2998    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
2999      int iColumn = column[j];
3000      if (whichColumn[iColumn]>=0) {
3001        if (type==1) {
3002          type=0;
3003        } else if (type>0) {
3004          assert (type>1000); 
3005          type = -(type-1000);
3006        }
3007      } else if (type==1) {
3008        type = 1001;
3009      } else if (type<0) {
3010        type --;
3011      } else if (type==0) {
3012        type = -1;
3013      } else {
3014        assert (type>1000);
3015        type++;
3016      }
3017    }
3018    markRow[iRow]=type;
3019    if (type<0&&type>-30&&false)
3020      printf("markrow on row %d is %d\n",iRow,markRow[iRow]);
3021  }
3022  int * bound = new int [numberColumns_+1];
3023  int * stack = new int [numberColumns_+1];
3024  int * flip = new int [numberColumns_+1];
3025  double * offset = new double[numberColumns_+1];
3026  double * size = new double [numberColumns_+1];
3027  double * rhsOffset = new double[numberRows_];
3028  int * build = new int[numberColumns_];
3029  int maxNumber=numberOutput;
3030  numJ=0;
3031  double minSize = rowLower_[knapsackRow];
3032  double maxSize = rowUpper_[knapsackRow];
3033  double knapsackOffset=0.0;
3034  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
3035    int iColumn = column[j];
3036    double lowerColumn=columnLower_[iColumn];
3037    double upperColumn=columnUpper_[iColumn];
3038    if (lowerColumn==upperColumn)
3039      continue;
3040    double gap = upperColumn-lowerColumn;
3041    if (gap>1.0e8)
3042      gap=1.0e8;
3043    assert (fabs(floor(gap+0.5)-gap)<1.0e-5);
3044    whichColumn[numJ]=iColumn;
3045    bound[numJ]=static_cast<int> (gap);
3046    if (elementByRow[j]>0.0) {
3047      flip[numJ]=1;
3048      offset[numJ]=lowerColumn;
3049      size[numJ++]=elementByRow[j];
3050    } else {
3051      flip[numJ]=-1;
3052      offset[numJ]=upperColumn;
3053      size[numJ++]=-elementByRow[j];
3054      lowerColumn = upperColumn;
3055    }
3056    knapsackOffset += elementByRow[j]*lowerColumn;
3057  } 
3058  int jRow;
3059  for (iRow=0;iRow<numberRows_;iRow++)
3060    whichRow[iRow]=iRow;
3061  ClpSimplex smallModel(this,numberRows_,whichRow,numJ,whichColumn,true,true,true);
3062  // modify rhs to allow for nonzero lower bounds
3063  //double * rowLower = smallModel.rowLower();
3064  //double * rowUpper = smallModel.rowUpper();
3065  //const double * columnLower = smallModel.columnLower();
3066  //const double * columnUpper = smallModel.columnUpper();
3067  const CoinPackedMatrix * matrix = smallModel.matrix();
3068  const double * element = matrix->getElements();
3069  const int * row = matrix->getIndices();
3070  const CoinBigIndex * columnStart = matrix->getVectorStarts();
3071  const int * columnLength = matrix->getVectorLengths();
3072  const double * objective = smallModel.objective();
3073  //double objectiveOffset=0.0;
3074  // would use for fixed?
3075  CoinZeroN(rhsOffset,numberRows_);
3076  double * rowActivity = smallModel.primalRowSolution();
3077  CoinZeroN(rowActivity,numberRows_);
3078  maxSize -= knapsackOffset;
3079  minSize -= knapsackOffset;
3080  // now generate
3081  int i;
3082  int iStack=numJ;
3083  for (i=0;i<numJ;i++) {
3084    stack[i]=0;
3085  }
3086  double tooMuch = 10.0*maxSize+10000;
3087  stack[numJ]=1;
3088  size[numJ]=tooMuch;
3089  bound[numJ]=0;
3090  double sum=tooMuch;
3091  // allow for all zero being OK
3092  stack[numJ-1]=-1;
3093  sum -= size[numJ-1];
3094  numberOutput=0;
3095  int nelCreate=0;
3096  /* typeRun is - 0 for initial sizes
3097                  1 for build
3098                  2 for reconstruct
3099  */
3100  int typeRun = buildObj ? 1 : 0;
3101  if (reConstruct>=0) {
3102    assert (buildRow&&buildElement);
3103    typeRun=2;
3104  }
3105  if (typeRun==1)
3106    buildStart[0]=0;
3107  while (iStack>=0) {
3108    if (sum>=minSize&&sum<=maxSize) {
3109      double checkSize =0.0;
3110      bool good=true;
3111      int nRow=0;
3112      double obj=0.0;
3113      CoinZeroN(rowActivity,nRow);
3114      for (iColumn=0;iColumn<numJ;iColumn++) {
3115        int iValue = stack[iColumn];
3116        if (iValue>bound[iColumn]) {
3117          good=false;
3118          break;
3119        } else {
3120          double realValue = offset[iColumn] + flip[iColumn]*iValue;
3121          if (realValue) {
3122            obj += objective[iColumn]*realValue;
3123            for (CoinBigIndex j=columnStart[iColumn];
3124                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
3125              double value = element[j]*realValue;
3126              int kRow = row[j];
3127              if (rowActivity[kRow]) {
3128                rowActivity[kRow] += value;
3129                if (!rowActivity[kRow])
3130                  rowActivity[kRow]=1.0e-100;
3131              } else {
3132                build[nRow++]=kRow;
3133                rowActivity[kRow]=value;
3134              }
3135            }
3136          }
3137        }
3138      }
3139      if (good) {
3140        for (jRow=0;jRow<nRow;jRow++) {
3141          int kRow=build[jRow];
3142          double value = rowActivity[kRow];
3143          if (value>high[kRow]||value<lo[kRow]) {
3144            good=false;
3145            break;
3146          }
3147        }
3148      }
3149      if (good) {
3150        if (typeRun==1) {
3151          buildObj[numberOutput]=obj;
3152          for (jRow=0;jRow<nRow;jRow++) {
3153            int kRow=build[jRow];
3154            double value = rowActivity[kRow];
3155            if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
3156              buildElement[nelCreate]=value;
3157              buildRow[nelCreate++]=kRow;
3158            }
3159          }
3160          buildStart[numberOutput+1]=nelCreate;
3161        } else if (!typeRun) {
3162          for (jRow=0;jRow<nRow;jRow++) {
3163            int kRow=build[jRow];
3164            double value = rowActivity[kRow];
3165            if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
3166              nelCreate++;
3167            }
3168          }
3169        }
3170        if (typeRun==2&&reConstruct==numberOutput) {
3171          // build and exit
3172          nelCreate=0;
3173          for (iColumn=0;iColumn<numJ;iColumn++) {
3174            int iValue = stack[iColumn];
3175            double realValue = offset[iColumn] + flip[iColumn]*iValue;
3176            if (realValue) {
3177              buildRow[nelCreate]=whichColumn[iColumn];
3178              buildElement[nelCreate++]=realValue;
3179            }
3180          }
3181          numberOutput=1;
3182          for (i=0;i<numJ;i++) {
3183            bound[i]=0;
3184          }
3185          break;
3186        }
3187        numberOutput++;
3188        if (numberOutput>maxNumber) {
3189          nelCreate=-numberOutput;
3190          numberOutput=-1;
3191          for (i=0;i<numJ;i++) {
3192            bound[i]=0;
3193          }
3194          break;
3195        } else if (typeRun==1&&numberOutput==maxNumber) {
3196          // On second run
3197          for (i=0;i<numJ;i++) {
3198            bound[i]=0;
3199          }
3200          break;
3201        } 
3202        for (int j=0;j<numJ;j++) {
3203          checkSize += stack[j]*size[j];
3204        }
3205        assert (fabs(sum-checkSize)<1.0e-3);
3206      }
3207      for (jRow=0;jRow<nRow;jRow++) {
3208        int kRow=build[jRow];
3209        rowActivity[kRow]=0.0;
3210      }
3211    }
3212    if (sum>maxSize||stack[iStack]>bound[iStack]) {
3213      sum -= size[iStack]*stack[iStack]; 
3214      stack[iStack--]=0;
3215      if (iStack>=0) {
3216        stack[iStack] ++;
3217        sum += size[iStack];
3218      }
3219    } else {
3220      // must be less
3221      // add to last possible
3222      iStack = numJ-1;
3223      sum += size[iStack];
3224      stack[iStack]++;
3225    }
3226  }
3227  //printf("%d will be created\n",numberOutput);
3228  delete [] whichColumn;
3229  delete [] whichRow;
3230  delete [] bound;
3231  delete [] stack;
3232  delete [] flip;
3233  delete [] size;
3234  delete [] offset;
3235  delete [] rhsOffset;
3236  delete [] build;
3237  delete [] markRow;
3238  delete [] lo;
3239  delete [] high;
3240  return nelCreate;
3241}
Note: See TracBrowser for help on using the repository browser.