source: branches/devel/Clp/src/ClpSimplexOther.cpp @ 989

Last change on this file since 989 was 989, checked in by forrest, 13 years ago

this may be a mistake

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