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

Last change on this file since 1366 was 1366, checked in by forrest, 11 years ago

try and improve quadratic and barrier

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