Logo Search packages:      
Sourcecode: rsymphony version File versions  Download package

ClpSimplexOther.cpp

/* $Id: ClpSimplexOther.cpp 1458 2009-11-05 12:34:07Z forrest $ */
// Copyright (C) 2004, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "CoinPragma.hpp"

#include <math.h>

#include "CoinHelperFunctions.hpp"
#include "ClpSimplexOther.hpp"
#include "ClpSimplexDual.hpp"
#include "ClpSimplexPrimal.hpp"
#include "ClpEventHandler.hpp"
#include "ClpHelperFunctions.hpp"
#include "ClpFactorization.hpp"
#include "ClpDualRowDantzig.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinBuild.hpp"
#include "CoinMpsIO.hpp"
#include "ClpMessage.hpp"
#include <cfloat>
#include <cassert>
#include <string>
#include <stdio.h>
#include <iostream>
/* Dual ranging.
   This computes increase/decrease in cost for each given variable and corresponding
   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns 
   and numberColumns.. for artificials/slacks.
   For non-basic variables the sequence number will be that of the non-basic variables.
   
   Up to user to provide correct length arrays.
   
*/
00036 void ClpSimplexOther::dualRanging(int numberCheck,const int * which,
                      double * costIncreased, int * sequenceIncreased,
                          double * costDecreased, int * sequenceDecreased,
                          double * valueIncrease, double * valueDecrease)
{
  rowArray_[1]->clear();
  columnArray_[1]->clear();
  // long enough for rows+columns
  assert(rowArray_[3]->capacity()>=numberRows_+numberColumns_);
  rowArray_[3]->clear();
  int * backPivot = rowArray_[3]->getIndices();
  int i;
  for ( i=0;i<numberRows_+numberColumns_;i++) {
    backPivot[i]=-1;
  }
  for (i=0;i<numberRows_;i++) {
    int iSequence = pivotVariable_[i];
    backPivot[iSequence]=i;
  }
  // dualTolerance may be zero if from CBC.  In fact use that fact
  bool inCBC = !dualTolerance_;
  if (inCBC)
    assert (integerType_);
  dualTolerance_ = dblParam_[ClpDualTolerance];
  double * arrayX = rowArray_[0]->denseVector();
  for ( i=0;i<numberCheck;i++) {
    rowArray_[0]->clear();
    //rowArray_[0]->checkClear();
    //rowArray_[1]->checkClear();
    //columnArray_[1]->checkClear();
    columnArray_[0]->clear();
    //columnArray_[0]->checkClear();
    int iSequence = which[i];
    if (iSequence<0) {
      costIncreased[i] = 0.0;
      sequenceIncreased[i] = -1;
      costDecreased[i] = 0.0;
      sequenceDecreased[i] = -1;
      continue;
    }
    double costIncrease=COIN_DBL_MAX;
    double costDecrease=COIN_DBL_MAX;
    int sequenceIncrease=-1;
    int sequenceDecrease=-1;
    if (valueIncrease) {
      assert (valueDecrease);
      valueIncrease[i]=iSequence<numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
      valueDecrease[i]=valueIncrease[i];
    }
    
    switch(getStatus(iSequence)) {
      
    case basic:
      {
      // non-trvial
      // Get pivot row
      int iRow=backPivot[iSequence];
      assert (iRow>=0);
      double plusOne=1.0;
        rowArray_[0]->createPacked(1,&iRow,&plusOne);
      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
      // put row of tableau in rowArray[0] and columnArray[0]
      matrix_->transposeTimes(this,-1.0,
                        rowArray_[0],columnArray_[1],columnArray_[0]);
        double alphaIncrease;
        double alphaDecrease;
      // do ratio test up and down
      checkDualRatios(rowArray_[0],columnArray_[0],costIncrease,sequenceIncrease,alphaIncrease,
                costDecrease,sequenceDecrease,alphaDecrease);
        if (!inCBC) {
        if (valueIncrease) {
          if (sequenceIncrease>=0)
            valueIncrease[i] = primalRanging1(sequenceIncrease,iSequence);
          if (sequenceDecrease>=0)
            valueDecrease[i] = primalRanging1(sequenceDecrease,iSequence);
        }
      } else { 
        int number = rowArray_[0]->getNumElements();
        double scale2=0.0;
        int j;
        for (j=0;j<number;j++) {
          scale2+=arrayX[j]*arrayX[j];
        }
        scale2 = 1.0/sqrt(scale2);
        valueIncrease[i]=scale2;
          if (sequenceIncrease>=0) {
            double djValue = dj_[sequenceIncrease];
            if (fabs(djValue)>10.0*dualTolerance_) {
              // we are going to use for cutoff so be exact
              costIncrease = fabs(djValue/alphaIncrease); 
              /* Not sure this is good idea as I don't think correct e.g.
                 suppose a continuous variable has dj slightly greater. */
              if(false&&sequenceIncrease<numberColumns_&&integerType_[sequenceIncrease]) {
                // can improve
                double movement = (columnScale_==NULL) ? 1.0 : 
                  rhsScale_*inverseColumnScale_[sequenceIncrease];
                costIncrease = CoinMax(fabs(djValue*movement),costIncrease);
              }
            } else {
              costIncrease=0.0;
            }
          }
          if (sequenceDecrease>=0) {
            double djValue = dj_[sequenceDecrease];
            if (fabs(djValue)>10.0*dualTolerance_) {
              // we are going to use for cutoff so be exact
              costDecrease = fabs(djValue/alphaDecrease); 
              if(sequenceDecrease<numberColumns_&&integerType_[sequenceDecrease]) {
                // can improve
                double movement = (columnScale_==NULL) ? 1.0 : 
                  rhsScale_*inverseColumnScale_[sequenceDecrease];
                costDecrease = CoinMax(fabs(djValue*movement),costDecrease);
              }
            } else {
              costDecrease=0.0;
            }
          }
        costIncrease *= scale2;
        costDecrease *= scale2;
        }
      }
      break;
    case isFixed:
      break;
    case isFree:
    case superBasic:
      costIncrease=0.0;
      costDecrease=0.0;
      sequenceIncrease=iSequence;
      sequenceDecrease=iSequence;
      break;
    case atUpperBound:
      costIncrease = CoinMax(0.0,-dj_[iSequence]);
      sequenceIncrease = iSequence;
      if (valueIncrease) 
      valueIncrease[i] = primalRanging1(iSequence,iSequence);
      break;
    case atLowerBound:
      costDecrease = CoinMax(0.0,dj_[iSequence]);
      sequenceDecrease = iSequence;
      if (valueIncrease) 
      valueDecrease[i] = primalRanging1(iSequence,iSequence);
      break;
    }
    double scaleFactor;
    if (rowScale_) {
      if (iSequence<numberColumns_) 
      scaleFactor = 1.0/(objectiveScale_*columnScale_[iSequence]);
      else
      scaleFactor = rowScale_[iSequence-numberColumns_]/objectiveScale_;
    } else {
      scaleFactor = 1.0/objectiveScale_;
    }
    if (costIncrease<1.0e30)
      costIncrease *= scaleFactor;
    if (costDecrease<1.0e30)
      costDecrease *= scaleFactor;
    if (optimizationDirection_==1.0) {
      costIncreased[i] = costIncrease;
      sequenceIncreased[i] = sequenceIncrease;
      costDecreased[i] = costDecrease;
      sequenceDecreased[i] = sequenceDecrease;
    } else if (optimizationDirection_==-1.0) {
      costIncreased[i] = costDecrease;
      sequenceIncreased[i] = sequenceDecrease;
      costDecreased[i] = costIncrease;
      sequenceDecreased[i] = sequenceIncrease;
      if (valueIncrease) {
      double temp = valueIncrease[i];
      valueIncrease[i]=valueDecrease[i];
      valueDecrease[i]=temp;
      }
    } else if (optimizationDirection_==0.0) {
      // !!!!!! ???
      costIncreased[i] = COIN_DBL_MAX;
      sequenceIncreased[i] = -1;
      costDecreased[i] = COIN_DBL_MAX;
      sequenceDecreased[i] = -1;
    } else {
      abort();
    }
  }
  rowArray_[0]->clear();
  //rowArray_[1]->clear();
  //columnArray_[1]->clear();
  columnArray_[0]->clear();
  //rowArray_[3]->clear();
  if (!optimizationDirection_)
    printf("*** ????? Ranging with zero optimization costs\n");
}
/* 
   Row array has row part of pivot row
   Column array has column part.
   This is used in dual ranging
*/
void
00232 ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
                                 CoinIndexedVector * columnArray,
                                 double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
                                 double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
{
  double acceptablePivot = 1.0e-9;
  double * work;
  int number;
  int * which;
  int iSection;

  double thetaDown = 1.0e31;
  double thetaUp = 1.0e31;
  int sequenceDown =-1;
  int sequenceUp = -1;
  double alphaDown=0.0;
  double alphaUp=0.0;

  int addSequence;

  for (iSection=0;iSection<2;iSection++) {

    int i;
    if (!iSection) {
      work = rowArray->denseVector();
      number = rowArray->getNumElements();
      which = rowArray->getIndices();
      addSequence = numberColumns_;
    } else {
      work = columnArray->denseVector();
      number = columnArray->getNumElements();
      which = columnArray->getIndices();
      addSequence = 0;
    }
    
    for (i=0;i<number;i++) {
      int iSequence = which[i];
      int iSequence2 = iSequence + addSequence;
      double alpha=work[i];
      if (fabs(alpha)<acceptablePivot)
      continue;
      double oldValue=dj_[iSequence2];

      switch(getStatus(iSequence2)) {
        
      case basic: 
      break;
      case ClpSimplex::isFixed:
      break;
      case isFree:
      case superBasic:
      // treat dj as if zero
      thetaDown=0.0;
      thetaUp=0.0;
      sequenceDown=iSequence2;
      sequenceUp=iSequence2;
      break;
      case atUpperBound:
      if (alpha>0.0) {
        // test up
        if (oldValue + thetaUp*alpha > dualTolerance_) {
          thetaUp = (dualTolerance_-oldValue)/alpha;
          sequenceUp = iSequence2;
            alphaUp=alpha;
        }
      } else {
        // test down
        if (oldValue - thetaDown*alpha > dualTolerance_) {
          thetaDown = -(dualTolerance_-oldValue)/alpha;
          sequenceDown = iSequence2;
            alphaDown=alpha;
        }
      }
      break;
      case atLowerBound:
      if (alpha<0.0) {
        // test up
        if (oldValue + thetaUp*alpha <- dualTolerance_) {
          thetaUp = -(dualTolerance_+oldValue)/alpha;
          sequenceUp = iSequence2;
            alphaUp=alpha;
        }
      } else {
        // test down
        if (oldValue - thetaDown*alpha < -dualTolerance_) {
          thetaDown = (dualTolerance_+oldValue)/alpha;
          sequenceDown = iSequence2;
            alphaDown=alpha;
        }
      }
      break;
      }
    }
  }
  if (sequenceUp>=0) {
    costIncrease = thetaUp;
    sequenceIncrease = sequenceUp;
    alphaIncrease = alphaUp;
  }
  if (sequenceDown>=0) {
    costDecrease = thetaDown;
    sequenceDecrease = sequenceDown;
    alphaDecrease = alphaDown;
  }
}
/** Primal ranging.
    This computes increase/decrease in value for each given variable and corresponding
    sequence numbers which would change basis.  Sequence numbers are 0..numberColumns 
    and numberColumns.. for artificials/slacks.
    For basic variables the sequence number will be that of the basic variables.
    
    Up to user to provide correct length arrays.
    
    When here - guaranteed optimal
*/
void 
00348 ClpSimplexOther::primalRanging(int numberCheck,const int * which,
              double * valueIncreased, int * sequenceIncreased,
              double * valueDecreased, int * sequenceDecreased)
{
  rowArray_[0]->clear();
  rowArray_[1]->clear();
  lowerIn_=-COIN_DBL_MAX;
  upperIn_=COIN_DBL_MAX;
  valueIn_ = 0.0;
  for ( int i=0;i<numberCheck;i++) {
    int iSequence = which[i];
    double valueIncrease=COIN_DBL_MAX;
    double valueDecrease=COIN_DBL_MAX;
    int sequenceIncrease=-1;
    int sequenceDecrease=-1;
    
    switch(getStatus(iSequence)) {
      
    case basic:
    case isFree:
    case superBasic:
      // Easy
      valueDecrease=CoinMax(0.0,upper_[iSequence]-solution_[iSequence]);
      valueIncrease=CoinMax(0.0,solution_[iSequence]-lower_[iSequence]);
      sequenceDecrease=iSequence;
      sequenceIncrease=iSequence;
      break;
    case isFixed:
    case atUpperBound:
    case atLowerBound:
      {
      // Non trivial
      // Other bound is ignored
      unpackPacked(rowArray_[1],iSequence);
      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
      // Get extra rows
      matrix_->extendUpdated(this,rowArray_[1],0);
      // do ratio test
      checkPrimalRatios(rowArray_[1],1);
      if (pivotRow_>=0) {
        valueIncrease = theta_;
        sequenceIncrease=pivotVariable_[pivotRow_];
      }
      checkPrimalRatios(rowArray_[1],-1);
      if (pivotRow_>=0) {
        valueDecrease = theta_;
        sequenceDecrease=pivotVariable_[pivotRow_];
      }
      rowArray_[1]->clear();
      }
      break;
    }
    double scaleFactor;
    if (rowScale_) {
      if (iSequence<numberColumns_) 
      scaleFactor = columnScale_[iSequence]/rhsScale_;
      else
      scaleFactor = 1.0/(rowScale_[iSequence-numberColumns_]*rhsScale_);
    } else {
      scaleFactor = 1.0/rhsScale_;
    }
    if (valueIncrease<1.0e30)
      valueIncrease *= scaleFactor;
    else
      valueIncrease = COIN_DBL_MAX;
    if (valueDecrease<1.0e30)
      valueDecrease *= scaleFactor;
    else
      valueDecrease = COIN_DBL_MAX;
    valueIncreased[i] = valueIncrease;
    sequenceIncreased[i] = sequenceIncrease;
    valueDecreased[i] = valueDecrease;
    sequenceDecreased[i] = sequenceDecrease;
  }
}
// Returns new value of whichOther when whichIn enters basis
double 
00425 ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
{
  rowArray_[0]->clear();
  rowArray_[1]->clear();
  int iSequence = whichIn;
  double newValue=solution_[whichOther];
  double alphaOther=0.0;
  Status status = getStatus(iSequence);
  assert (status==atLowerBound||status==atUpperBound);
  int wayIn = (status==atLowerBound) ? 1 : -1;
  
  switch(getStatus(iSequence)) {
    
  case basic:
  case isFree:
  case superBasic:
    assert (whichIn==whichOther);
    // Easy
    newValue = wayIn>0 ? upper_[iSequence] : lower_[iSequence];
    break;
  case isFixed:
  case atUpperBound:
  case atLowerBound:
    // Non trivial
    {
      // Other bound is ignored
      unpackPacked(rowArray_[1],iSequence);
      factorization_->updateColumn(rowArray_[2],rowArray_[1]);
      // Get extra rows
      matrix_->extendUpdated(this,rowArray_[1],0);
      // do ratio test
      double acceptablePivot=1.0e-7;
      double * work=rowArray_[1]->denseVector();
      int number=rowArray_[1]->getNumElements();
      int * which=rowArray_[1]->getIndices();
      
      // we may need to swap sign
      double way = wayIn;
      double theta = 1.0e30;
      for (int iIndex=0;iIndex<number;iIndex++) {

      int iRow = which[iIndex];
      double alpha = work[iIndex]*way;
      int iPivot=pivotVariable_[iRow];
      if (iPivot==whichOther) {
        alphaOther=alpha;
        continue;
      }
      double oldValue = solution_[iPivot];
      if (fabs(alpha)>acceptablePivot) {
        if (alpha>0.0) {
          // basic variable going towards lower bound
          double bound = lower_[iPivot];
          oldValue -= bound;
          if (oldValue-theta*alpha<0.0) {
            theta = CoinMax(0.0,oldValue/alpha);
          }
        } else {
          // basic variable going towards upper bound
          double bound = upper_[iPivot];
          oldValue = oldValue-bound;
          if (oldValue-theta*alpha>0.0) {
            theta = CoinMax(0.0,oldValue/alpha);
          }
        }
      }
      }
      if (whichIn!=whichOther) {
      if (theta<1.0e30)
        newValue -= theta*alphaOther;
      else
        newValue = alphaOther>0.0 ? -1.0e30 : 1.0e30;
      } else {
      newValue += theta*wayIn;
      }
    }
    rowArray_[1]->clear();
    break;
  }
  double scaleFactor;
  if (rowScale_) {
    if (whichOther<numberColumns_) 
      scaleFactor = columnScale_[whichOther]/rhsScale_;
    else
      scaleFactor = 1.0/(rowScale_[whichOther-numberColumns_]*rhsScale_);
  } else {
    scaleFactor = 1.0/rhsScale_;
  }
  if (newValue<1.0e29)
    if (newValue>-1.0e29)
      newValue *= scaleFactor;
    else
      newValue = -COIN_DBL_MAX;
  else
    newValue = COIN_DBL_MAX;
  return newValue;
}
/* 
   Row array has pivot column
   This is used in primal ranging
*/
void 
00527 ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
                   int direction)
{
  // sequence stays as row number until end
  pivotRow_=-1;
  double acceptablePivot=1.0e-7;
  double * work=rowArray->denseVector();
  int number=rowArray->getNumElements();
  int * which=rowArray->getIndices();

  // we need to swap sign if going down
  double way = direction;
  theta_ = 1.0e30;
  for (int iIndex=0;iIndex<number;iIndex++) {
    
    int iRow = which[iIndex];
    double alpha = work[iIndex]*way;
    int iPivot=pivotVariable_[iRow];
    double oldValue = solution_[iPivot];
    if (fabs(alpha)>acceptablePivot) {
      if (alpha>0.0) {
      // basic variable going towards lower bound
      double bound = lower_[iPivot];
      oldValue -= bound;
      if (oldValue-theta_*alpha<0.0) {
        pivotRow_ = iRow;
        theta_ = CoinMax(0.0,oldValue/alpha);
      }
      } else {
      // basic variable going towards upper bound
      double bound = upper_[iPivot];
      oldValue = oldValue-bound;
      if (oldValue-theta_*alpha>0.0) {
        pivotRow_ = iRow;
        theta_ = CoinMax(0.0,oldValue/alpha);
      }
      }
    }
  }
}
/* Write the basis in MPS format to the specified file.
   If writeValues true writes values of structurals
   (and adds VALUES to end of NAME card)
   
   Row and column names may be null.
   formatType is
   <ul>
   <li> 0 - normal
   <li> 1 - extra accuracy 
   <li> 2 - IEEE hex (later)
   </ul>
   
   Returns non-zero on I/O error

   This is based on code contributed by Thorsten Koch
*/
int 
00584 ClpSimplexOther::writeBasis(const char *filename,
                      bool writeValues,
                      int formatType) const
{
  formatType=CoinMax(0,formatType);
  formatType=CoinMin(2,formatType);
  if (!writeValues)
    formatType=0;
  // See if INTEL if IEEE
  if (formatType==2) {
    // test intel here and add 1 if not intel
    double value=1.0;
    char x[8];
    memcpy(x,&value,8);
    if (x[0]==63) {
      formatType ++; // not intel
    } else {
      assert (x[0]==0);
    }
  }
  
  char number[20];
  FILE * fp = fopen(filename,"w");
  if (!fp)
    return -1;
   
  // NAME card

  if (strcmp(strParam_[ClpProbName].c_str(),"")==0) {
    fprintf(fp, "NAME          BLANK      ");
  } else {
    fprintf(fp, "NAME          %s       ",strParam_[ClpProbName].c_str());
  }
  if (formatType>=2)
    fprintf(fp,"FREEIEEE");
  else if (writeValues)
    fprintf(fp,"VALUES");
  // finish off name 
  fprintf(fp,"\n");
  int iRow=0;
  for(int iColumn =0; iColumn < numberColumns_; iColumn++) {
    bool printit=false;
    if( getColumnStatus(iColumn) == ClpSimplex::basic) {
      printit=true;
      // Find non basic row
      for(; iRow < numberRows_; iRow++) {
      if (getRowStatus(iRow) != ClpSimplex::basic) 
        break;
      }
      if (lengthNames_) {
      if (iRow!=numberRows_) {
        fprintf(fp, " %s %-8s       %s",
              getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
              columnNames_[iColumn].c_str(),
              rowNames_[iRow].c_str());
        iRow++;
      } else {
        // Allow for too many basics!
        fprintf(fp, " BS %-8s       ",
              columnNames_[iColumn].c_str());
        // Dummy row name if values
        if (writeValues)
          fprintf(fp,"      _dummy_");
      }
      } else {
      // no names
      if (iRow!=numberRows_) {
        fprintf(fp, " %s C%7.7d     R%7.7d",
              getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
              iColumn,iRow);
        iRow++;
      } else {
        // Allow for too many basics!
        fprintf(fp, " BS C%7.7d",iColumn);
        // Dummy row name if values
        if (writeValues)
          fprintf(fp,"      _dummy_");
      }
      }
    } else  {
      if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
      printit=true;
      if (lengthNames_) 
        fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
      else 
        fprintf(fp, " UL C%7.7d", iColumn);
      // Dummy row name if values
      if (writeValues)
        fprintf(fp,"      _dummy_");
      }
    }
    if (printit&&writeValues) {
      // add value
      CoinConvertDouble(0,formatType,columnActivity_[iColumn],number);
      fprintf(fp,"     %s",number);
    }
    if (printit)
      fprintf(fp,"\n");
  }
  fprintf(fp, "ENDATA\n");
  fclose(fp);
  return 0;
}
// Read a basis from the given filename
int 
00689 ClpSimplexOther::readBasis(const char *fileName)
{
  int status=0;
  bool canOpen=false;
  if (!strcmp(fileName,"-")||!strcmp(fileName,"stdin")) {
    // stdin
    canOpen=true;
  } else {
    FILE *fp=fopen(fileName,"r");
    if (fp) {
      // can open - lets go for it
      fclose(fp);
      canOpen=true;
    } else {
      handler_->message(CLP_UNABLE_OPEN,messages_)
      <<fileName<<CoinMessageEol;
      return -1;
    }
  }
  CoinMpsIO m;
  m.passInMessageHandler(handler_);
  *m.messagesPointer()=coinMessages();
  bool savePrefix =m.messageHandler()->prefix();
  m.messageHandler()->setPrefix(handler_->prefix());
  status=m.readBasis(fileName,"",columnActivity_,status_+numberColumns_,
                 status_,
                 columnNames_,numberColumns_,
                 rowNames_,numberRows_);
  m.messageHandler()->setPrefix(savePrefix);
  if (status>=0) {
    if (!status) {
      // set values
      int iColumn,iRow;
      for (iRow=0;iRow<numberRows_;iRow++) {
      if (getRowStatus(iRow)==atLowerBound)
        rowActivity_[iRow]=rowLower_[iRow];
      else if (getRowStatus(iRow)==atUpperBound)
        rowActivity_[iRow]=rowUpper_[iRow];
      }
      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
      if (getColumnStatus(iColumn)==atLowerBound)
        columnActivity_[iColumn]=columnLower_[iColumn];
      else if (getColumnStatus(iColumn)==atUpperBound)
        columnActivity_[iColumn]=columnUpper_[iColumn];
      }
    } else {
      memset(rowActivity_,0,numberRows_*sizeof(double));
      matrix_->times(-1.0,columnActivity_,rowActivity_);
    }
  } else {
    // errors
    handler_->message(CLP_IMPORT_ERRORS,messages_)
      <<status<<fileName<<CoinMessageEol;
  }
  return status;
}
/* Creates dual of a problem if looks plausible
   (defaults will always create model)
   fractionRowRanges is fraction of rows allowed to have ranges 
   fractionColumnRanges is fraction of columns allowed to have ranges 
*/
ClpSimplex * 
00751 ClpSimplexOther::dualOfModel(double fractionRowRanges,double fractionColumnRanges) const
{
  const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
  bool changed=false;
  int numberChanged=0;
  int iColumn;
  // check if we need to change bounds to rows
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    if (columnUpper_[iColumn]<1.0e20&&
        columnLower_[iColumn]>-1.0e20) {
      changed=true;
      numberChanged++;
    }
  }
  int iRow;
  int numberExtraRows=0;
  if (numberChanged<=fractionColumnRanges*numberColumns_) {
    for (iRow=0;iRow<numberRows_;iRow++) {
      if (rowLower_[iRow]>-1.0e20&&
        rowUpper_[iRow]<1.0e20) {
      if (rowUpper_[iRow]!=rowLower_[iRow])
        numberExtraRows++;
      }
    }
    if (numberExtraRows>fractionRowRanges*numberRows_)
      return NULL;
  } else {
    return NULL;
  }
  if (changed) {
    ClpSimplex * model3 = new ClpSimplex(*model2);
    CoinBuild build;
    double one=1.0;
    int numberColumns = model3->numberColumns();
    const double * columnLower = model3->columnLower();
    const double * columnUpper = model3->columnUpper();
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (columnUpper[iColumn]<1.0e20&&
          columnLower[iColumn]>-1.0e20) {
        if (fabs(columnLower[iColumn])<fabs(columnUpper[iColumn])) {
          double value = columnUpper[iColumn];
          model3->setColumnUpper(iColumn,COIN_DBL_MAX);
          build.addRow(1,&iColumn,&one,-COIN_DBL_MAX,value);
        } else {
          double value = columnLower[iColumn];
          model3->setColumnLower(iColumn,-COIN_DBL_MAX);
          build.addRow(1,&iColumn,&one,value,COIN_DBL_MAX);
        }
      }
    }
    model3->addRows(build);
    model2=model3;
  }
  int numberColumns = model2->numberColumns();
  const double * columnLower = model2->columnLower();
  const double * columnUpper = model2->columnUpper();
  int numberRows = model2->numberRows();
  double * rowLower = CoinCopyOfArray(model2->rowLower(),numberRows);
  double * rowUpper = CoinCopyOfArray(model2->rowUpper(),numberRows);

  const double * objective = model2->objective();
  CoinPackedMatrix * matrix = model2->matrix();
  // get transpose
  CoinPackedMatrix rowCopy = *matrix;
  const int * row = matrix->getIndices();
  const int * columnLength = matrix->getVectorLengths();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const double * elementByColumn = matrix->getElements();
  double objOffset=0.0;
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    double offset=0.0;
    double objValue =optimizationDirection_*objective[iColumn];
    if (columnUpper[iColumn]>1.0e20) {
      if (columnLower[iColumn]>-1.0e20)
        offset=columnLower[iColumn];
    } else if (columnLower[iColumn]<-1.0e20) {
      offset=columnUpper[iColumn];
    } else {
      // taken care of before
      abort();
    }
    if (offset) {
      objOffset += offset*objValue;
      for (CoinBigIndex j=columnStart[iColumn];
           j<columnStart[iColumn]+columnLength[iColumn];j++) {
        int iRow = row[j];
        if (rowLower[iRow]>-1.0e20)
          rowLower[iRow] -= offset*elementByColumn[j];
        if (rowUpper[iRow]<1.0e20)
          rowUpper[iRow] -= offset*elementByColumn[j];
      }
    }
  }
  int * which = new int[numberRows+numberExtraRows];
  rowCopy.reverseOrdering();
  rowCopy.transpose();
  double * fromRowsLower = new double[numberRows+numberExtraRows];
  double * fromRowsUpper = new double[numberRows+numberExtraRows];
  double * newObjective = new double[numberRows+numberExtraRows];
  double * fromColumnsLower = new double[numberColumns];
  double * fromColumnsUpper = new double[numberColumns];
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    double objValue =optimizationDirection_*objective[iColumn];
    // Offset is already in
    if (columnUpper[iColumn]>1.0e20) {
      if (columnLower[iColumn]>-1.0e20) {
        fromColumnsLower[iColumn]=-COIN_DBL_MAX;
        fromColumnsUpper[iColumn]=objValue;
      } else {
        // free
        fromColumnsLower[iColumn]=objValue;
        fromColumnsUpper[iColumn]=objValue;
      }
    } else if (columnLower[iColumn]<-1.0e20) {
      fromColumnsLower[iColumn]=objValue;
      fromColumnsUpper[iColumn]=COIN_DBL_MAX;
    } else {
      abort();
    }
  }
  int kRow=0;
  int kExtraRow=numberRows;
  for (iRow=0;iRow<numberRows;iRow++) {
    if (rowLower[iRow]<-1.0e20) {
      assert (rowUpper[iRow]<1.0e20);
      newObjective[kRow]=-rowUpper[iRow];
      fromRowsLower[kRow]=-COIN_DBL_MAX;
      fromRowsUpper[kRow]=0.0;
      which[kRow]=iRow;
      kRow++;
    } else if (rowUpper[iRow]>1.0e20) {
      newObjective[kRow]=-rowLower[iRow];
      fromRowsLower[kRow]=0.0;
      fromRowsUpper[kRow]=COIN_DBL_MAX;
      which[kRow]=iRow;
      kRow++;
    } else {
      if (rowUpper[iRow]==rowLower[iRow]) {
        newObjective[kRow]=-rowLower[iRow];
        fromRowsLower[kRow]=-COIN_DBL_MAX;;
        fromRowsUpper[kRow]=COIN_DBL_MAX;
        which[kRow]=iRow;
        kRow++;
      } else {
        // range
        newObjective[kRow]=-rowUpper[iRow];
        fromRowsLower[kRow]=-COIN_DBL_MAX;
        fromRowsUpper[kRow]=0.0;
        which[kRow]=iRow;
        kRow++;
        newObjective[kExtraRow]=-rowLower[iRow];
        fromRowsLower[kExtraRow]=0.0;
        fromRowsUpper[kExtraRow]=COIN_DBL_MAX;
        which[kExtraRow]=iRow;
        kExtraRow++;
      }
    }
  }
  if (numberExtraRows) {
    CoinPackedMatrix newCopy;
    newCopy.setExtraGap(0.0);
    newCopy.setExtraMajor(0.0);
    newCopy.submatrixOfWithDuplicates(rowCopy,kExtraRow,which);
    rowCopy=newCopy;
  }
  ClpSimplex * modelDual = new ClpSimplex();
  modelDual->loadProblem(rowCopy,fromRowsLower,fromRowsUpper,newObjective,
                        fromColumnsLower,fromColumnsUpper);
  modelDual->setObjectiveOffset(objOffset);
  modelDual->setDualBound(model2->dualBound());
  modelDual->setInfeasibilityCost(model2->infeasibilityCost());
  modelDual->setDualTolerance(model2->dualTolerance());
  modelDual->setPrimalTolerance(model2->primalTolerance());
  modelDual->setPerturbation(model2->perturbation());
  modelDual->setSpecialOptions(model2->specialOptions());
  modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
  delete [] fromRowsLower;
  delete [] fromRowsUpper;
  delete [] fromColumnsLower;
  delete [] fromColumnsUpper;
  delete [] newObjective;
  delete [] which;
  delete [] rowLower;
  delete [] rowUpper;
  if (changed)
    delete model2;
  modelDual->createStatus();
  return modelDual;
}
// Restores solution from dualized problem
int
00942 ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
{
  int returnCode=0;;
  createStatus();
  // Number of rows in dual problem was original number of columns
  assert (numberColumns_==dualProblem->numberRows());
  // If slack on d-row basic then column at bound otherwise column basic
  // If d-column basic then rhs tight
  int numberBasic=0;
  int iRow,iColumn=0;
  // Get number of extra rows from ranges
  int numberExtraRows=0;
  for (iRow=0;iRow<numberRows_;iRow++) {
    if (rowLower_[iRow]>-1.0e20&&
        rowUpper_[iRow]<1.0e20) {
      if (rowUpper_[iRow]!=rowLower_[iRow])
         numberExtraRows++;
    }
  }
  const double * objective = this->objective();
  const double * dualDual = dualProblem->dualRowSolution();
  const double * dualDj = dualProblem->dualColumnSolution();
  const double * dualSol = dualProblem->primalColumnSolution();
  const double * dualActs = dualProblem->primalRowSolution();
#if 0
  ClpSimplex thisCopy=*this;
  thisCopy.dual(); // for testing
  const double * primalDual = thisCopy.dualRowSolution();
  const double * primalDj = thisCopy.dualColumnSolution();
  const double * primalSol = thisCopy.primalColumnSolution();
  const double * primalActs = thisCopy.primalRowSolution();
  char ss[]={'F','B','U','L','S','F'};
  printf ("Dual problem row info %d rows\n",dualProblem->numberRows());
  for (iRow=0;iRow<dualProblem->numberRows();iRow++)
    printf("%d at %c primal %g dual %g\n",
           iRow,ss[dualProblem->getRowStatus(iRow)],
           dualActs[iRow],dualDual[iRow]);
  printf ("Dual problem column info %d columns\n",dualProblem->numberColumns());
  for (iColumn=0;iColumn<dualProblem->numberColumns();iColumn++)
    printf("%d at %c primal %g dual %g\n",
           iColumn,ss[dualProblem->getColumnStatus(iColumn)],
           dualSol[iColumn],dualDj[iColumn]);
  printf ("Primal problem row info %d rows\n",thisCopy.numberRows());
  for (iRow=0;iRow<thisCopy.numberRows();iRow++)
    printf("%d at %c primal %g dual %g\n",
           iRow,ss[thisCopy.getRowStatus(iRow)],
           primalActs[iRow],primalDual[iRow]);
  printf ("Primal problem column info %d columns\n",thisCopy.numberColumns());
  for (iColumn=0;iColumn<thisCopy.numberColumns();iColumn++)
    printf("%d at %c primal %g dual %g\n",
           iColumn,ss[thisCopy.getColumnStatus(iColumn)],
           primalSol[iColumn],primalDj[iColumn]);
#endif
  // position at bound information
  int jColumn=numberRows_;
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    double objValue =optimizationDirection_*objective[iColumn];
    Status status = dualProblem->getRowStatus(iColumn);
    double otherValue = COIN_DBL_MAX;
    if (columnUpper_[iColumn]<1.0e20&&
        columnLower_[iColumn]>-1.0e20) {
        if (fabs(columnLower_[iColumn])<fabs(columnUpper_[iColumn])) {
          otherValue = columnUpper_[iColumn]+dualDj[jColumn];
        } else {
          otherValue = columnLower_[iColumn]+dualDj[jColumn];
        }
        jColumn++;
    }
    if (status==basic) {
      // column is at bound
      if (otherValue==COIN_DBL_MAX) {
        reducedCost_[iColumn]=objValue-dualActs[iColumn];
        if (columnUpper_[iColumn]>1.0e20) {
        if (columnLower_[iColumn]>-1.0e20) {
          if (columnUpper_[iColumn]>columnLower_[iColumn])
            setColumnStatus(iColumn,atLowerBound);
          else
            setColumnStatus(iColumn,isFixed);
          columnActivity_[iColumn]=columnLower_[iColumn];
        } else {
          // free
          setColumnStatus(iColumn,isFree);
          columnActivity_[iColumn]=0.0;
        }
        } else {
          setColumnStatus(iColumn,atUpperBound);
          columnActivity_[iColumn]=columnUpper_[iColumn];
        }
      } else {
        reducedCost_[iColumn]=objValue-dualActs[iColumn];
        //printf("other dual sol %g\n",otherValue);
        if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
        if (columnUpper_[iColumn]>columnLower_[iColumn])
          setColumnStatus(iColumn,atLowerBound);
        else
          setColumnStatus(iColumn,isFixed);
          columnActivity_[iColumn]=columnLower_[iColumn];
        } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
        if (columnUpper_[iColumn]>columnLower_[iColumn])
          setColumnStatus(iColumn,atUpperBound);
        else
          setColumnStatus(iColumn,isFixed);
          columnActivity_[iColumn]=columnUpper_[iColumn];
        } else {
          abort();
        }
      }
    } else {
      if (otherValue==COIN_DBL_MAX) {
        // column basic
        setColumnStatus(iColumn,basic);
        numberBasic++;
      if (columnLower_[iColumn]>-1.0e20) {
        columnActivity_[iColumn]=-dualDual[iColumn] + columnLower_[iColumn];
      } else if (columnUpper_[iColumn]<1.0e20) {
        columnActivity_[iColumn]=-dualDual[iColumn] + columnUpper_[iColumn];
      } else {
        columnActivity_[iColumn]=-dualDual[iColumn];
      }
        reducedCost_[iColumn]=0.0;
      } else {
        // may be at other bound
        //printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
        if (dualProblem->getColumnStatus(jColumn-1)!=basic) {
          // column basic
          setColumnStatus(iColumn,basic);
          numberBasic++;
        //printf("Col %d otherV %g dualDual %g\n",iColumn,
        // otherValue,dualDual[iColumn]);
          columnActivity_[iColumn]=-dualDual[iColumn];
          columnActivity_[iColumn]=otherValue;
          reducedCost_[iColumn]=0.0;
        } else {
          reducedCost_[iColumn]=objValue-dualActs[iColumn];
          if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
          if (columnUpper_[iColumn]>columnLower_[iColumn])
            setColumnStatus(iColumn,atLowerBound);
          else
            setColumnStatus(iColumn,isFixed);
            columnActivity_[iColumn]=columnLower_[iColumn];
          } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
          if (columnUpper_[iColumn]>columnLower_[iColumn])
            setColumnStatus(iColumn,atUpperBound);
          else
            setColumnStatus(iColumn,isFixed);
            columnActivity_[iColumn]=columnUpper_[iColumn];
          } else {
            abort();
          }
        }
      }
    }
  }
  // now rows
  int kExtraRow=jColumn;
  int numberRanges=0;
  for (iRow=0;iRow<numberRows_;iRow++) {
    Status status = dualProblem->getColumnStatus(iRow);
    if (status==basic) {
      // row is at bound
      dual_[iRow]=dualSol[iRow];;
    } else {
      // row basic
      setRowStatus(iRow,basic);
      numberBasic++;
      dual_[iRow]=0.0;
    }
    if (rowLower_[iRow]<-1.0e20) {
      if (status==basic) {
        rowActivity_[iRow]=rowUpper_[iRow];
        setRowStatus(iRow,atUpperBound);
      } else {
      assert (dualDj[iRow]<1.0e-5);
        rowActivity_[iRow]=rowUpper_[iRow]+dualDj[iRow];
      }        
    } else if (rowUpper_[iRow]>1.0e20) {
      if (status==basic) {
        rowActivity_[iRow]=rowLower_[iRow];
        setRowStatus(iRow,atLowerBound);
      } else {
        rowActivity_[iRow]=rowLower_[iRow]+dualDj[iRow];
      assert (dualDj[iRow]>-1.0e-5);
      }        
    } else {
      if (rowUpper_[iRow]==rowLower_[iRow]) {
        rowActivity_[iRow]=rowLower_[iRow];
        if (status==basic) {
          setRowStatus(iRow,isFixed);
        }        
      } else {
        // range
        numberRanges++;
      Status statusL = dualProblem->getColumnStatus(kExtraRow);
      //printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
      //     iRow,status,kExtraRow,statusL, dualSol[iRow],
      //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
      if (status==basic) {
        assert (statusL!=basic);
        rowActivity_[iRow]=rowUpper_[iRow];
        setRowStatus(iRow,atUpperBound);
      } else if (statusL==basic) {
        numberBasic--; // already counted
        rowActivity_[iRow]=rowLower_[iRow];
        setRowStatus(iRow,atLowerBound);
        dual_[iRow]=dualSol[kExtraRow];;
      } else {
        rowActivity_[iRow]=rowLower_[iRow]-dualDj[iRow];
        assert (dualDj[iRow]<1.0e-5);
        // row basic
        //setRowStatus(iRow,basic);
        //numberBasic++;
        dual_[iRow]=0.0;
      }
      kExtraRow++;
      }
    }
  }
  if (numberBasic!=numberRows_) {
    printf("Bad basis - ranges - coding needed\n");
    assert (numberRanges);
    abort();
  }
  if (optimizationDirection_<0.0) {
    for (iRow=0;iRow<numberRows_;iRow++) {
      dual_[iRow]=-dual_[iRow];
    }
  }
  // redo row activities
  memset(rowActivity_,0,numberRows_*sizeof(double));
  matrix_->times(1.0,columnActivity_,rowActivity_);
  // redo reduced costs
  memcpy(reducedCost_,this->objective(),numberColumns_*sizeof(double));
  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
  checkSolutionInternal();
  if (sumDualInfeasibilities_>1.0e-5||sumPrimalInfeasibilities_>1.0e-5) {
    returnCode=1;
#ifdef CLP_INVESTIGATE
    printf("There are %d dual infeasibilities summing to %g ",
         numberDualInfeasibilities_,sumDualInfeasibilities_);
    printf("and %d primal infeasibilities summing to %g\n",
         numberPrimalInfeasibilities_,sumPrimalInfeasibilities_);
#endif
  }
  // Below will go to ..DEBUG later
#if 1 //ndef NDEBUG
  // Check if correct
  double * columnActivity = CoinCopyOfArray(columnActivity_,numberColumns_);
  double * rowActivity = CoinCopyOfArray(rowActivity_,numberRows_);
  double * reducedCost = CoinCopyOfArray(reducedCost_,numberColumns_);
  double * dual = CoinCopyOfArray(dual_,numberRows_);
  this->dual(); //primal();
  CoinRelFltEq eq(1.0e-5);
  for (iRow=0;iRow<numberRows_;iRow++) {
    assert(eq(dual[iRow],dual_[iRow]));
  }
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    assert(eq(columnActivity[iColumn],columnActivity_[iColumn]));
  }
  for (iRow=0;iRow<numberRows_;iRow++) {
    assert(eq(rowActivity[iRow],rowActivity_[iRow]));
  }
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    assert(eq(reducedCost[iColumn],reducedCost_[iColumn]));
  }
  delete [] columnActivity;
  delete [] rowActivity;
  delete [] reducedCost;
  delete [] dual;
#endif
  return returnCode;
}
/* Does very cursory presolve.
   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
*/
ClpSimplex * 
01217 ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
                        int & nBound, bool moreBounds, bool tightenBounds)
{
  //#define CHECK_STATUS
#ifdef CHECK_STATUS
  {
    int n=0;
    int i;
    for (i=0;i<numberColumns_;i++)
      if (getColumnStatus(i)==ClpSimplex::basic)
        n++;
    for (i=0;i<numberRows_;i++)
      if (getRowStatus(i)==ClpSimplex::basic)
        n++;
    assert (n==numberRows_);
  }
#endif
  
  const double * element = matrix_->getElements();
  const int * row = matrix_->getIndices();
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * columnLength = matrix_->getVectorLengths();

  CoinZeroN(rhs,numberRows_);
  int iColumn;
  int iRow;
  CoinZeroN(whichRow,numberRows_);
  int * backColumn = whichColumn+numberColumns_;
  int numberRows2=0;
  int numberColumns2=0;
  double offset=0.0;
  const double * objective = this->objective();
  double * solution = columnActivity_;
  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    double lower = columnLower_[iColumn];
    double upper = columnUpper_[iColumn];
    if (upper>lower||getColumnStatus(iColumn)==ClpSimplex::basic) {
      backColumn[iColumn]=numberColumns2;
      whichColumn[numberColumns2++]=iColumn;
      for (CoinBigIndex j = columnStart[iColumn];
           j<columnStart[iColumn]+columnLength[iColumn];j++) {
        int iRow = row[j];
        int n=whichRow[iRow];
        if (n==0&&element[j])
          whichRow[iRow]=-iColumn-1;
        else if (n<0) 
          whichRow[iRow]=2;
      }
    } else {
      // fixed
      backColumn[iColumn]=-1;
      solution[iColumn]=upper;
      if (upper) {
        offset += objective[iColumn]*upper;
        for (CoinBigIndex j = columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
          int iRow = row[j];
          double value = element[j];
          rhs[iRow] += upper*value;
        }
      }
    }
  }
  int returnCode=0;
  double tolerance = primalTolerance();
  nBound=2*numberRows_;
  for (iRow=0;iRow<numberRows_;iRow++) {
    int n=whichRow[iRow];
    if (n>0) {
      whichRow[numberRows2++]=iRow;
    } else if (n<0) {
      //whichRow[numberRows2++]=iRow;
      //continue;
      // Can only do in certain circumstances as we don't know current value
      if (rowLower_[iRow]==rowUpper_[iRow]||getRowStatus(iRow)==ClpSimplex::basic) {
        // save row and column for bound
        whichRow[--nBound]=iRow;
        whichRow[nBound+numberRows_]=-n-1;
      } else if (moreBounds) {
        // save row and column for bound
        whichRow[--nBound]=iRow;
        whichRow[nBound+numberRows_]=-n-1;
      } else {
        whichRow[numberRows2++]=iRow;
      }
    } else {
      // empty
      double rhsValue = rhs[iRow];
      if (rhsValue<rowLower_[iRow]-tolerance||rhsValue>rowUpper_[iRow]+tolerance) {
        returnCode=1; // infeasible
      }
    }
  }
  ClpSimplex * small=NULL;
  if (!returnCode) {
    small = new ClpSimplex(this,numberRows2,whichRow,
                     numberColumns2,whichColumn,true,false);
#if 0
    ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
    if (rowCopy) {
      assert(!small->rowCopy());
      small->setNewRowCopy(new ClpPackedMatrix(*rowCopy,numberRows2,whichRow,
                                    numberColumns2,whichColumn));
    }
#endif
    // Set some stuff
    small->setDualBound(dualBound_);
    small->setInfeasibilityCost(infeasibilityCost_);
    small->setSpecialOptions(specialOptions_);
    small->setPerturbation(perturbation_);
    small->defaultFactorizationFrequency();
    small->setAlphaAccuracy(alphaAccuracy_);
    // If no rows left then no tightening!
    if (!numberRows2||!numberColumns2) 
      tightenBounds=false;

    int numberElements=getNumElements();
    int numberElements2=small->getNumElements();
    small->setObjectiveOffset(objectiveOffset()-offset);
    handler_->message(CLP_CRUNCH_STATS,messages_)
      <<numberRows2<< -(numberRows_ - numberRows2)
      <<numberColumns2<< -(numberColumns_ - numberColumns2)
      <<numberElements2<< -(numberElements - numberElements2)
      <<CoinMessageEol;
    // And set objective value to match
    small->setObjectiveValue(this->objectiveValue());
    double * rowLower2 = small->rowLower();
    double * rowUpper2 = small->rowUpper();
    int jRow;
    for (jRow=0;jRow<numberRows2;jRow++) {
      iRow = whichRow[jRow];
      if (rowLower2[jRow]>-1.0e20)
        rowLower2[jRow] -= rhs[iRow];
      if (rowUpper2[jRow]<1.0e20)
        rowUpper2[jRow] -= rhs[iRow];
    }
    // and bounds
    double * columnLower2 = small->columnLower();
    double * columnUpper2 = small->columnUpper();
    const char * integerInformation = integerType_;
    for (jRow=nBound;jRow<2*numberRows_;jRow++) {
      iRow = whichRow[jRow];
      iColumn = whichRow[jRow+numberRows_];
      double lowerRow = rowLower_[iRow];
      if (lowerRow>-1.0e20)
        lowerRow -= rhs[iRow];
      double upperRow = rowUpper_[iRow];
      if (upperRow<1.0e20)
        upperRow -= rhs[iRow];
      int jColumn = backColumn[iColumn];
      double lower = columnLower2[jColumn];
      double upper = columnUpper2[jColumn];
      double value=0.0;
      for (CoinBigIndex j = columnStart[iColumn];
           j<columnStart[iColumn]+columnLength[iColumn];j++) {
        if (iRow==row[j]) {
          value=element[j];
          break;
        }
      }
      assert (value);
      // convert rowLower and Upper to implied bounds on column
      double newLower=-COIN_DBL_MAX;
      double newUpper=COIN_DBL_MAX;
      if (value>0.0) {
        if (lowerRow>-1.0e20)
          newLower = lowerRow/value;
        if (upperRow<1.0e20)
          newUpper = upperRow/value;
      } else {
        if (upperRow<1.0e20)
          newLower = upperRow/value;
        if (lowerRow>-1.0e20)
          newUpper = lowerRow/value;
      }
      if (integerInformation&&integerInformation[iColumn]) {
        if (newLower-floor(newLower)<10.0*tolerance) 
          newLower=floor(newLower);
        else
          newLower=ceil(newLower);
        if (ceil(newUpper)-newUpper<10.0*tolerance) 
          newUpper=ceil(newUpper);
        else
          newUpper=floor(newUpper);
      }
      newLower = CoinMax(lower,newLower);
      newUpper = CoinMin(upper,newUpper);
      if (newLower>newUpper+tolerance) {
        //printf("XXYY inf on bound\n");
        returnCode=1;
      }
      columnLower2[jColumn]=newLower;
      columnUpper2[jColumn]=CoinMax(newLower,newUpper);
      if (getRowStatus(iRow)!=ClpSimplex::basic) {
        if (getColumnStatus(iColumn)==ClpSimplex::basic) {
          if (columnLower2[jColumn]==columnUpper2[jColumn]) {
            // can only get here if will be fixed
            small->setColumnStatus(jColumn,ClpSimplex::isFixed);
          } else {
            // solution is valid
            if (fabs(columnActivity_[iColumn]-columnLower2[jColumn])<
                fabs(columnActivity_[iColumn]-columnUpper2[jColumn]))
              small->setColumnStatus(jColumn,ClpSimplex::atLowerBound);
            else
              small->setColumnStatus(jColumn,ClpSimplex::atUpperBound);
          }
        } else {
          //printf("what now neither basic\n");
        }
      }
    }
    if (returnCode) {
      delete small;
      small=NULL;
    } else if (tightenBounds&&integerInformation) {
      // See if we can tighten any bounds
      // use rhs for upper and small duals for lower
      double * up = rhs;
      double * lo = small->dualRowSolution();
      const double * element = small->clpMatrix()->getElements();
      const int * row = small->clpMatrix()->getIndices();
      const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
      //const int * columnLength = small->clpMatrix()->getVectorLengths();
      CoinZeroN(lo,numberRows2);
      CoinZeroN(up,numberRows2);
      for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
        double upper=columnUpper2[iColumn];
        double lower=columnLower2[iColumn];
        //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
          int iRow=row[j];
          double value = element[j];
          if (value>0.0) {
            if (upper<1.0e20)
              up[iRow] += upper*value;
            else
              up[iRow] = COIN_DBL_MAX;
            if (lower>-1.0e20)
              lo[iRow] += lower*value;
            else
              lo[iRow] = -COIN_DBL_MAX;
          } else {
            if (upper<1.0e20)
              lo[iRow] += upper*value;
            else
              lo[iRow] = -COIN_DBL_MAX;
            if (lower>-1.0e20)
              up[iRow] += lower*value;
            else
              up[iRow] = COIN_DBL_MAX;
          }
        }
      }
      double * rowLower2 = small->rowLower();
      double * rowUpper2 = small->rowUpper();
      bool feasible=true;
      // make safer
      for (int iRow=0;iRow<numberRows2;iRow++) {
        double lower = lo[iRow];
        if (lower>rowUpper2[iRow]+tolerance) {
          feasible=false;
          break;
        } else {
          lo[iRow] = CoinMin(lower-rowUpper2[iRow],0.0)-tolerance;
        }
        double upper = up[iRow];
        if (upper<rowLower2[iRow]-tolerance) {
          feasible=false;
          break;
        } else {
          up[iRow] = CoinMax(upper-rowLower2[iRow],0.0)+tolerance;
        }
      }
      if (!feasible) {
        delete small;
        small=NULL;
      } else {
        // and tighten
        for (int iColumn=0;iColumn<numberColumns2;iColumn++) {
          if (integerInformation[whichColumn[iColumn]]) {
            double upper=columnUpper2[iColumn];
            double lower=columnLower2[iColumn];
            double newUpper = upper;
            double newLower = lower;
            double difference = upper-lower;
            if (lower>-1000.0&&upper<1000.0) {
              for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
                int iRow=row[j];
                double value = element[j];
                if (value>0.0) {
                  double upWithOut = up[iRow] - value*difference;
                  if (upWithOut<0.0) {
                    newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
                  }
                  double lowWithOut = lo[iRow] + value*difference;
                  if (lowWithOut>0.0) {
                    newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
                  }
                } else {
                  double upWithOut = up[iRow] + value*difference;
                  if (upWithOut<0.0) {
                    newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
                  }
                  double lowWithOut = lo[iRow] - value*difference;
                  if (lowWithOut>0.0) {
                    newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
                  }
                }
              }
              if (newLower>lower||newUpper<upper) {
                if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
                  newUpper = floor(newUpper);
                else
                  newUpper = floor(newUpper+0.5);
                if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
                  newLower = ceil(newLower);
                else
                  newLower = ceil(newLower-0.5);
                // change may be too small - check
                if (newLower>lower||newUpper<upper) {
                  if (newUpper>=newLower) {
                    // Could also tighten in this
                    //printf("%d bounds %g %g tightened to %g %g\n",
                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
                    //     newLower,newUpper);
#if 1
                    columnUpper2[iColumn]=newUpper;
                    columnLower2[iColumn]=newLower;
                    columnUpper_[whichColumn[iColumn]]=newUpper;
                    columnLower_[whichColumn[iColumn]]=newLower;
#endif
                    // and adjust bounds on rows
                    newUpper -= upper;
                    newLower -= lower;
                    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
                      int iRow=row[j];
                      double value = element[j];
                      if (value>0.0) {
                        up[iRow] += newUpper*value;
                        lo[iRow] += newLower*value;
                      } else {
                        lo[iRow] += newUpper*value;
                        up[iRow] += newLower*value;
                      }
                    }
                  } else {
                    // infeasible
                    //printf("%d bounds infeasible %g %g tightened to %g %g\n",
                    //     iColumn,columnLower2[iColumn],columnUpper2[iColumn],
                    //     newLower,newUpper);
#if 1
                    delete small;
                    small=NULL;
                    break;
#endif
                  }
                }
              }
            }
          }
        }
      }
    }
  }
#if 0
  if (small) {
    static int which=0;
    which++;
    char xxxx[20];
    sprintf(xxxx,"bad%d.mps",which);
    small->writeMps(xxxx,0,1);
    sprintf(xxxx,"largebad%d.mps",which);
    writeMps(xxxx,0,1);
    printf("bad%d %x old size %d %d new %d %d\n",which,small,
         numberRows_,numberColumns_,small->numberRows(),small->numberColumns());
#if 0
    for (int i=0;i<numberColumns_;i++) 
      printf("Bound %d %g %g\n",i,columnLower_[i],columnUpper_[i]);
    for (int i=0;i<numberRows_;i++) 
      printf("Row bound %d %g %g\n",i,rowLower_[i],rowUpper_[i]);
#endif
  }
#endif
#ifdef CHECK_STATUS
  {
    int n=0;
    int i;
    for (i=0;i<small->numberColumns();i++)
      if (small->getColumnStatus(i)==ClpSimplex::basic)
        n++;
    for (i=0;i<small->numberRows();i++)
      if (small->getRowStatus(i)==ClpSimplex::basic)
        n++;
    assert (n==small->numberRows());
  }
#endif
  return small;
}
/* After very cursory presolve.
   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
*/
void 
01619 ClpSimplexOther::afterCrunch(const ClpSimplex & small,
                             const int * whichRow, 
                             const int * whichColumn, int nBound)
{
#ifndef NDEBUG
  for (int i=0;i<small.numberRows();i++)
    assert (whichRow[i]>=0&&whichRow[i]<numberRows_);
  for (int i=0;i<small.numberColumns();i++)
    assert (whichColumn[i]>=0&&whichColumn[i]<numberColumns_);
#endif
  getbackSolution(small,whichRow,whichColumn);
  // and deal with status for bounds
  const double * element = matrix_->getElements();
  const int * row = matrix_->getIndices();
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * columnLength = matrix_->getVectorLengths();
  double tolerance = primalTolerance();
  double djTolerance = dualTolerance();
  for (int jRow=nBound;jRow<2*numberRows_;jRow++) {
    int iRow = whichRow[jRow];
    int iColumn = whichRow[jRow+numberRows_];
    if (getColumnStatus(iColumn)!=ClpSimplex::basic) {
      double lower = columnLower_[iColumn];
      double upper = columnUpper_[iColumn];
      double value = columnActivity_[iColumn];
      double djValue = reducedCost_[iColumn];
      dual_[iRow]=0.0;
      if (upper>lower) {
        if (value<lower+tolerance&&djValue>-djTolerance) {
          setColumnStatus(iColumn,ClpSimplex::atLowerBound);
          setRowStatus(iRow,ClpSimplex::basic);
        } else if (value>upper-tolerance&&djValue<djTolerance) {
          setColumnStatus(iColumn,ClpSimplex::atUpperBound);
          setRowStatus(iRow,ClpSimplex::basic);
        } else {
          // has to be basic
          setColumnStatus(iColumn,ClpSimplex::basic);
          reducedCost_[iColumn] = 0.0;
          double value=0.0;
          for (CoinBigIndex j = columnStart[iColumn];
               j<columnStart[iColumn]+columnLength[iColumn];j++) {
            if (iRow==row[j]) {
              value=element[j];
              break;
            }
          }
          dual_[iRow]=djValue/value;
          if (rowUpper_[iRow]>rowLower_[iRow]) {
            if (fabs(rowActivity_[iRow]-rowLower_[iRow])<
                fabs(rowActivity_[iRow]-rowUpper_[iRow]))
              setRowStatus(iRow,ClpSimplex::atLowerBound);
            else
              setRowStatus(iRow,ClpSimplex::atUpperBound);
          } else {
            setRowStatus(iRow,ClpSimplex::isFixed);
          }
        }
      } else {
        // row can always be basic
        setRowStatus(iRow,ClpSimplex::basic);
      }
    } else {
      // row can always be basic
      setRowStatus(iRow,ClpSimplex::basic);
    }
  }
  //#ifndef NDEBUG
#if 0
  if  (small.status()==0) {
    int n=0;
    int i;
    for (i=0;i<numberColumns;i++)
      if (getColumnStatus(i)==ClpSimplex::basic)
        n++;
    for (i=0;i<numberRows;i++)
      if (getRowStatus(i)==ClpSimplex::basic)
        n++;
    assert (n==numberRows);
  }
#endif
}
/* Tightens integer bounds - returns number tightened or -1 if infeasible
 */
int 
01703 ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
{
  // See if we can tighten any bounds
  // use rhs for upper and small duals for lower
  double * up = rhsSpace;
  double * lo = dual_;
  const double * element = matrix_->getElements();
  const int * row = matrix_->getIndices();
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * columnLength = matrix_->getVectorLengths();
  CoinZeroN(lo,numberRows_);
  CoinZeroN(up,numberRows_);
  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    double upper=columnUpper_[iColumn];
    double lower=columnLower_[iColumn];
    //assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
    for (CoinBigIndex j=columnStart[iColumn];
         j<columnStart[iColumn]+columnLength[iColumn];j++) {
      int iRow=row[j];
      double value = element[j];
      if (value>0.0) {
        if (upper<1.0e20)
          up[iRow] += upper*value;
        else
          up[iRow] = COIN_DBL_MAX;
        if (lower>-1.0e20)
          lo[iRow] += lower*value;
        else
          lo[iRow] = -COIN_DBL_MAX;
      } else {
        if (upper<1.0e20)
          lo[iRow] += upper*value;
        else
          lo[iRow] = -COIN_DBL_MAX;
        if (lower>-1.0e20)
          up[iRow] += lower*value;
        else
          up[iRow] = COIN_DBL_MAX;
      }
    }
  }
  bool feasible=true;
  // make safer
  double tolerance = primalTolerance();
  for (int iRow=0;iRow<numberRows_;iRow++) {
    double lower = lo[iRow];
    if (lower>rowUpper_[iRow]+tolerance) {
      feasible=false;
      break;
    } else {
      lo[iRow] = CoinMin(lower-rowUpper_[iRow],0.0)-tolerance;
    }
    double upper = up[iRow];
    if (upper<rowLower_[iRow]-tolerance) {
      feasible=false;
      break;
    } else {
      up[iRow] = CoinMax(upper-rowLower_[iRow],0.0)+tolerance;
    }
  }
  int numberTightened=0;
  if (!feasible) {
    return -1;
  } else if (integerType_) {
    // and tighten
    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
      if (integerType_[iColumn]) {
        double upper=columnUpper_[iColumn];
        double lower=columnLower_[iColumn];
        double newUpper = upper;
        double newLower = lower;
        double difference = upper-lower;
        if (lower>-1000.0&&upper<1000.0) {
          for (CoinBigIndex j=columnStart[iColumn];
               j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow=row[j];
            double value = element[j];
            if (value>0.0) {
              double upWithOut = up[iRow] - value*difference;
              if (upWithOut<0.0) {
                newLower = CoinMax(newLower,lower-(upWithOut+tolerance)/value);
              }
              double lowWithOut = lo[iRow] + value*difference;
              if (lowWithOut>0.0) {
                newUpper = CoinMin(newUpper,upper-(lowWithOut-tolerance)/value);
              }
            } else {
              double upWithOut = up[iRow] + value*difference;
              if (upWithOut<0.0) {
                newUpper = CoinMin(newUpper,upper-(upWithOut+tolerance)/value);
              }
              double lowWithOut = lo[iRow] - value*difference;
              if (lowWithOut>0.0) {
                newLower = CoinMax(newLower,lower-(lowWithOut-tolerance)/value);
              }
            }
          }
          if (newLower>lower||newUpper<upper) {
            if (fabs(newUpper-floor(newUpper+0.5))>1.0e-6)
              newUpper = floor(newUpper);
            else
              newUpper = floor(newUpper+0.5);
            if (fabs(newLower-ceil(newLower-0.5))>1.0e-6)
              newLower = ceil(newLower);
            else
              newLower = ceil(newLower-0.5);
            // change may be too small - check
            if (newLower>lower||newUpper<upper) {
              if (newUpper>=newLower) {
                numberTightened++;
                //printf("%d bounds %g %g tightened to %g %g\n",
                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
                //     newLower,newUpper);
                columnUpper_[iColumn]=newUpper;
                columnLower_[iColumn]=newLower;
                // and adjust bounds on rows
                newUpper -= upper;
                newLower -= lower;
                for (CoinBigIndex j=columnStart[iColumn];
                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
                  int iRow=row[j];
                  double value = element[j];
                  if (value>0.0) {
                    up[iRow] += newUpper*value;
                    lo[iRow] += newLower*value;
                  } else {
                    lo[iRow] += newUpper*value;
                    up[iRow] += newLower*value;
                  }
                }
              } else {
                // infeasible
                //printf("%d bounds infeasible %g %g tightened to %g %g\n",
                //     iColumn,columnLower_[iColumn],columnUpper_[iColumn],
                //     newLower,newUpper);
                return -1;
              }
            }
          }
        }
      }
    }
  }
  return numberTightened;
}
/* Parametrics
   This is an initial slow version.
   The code uses current bounds + theta * change (if change array not NULL)
   and similarly for objective.
   It starts at startingTheta and returns ending theta in endingTheta.
   If reportIncrement 0.0 it will report on any movement
   If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
   If it can not reach input endingTheta return code will be 1 for infeasible,
   2 for unbounded, if error on ranges -1,  otherwise 0.
   Normal report is just theta and objective but
   if event handler exists it may do more
   On exit endingTheta is maximum reached (can be used for next startingTheta)
*/
int 
01862 ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,double reportIncrement,
                             const double * changeLowerBound, const double * changeUpperBound,
                             const double * changeLowerRhs, const double * changeUpperRhs,
                             const double * changeObjective)
{
  bool needToDoSomething=true;
  bool canTryQuick = (reportIncrement) ? true : false;
  // Save copy of model
  ClpSimplex copyModel = *this;
  int savePerturbation = perturbation_;
  perturbation_=102; // switch off
  while (needToDoSomething) {
    needToDoSomething=false;
    algorithm_ = -1;
    
    // save data
    ClpDataSave data = saveData();
    int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0,NULL,0);
    int iRow,iColumn;
    double * chgUpper=NULL;
    double * chgLower=NULL;
    double * chgObjective=NULL;
    
    // Dantzig (as will not be used) (out later)
    ClpDualRowPivot * savePivot = dualRowPivot_;
    dualRowPivot_ = new ClpDualRowDantzig();
    
    if (!returnCode) {
      // Find theta when bounds will cross over and create arrays
      int numberTotal = numberRows_+numberColumns_;
      chgLower = new double[numberTotal];
      memset(chgLower,0,numberTotal*sizeof(double));
      chgUpper = new double[numberTotal];
      memset(chgUpper,0,numberTotal*sizeof(double));
      chgObjective = new double[numberTotal];
      memset(chgObjective,0,numberTotal*sizeof(double));
      assert (!rowScale_);
      double maxTheta=1.0e50;
      if (changeLowerRhs||changeUpperRhs) {
        for (iRow=0;iRow<numberRows_;iRow++) {
          double lower = rowLower_[iRow];
          double upper = rowUpper_[iRow];
          if (lower>upper) {
            maxTheta=-1.0;
            break;
          }
          double changeLower = (changeLowerRhs) ? changeLowerRhs[iRow] : 0.0;
          double changeUpper = (changeUpperRhs) ? changeUpperRhs[iRow] : 0.0;
          if (lower>-1.0e20&&upper<1.0e20) {
            if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
              maxTheta = (upper-lower)/(changeLower-changeUpper);
            }
          }
          if (lower>-1.0e20) {
            lower_[numberColumns_+iRow] += startingTheta*changeLower;
            chgLower[numberColumns_+iRow]=changeLower;
          }
          if (upper<1.0e20) {
            upper_[numberColumns_+iRow] += startingTheta*changeUpper;
            chgUpper[numberColumns_+iRow]=changeUpper;
          }
        }
      }
      if (maxTheta>0.0) {
        if (changeLowerBound||changeUpperBound) {
          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
            double lower = columnLower_[iColumn];
            double upper = columnUpper_[iColumn];
            if (lower>upper) {
              maxTheta=-1.0;
              break;
            }
            double changeLower = (changeLowerBound) ? changeLowerBound[iColumn] : 0.0;
            double changeUpper = (changeUpperBound) ? changeUpperBound[iColumn] : 0.0;
            if (lower>-1.0e20&&upper<1.0e20) {
              if (lower+maxTheta*changeLower>upper+maxTheta*changeUpper) {
                maxTheta = (upper-lower)/(changeLower-changeUpper);
              }
            }
            if (lower>-1.0e20) {
              lower_[iColumn] += startingTheta*changeLower;
              chgLower[iColumn]=changeLower;
            }
            if (upper<1.0e20) {
              upper_[iColumn] += startingTheta*changeUpper;
              chgUpper[iColumn]=changeUpper;
            }
          }
        }
        if (maxTheta==1.0e50)
          maxTheta = COIN_DBL_MAX;
      }
      if (maxTheta<0.0) {
        // bad ranges or initial
        returnCode = -1;
      }
      endingTheta = CoinMin(endingTheta,maxTheta);
      if (endingTheta<startingTheta) {
        // bad initial
        returnCode = -2;
      }
    }
    double saveEndingTheta=endingTheta;
    if (!returnCode) {
      if (changeObjective) {
        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
          chgObjective[iColumn] = changeObjective[iColumn];
          cost_[iColumn] += startingTheta*changeObjective[iColumn];
        }
      }
      double * saveDuals=NULL;
      reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0,saveDuals,-1,data);
      assert (!problemStatus_);
      // Now do parametrics
      printf("at starting theta of %g objective value is %g\n",startingTheta,
             objectiveValue());
      while (!returnCode) {
        assert (reportIncrement);
        returnCode = parametricsLoop(startingTheta,endingTheta,reportIncrement,
                                     chgLower,chgUpper,chgObjective,data,
                                     canTryQuick);
        if (!returnCode) {
          //double change = endingTheta-startingTheta;
          startingTheta=endingTheta;
          endingTheta = saveEndingTheta;
          //for (int i=0;i<numberTotal;i++) {
          //lower_[i] += change*chgLower[i];
          //upper_[i] += change*chgUpper[i];
          //cost_[i] += change*chgObjective[i];
          //}
          printf("at theta of %g objective value is %g\n",startingTheta,
                 objectiveValue());
          if (startingTheta>=endingTheta)
            break;
        } else if (returnCode==-1) {
          // trouble - do external solve
          needToDoSomething=true;
        } else {
          abort();
        }
      }
    }
    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
    
    delete dualRowPivot_;
    dualRowPivot_ = savePivot;
    // Restore any saved stuff
    restoreData(data);
    if (needToDoSomething) {
      double saveStartingTheta=startingTheta; // known to be feasible
      int cleanedUp=1;
      while (cleanedUp) {
        // tweak
        if (cleanedUp==1) {
          if (!reportIncrement)
            startingTheta = CoinMin(startingTheta+1.0e-5,saveEndingTheta);
          else
            startingTheta = CoinMin(startingTheta+reportIncrement,saveEndingTheta);
        } else {
          // restoring to go slowly
          startingTheta=saveStartingTheta;
        }
        // only works if not scaled
        int i;
        const double * obj1 = objective();
        double * obj2 = copyModel.objective();
        const double * lower1 = columnLower_;
        double * lower2 = copyModel.columnLower();
        const double * upper1 = columnUpper_;
        double * upper2 = copyModel.columnUpper();
        for (i=0;i<numberColumns_;i++) {
          obj2[i] = obj1[i] + startingTheta*chgObjective[i];
          lower2[i] = lower1[i] + startingTheta*chgLower[i];
          upper2[i] = upper1[i] + startingTheta*chgUpper[i];
        }
        lower1 = rowLower_;
        lower2 = copyModel.rowLower();
        upper1 = rowUpper_;
        upper2 = copyModel.rowUpper();
        for (i=0;i<numberRows_;i++) {
          lower2[i] = lower1[i] + startingTheta*chgLower[i+numberColumns_];
          upper2[i] = upper1[i] + startingTheta*chgUpper[i+numberColumns_];
        }
        copyModel.dual();
        if (copyModel.problemStatus()) {
          printf("Can not get to theta of %g\n",startingTheta);
          canTryQuick=false; // do slowly to get exact amount
          // back to last known good
          if (cleanedUp==1)
            cleanedUp=2;
          else
            abort();
        } else {
          // and move stuff back
          int numberTotal = numberRows_+numberColumns_;
          CoinMemcpyN(copyModel.statusArray(),numberTotal,status_);
          CoinMemcpyN(copyModel.primalColumnSolution(),numberColumns_,columnActivity_);
          CoinMemcpyN(copyModel.primalRowSolution(),numberRows_,rowActivity_);
          cleanedUp=0;
        }
      }
    }
    delete [] chgLower;
    delete [] chgUpper;
    delete [] chgObjective;
  }
  perturbation_ = savePerturbation;
  return problemStatus_;
}
int 
02072 ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta,double reportIncrement,
                                 const double * changeLower, const double * changeUpper,
                                 const double * changeObjective, ClpDataSave & data,
                                 bool canTryQuick)
{
  // stuff is already at starting
  // For this crude version just try and go to end
  double change=0.0;
  if (reportIncrement&&canTryQuick) { 
    endingTheta = CoinMin(endingTheta,startingTheta+reportIncrement);
    change = endingTheta-startingTheta;
  }
  int numberTotal = numberRows_+numberColumns_;
  int i;
  for ( i=0;i<numberTotal;i++) {
    lower_[i] += change*changeLower[i];
    upper_[i] += change*changeUpper[i];
    switch(getStatus(i)) {
      
    case basic:
    case isFree:
    case superBasic:
      break;
    case isFixed:
    case atUpperBound:
      solution_[i]=upper_[i];
      break;
    case atLowerBound:
      solution_[i]=lower_[i];
      break;
    }
    cost_[i] += change*changeObjective[i];
  }
  problemStatus_=-1;
  
  // This says whether to restore things etc
  // startup will have factorized so can skip
  int factorType=0;
  // Start check for cycles
  progress_.startCheck();
  // Say change made on first iteration
  changeMade_=1;
  /*
    Status of problem:
    0 - optimal
    1 - infeasible
    2 - unbounded
    -1 - iterating
    -2 - factorization wanted
    -3 - redo checking without factorization
    -4 - looks infeasible
  */
  while (problemStatus_<0) {
    int iRow, iColumn;
    // clear
    for (iRow=0;iRow<4;iRow++) {
      rowArray_[iRow]->clear();
    }    
    
    for (iColumn=0;iColumn<2;iColumn++) {
      columnArray_[iColumn]->clear();
    }    
    
    // give matrix (and model costs and bounds a chance to be
    // refreshed (normally null)
    matrix_->refresh(this);
    // may factorize, checks if problem finished
    statusOfProblemInParametrics(factorType,data);
    // Say good factorization
    factorType=1;
    if (data.sparseThreshold_) {
      // use default at present
      factorization_->sparseThreshold(0);
      factorization_->goSparse();
    }
    
    // exit if victory declared
    if (problemStatus_>=0)
      break;
    
    // test for maximum iterations
    if (hitMaximumIterations()) {
      problemStatus_=3;
      break;
    }
    // Check event
    {
      int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
      if (status>=0) {
        problemStatus_=5;
        secondaryStatus_=ClpEventHandler::endOfFactorization;
        break;
      }
    }
    // Do iterations
    if (canTryQuick) {
      double * saveDuals=NULL;
      reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals,0);
    } else {
      whileIterating(startingTheta,  endingTheta, reportIncrement,
                     changeLower, changeUpper,
                     changeObjective);
    }
  }
  if (!problemStatus_) {
    theta_=change+startingTheta;
    eventHandler_->event(ClpEventHandler::theta);
    return 0;
  } else if (problemStatus_==10) {
    return -1;
  } else {
    return problemStatus_;
  }
}
/* Checks if finished.  Updates status */
void 
02188 ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
{
  if (type==2) {
    // trouble - go to recovery
    problemStatus_=10;
    return;
  }
  if (problemStatus_>-3||factorization_->pivots()) {
    // factorize
    // later on we will need to recover from singularities
    // also we could skip if first time
    if (type) {
      // is factorization okay?
      if (internalFactorize(1)) {
        // trouble - go to recovery
        problemStatus_=10;
        return;
      }
    }
    if (problemStatus_!=-4||factorization_->pivots()>10)
      problemStatus_=-3;
  }
  // at this stage status is -3 or -4 if looks infeasible
  // get primal and dual solutions
  gutsOfSolution(NULL,NULL);
  double realDualInfeasibilities=sumDualInfeasibilities_;
  // If bad accuracy treat as singular
  if ((largestPrimalError_>1.0e15||largestDualError_>1.0e15)&&numberIterations_) {
    // trouble - go to recovery
    problemStatus_=10;
    return;
  } else if (largestPrimalError_<1.0e-7&&largestDualError_<1.0e-7) {
    // Can reduce tolerance
    double newTolerance = CoinMax(0.99*factorization_->pivotTolerance(),saveData.pivotTolerance_);
    factorization_->pivotTolerance(newTolerance);
  } 
  // Check if looping
  int loop;
  if (type!=2) 
    loop = progress_.looping();
  else
    loop=-1;
  if (loop>=0) {
    problemStatus_ = loop; //exit if in loop
    if (!problemStatus_) {
      // declaring victory
      numberPrimalInfeasibilities_ = 0;
      sumPrimalInfeasibilities_ = 0.0;
    } else {
      problemStatus_ = 10; // instead - try other algorithm
    }
    return;
  } else if (loop<-1) {
    // something may have changed
    gutsOfSolution(NULL,NULL);
  }
  progressFlag_ = 0; //reset progress flag
  if (handler_->detail(CLP_SIMPLEX_STATUS,messages_)<100) {
    handler_->message(CLP_SIMPLEX_STATUS,messages_)
      <<numberIterations_<<objectiveValue();
    handler_->printing(sumPrimalInfeasibilities_>0.0)
      <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
    handler_->printing(sumDualInfeasibilities_>0.0)
      <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
    handler_->printing(numberDualInfeasibilitiesWithoutFree_
                       <numberDualInfeasibilities_)
                         <<numberDualInfeasibilitiesWithoutFree_;
    handler_->message()<<CoinMessageEol;
  }
  /* If we are primal feasible and any dual infeasibilities are on
     free variables then it is better to go to primal */
  if (!numberPrimalInfeasibilities_&&!numberDualInfeasibilitiesWithoutFree_&&
      numberDualInfeasibilities_) {
    problemStatus_=10;
    return;
  }
  
  // check optimal
  // give code benefit of doubt
  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
    // say optimal (with these bounds etc)
    numberDualInfeasibilities_ = 0;
    sumDualInfeasibilities_ = 0.0;
    numberPrimalInfeasibilities_ = 0;
    sumPrimalInfeasibilities_ = 0.0;
  }
  if (dualFeasible()||problemStatus_==-4) {
    progress_.modifyObjective(objectiveValue_
                         -sumDualInfeasibilities_*dualBound_);
  }
  if (numberPrimalInfeasibilities_) {
    if (problemStatus_==-4||problemStatus_==-5) {
      problemStatus_=1; // infeasible
    }
  } else if (numberDualInfeasibilities_) {
    // clean up
    problemStatus_=10;
  } else {
    problemStatus_=0;
  }
  lastGoodIteration_ = numberIterations_;
  if (problemStatus_<0) {
    sumDualInfeasibilities_=realDualInfeasibilities; // back to say be careful
    if (sumDualInfeasibilities_)
      numberDualInfeasibilities_=1;
  }
  // Allow matrices to be sorted etc
  int fake=-999; // signal sort
  matrix_->correctSequence(this,fake,fake);
}
/* This has the flow between re-factorizations
   Reasons to come out:
   -1 iterations etc
   -2 inaccuracy 
   -3 slight inaccuracy (and done iterations)
   +0 looks optimal (might be unbounded - but we will investigate)
   +1 looks infeasible
   +3 max iterations 
*/
int 
02309 ClpSimplexOther::whileIterating(double startingTheta, double & endingTheta,double /*reportIncrement*/,
                                const double * changeLower, const double * changeUpper,
                                const double * changeObjective)
{
  {
    int i;
    for (i=0;i<4;i++) {
      rowArray_[i]->clear();
    }    
    for (i=0;i<2;i++) {
      columnArray_[i]->clear();
    }    
  }      
  // if can't trust much and long way from optimal then relax
  if (largestPrimalError_>10.0)
    factorization_->relaxAccuracyCheck(CoinMin(1.0e2,largestPrimalError_/10.0));
  else
    factorization_->relaxAccuracyCheck(1.0);
  // status stays at -1 while iterating, >=0 finished, -2 to invert
  // status -3 to go to top without an invert
  int returnCode = -1;
  double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
  double useTheta = startingTheta;
  double * primalChange = new double[numberRows_];
  double * dualChange = new double[numberColumns_];
  int numberTotal = numberColumns_+numberRows_;
  int iSequence;
  // See if bounds
  int type=0;
  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    if (changeLower[iSequence]||changeUpper[iSequence]) {
      type=1;
      break;
    }
  }
  // See if objective
  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    if (changeObjective[iSequence]) {
      type |= 2;
      break;
    }
  }
  assert (type);
  while (problemStatus_==-1) {
    double increaseTheta = CoinMin(endingTheta-useTheta,1.0e50);
    
    // Get theta for bounds - we know can't crossover
    int pivotType = nextTheta(type,increaseTheta,primalChange,dualChange,
                              changeLower,changeUpper,changeObjective);
    if (pivotType)
      abort();
    // choose row to go out
    // dualRow will go to virtual row pivot choice algorithm
    reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
    if (pivotRow_>=0) {
      // we found a pivot row
      if (handler_->detail(CLP_SIMPLEX_PIVOTROW,messages_)<100) {
        handler_->message(CLP_SIMPLEX_PIVOTROW,messages_)
          <<pivotRow_
          <<CoinMessageEol;
      }
      // check accuracy of weights
      dualRowPivot_->checkAccuracy();
      // Get good size for pivot
      // Allow first few iterations to take tiny
      double acceptablePivot=1.0e-9;
      if (numberIterations_>100)
        acceptablePivot=1.0e-8;
      if (factorization_->pivots()>10||
        (factorization_->pivots()&&saveSumDual))
      acceptablePivot=1.0e-5; // if we have iterated be more strict
      else if (factorization_->pivots()>5)
      acceptablePivot=1.0e-6; // if we have iterated be slightly more strict
      else if (factorization_->pivots())
        acceptablePivot=1.0e-8; // relax
      double bestPossiblePivot=1.0;
      // get sign for finding row of tableau
      // normal iteration
      // create as packed
      double direction=directionOut_;
      rowArray_[0]->createPacked(1,&pivotRow_,&direction);
      factorization_->updateColumnTranspose(rowArray_[1],rowArray_[0]);
      // put row of tableau in rowArray[0] and columnArray[0]
      matrix_->transposeTimes(this,-1.0,
                        rowArray_[0],rowArray_[3],columnArray_[0]);
      // do ratio test for normal iteration
      bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)->dualColumn(rowArray_[0],
                                                                columnArray_[0],columnArray_[1],
                                                                rowArray_[3],acceptablePivot,NULL);
      if (sequenceIn_>=0) {
      // normal iteration
      // update the incoming column
      double btranAlpha = -alpha_*directionOut_; // for check
      unpackPacked(rowArray_[1]);
      // moved into updateWeights factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
      // and update dual weights (can do in parallel - with extra array)
      alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
                                    rowArray_[2],
                                    rowArray_[3],
                                    rowArray_[1]);
      // see if update stable
#ifdef CLP_DEBUG
      if ((handler_->logLevel()&32))
        printf("btran alpha %g, ftran alpha %g\n",btranAlpha,alpha_);
#endif
      double checkValue=1.0e-7;
      // if can't trust much and long way from optimal then relax
      if (largestPrimalError_>10.0)
        checkValue = CoinMin(1.0e-4,1.0e-8*largestPrimalError_);
      if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
          fabs(btranAlpha-alpha_)>checkValue*(1.0+fabs(alpha_))) {
        handler_->message(CLP_DUAL_CHECK,messages_)
          <<btranAlpha
          <<alpha_
          <<CoinMessageEol;
        if (factorization_->pivots()) {
          dualRowPivot_->unrollWeights();
          problemStatus_=-2; // factorize now
          rowArray_[0]->clear();
          rowArray_[1]->clear();
          columnArray_[0]->clear();
          returnCode=-2;
          break;
        } else {
          // take on more relaxed criterion
            double test;
          if (fabs(btranAlpha)<1.0e-8||fabs(alpha_)<1.0e-8)
              test = 1.0e-1*fabs(alpha_);
            else
              test = 1.0e-4*(1.0+fabs(alpha_));
          if (fabs(btranAlpha)<1.0e-12||fabs(alpha_)<1.0e-12||
            fabs(btranAlpha-alpha_)>test) {
            dualRowPivot_->unrollWeights();
            // need to reject something
            char x = isColumn(sequenceOut_) ? 'C' :'R';
            handler_->message(CLP_SIMPLEX_FLAG,messages_)
            <<x<<sequenceWithin(sequenceOut_)
            <<CoinMessageEol;
            setFlagged(sequenceOut_);
            progress_.clearBadTimes();
            lastBadIteration_ = numberIterations_; // say be more cautious
            rowArray_[0]->clear();
            rowArray_[1]->clear();
            columnArray_[0]->clear();
              if (fabs(alpha_)<1.0e-10&&fabs(btranAlpha)<1.0e-8&&numberIterations_>100) {
                //printf("I think should declare infeasible\n");
                problemStatus_=1;
                returnCode=1;
                break;
              }
            continue;
          }
        }
      }
      // update duals BEFORE replaceColumn so can do updateColumn
      double objectiveChange=0.0;
      // do duals first as variables may flip bounds
      // rowArray_[0] and columnArray_[0] may have flips
      // so use rowArray_[3] for work array from here on
      int nswapped = 0;
      //rowArray_[0]->cleanAndPackSafe(1.0e-60);
      //columnArray_[0]->cleanAndPackSafe(1.0e-60);
        nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],
                                     rowArray_[2],theta_,
                                     objectiveChange,false);

      // which will change basic solution
      if (nswapped) {
        factorization_->updateColumn(rowArray_[3],rowArray_[2]);
        dualRowPivot_->updatePrimalSolution(rowArray_[2],
                                    1.0,objectiveChange);
        // recompute dualOut_
        valueOut_ = solution_[sequenceOut_];
        if (directionOut_<0) {
          dualOut_ = valueOut_ - upperOut_;
        } else {
          dualOut_ = lowerOut_ - valueOut_;
        }
      }
      // amount primal will move
      double movement = -dualOut_*directionOut_/alpha_;
      // so objective should increase by fabs(dj)*movement
      // but we already have objective change - so check will be good
      if (objectiveChange+fabs(movement*dualIn_)<-1.0e-5) {
#ifdef CLP_DEBUG
        if (handler_->logLevel()&32)
          printf("movement %g, swap change %g, rest %g  * %g\n",
               objectiveChange+fabs(movement*dualIn_),
               objectiveChange,movement,dualIn_);
#endif
        if(factorization_->pivots()) {
          // going backwards - factorize
          dualRowPivot_->unrollWeights();
          problemStatus_=-2; // factorize now
          returnCode=-2;
          break;
        }
      }
      CoinAssert(fabs(dualOut_)<1.0e50);
      // if stable replace in basis
      int updateStatus = factorization_->replaceColumn(this,
                                           rowArray_[2],
                                           rowArray_[1],
                                           pivotRow_,
                                           alpha_);
      // if no pivots, bad update but reasonable alpha - take and invert
      if (updateStatus==2&&
               !factorization_->pivots()&&fabs(alpha_)>1.0e-5)
        updateStatus=4;
      if (updateStatus==1||updateStatus==4) {
        // slight error
        if (factorization_->pivots()>5||updateStatus==4) {
          problemStatus_=-2; // factorize now
          returnCode=-3;
        }
      } else if (updateStatus==2) {
        // major error
        dualRowPivot_->unrollWeights();
        // later we may need to unwind more e.g. fake bounds
        if (factorization_->pivots()) {
          problemStatus_=-2; // factorize now
          returnCode=-2;
          break;
        } else {
          // need to reject something
          char x = isColumn(sequenceOut_) ? 'C' :'R';
          handler_->message(CLP_SIMPLEX_FLAG,messages_)
            <<x<<sequenceWithin(sequenceOut_)
            <<CoinMessageEol;
          setFlagged(sequenceOut_);
          progress_.clearBadTimes();
          lastBadIteration_ = numberIterations_; // say be more cautious
          rowArray_[0]->clear();
          rowArray_[1]->clear();
          columnArray_[0]->clear();
          // make sure dual feasible
          // look at all rows and columns
          double objectiveChange=0.0;
          reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
                        0.0,objectiveChange,true);
          continue;
        }
      } else if (updateStatus==3) {
        // out of memory
        // increase space if not many iterations
        if (factorization_->pivots()<
            0.5*factorization_->maximumPivots()&&
            factorization_->pivots()<200)
          factorization_->areaFactor(
                               factorization_->areaFactor() * 1.1);
        problemStatus_=-2; // factorize now
      } else if (updateStatus==5) {
        problemStatus_=-2; // factorize now
      }
      // update primal solution
      if (theta_<0.0) {
#ifdef CLP_DEBUG
        if (handler_->logLevel()&32)
          printf("negative theta %g\n",theta_);
#endif
        theta_=0.0;
      }
      // do actual flips
      reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0],columnArray_[0]);
      //rowArray_[1]->expand();
      dualRowPivot_->updatePrimalSolution(rowArray_[1],
                                  movement,
                                  objectiveChange);
      // modify dualout
      dualOut_ /= alpha_;
      dualOut_ *= -directionOut_;
      //setStatus(sequenceIn_,basic);
      dj_[sequenceIn_]=0.0;
      double oldValue=valueIn_;
      if (directionIn_==-1) {
        // as if from upper bound
        valueIn_ = upperIn_+dualOut_;
      } else {
        // as if from lower bound
        valueIn_ = lowerIn_+dualOut_;
      }
      objectiveChange += cost_[sequenceIn_]*(valueIn_-oldValue);
      // outgoing
      // set dj to zero unless values pass
      if (directionOut_>0) {
        valueOut_ = lowerOut_;
          dj_[sequenceOut_] = theta_;
      } else {
        valueOut_ = upperOut_;
          dj_[sequenceOut_] = -theta_;
      }
      solution_[sequenceOut_]=valueOut_;
      int whatNext=housekeeping(objectiveChange);
      // and set bounds correctly
      reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_); 
      reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
      if (whatNext==1) {
        problemStatus_ =-2; // refactorize
      } else if (whatNext==2) {
        // maximum iterations or equivalent
        problemStatus_= 3;
        returnCode=3;
        break;
      }
      // Check event
      {
        int status = eventHandler_->event(ClpEventHandler::endOfIteration);
        if (status>=0) {
          problemStatus_=5;
          secondaryStatus_=ClpEventHandler::endOfIteration;
          returnCode=4;
          break;
        }
      }
      } else {
      // no incoming column is valid
        pivotRow_=-1;
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32)
        printf("** no column pivot\n");
#endif
      if (factorization_->pivots()<5) {
          // If not in branch and bound etc save ray
          if ((specialOptions_&(1024|4096))==0) {
          // create ray anyway
          delete [] ray_;
          ray_ = new double [ numberRows_];
          rowArray_[0]->expand(); // in case packed
          ClpDisjointCopyN(rowArray_[0]->denseVector(),numberRows_,ray_);
          }
        // If we have just factorized and infeasibility reasonable say infeas
        if (((specialOptions_&4096)!=0||bestPossiblePivot<1.0e-11)&&dualBound_>1.0e8) {
          if (valueOut_>upperOut_+1.0e-3||valueOut_<lowerOut_-1.0e-3
            || (specialOptions_&64)==0) {
            // say infeasible
            problemStatus_=1;
              // unless primal feasible!!!!
              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
              //     numberDualInfeasibilities_,sumDualInfeasibilities_);
              if (numberDualInfeasibilities_)
                problemStatus_=10;
            rowArray_[0]->clear();
            columnArray_[0]->clear();
            returnCode=1;
            break;
          }
        }
        // If special option set - put off as long as possible
        if ((specialOptions_&64)==0) {
          problemStatus_=-4; //say looks infeasible
        } else {
          // flag
          char x = isColumn(sequenceOut_) ? 'C' :'R';
          handler_->message(CLP_SIMPLEX_FLAG,messages_)
            <<x<<sequenceWithin(sequenceOut_)
            <<CoinMessageEol;
          setFlagged(sequenceOut_);
          if (!factorization_->pivots()) {
            rowArray_[0]->clear();
            columnArray_[0]->clear();
            continue;
          }
        }
      }
      rowArray_[0]->clear();
      columnArray_[0]->clear();
      returnCode=1;
      break;
      }
    } else {
      // no pivot row
#ifdef CLP_DEBUG
      if (handler_->logLevel()&32)
      printf("** no row pivot\n");
#endif
      int numberPivots = factorization_->pivots();
      bool specialCase;
      int useNumberFake;
      returnCode=0;
      if (numberPivots<20&&
        (specialOptions_&2048)!=0&&!numberChanged_&&perturbation_>=100
        &&dualBound_>1.0e8) {
      specialCase=true;
      // as dual bound high - should be okay
      useNumberFake=0;
      } else {
      specialCase=false;
      useNumberFake=numberFake_;
      }
      if (!numberPivots||specialCase) {
      // may have crept through - so may be optimal
      // check any flagged variables
      int iRow;
      for (iRow=0;iRow<numberRows_;iRow++) {
        int iPivot=pivotVariable_[iRow];
        if (flagged(iPivot))
          break;
      }
      if (iRow<numberRows_&&numberPivots) {
        // try factorization
        returnCode=-2;
      }
        
      if (useNumberFake||numberDualInfeasibilities_) {
        // may be dual infeasible
        problemStatus_=-5;
      } else {
        if (iRow<numberRows_) {
          problemStatus_=-5;
        } else {
            if (numberPivots) {
              // objective may be wrong
              objectiveValue_ = innerProduct(cost_,
                                                                        numberColumns_+numberRows_,
                                                                        solution_);
              objectiveValue_ += objective_->nonlinearOffset();
              objectiveValue_ /= (objectiveScale_*rhsScale_);
              if ((specialOptions_&16384)==0) {
                // and dual_ may be wrong (i.e. for fixed or basic)
                CoinIndexedVector * arrayVector = rowArray_[1];
                arrayVector->clear();
                int iRow;
                double * array = arrayVector->denseVector();
                /* Use dual_ instead of array
                   Even though dual_ is only numberRows_ long this is
                   okay as gets permuted to longer rowArray_[2]
                */
                arrayVector->setDenseVector(dual_);
                int * index = arrayVector->getIndices();
                int number=0;
                for (iRow=0;iRow<numberRows_;iRow++) {
                  int iPivot=pivotVariable_[iRow];
                  double value = cost_[iPivot];
                  dual_[iRow]=value;
                  if (value) {
                    index[number++]=iRow;
                  }
                }
                arrayVector->setNumElements(number);
                // Extended duals before "updateTranspose"
                matrix_->dualExpanded(this,arrayVector,NULL,0);
                // Btran basic costs
                rowArray_[2]->clear();
                factorization_->updateColumnTranspose(rowArray_[2],arrayVector);
                // and return vector
                arrayVector->setDenseVector(array);
              }
            }
          problemStatus_=0;
          sumPrimalInfeasibilities_=0.0;
            if ((specialOptions_&(1024+16384))!=0) {
              CoinIndexedVector * arrayVector = rowArray_[1];
              arrayVector->clear();
              double * rhs = arrayVector->denseVector();
              times(1.0,solution_,rhs);
              bool bad2=false;
              int i;
              for ( i=0;i<numberRows_;i++) {
                if (rhs[i]<rowLowerWork_[i]-primalTolerance_||
                    rhs[i]>rowUpperWork_[i]+primalTolerance_) {
                  bad2=true;
                } else if (fabs(rhs[i]-rowActivityWork_[i])>1.0e-3) {
                }
                rhs[i]=0.0;
              }
              for ( i=0;i<numberColumns_;i++) {
                if (solution_[i]<columnLowerWork_[i]-primalTolerance_||
                    solution_[i]>columnUpperWork_[i]+primalTolerance_) {
                  bad2=true;
                }
              }
              if (bad2) {
                problemStatus_=-3;
                returnCode=-2;
                // Force to re-factorize early next time
                int numberPivots = factorization_->pivots();
                forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
              }
            }
        }
      }
      } else {
      problemStatus_=-3;
        returnCode=-2;
      // Force to re-factorize early next time
      int numberPivots = factorization_->pivots();
      forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
      }
      break;
    }
  }
  delete [] primalChange;
  delete [] dualChange;
  return returnCode;
}
// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
int 
02806 ClpSimplexOther::nextTheta(int type, double maxTheta, double * primalChange, double * /*dualChange*/,
                           const double * changeLower, const double * changeUpper,
                           const double * /*changeObjective*/)
{
  int numberTotal = numberColumns_+numberRows_;
  int iSequence;
  int iRow;
  theta_=maxTheta;
  bool toLower=false;
  if ((type&1)!=0) {
    // get change 
    for (iSequence=0;iSequence<numberTotal;iSequence++) {
      primalChange[iSequence]=0.0;
      switch(getStatus(iSequence)) {
        
      case basic:
      case isFree:
      case superBasic:
        break;
      case isFixed:
      case atUpperBound:
        primalChange[iSequence]=changeUpper[iSequence];
        break;
      case atLowerBound:
        primalChange[iSequence]=changeLower[iSequence];
        break;
      }
    }
    // use array
    double * array = rowArray_[1]->denseVector();
    times(1.0,primalChange,array);
    int * index = rowArray_[1]->getIndices();
    int number=0;
    for (iRow=0;iRow<numberRows_;iRow++) {
      double value = array[iRow];
      if (value) {
      array[iRow]=value;
      index[number++]=iRow;
      }
    }
    // ftran it
    rowArray_[1]->setNumElements(number);
    factorization_->updateColumn(rowArray_[0],rowArray_[1]);
    number=rowArray_[1]->getNumElements();
    pivotRow_=-1;
    for (iRow=0;iRow<number;iRow++) {
      int iPivot = index[iRow];
      iSequence = pivotVariable_[iPivot];
      // solution value will be sol - theta*alpha
      // bounds will be bounds + change *theta
      double currentSolution = solution_[iSequence];
      double currentLower = lower_[iSequence];
      double currentUpper = upper_[iSequence];
      double alpha = array[iPivot];
      assert (currentSolution>=currentLower-primalTolerance_);
      assert (currentSolution<=currentUpper+primalTolerance_);
      double thetaCoefficient;
      double hitsLower = COIN_DBL_MAX;
      thetaCoefficient = changeLower[iSequence]+alpha;
      if (fabs(thetaCoefficient)>1.0e-8)
        hitsLower = (currentSolution-currentLower)/thetaCoefficient;
      if (hitsLower<0.0) {
        // does not hit - but should we check further
        hitsLower=COIN_DBL_MAX;
      }
      double hitsUpper = COIN_DBL_MAX;
      thetaCoefficient = changeUpper[iSequence]+alpha;
      if (fabs(thetaCoefficient)>1.0e-8)
        hitsUpper = (currentSolution-currentUpper)/thetaCoefficient;
      if (hitsUpper<0.0) {
        // does not hit - but should we check further
        hitsUpper=COIN_DBL_MAX;
      }
      if (CoinMin(hitsLower,hitsUpper)<theta_) {
        theta_ = CoinMin(hitsLower,hitsUpper);
        toLower = hitsLower<hitsUpper;
        pivotRow_=iPivot;
      }
    }
  }
  if ((type&2)!=0) {
    abort();
  }
  if (pivotRow_>=0) {
    sequenceOut_ = pivotVariable_[pivotRow_];
    valueOut_ = solution_[sequenceOut_];
    lowerOut_ = lower_[sequenceOut_];
    upperOut_ = upper_[sequenceOut_];
    if (!toLower) {
      directionOut_ = -1;
      dualOut_ = valueOut_ - upperOut_;
    } else if (valueOut_<lowerOut_) {
      directionOut_ = 1;
      dualOut_ = lowerOut_ - valueOut_;
    }
    return 0;
  } else {
    return -1;
  }
}
/* Expands out all possible combinations for a knapsack
   If buildObj NULL then just computes space needed - returns number elements
   On entry numberOutput is maximum allowed, on exit it is number needed or
   -1 (as will be number elements) if maximum exceeded.  numberOutput will have at
   least space to return values which reconstruct input.
   Rows returned will be original rows but no entries will be returned for
   any rows all of whose entries are in knapsack.  So up to user to allow for this.
   If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
   in expanded knapsack.  Values in buildRow and buildElement;
*/
int 
02917 ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
                        double * buildObj, CoinBigIndex * buildStart,
                        int * buildRow, double * buildElement,int reConstruct) const
{
  int iRow;
  int iColumn;
  // Get column copy
  CoinPackedMatrix * columnCopy = matrix();
  // Get a row copy in standard format
  CoinPackedMatrix matrixByRow;
  matrixByRow.reverseOrderedCopyOf(*columnCopy);
  const double * elementByRow = matrixByRow.getElements();
  const int * column = matrixByRow.getIndices();
  const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
  const int * rowLength = matrixByRow.getVectorLengths();
  CoinBigIndex j;
  int * whichColumn = new int [numberColumns_];
  int * whichRow = new int [numberRows_];
  int numJ=0;
  // Get what other columns can compensate for
  double * lo = new double [numberRows_];
  double * high = new double [numberRows_];
  {
    // Use to get tight column bounds
    ClpSimplex tempModel(*this);
    tempModel.tightenPrimalBounds(0.0,0,true);
    // Now another model without knapsacks
    int nCol=0;
    for (iRow=0;iRow<numberRows_;iRow++) {
      whichRow[iRow]=iRow;
    }
    for (iColumn=0;iColumn<numberColumns_;iColumn++)
      whichColumn[iColumn]=-1;
    for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
      int iColumn = column[j];
      if (columnUpper_[iColumn]>columnLower_[iColumn]) {
      whichColumn[iColumn]=0;
      } else {
      assert (!columnLower_[iColumn]); // fix later
      }
    }
    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
      if (whichColumn[iColumn]<0)
      whichColumn[nCol++]=iColumn;
    }
    ClpSimplex tempModel2(&tempModel,numberRows_,whichRow,nCol,whichColumn,false,false,false);
    // Row copy
    CoinPackedMatrix matrixByRow;
    matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
    const double * elementByRow = matrixByRow.getElements();
    const int * column = matrixByRow.getIndices();
    const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
    const int * rowLength = matrixByRow.getVectorLengths();
    const double * columnLower = tempModel2.getColLower();
    const double * columnUpper = tempModel2.getColUpper();
    for (iRow = 0; iRow < numberRows_; iRow++) {
      lo[iRow]=-COIN_DBL_MAX;
      high[iRow]=COIN_DBL_MAX;
      if (rowLower_[iRow]>-1.0e20||rowUpper_[iRow]<1.0e20) {

      // possible row
      int infiniteUpper = 0;
      int infiniteLower = 0;
      double maximumUp = 0.0;
      double maximumDown = 0.0;
      CoinBigIndex rStart = rowStart[iRow];
      CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
      CoinBigIndex j;
      // Compute possible lower and upper ranges

      for (j = rStart; j < rEnd; ++j) {
        double value=elementByRow[j];
        iColumn = column[j];
        if (value > 0.0) {
          if (columnUpper[iColumn] >= 1.0e20) {
            ++infiniteUpper;
          } else {
            maximumUp += columnUpper[iColumn] * value;
          }
          if (columnLower[iColumn] <= -1.0e20) {
            ++infiniteLower;
          } else {
            maximumDown += columnLower[iColumn] * value;
          }
        } else if (value<0.0) {
          if (columnUpper[iColumn] >= 1.0e20) {
            ++infiniteLower;
          } else {
            maximumDown += columnUpper[iColumn] * value;
          }
          if (columnLower[iColumn] <= -1.0e20) {
            ++infiniteUpper;
          } else {
            maximumUp += columnLower[iColumn] * value;
          }
        }
      }
      // Build in a margin of error
      maximumUp += 1.0e-8*fabs(maximumUp)+1.0e-7;
      maximumDown -= 1.0e-8*fabs(maximumDown)+1.0e-7;
      // we want to save effective rhs
      double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
      double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
      if (up==COIN_DBL_MAX||rowLower_[iRow]==-COIN_DBL_MAX) {
        // However low we go it doesn't matter
        lo[iRow]=-COIN_DBL_MAX;
      } else {
        // If we go below this then can not be feasible
        lo[iRow]=rowLower_[iRow]-up;
      }
      if (down==-COIN_DBL_MAX||rowUpper_[iRow]==COIN_DBL_MAX) {
        // However high we go it doesn't matter
        high[iRow]=COIN_DBL_MAX;
      } else {
        // If we go above this then can not be feasible
        high[iRow]=rowUpper_[iRow]-down;
      }
      }
    }
  }
  numJ =0;
  for (iColumn=0;iColumn<numberColumns_;iColumn++)
    whichColumn[iColumn]=-1;
  int * markRow = new int [numberRows_];
  for (iRow=0;iRow<numberRows_;iRow++) 
    markRow[iRow]=1;
  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
    int iColumn = column[j];
    if (columnUpper_[iColumn]>columnLower_[iColumn]) {
      whichColumn[iColumn]=numJ;
      numJ++;
    }
  }
  /* mark rows
     -n in knapsack and n other variables
     1 no entries
     n+1000 not involved in knapsack but n entries
     0 only in knapsack 
  */
  for (iRow=0;iRow<numberRows_;iRow++) {
    int type=1;
    for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
      int iColumn = column[j];
      if (whichColumn[iColumn]>=0) {
      if (type==1) {
        type=0;
      } else if (type>0) {
        assert (type>1000); 
        type = -(type-1000);
      }
      } else if (type==1) {
      type = 1001;
      } else if (type<0) {
      type --;
      } else if (type==0) {
      type = -1;
      } else {
      assert (type>1000);
      type++;
      }
    }
    markRow[iRow]=type;
    if (type<0&&type>-30&&false)
      printf("markrow on row %d is %d\n",iRow,markRow[iRow]);
  }
  int * bound = new int [numberColumns_+1];
  int * stack = new int [numberColumns_+1];
  int * flip = new int [numberColumns_+1];
  double * offset = new double[numberColumns_+1];
  double * size = new double [numberColumns_+1];
  double * rhsOffset = new double[numberRows_];
  int * build = new int[numberColumns_];
  int maxNumber=numberOutput;
  numJ=0;
  double minSize = rowLower_[knapsackRow];
  double maxSize = rowUpper_[knapsackRow];
  double knapsackOffset=0.0;
  for (j=rowStart[knapsackRow];j<rowStart[knapsackRow]+rowLength[knapsackRow];j++) {
    int iColumn = column[j];
    double lowerColumn=columnLower_[iColumn];
    double upperColumn=columnUpper_[iColumn];
    if (lowerColumn==upperColumn)
      continue;
    double gap = upperColumn-lowerColumn;
    if (gap>1.0e8)
      gap=1.0e8;
    assert (fabs(floor(gap+0.5)-gap)<1.0e-5);
    whichColumn[numJ]=iColumn;
    bound[numJ]=static_cast<int> (gap);
    if (elementByRow[j]>0.0) {
      flip[numJ]=1;
      offset[numJ]=lowerColumn;
      size[numJ++]=elementByRow[j];
    } else {
      flip[numJ]=-1;
      offset[numJ]=upperColumn;
      size[numJ++]=-elementByRow[j];
      lowerColumn = upperColumn;
    }
    knapsackOffset += elementByRow[j]*lowerColumn;
  }  
  int jRow;
  for (iRow=0;iRow<numberRows_;iRow++)
    whichRow[iRow]=iRow;
  ClpSimplex smallModel(this,numberRows_,whichRow,numJ,whichColumn,true,true,true);
  // modify rhs to allow for nonzero lower bounds
  //double * rowLower = smallModel.rowLower();
  //double * rowUpper = smallModel.rowUpper();
  //const double * columnLower = smallModel.columnLower();
  //const double * columnUpper = smallModel.columnUpper();
  const CoinPackedMatrix * matrix = smallModel.matrix();
  const double * element = matrix->getElements();
  const int * row = matrix->getIndices();
  const CoinBigIndex * columnStart = matrix->getVectorStarts();
  const int * columnLength = matrix->getVectorLengths();
  const double * objective = smallModel.objective();
  //double objectiveOffset=0.0;
  // would use for fixed?
  CoinZeroN(rhsOffset,numberRows_);
  double * rowActivity = smallModel.primalRowSolution();
  CoinZeroN(rowActivity,numberRows_);
  maxSize -= knapsackOffset;
  minSize -= knapsackOffset;
  // now generate
  int i;
  int iStack=numJ;
  for (i=0;i<numJ;i++) {
    stack[i]=0;
  }
  double tooMuch = 10.0*maxSize+10000;
  stack[numJ]=1;
  size[numJ]=tooMuch;
  bound[numJ]=0;
  double sum=tooMuch;
  // allow for all zero being OK
  stack[numJ-1]=-1;
  sum -= size[numJ-1];
  numberOutput=0;
  int nelCreate=0;
  /* typeRun is - 0 for initial sizes
                  1 for build
              2 for reconstruct
  */
  int typeRun = buildObj ? 1 : 0;
  if (reConstruct>=0) {
    assert (buildRow&&buildElement);
    typeRun=2;
  }
  if (typeRun==1)
    buildStart[0]=0;
  while (iStack>=0) {
    if (sum>=minSize&&sum<=maxSize) {
      double checkSize =0.0;
      bool good=true;
      int nRow=0;
      double obj=0.0;
      CoinZeroN(rowActivity,numberRows_);
      for (iColumn=0;iColumn<numJ;iColumn++) {
      int iValue = stack[iColumn];
      if (iValue>bound[iColumn]) {
        good=false;
        break;
      } else {
        double realValue = offset[iColumn] + flip[iColumn]*iValue;
        if (realValue) {
          obj += objective[iColumn]*realValue;
          for (CoinBigIndex j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            double value = element[j]*realValue;
            int kRow = row[j];
            if (rowActivity[kRow]) {
            rowActivity[kRow] += value;
            if (!rowActivity[kRow])
              rowActivity[kRow]=1.0e-100;
            } else {
            build[nRow++]=kRow;
            rowActivity[kRow]=value;
            }
          }
        }
      }
      }
      if (good) {
      for (jRow=0;jRow<nRow;jRow++) {
        int kRow=build[jRow];
        double value = rowActivity[kRow];
        if (value>high[kRow]||value<lo[kRow]) {
          good=false;
          break;
        }
      }
      }
      if (good) {
      if (typeRun==1) {
        buildObj[numberOutput]=obj;
        for (jRow=0;jRow<nRow;jRow++) {
          int kRow=build[jRow];
          double value = rowActivity[kRow];
          if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
            buildElement[nelCreate]=value;
            buildRow[nelCreate++]=kRow;
          }
        }
        buildStart[numberOutput+1]=nelCreate;
      } else if (!typeRun) {
        for (jRow=0;jRow<nRow;jRow++) {
          int kRow=build[jRow];
          double value = rowActivity[kRow];
          if (markRow[kRow]<0&&fabs(value)>1.0e-13) {
            nelCreate++;
          }
        }
      }
      if (typeRun==2&&reConstruct==numberOutput) {
        // build and exit
        nelCreate=0;
        for (iColumn=0;iColumn<numJ;iColumn++) {
          int iValue = stack[iColumn];
          double realValue = offset[iColumn] + flip[iColumn]*iValue;
          if (realValue) {
            buildRow[nelCreate]=whichColumn[iColumn];
            buildElement[nelCreate++]=realValue;
          }
        }
        numberOutput=1;
        for (i=0;i<numJ;i++) {
          bound[i]=0;
        }
        break;
      }
      numberOutput++;
      if (numberOutput>maxNumber) {
        nelCreate=-numberOutput;
        numberOutput=-1;
        for (i=0;i<numJ;i++) {
          bound[i]=0;
        }
        break;
      } else if (typeRun==1&&numberOutput==maxNumber) {
        // On second run 
        for (i=0;i<numJ;i++) {
          bound[i]=0;
        }
        break;
      } 
      for (int j=0;j<numJ;j++) {
        checkSize += stack[j]*size[j];
      }
      assert (fabs(sum-checkSize)<1.0e-3);
      }
      for (jRow=0;jRow<nRow;jRow++) {
      int kRow=build[jRow];
      rowActivity[kRow]=0.0;
      }
    }
    if (sum>maxSize||stack[iStack]>bound[iStack]) {
      sum -= size[iStack]*stack[iStack]; 
      stack[iStack--]=0;
      if (iStack>=0) {
      stack[iStack] ++;
      sum += size[iStack];
      }
    } else {
      // must be less
      // add to last possible
      iStack = numJ-1;
      sum += size[iStack];
      stack[iStack]++;
    }
  }
  //printf("%d will be created\n",numberOutput);
  delete [] whichColumn;
  delete [] whichRow;
  delete [] bound;
  delete [] stack;
  delete [] flip;
  delete [] size;
  delete [] offset;
  delete [] rhsOffset;
  delete [] build;
  delete [] markRow;
  delete [] lo;
  delete [] high;
  return nelCreate;
}
// Quick try at cleaning up duals if postsolve gets wrong
void 
03304 ClpSimplexOther::cleanupAfterPostsolve()
{
  // First mark singleton equality rows
  char * mark = new char [ numberRows_];
  memset(mark,0,numberRows_);
  const int * row = matrix_->getIndices();
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * columnLength = matrix_->getVectorLengths(); 
  const double * element = matrix_->getElements();
  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    for (CoinBigIndex j =columnStart[iColumn];
       j<columnStart[iColumn]+columnLength[iColumn];j++) {
      int iRow = row[j];
      if (mark[iRow])
      mark[iRow]=2;
      else
      mark[iRow]=1;
    }
  }
  // for now just == rows
  for (int iRow=0;iRow<numberRows_;iRow++) {
    if (rowUpper_[iRow]>rowLower_[iRow])
      mark[iRow]=3;
  }
  double dualTolerance=dblParam_[ClpDualTolerance];
  double primalTolerance=dblParam_[ClpPrimalTolerance];
  int numberCleaned=0;
  double maxmin = optimizationDirection_;
  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    double dualValue = reducedCost_[iColumn]*maxmin;
    double primalValue = columnActivity_[iColumn];
    double lower = columnLower_[iColumn];
    double upper = columnUpper_[iColumn];
    int way=0;
    switch(getColumnStatus(iColumn)) {
      
    case basic:
      // dual should be zero
      if (dualValue>dualTolerance) {
      way=-1;
      } else if (dualValue<-dualTolerance) {
      way=1;
      }
      break;
    case ClpSimplex::isFixed:
      break;
    case atUpperBound:
      // dual should not be positive
      if (dualValue>dualTolerance) {
      way=-1;
      }
      break;
    case atLowerBound:
      // dual should not be negative
      if (dualValue<-dualTolerance) {
      way=1;
      }
      break;
    case superBasic:
    case isFree:
      if (primalValue<upper-primalTolerance) {
      // dual should not be negative
      if (dualValue<-dualTolerance) {
        way=1;
      }
      }
      if (primalValue>lower+primalTolerance) {
      // dual should not be positive
      if (dualValue>dualTolerance) {
        way=-1;
      }
      }
      break;
    }
    if (way) {
      // see if can find singleton row
      for (CoinBigIndex j =columnStart[iColumn];
         j<columnStart[iColumn]+columnLength[iColumn];j++) {
      int iRow = row[j];
      if (mark[iRow]==1) {
        double value=element[j];
        // dj - addDual*value == 0.0
        double addDual = dualValue/value;
        dual_[iRow] += addDual;
        reducedCost_[iColumn]=0.0;
        numberCleaned++;
        break;
      }
      }
    }
  }
  delete [] mark;
#ifdef CLP_INVESTIGATE
  printf("cleanupAfterPostsolve cleaned up %d columns\n",numberCleaned);
#endif
  // Redo
  memcpy(reducedCost_,this->objective(),numberColumns_*sizeof(double));
  matrix_->transposeTimes(-1.0,dual_,reducedCost_);
  checkSolutionInternal();
}

Generated by  Doxygen 1.6.0   Back to index