// 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.
   
*/
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];
  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];
    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 (valueIncrease) {
	  if (sequenceIncrease>=0)
	    valueIncrease[i] = primalRanging1(sequenceIncrease,iSequence);
	  if (sequenceDecrease>=0)
	    valueDecrease[i] = primalRanging1(sequenceDecrease,iSequence);
	}
        if (inCBC) { 
          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_/columnScale_[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_/columnScale_[sequenceDecrease];
                costDecrease = CoinMax(fabs(djValue*movement),costDecrease);
              }
            } else {
              costDecrease=0.0;
            }
          }
        }
      }
      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 (!auxiliaryModel_) {
      if (rowScale_) {
        if (iSequence<numberColumns_) 
          scaleFactor = 1.0/(objectiveScale_*columnScale_[iSequence]);
        else
          scaleFactor = rowScale_[iSequence-numberColumns_]/objectiveScale_;
      } else {
        scaleFactor = 1.0/objectiveScale_;
      }
    } else {
      if (auxiliaryModel_->rowScale()) {
        if (iSequence<numberColumns_) 
          scaleFactor = 1.0/(objectiveScale_*auxiliaryModel_->columnScale()[iSequence]);
        else
          scaleFactor = auxiliaryModel_->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
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 
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 
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 
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 
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,"IEEE");
  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 
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 * 
ClpSimplexOther::dualOfModel(double fractionRowRanges,double fractionColumnRanges) const
{
  const ClpSimplex * model2 = (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;
  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[kRow]=-rowLower[iRow];
        fromRowsLower[kRow]=0.0;
        fromRowsUpper[kRow]=COIN_DBL_MAX;
        which[kRow]=iRow;
        kRow++;
      }
    }
  }
  if (numberExtraRows) {
    CoinPackedMatrix newCopy;
    newCopy.setExtraGap(0.0);
    newCopy.setExtraMajor(0.0);
    newCopy.submatrixOfWithDuplicates(rowCopy,kRow,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
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
  const double * primalDual = this->dualRowSolution();
  const double * primalDj = this->dualColumnSolution();
  const double * primalSol = this->primalColumnSolution();
  const double * primalActs = this->primalRowSolution();
  char ss[]={'F','B','U','L','S','F'};
  dual(); // for testing
  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",this->numberRows());
  for (iRow=0;iRow<this->numberRows();iRow++)
    printf("%d at %c primal %g dual %g\n",
           iRow,ss[this->getRowStatus(iRow)],
           primalActs[iRow],primalDual[iRow]);
  printf ("Primal problem column info %d columns\n",this->numberColumns());
  for (iColumn=0;iColumn<this->numberColumns();iColumn++)
    printf("%d at %c primal %g dual %g\n",
           iColumn,ss[this->getColumnStatus(iColumn)],
           primalSol[iColumn],primalDj[iColumn]);
#endif
  // position at bound information
  int jColumn=numberRows_+numberExtraRows;
  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) {
          setColumnStatus(iColumn,atLowerBound);
          columnActivity_[iColumn]=columnLower_[iColumn];
        } 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) {
          setColumnStatus(iColumn,atLowerBound);
          columnActivity_[iColumn]=columnLower_[iColumn];
        } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
          setColumnStatus(iColumn,atUpperBound);
          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];
	}
        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++;
          columnActivity_[iColumn]=-dualDual[iColumn];
          reducedCost_[iColumn]=0.0;
        } else {
          reducedCost_[iColumn]=objValue-dualActs[iColumn];
          if (fabs(otherValue-columnLower_[iColumn])<1.0e-5) {
            setColumnStatus(iColumn,atLowerBound);
            columnActivity_[iColumn]=columnLower_[iColumn];
          } else if (fabs(otherValue-columnUpper_[iColumn])<1.0e-5) {
            setColumnStatus(iColumn,atUpperBound);
            columnActivity_[iColumn]=columnUpper_[iColumn];
          } else {
            abort();
          }
        }
      }
    }
  }
  // now rows
  int kRow=0;
  int numberRanges=0;
  for (iRow=0;iRow<numberRows_;iRow++) {
    Status status = dualProblem->getColumnStatus(kRow);
    if (status==basic) {
      // row is at bound
      dual_[iRow]=dualSol[kRow];;
    } 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 {
        rowActivity_[iRow]=rowUpper_[iRow]+dualSol[kRow];
      }        
      kRow++;
    } else if (rowUpper_[iRow]>1.0e20) {
      if (status==basic) {
        rowActivity_[iRow]=rowLower_[iRow];
        setRowStatus(iRow,atLowerBound);
      } else {
        rowActivity_[iRow]=rowLower_[iRow]+dualSol[kRow];
      }        
      kRow++;
    } else {
      if (rowUpper_[iRow]==rowLower_[iRow]) {
        rowActivity_[iRow]=rowLower_[iRow];
        if (status==basic) {
          setRowStatus(iRow,atLowerBound);
        }        
        kRow++;
      } else {
        // range
        numberRanges++;
	Status statusL = dualProblem->getColumnStatus(kRow+1);
	if (status==basic) {
	  assert (statusL!=basic);
	  rowActivity_[iRow]=rowUpper_[iRow];
	  setRowStatus(iRow,atUpperBound);
	} else if (statusL==basic) {
	  rowActivity_[iRow]=rowLower_[iRow];
	  setRowStatus(iRow,atLowerBound);
	} else {
	  rowActivity_[iRow]=rowLower_[iRow]+dualSol[kRow];
	  // row basic
	  setRowStatus(iRow,basic);
	  numberBasic++;
	  dual_[iRow]=0.0;
	}
        kRow++;
        kRow++;
      }
    }
  }
  if (numberRanges) {
    printf("%d ranges - coding needed\n",numberRanges);
    returnCode=1;
  }
  if (numberBasic!=numberRows_) {
    printf("Bad basis - ranges?\n");
    assert (numberRanges);
  }
  if (optimizationDirection_<0.0) {
    for (iRow=0;iRow<numberRows_;iRow++) {
      dual_[iRow]=-dual_[iRow];
    }
    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
      reducedCost_[iColumn]=-reducedCost_[iColumn];
    }
  }
  // redo row activities
  memset(rowActivity_,0,numberRows_*sizeof(double));
  matrix_->times(-1.0,columnActivity_,rowActivity_);
  checkSolutionInternal();
  return 1; //temp
  return returnCode;
}
/* Does very cursory presolve.
   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
*/
ClpSimplex * 
ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
                        int & nBound, bool moreBounds, bool tightenBounds)
{
#if 0
  //#ifndef NDEBUG
  {
    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);
    // 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
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  return small;
}
/* After very cursory presolve.
   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
*/
void 
ClpSimplexOther::afterCrunch(const ClpSimplex & small,
                             const int * whichRow, 
                             const int * whichColumn, int nBound)
{
  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 
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 
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 = ((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;
      ((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();
        }
      }
    }
    ((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_;
          memcpy(status_,copyModel.statusArray(),numberTotal);
          memcpy(columnActivity_,copyModel.primalColumnSolution(),numberColumns_*sizeof(double));
          memcpy(rowActivity_,copyModel.primalRowSolution(),numberRows_*sizeof(double));
          cleanedUp=0;
        }
      }
    }
    delete [] chgLower;
    delete [] chgUpper;
    delete [] chgObjective;
  }
  perturbation_ = savePerturbation;
  return problemStatus_;
}
int 
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;
      ((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 
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 
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
    ((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 = ((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]);
	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_[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 = ((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;
	    ((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
	((ClpSimplexDual *) this)->flipBounds(rowArray_[0],columnArray_[0],theta_);
	//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
	((ClpSimplexDual *) this)->originalBound(sequenceIn_); 
	((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 
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;
  }
}