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

Last change on this file since 1304 was 1286, checked in by forrest, 12 years ago

changes for factorization and aux region

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