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

Last change on this file since 1402 was 1402, checked in by forrest, 10 years ago

get rid of compiler warnings

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