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

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

ranging

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