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

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

changes to simplex and lots of stuff and start Mumps cholesky

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