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

ClpDynamicMatrix.cpp

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


#include <cstdio>

#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "ClpDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
00028 ClpDynamicMatrix::ClpDynamicMatrix () 
  : ClpPackedMatrix(),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    savedBestGubDual_(0.0),
    savedBestSet_(0),
    backToPivotRow_(NULL),
    keyVariable_(NULL),
    toIndex_(NULL),
    fromIndex_(NULL), 
    numberSets_(0),
    numberActiveSets_(0),
    objectiveOffset_(0.0),
    lowerSet_(NULL),
    upperSet_(NULL),
    status_(NULL),
    model_(NULL),
    firstAvailable_(0),
    firstAvailableBefore_(0),
    firstDynamic_(0),
    lastDynamic_(0),
    numberStaticRows_(0),
    numberElements_(0),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    noCheck_(-1),
    infeasibilityWeight_(0.0),
    numberGubColumns_(0),
    maximumGubColumns_(0),
    maximumElements_(0),
    startSet_(NULL),
    next_(NULL),
    startColumn_(NULL),
    row_(NULL),
    element_(NULL),
    cost_(NULL),
    id_(NULL),
    dynamicStatus_(NULL),
    columnLower_(NULL),
    columnUpper_(NULL)
{
  setType(15);
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
00077 ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs) 
  : ClpPackedMatrix(rhs)
{  
  objectiveOffset_ = rhs.objectiveOffset_;
  numberSets_ = rhs.numberSets_;
  numberActiveSets_ = rhs.numberActiveSets_;
  firstAvailable_ = rhs.firstAvailable_;
  firstAvailableBefore_ = rhs.firstAvailableBefore_;
  firstDynamic_ = rhs.firstDynamic_;
  lastDynamic_ = rhs.lastDynamic_;
  numberStaticRows_ = rhs.numberStaticRows_;
  numberElements_ = rhs.numberElements_;
  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_);
  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
  toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_);
  lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
  upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
  status_ = ClpCopyOfArray(rhs.status_,numberSets_);
  model_ = rhs.model_;
  sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
  savedBestGubDual_=rhs.savedBestGubDual_;
  savedBestSet_=rhs.savedBestSet_;
  noCheck_ = rhs.noCheck_;
  infeasibilityWeight_=rhs.infeasibilityWeight_;
  // Now secondary data
  numberGubColumns_ = rhs.numberGubColumns_;
  maximumGubColumns_=rhs.maximumGubColumns_;
  maximumElements_ = rhs.maximumElements_;
  startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_);
  next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_);
  startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1);
  row_ = ClpCopyOfArray(rhs.row_,maximumElements_);
  element_ = ClpCopyOfArray(rhs.element_,maximumElements_);
  cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_);
  id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
  columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_);
  columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_);
  dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_);
}

/* This is the real constructor*/
00124 ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * model, int numberSets,
                           int numberGubColumns, const int * starts,
                           const double * lower, const double * upper,
                           const CoinBigIndex * startColumn, const int * row,
                           const double * element, const double * cost,
                           const double * columnLower, const double * columnUpper,
                           const unsigned char * status,
                           const unsigned char * dynamicStatus)
  : ClpPackedMatrix()
{
  setType(15);
  objectiveOffset_ = model->objectiveOffset();
  model_ = model;
  numberSets_ = numberSets;
  numberGubColumns_ = numberGubColumns;
  maximumGubColumns_=numberGubColumns_;
  if (numberGubColumns_)
    maximumElements_ = startColumn[numberGubColumns_];
  else 
    maximumElements_ = 0;
  startSet_ = new int [numberSets_];
  next_ = new int [maximumGubColumns_];
  // fill in startSet and next
  int iSet;
  if (numberGubColumns_) {
    for (iSet=0;iSet<numberSets_;iSet++) {
      int first = starts[iSet];
      int last = starts[iSet+1]-1;
      startSet_[iSet]=first;
      for (int i=first;i<last;i++)
      next_[i]=i+1;
      next_[last]=-iSet-1;
    }
  }
  int numberColumns = model->numberColumns();
  int numberRows = model->numberRows();
  numberStaticRows_ = numberRows;
  savedBestGubDual_=0.0;
  savedBestSet_=0;
  // Number of columns needed
  int frequency = model->factorizationFrequency();
  int numberGubInSmall = numberRows + frequency + CoinMin(frequency,numberSets_)+4;
  // for small problems this could be too big
  //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
  int numberNeeded = numberGubInSmall + numberColumns;
  firstAvailable_ = numberColumns;
  firstAvailableBefore_ = firstAvailable_;
  firstDynamic_ = numberColumns;
  lastDynamic_ = numberNeeded;
  startColumn_ = ClpCopyOfArray(startColumn,numberGubColumns_+1);
  if (!numberGubColumns_) {
    startColumn_ = new CoinBigIndex [1];
    startColumn_[0]=0;
  }
  CoinBigIndex numberElements = startColumn_[numberGubColumns_];
  row_ = ClpCopyOfArray(row,numberElements);
  element_ = new double[numberElements];
  CoinBigIndex i;
  for (i=0;i<numberElements;i++)
    element_[i]=element[i];
  cost_ = new double[numberGubColumns_];
  for (i=0;i<numberGubColumns_;i++) {
    cost_[i]=cost[i];
    // I don't think I need sorted but ...
    CoinSort_2(row_+startColumn_[i],row_+startColumn_[i+1],element_+startColumn_[i]);
  }
  if (columnLower) {
    columnLower_ = new double[numberGubColumns_];
    for (i=0;i<numberGubColumns_;i++) 
      columnLower_[i]=columnLower[i];
  } else {
    columnLower_=NULL;
  }
  if (columnUpper) {
    columnUpper_ = new double[numberGubColumns_];
    for (i=0;i<numberGubColumns_;i++) 
      columnUpper_[i]=columnUpper[i];
  } else {
    columnUpper_=NULL;
  }
  lowerSet_ = new double[numberSets_];
  for (i=0;i<numberSets_;i++) {
    if (lower[i]>-1.0e20)
      lowerSet_[i]=lower[i];
    else
      lowerSet_[i] = -1.0e30;
  }
  upperSet_ = new double[numberSets_];
  for (i=0;i<numberSets_;i++) {
    if (upper[i]<1.0e20)
      upperSet_[i]=upper[i];
    else
      upperSet_[i] = 1.0e30;
  }
  id_ = new int[numberGubInSmall]; 
  for (i=0;i<numberGubInSmall;i++)
    id_[i]=-1;
  ClpPackedMatrix* originalMatrixA =
    dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
  assert (originalMatrixA);
  CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
  originalMatrixA->setMatrixNull(); // so can be deleted safely
  // guess how much space needed
  double guess = numberElements;
  guess /= static_cast<double> (numberColumns);
  guess *= 2*numberGubInSmall;
  numberElements_ = static_cast<int> (guess);
  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
  matrix_ = originalMatrix;
  flags_ &= ~1;
  // resize model (matrix stays same)
  int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(frequency,numberRows))+1;
  model->resize(newRowSize,numberNeeded);
  for (i=numberRows;i<newRowSize;i++)
    model->setRowStatus(i,ClpSimplex::basic);
  if (columnUpper_) {
    // set all upper bounds so we have enough space
    double * columnUpper = model->columnUpper();
    for(i=firstDynamic_;i<lastDynamic_;i++)
      columnUpper[i]=1.0e10;
  }
  // resize matrix
  // extra 1 is so can keep number of elements handy
  originalMatrix->reserve(numberNeeded,numberElements_,true);
  originalMatrix->reserve(numberNeeded+1,numberElements_,false);
  originalMatrix->getMutableVectorStarts()[numberColumns]=originalMatrix->getNumElements();
  originalMatrix->setDimensions(newRowSize,-1);
  numberActiveColumns_=firstDynamic_;
  // redo number of columns
  numberColumns = matrix_->getNumCols();
  backToPivotRow_ = new int[numberNeeded];
  keyVariable_ = new int[numberSets_];
  if (status) {
    status_ = ClpCopyOfArray(status,numberSets_);
    assert (dynamicStatus);
    dynamicStatus_ = ClpCopyOfArray(dynamicStatus,numberGubColumns_);
  } else {
    status_= new unsigned char [numberSets_];
    memset(status_,0,numberSets_);
    int i;
    for (i=0;i<numberSets_;i++) {
      // make slack key
      setStatus(i,ClpSimplex::basic);
    }
    dynamicStatus_ = new unsigned char [numberGubColumns_];
    memset(dynamicStatus_,0,numberGubColumns_); // for clarity
    for (i=0;i<numberGubColumns_;i++)
      setDynamicStatus(i,atLowerBound);
  }
  toIndex_ = new int[numberSets_];
  for (iSet=0;iSet<numberSets_;iSet++) 
    toIndex_[iSet]=-1;
  fromIndex_ = new int [newRowSize-numberStaticRows_+1];
  numberActiveSets_=0;
  rhsOffset_=NULL;
  if (numberGubColumns_) {
    if (!status) {
      gubCrash();
    } else {
      initialProblem();
    }
  }
  noCheck_ = -1;
  infeasibilityWeight_=0.0;
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
00293 ClpDynamicMatrix::~ClpDynamicMatrix ()
{
  delete [] backToPivotRow_;
  delete [] keyVariable_;
  delete [] toIndex_;
  delete [] fromIndex_;
  delete [] lowerSet_;
  delete [] upperSet_;
  delete [] status_;
  delete [] startSet_;
  delete [] next_;
  delete [] startColumn_;
  delete [] row_;
  delete [] element_;
  delete [] cost_;
  delete [] id_;
  delete [] dynamicStatus_;
  delete [] columnLower_;
  delete [] columnUpper_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpDynamicMatrix &
ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs)
{
  if (this != &rhs) {
    ClpPackedMatrix::operator=(rhs);
    delete [] backToPivotRow_;
    delete [] keyVariable_;
    delete [] toIndex_;
    delete [] fromIndex_;
    delete [] lowerSet_;
    delete [] upperSet_;
    delete [] status_;
    delete [] startSet_;
    delete [] next_;
    delete [] startColumn_;
    delete [] row_;
    delete [] element_;
    delete [] cost_;
    delete [] id_;
    delete [] dynamicStatus_;
    delete [] columnLower_;
    delete [] columnUpper_;
    objectiveOffset_ = rhs.objectiveOffset_;
    numberSets_ = rhs.numberSets_;
    numberActiveSets_ = rhs.numberActiveSets_;
    firstAvailable_ = rhs.firstAvailable_;
    firstAvailableBefore_ = rhs.firstAvailableBefore_;
    firstDynamic_ = rhs.firstDynamic_;
    lastDynamic_ = rhs.lastDynamic_;
    numberStaticRows_ = rhs.numberStaticRows_;
    numberElements_ = rhs.numberElements_;
    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_);
    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
    toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_);
    lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
    upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
    status_ = ClpCopyOfArray(rhs.status_,numberSets_);
    model_ = rhs.model_;
    sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    savedBestGubDual_=rhs.savedBestGubDual_;
    savedBestSet_=rhs.savedBestSet_;
    noCheck_ = rhs.noCheck_;
    infeasibilityWeight_=rhs.infeasibilityWeight_;
    // Now secondary data
    numberGubColumns_ = rhs.numberGubColumns_;
    maximumGubColumns_=rhs.maximumGubColumns_;
    maximumElements_ = rhs.maximumElements_;
    startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_);
    next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_);
    startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1);
    row_ = ClpCopyOfArray(rhs.row_,maximumElements_);
    element_ = ClpCopyOfArray(rhs.element_,maximumElements_);
    cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_);
    id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
    columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_);
    columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_);
    dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_);
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
00386 ClpMatrixBase * ClpDynamicMatrix::clone() const
{
  return new ClpDynamicMatrix(*this);
}
// Partial pricing 
void 
00392 ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
                        int & bestSequence, int & numberWanted)
{
  numberWanted=currentWanted_;
  assert(!model->rowScale());
  if (numberSets_) {
    // Do packed part before gub
    // always???
    //printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
  } else {
    // no gub
    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
    return;
  }
  if (numberWanted>0) {
    // and do some proportion of full set
    int startG2 = static_cast<int> (startFraction*numberSets_);
    int endG2 = static_cast<int> (endFraction*numberSets_+0.1);
    endG2 = CoinMin(endG2,numberSets_);
    //printf("gub price - set start %d end %d\n",
    //   startG2,endG2);
    double tolerance=model->currentDualTolerance();
    double * reducedCost = model->djRegion();
    const double * duals = model->dualRowSolution();
    double bestDj;
    int numberRows = model->numberRows();
    int slackOffset = lastDynamic_+numberRows;
    int structuralOffset = slackOffset+numberSets_;
    // If nothing found yet can go all the way to end
    int endAll = endG2;
    if (bestSequence<0&&!startG2)
      endAll = numberSets_;
    if (bestSequence>=0) {
      if (bestSequence!=savedBestSequence_)
      bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
      else
      bestDj = savedBestDj_;
    } else {
      bestDj=tolerance;
    }
    int saveSequence = bestSequence;
    double djMod=0.0;
    double bestDjMod=0.0;
    //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
    //     startG2,endG2,numberWanted);
    int bestSet=-1;
#if 0    
    // make sure first available is clean (in case last iteration rejected)
    cost[firstAvailable_]=0.0;
    length[firstAvailable_]=0;
    model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
    model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
    {
      for (int i=firstAvailable_;i<lastDynamic_;i++)
      assert(!cost[i]);
    }
#endif
    int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 
    int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
    for (int iSet = startG2;iSet<endAll;iSet++) {
      if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
      // give up
      numberWanted=0;
      break;
      } else if (iSet==endG2&&bestSequence>=0) {
      break;
      }
      int gubRow = toIndex_[iSet];
      if (gubRow>=0) {
      djMod = duals[gubRow+numberStaticRows_]; // have I got sign right?
      } else {
      int iBasic = keyVariable_[iSet];
      if (iBasic>=maximumGubColumns_) {
        djMod = 0.0; // set not in
      } else {
        // get dj without 
        djMod=0.0;
        for (CoinBigIndex j=startColumn_[iBasic];
             j<startColumn_[iBasic+1];j++) {
          int jRow=row_[j];
          djMod -= duals[jRow]*element_[j];
        }
        djMod += cost_[iBasic];
        // See if gub slack possible - dj is djMod
        if (getStatus(iSet)==ClpSimplex::atLowerBound) {
          double value = -djMod;
          if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
            // check flagged variable and correct dj
            if (!flagged(iSet)) {
              bestDj=value;
              bestSequence = slackOffset+iSet;
              bestDjMod = djMod;
              bestSet = iSet;
            } else {
              // just to make sure we don't exit before got something
              numberWanted++;
              abort();
            }
            }
          }
        } else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
          double value = djMod;
          if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
            // check flagged variable and correct dj
            if (!flagged(iSet)) {
              bestDj=value;
              bestSequence = slackOffset+iSet;
              bestDjMod = djMod;
              bestSet = iSet;
            } else {
              // just to make sure we don't exit before got something
              numberWanted++;
              abort();
            }
            }
          }
        }
      }
      }
      int iSequence= startSet_[iSet];
      while (iSequence>=0) {
      DynamicStatus status = getDynamicStatus(iSequence);
      if (status==atLowerBound||status==atUpperBound) {
        double value=cost_[iSequence]-djMod;
        for (CoinBigIndex j=startColumn_[iSequence];
             j<startColumn_[iSequence+1];j++) {
          int jRow=row_[j];
          value -= duals[jRow]*element_[j];
        }
        // change sign if at lower bound
        if (status==atLowerBound)
          value = -value;
        if (value>tolerance) {
          numberWanted--;
          if (value>bestDj) {
            // check flagged variable and correct dj
            if (!flagged(iSequence)) {
            bestDj=value;
            bestSequence = structuralOffset+iSequence;
            bestDjMod = djMod;
            bestSet = iSet;
            } else {
            // just to make sure we don't exit before got something
            numberWanted++;
            }
          }
        }
      }
      iSequence = next_[iSequence]; //onto next in set
      }
      if (numberWanted<=0) {
      numberWanted=0;
      break;
      }
    }
    if (bestSequence!=saveSequence) {
      savedBestGubDual_ = bestDjMod;
      savedBestDj_ = bestDj;
      savedBestSequence_ = bestSequence;
      savedBestSet_ = bestSet;
    }
    // See if may be finished
    if (!startG2&&bestSequence<0)
      infeasibilityWeight_=model_->infeasibilityCost();
    else if (bestSequence>=0)
      infeasibilityWeight_=-1.0;
  }
  currentWanted_=numberWanted;
}
/* Returns effective RHS if it is being used.  This is used for long problems
   or big gub or anywhere where going through full columns is
   expensive.  This may re-compute */
double * 
00570 ClpDynamicMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,
                  bool
#ifdef CLP_DEBUG
 check
#endif
)
{
  //forceRefresh=true;
  if (!model_->numberIterations())
    forceRefresh=true;
  //check=false;
#ifdef CLP_DEBUG
  double * saveE=NULL;
  if (rhsOffset_&&check) {
    int numberRows = model->numberRows();
    saveE = new double[numberRows];
  }
#endif
  if (rhsOffset_) {
#ifdef CLP_DEBUG
    if (check) {
      // no need - but check anyway
      int numberRows = model->numberRows();
      double * rhs = new double[numberRows];
      int iRow;
      int iSet;
      CoinZeroN(rhs,numberRows);
      // do ones at bounds before gub
      const double * smallSolution = model->solutionRegion();
      const double * element =matrix_->getElements();
      const int * row = matrix_->getIndices();
      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
      const int * length = matrix_->getVectorLengths();
      int iColumn;
      double objectiveOffset=0.0;
      for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
      if (model->getStatus(iColumn)!=ClpSimplex::basic) {
        double value = smallSolution[iColumn];
        for (CoinBigIndex j=startColumn[iColumn];
             j<startColumn[iColumn]+length[iColumn];j++) {
          int jRow=row[j];
          rhs[jRow] -= value*element[j];
        }
      }
      }
      if (columnLower_||columnUpper_) {
      double * solution = new double [numberGubColumns_];
      for (iSet=0;iSet<numberSets_;iSet++) {
        int j= startSet_[iSet];
        while (j>=0) {
          double value=0.0;
          if (getDynamicStatus(j)!=inSmall) {
            if (getDynamicStatus(j)==atLowerBound) {
            if (columnLower_) 
              value= columnLower_[j];
            } else if (getDynamicStatus(j)==atUpperBound) { 
            value= columnUpper_[j];
            } else if (getDynamicStatus(j)==soloKey) {
            value =keyValue(iSet);
            }
            objectiveOffset += value*cost_[j];
          } 
          solution[j]=value;
          j = next_[j]; //onto next in set
        }
      }
      // ones in gub and in small problem
      for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
        if (model_->getStatus(iColumn)!=ClpSimplex::basic) {
          int jFull = id_[iColumn-firstDynamic_];
          solution[jFull]=smallSolution[iColumn];
        }
      }
      for (iSet=0;iSet<numberSets_;iSet++) {
        int kRow = toIndex_[iSet];
        if (kRow>=0)
          kRow += numberStaticRows_;
        int j= startSet_[iSet];
        while (j>=0) {
          double value = solution[j];
          if (value) {
            for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) {
            int iRow = row_[k];
            rhs[iRow] -= element_[k]*value;
            }
            if (kRow>=0)
            rhs[kRow] -= value;
          }
          j = next_[j]; //onto next in set
        }
      }
      delete [] solution;
      } else {
      // bounds
      ClpSimplex::Status iStatus;
      for (int iSet=0;iSet<numberSets_;iSet++) {
        int kRow = toIndex_[iSet];
        if (kRow<0) {
          int iColumn = keyVariable_[iSet];
          if (iColumn<maximumGubColumns_) {
            // key is not treated as basic
            double b=0.0;
            // key is structural - where is slack
            iStatus = getStatus(iSet);
            assert (iStatus!=ClpSimplex::basic);
            if (iStatus==ClpSimplex::atLowerBound)
            b=lowerSet_[iSet];
            else
            b=upperSet_[iSet];
            if (b) {
            objectiveOffset += b*cost_[iColumn];
            for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
              int iRow = row_[j];
              rhs[iRow] -= element_[j]*b;
            }
            }
          }
        }
      }
      }
      if (fabs(model->objectiveOffset()-objectiveOffset_-objectiveOffset)>1.0e-1)
      printf("old offset %g, true %g\n",model->objectiveOffset(),
             objectiveOffset_-objectiveOffset);
      for (iRow=0;iRow<numberRows;iRow++) {
      if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
        printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
      }
      CoinMemcpyN(rhs,numberRows,saveE);
      delete [] rhs;
    }
#endif
    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
                   lastRefresh_+refreshFrequency_)) {
      int numberRows = model->numberRows();
      int iSet;
      CoinZeroN(rhsOffset_,numberRows);
      // do ones at bounds before gub
      const double * smallSolution = model->solutionRegion();
      const double * element =matrix_->getElements();
      const int * row = matrix_->getIndices();
      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
      const int * length = matrix_->getVectorLengths();
      int iColumn;
      double objectiveOffset=0.0;
      for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
      if (model->getStatus(iColumn)!=ClpSimplex::basic) {
        double value = smallSolution[iColumn];
        for (CoinBigIndex j=startColumn[iColumn];
             j<startColumn[iColumn]+length[iColumn];j++) {
          int jRow=row[j];
          rhsOffset_[jRow] -= value*element[j];
        }
      }
      }
      if (columnLower_||columnUpper_) {
      double * solution = new double [numberGubColumns_];
      for (iSet=0;iSet<numberSets_;iSet++) {
        int j= startSet_[iSet];
        while (j>=0) {
          double value=0.0;
          if (getDynamicStatus(j)!=inSmall) {
            if (getDynamicStatus(j)==atLowerBound) {
            if (columnLower_) 
              value= columnLower_[j];
            } else if (getDynamicStatus(j)==atUpperBound) { 
            value= columnUpper_[j];
            } else if (getDynamicStatus(j)==soloKey) {
            value =keyValue(iSet);
            }
            objectiveOffset += value*cost_[j];
          } 
          solution[j]=value;
          j = next_[j]; //onto next in set
        }
      }
      // ones in gub and in small problem
      for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
        if (model_->getStatus(iColumn)!=ClpSimplex::basic) {
          int jFull = id_[iColumn-firstDynamic_];
          solution[jFull]=smallSolution[iColumn];
        }
      }
      for (iSet=0;iSet<numberSets_;iSet++) {
        int kRow = toIndex_[iSet];
        if (kRow>=0)
          kRow += numberStaticRows_;
        int j= startSet_[iSet];
        while (j>=0) {
          double value = solution[j];
          if (value) {
            for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) {
            int iRow = row_[k];
            rhsOffset_[iRow] -= element_[k]*value;
            }
            if (kRow>=0)
            rhsOffset_[kRow] -= value;
          }
          j = next_[j]; //onto next in set
        }
      }
      delete [] solution;
      } else {
      // bounds
      ClpSimplex::Status iStatus;
      for (int iSet=0;iSet<numberSets_;iSet++) {
        int kRow = toIndex_[iSet];
        if (kRow<0) {
          int iColumn = keyVariable_[iSet];
          if (iColumn<maximumGubColumns_) {
            // key is not treated as basic
            double b=0.0;
            // key is structural - where is slack
            iStatus = getStatus(iSet);
            assert (iStatus!=ClpSimplex::basic);
            if (iStatus==ClpSimplex::atLowerBound)
            b=lowerSet_[iSet];
            else
            b=upperSet_[iSet];
            if (b) {
            objectiveOffset += b*cost_[iColumn];
            for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
              int iRow = row_[j];
              rhsOffset_[iRow] -= element_[j]*b;
            }
            }
          }
        }
      }
      }
      model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
#ifdef CLP_DEBUG
      if (saveE) {
      int iRow;
      for (iRow=0;iRow<numberRows;iRow++) {
        if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3)
          printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]);
      }
      delete [] saveE;
      }
#endif
      lastRefresh_ = model->numberIterations();
    }
  }
  return rhsOffset_;
}
/*
  update information for a pivot (and effective rhs)
*/
int 
00819 ClpDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
{
  
  // now update working model
  int sequenceIn = model->sequenceIn();
  int sequenceOut = model->sequenceOut();
  int numberColumns = model->numberColumns();
  if (sequenceIn!=sequenceOut&&sequenceIn<numberColumns) 
    backToPivotRow_[sequenceIn]=model->pivotRow();
  if (sequenceIn>=firstDynamic_&&sequenceIn<numberColumns) {
    int bigSequence = id_[sequenceIn-firstDynamic_];
    if (getDynamicStatus(bigSequence)!=inSmall) {
      firstAvailable_++;
      setDynamicStatus(bigSequence,inSmall);
    }
  }
  // make sure slack is synchronized
  if (sequenceIn>=numberColumns+numberStaticRows_) {
    int iDynamic = sequenceIn-numberColumns-numberStaticRows_;
    int iSet = fromIndex_[iDynamic];
    setStatus(iSet,model->getStatus(sequenceIn));
  }
  if (sequenceOut>=numberColumns+numberStaticRows_) {
    int iDynamic = sequenceOut-numberColumns-numberStaticRows_;
    int iSet = fromIndex_[iDynamic];
    // out may have gone through barrier - so check
    double valueOut = model->lowerRegion()[sequenceOut];
    if (fabs(valueOut -static_cast<double> (lowerSet_[iSet]))<
      fabs(valueOut -static_cast<double> (upperSet_[iSet])))
      setStatus(iSet,ClpSimplex::atLowerBound);
    else
      setStatus(iSet,ClpSimplex::atUpperBound);
    if (lowerSet_[iSet]==upperSet_[iSet])
      setStatus(iSet,ClpSimplex::isFixed);
    if (getStatus(iSet)!=model->getStatus(sequenceOut))
      printf("** set %d status %d, var status %d\n",iSet,
           getStatus(iSet),model->getStatus(sequenceOut));
  }
  ClpMatrixBase::updatePivot(model,oldInValue,oldOutValue);
#ifdef CLP_DEBUG
  char * inSmall = new char [numberGubColumns_];
  memset(inSmall,0,numberGubColumns_);
  const double * solution = model->solutionRegion();
  for (int i=0;i<numberGubColumns_;i++)
    if (getDynamicStatus(i)==ClpDynamicMatrix::inSmall)
      inSmall[i]=1;
  for (int i=firstDynamic_;i<firstAvailable_;i++) {
    int k=id_[i-firstDynamic_];
    inSmall[k]=0;
    //if (k>=23289&&k<23357&&solution[i])
    //printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
  }
  for (int i=0;i<numberGubColumns_;i++)
    assert (!inSmall[i]);
  delete [] inSmall;
#ifndef NDEBUG
  for (int i=0;i<numberActiveSets_;i++) {
    int iSet = fromIndex_[i];
    assert (toIndex_[iSet]==i);
    //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
    //printf("*** set %d status %d, var status %d\n",iSet,
    //     getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
    //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
    //if (iSet==1035) {
    //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
    //     model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
    //}
  }
#endif
#endif
  return 0;
}
/*
     utility dual function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void 
00897 ClpDynamicMatrix::dualExpanded(ClpSimplex * model,
                         CoinIndexedVector * /*array*/,
                         double * /*other*/,int mode)
{
  switch (mode) {
    // modify costs before transposeUpdate
  case 0:
    break;
    // create duals for key variables (without check on dual infeasible)
  case 1:
    break;
    // as 1 but check slacks and compute djs
  case 2:
    {
      // do pivots
      int * pivotVariable = model->pivotVariable();
      int numberRows=numberStaticRows_+numberActiveSets_;
      int numberColumns = model->numberColumns();
      for (int iRow=0;iRow<numberRows;iRow++) {
      int iPivot=pivotVariable[iRow];
      if (iPivot<numberColumns) 
        backToPivotRow_[iPivot]=iRow;
      }
      if (noCheck_>=0) {
      if (infeasibilityWeight_!=model_->infeasibilityCost()) {
        // don't bother checking
        sumDualInfeasibilities_=100.0;
        numberDualInfeasibilities_=1;
        sumOfRelaxedDualInfeasibilities_ =100.0;
        return;
      }
      }
      // In theory we should subtract out ones we have done but ....
      // If key slack then dual 0.0
      // If not then slack could be dual infeasible
      // dj for key is zero so that defines dual on set
      int i;
      double * dual = model->dualRowSolution();
      double dualTolerance = model->dualTolerance();
      double relaxedTolerance=dualTolerance;
      // we can't really trust infeasibilities if there is dual error
      double error = CoinMin(1.0e-2,model->largestDualError());
      // allow tolerance at least slightly bigger than standard
      relaxedTolerance = relaxedTolerance +  error;
      // but we will be using difference
      relaxedTolerance -= dualTolerance;
      sumDualInfeasibilities_=0.0;
      numberDualInfeasibilities_=0;
      sumOfRelaxedDualInfeasibilities_ =0.0;
      for (i=0;i<numberSets_;i++) {
      double value=0.0;
      int gubRow = toIndex_[i];
      if (gubRow<0) {
        int kColumn = keyVariable_[i];
        if (kColumn<maximumGubColumns_) {
          // dj without set
          value = cost_[kColumn];
          for (CoinBigIndex j=startColumn_[kColumn];
             j<startColumn_[kColumn+1];j++) {
            int iRow = row_[j];
            value -= dual[iRow]*element_[j];
          }
          double infeasibility = 0.0;
          if (getStatus(i)==ClpSimplex::atLowerBound) {
            if (-value>dualTolerance) 
            infeasibility=-value-dualTolerance;
          } else if (getStatus(i)==ClpSimplex::atUpperBound) {
            if (value>dualTolerance) 
            infeasibility=value-dualTolerance;
          }
          if (infeasibility>0.0) {
            sumDualInfeasibilities_ += infeasibility;
            if (infeasibility>relaxedTolerance) 
            sumOfRelaxedDualInfeasibilities_ += infeasibility;
            numberDualInfeasibilities_ ++;
          }
        }
      } else {
        value = dual[gubRow+numberStaticRows_];
      }
      // Now subtract out from all 
      int k= startSet_[i];
      while (k>=0) {
        if (getDynamicStatus(k)!=inSmall) {
          double djValue = cost_[k]-value;
          for (CoinBigIndex j=startColumn_[k];
             j<startColumn_[k+1];j++) {
            int iRow = row_[j];
            djValue -= dual[iRow]*element_[j];
          }
          double infeasibility=0.0;
          if (getDynamicStatus(k)==atLowerBound) {
            if (djValue<-dualTolerance) 
            infeasibility=-djValue-dualTolerance;
          } else if (getDynamicStatus(k)==atUpperBound) {
            // at upper bound
            if (djValue>dualTolerance) 
            infeasibility=djValue-dualTolerance;
          }
          if (infeasibility>0.0) {
            sumDualInfeasibilities_ += infeasibility;
            if (infeasibility>relaxedTolerance) 
            sumOfRelaxedDualInfeasibilities_ += infeasibility;
            numberDualInfeasibilities_ ++;
          }
        }
        k = next_[k]; //onto next in set
      }
      }
    }
    infeasibilityWeight_=-1.0;
    break;
    // Report on infeasibilities of key variables
  case 3:
    {
      model->setSumDualInfeasibilities(model->sumDualInfeasibilities()+
                               sumDualInfeasibilities_);
      model->setNumberDualInfeasibilities(model->numberDualInfeasibilities()+
                               numberDualInfeasibilities_);
      model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities()+
                               sumOfRelaxedDualInfeasibilities_);
    }
    break;
    // modify costs before transposeUpdate for partial pricing
  case 4:
    break;
  }
}
/*
     general utility function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
int
01031 ClpDynamicMatrix::generalExpanded(ClpSimplex * model,int mode,int &number)
{
  int returnCode=0;
  switch (mode) {
    // Fill in pivotVariable
  case 0:
    {
      // If no effective rhs - form it
      if (!rhsOffset_) {
      rhsOffset_ = new double[model->numberRows()];
      rhsOffset(model,true);
      }
      int i;
      int numberBasic=number;
      int numberColumns = model->numberColumns();
      // Use different array so can build from true pivotVariable_
      //int * pivotVariable = model->pivotVariable();
      int * pivotVariable = model->rowArray(0)->getIndices();
      for (i=0;i<numberColumns;i++) {
      if (model->getColumnStatus(i) == ClpSimplex::basic) 
        pivotVariable[numberBasic++]=i;
      }
      number = numberBasic;
    }
    break;
    // Do initial extra rows + maximum basic
  case 2:
    {
      number = model->numberRows();
    }
    break;
    // Before normal replaceColumn
  case 3:
    {
      if (numberActiveSets_+numberStaticRows_==model_->numberRows()) {
      // no space - re-factorize
      returnCode=4;
      number=-1; // say no need for normal replaceColumn
      }
    }
    break;
    // To see if can dual or primal
  case 4:
    {
      returnCode= 1;
    }
    break;
    // save status
  case 5:
    {
    }
    break;
    // restore status
  case 6:
    {
      // maybe mismatch - redo all in small model
      //for (int i=firstDynamic_;i<firstAvailable_;i++) {
      //int jColumn = id_[i-firstDynamic_];
      //setDynamicStatus(jColumn,inSmall);
      //}
    }
    break;
    // unflag all variables
  case 8:
    {
      for (int i=0;i<numberGubColumns_;i++) {
      if (flagged(i)) {
        unsetFlagged(i);
        returnCode++;
      }
      }
    }
    break;
    // redo costs in primal
  case 9:
    {
      double * cost = model->costRegion();
      double * solution = model->solutionRegion();
      double * columnLower = model->lowerRegion();
      double * columnUpper = model->upperRegion();
      int i;
      bool doCosts = (number&4)!=0;
      bool doBounds = (number&1)!=0;
      for ( i=firstDynamic_;i<firstAvailable_;i++) {
      int jColumn = id_[i-firstDynamic_];
      if (doBounds) {
        if (!columnLower_&&!columnUpper_) {
          columnLower[i] = 0.0;
          columnUpper[i] = COIN_DBL_MAX;
        }  else {
          if (columnLower_) 
            columnLower[i] = columnLower_[jColumn];
          else
            columnLower[i] = 0.0;
          if (columnUpper_)
            columnUpper[i] = columnUpper_[jColumn];
          else
            columnUpper[i] = COIN_DBL_MAX;
        }
      }
      if (doCosts) {
        cost[i]=cost_[jColumn];
        // Original bounds
        if (model->nonLinearCost())
          model->nonLinearCost()->setOne(i,solution[i],
                                 this->columnLower(jColumn),
                                 this->columnUpper(jColumn),cost_[jColumn]);
      }
      }
      // and active sets
      for (i=0;i<numberActiveSets_;i++ ) {
      int iSet = fromIndex_[i];
      int iSequence=lastDynamic_+numberStaticRows_+i;
      if (doBounds) {
        if (lowerSet_[iSet]>-1.0e20)
          columnLower[iSequence] = lowerSet_[iSet];
        else
          columnLower[iSequence]=-COIN_DBL_MAX;
        if (upperSet_[iSet]<1.0e20)
          columnUpper[iSequence] = upperSet_[iSet];
        else
          columnUpper[iSequence]=COIN_DBL_MAX;
      }
      if (doCosts) {
        if (model->nonLinearCost()) {
          double trueLower;
          if (lowerSet_[iSet]>-1.0e20)
            trueLower = lowerSet_[iSet];
          else
            trueLower=-COIN_DBL_MAX;
          double trueUpper;
          if (upperSet_[iSet]<1.0e20)
            trueUpper = upperSet_[iSet];
          else
            trueUpper=COIN_DBL_MAX;
          model->nonLinearCost()->setOne(iSequence,solution[iSequence],
                                 trueLower,trueUpper,0.0);
        }
      }
      }
    }
    break;
    // return 1 if there may be changing bounds on variable (column generation)
  case 10:
    {
      // return 1 as bounds on rhs will change
      returnCode = 1;
    }
    break;
    // make sure set is clean
  case 7:
    {
      // first flag
      if (number>=firstDynamic_&&number<lastDynamic_) {
      assert (number==model->sequenceIn());
      // id will be sitting at firstAvailable
      int sequence = id_[firstAvailable_-firstDynamic_];
      setFlagged(sequence);
      model->clearFlagged(firstAvailable_);
      }
    }
  case 11:
    {
      int sequenceIn = model->sequenceIn();
      if (sequenceIn>=firstDynamic_&&sequenceIn<lastDynamic_) {
      assert (number==model->sequenceIn());
      // take out variable (but leave key)
      double * cost = model->costRegion();
      double * columnLower = model->lowerRegion();
      double * columnUpper = model->upperRegion();
      double * solution = model->solutionRegion();
      int * length = matrix_->getMutableVectorLengths();
#ifndef NDEBUG
      if (length[sequenceIn]) {
        int * row = matrix_->getMutableIndices();
        CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
        int iRow = row[startColumn[sequenceIn]+length[sequenceIn]-1];
        int which = iRow-numberStaticRows_;
        assert (which>=0);
        int iSet = fromIndex_[which];
        assert (toIndex_[iSet]==which);
      }
#endif
      // no need firstAvailable_--;
      solution[firstAvailable_]=0.0;
      cost[firstAvailable_]=0.0;
      length[firstAvailable_]=0;
      model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
      model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
      columnLower[firstAvailable_]=0.0;
      columnUpper[firstAvailable_]=COIN_DBL_MAX;

      // not really in small problem
      int iBig=id_[sequenceIn-firstDynamic_];
      if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound) {
        setDynamicStatus(iBig,atLowerBound);
        if (columnLower_)
          modifyOffset(sequenceIn,columnLower_[iBig]);
      } else {
        setDynamicStatus(iBig,atUpperBound);
        modifyOffset(sequenceIn,columnUpper_[iBig]);
      }
      }
    }
    break;
  default:
    break;
  }
  return returnCode;
}
/* Purely for column generation and similar ideas.  Allows
   matrix and any bounds or costs to be updated (sensibly).
   Returns non-zero if any changes.
*/
int 
01246 ClpDynamicMatrix::refresh(ClpSimplex * model)
{
  // If at beginning then create initial problem
  if (firstDynamic_==firstAvailable_) {
    initialProblem();
    return 1;
  } else if (!model->nonLinearCost()) {
    // will be same as last time
    return 1;
  }
  // lookup array
  int * active = new int [numberActiveSets_];
  CoinZeroN(active,numberActiveSets_);
  int iColumn;
  int numberColumns = model->numberColumns();
  double * element =  matrix_->getMutableElements();
  int * row = matrix_->getMutableIndices();
  CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
  int * length = matrix_->getMutableVectorLengths();
  double * cost = model->costRegion();
  double * columnLower = model->lowerRegion();
  double * columnUpper = model->upperRegion();
  CoinBigIndex numberElements=startColumn[firstDynamic_];
  // first just do lookup and basic stuff
  int currentNumber = firstAvailable_;
  firstAvailable_ = firstDynamic_;
  double objectiveChange=0.0;
  double * solution = model->solutionRegion();
  int currentNumberActiveSets=numberActiveSets_;
  for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
    int iRow = row[startColumn[iColumn]+length[iColumn]-1];
    int which = iRow-numberStaticRows_;
    assert (which>=0);
    int iSet = fromIndex_[which];
    assert (toIndex_[iSet]==which);
    if (model->getStatus(iColumn)==ClpSimplex::basic) {
      active[which]++;
      // may as well make key
      keyVariable_[iSet]=id_[iColumn-firstDynamic_];
    }
  }
  int i;
  numberActiveSets_=0;
  int numberDeleted=0;
  for (i=0;i<currentNumberActiveSets;i++) {
    int iRow = i+numberStaticRows_;
    int numberActive = active[i];
    int iSet = fromIndex_[i];
    if (model->getRowStatus(iRow)==ClpSimplex::basic) {
      numberActive++;
      // may as well make key
      keyVariable_[iSet]=maximumGubColumns_+iSet;
    }
    if (numberActive>1) {
      // keep
      active[i]=numberActiveSets_;
      numberActiveSets_++;
    } else {
      active[i]=-1;
    }
  }

  for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
    int iRow = row[startColumn[iColumn]+length[iColumn]-1];
    int which = iRow-numberStaticRows_;
    int jColumn = id_[iColumn-firstDynamic_];
    if (model->getStatus(iColumn)==ClpSimplex::atLowerBound||
      model->getStatus(iColumn)==ClpSimplex::atUpperBound) {
      double value = solution[iColumn];
      bool toLowerBound=true;
      if (columnUpper_) {
      if (!columnLower_) {
        if (fabs(value-columnUpper_[jColumn])<fabs(value)) 
        toLowerBound=false;
      } else if (fabs(value-columnUpper_[jColumn])<fabs(value-columnLower_[iColumn])) {
        toLowerBound=false;
      }
      }
      if (toLowerBound) {
      // throw out to lower bound
      if (columnLower_) {
        setDynamicStatus(jColumn,atLowerBound);
        // treat solution as if exactly at a bound
        double value = columnLower[iColumn];
        objectiveChange += cost[iColumn]*value;
      } else {
        setDynamicStatus(jColumn,atLowerBound);
      }
      } else {
      // throw out to upper bound
      setDynamicStatus(jColumn,atUpperBound);
      double value = columnUpper[iColumn];
      objectiveChange += cost[iColumn]*value;
      }
      numberDeleted++;
    } else {
      assert(model->getStatus(iColumn)!=ClpSimplex::superBasic); // deal with later
      int iPut = active[which];
      if (iPut>=0) {
      // move
      id_[firstAvailable_-firstDynamic_] = jColumn;
      int numberThis = startColumn_[jColumn+1]-startColumn_[jColumn]+1;
      length[firstAvailable_]=numberThis;
      cost[firstAvailable_]=cost[iColumn];
      columnLower[firstAvailable_]=columnLower[iColumn];
      columnUpper[firstAvailable_]=columnUpper[iColumn];
      model->nonLinearCost()->setOne(firstAvailable_,solution[iColumn],0.0,COIN_DBL_MAX,
                               cost_[jColumn]);
      CoinBigIndex base = startColumn_[jColumn];
      for (int j=0;j<numberThis-1;j++) {
        row[numberElements]=row_[base+j];
        element[numberElements++]=element_[base+j];
      }
      row[numberElements]=iPut+numberStaticRows_;
      element[numberElements++]=1.0;
      model->setStatus(firstAvailable_,model->getStatus(iColumn));
      solution[firstAvailable_] = solution[iColumn];
      firstAvailable_++;
      startColumn[firstAvailable_]=numberElements;
      }
    }
  }
  // zero solution
  CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
  // zero cost
  CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
  // zero lengths
  CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
  for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
    model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
    model->setStatus(iColumn,ClpSimplex::atLowerBound);
    columnLower[iColumn]=0.0;
    columnUpper[iColumn]=COIN_DBL_MAX;
  }
  // move rhs and set rest to infinite
  numberActiveSets_=0;
  for (i=0;i<currentNumberActiveSets;i++) {
    int iSet = fromIndex_[i];
    assert (toIndex_[iSet]==i);
    int iRow = i+numberStaticRows_;
    int iPut=active[i];
    if (iPut>=0) {
      int in = iRow+numberColumns;
      int out = iPut+numberColumns+numberStaticRows_;
      solution[out]=solution[in];
      columnLower[out]=columnLower[in];
      columnUpper[out]=columnUpper[in];
      cost[out]=cost[in];
      double trueLower;
      if (lowerSet_[iSet]>-1.0e20)
      trueLower = lowerSet_[iSet];
      else
      trueLower=-COIN_DBL_MAX;
      double trueUpper;
      if (upperSet_[iSet]<1.0e20)
      trueUpper = upperSet_[iSet];
      else
      trueUpper=COIN_DBL_MAX;
      model->nonLinearCost()->setOne(out,solution[out],
                             trueLower,trueUpper,0.0);
      model->setStatus(out,model->getStatus(in));
      toIndex_[iSet]=numberActiveSets_;
      rhsOffset_[numberActiveSets_+numberStaticRows_]= rhsOffset_[i+numberStaticRows_];
      fromIndex_[numberActiveSets_++]=fromIndex_[i];
    } else {
      // adjust offset
      // put out as key
      int jColumn = keyVariable_[iSet];
      toIndex_[iSet]=-1;
      if (jColumn<maximumGubColumns_) {
      setDynamicStatus(jColumn,soloKey);
      double value = keyValue(iSet);
      objectiveChange += cost_[jColumn]*value;
      modifyOffset(jColumn,-value);
      }
    }
  }
  model->setObjectiveOffset(model->objectiveOffset()-objectiveChange);
  delete [] active;
  for (i=numberActiveSets_;i<currentNumberActiveSets;i++) {
    int iSequence = i+numberStaticRows_+numberColumns;
    solution[iSequence]=0.0;
    columnLower[iSequence]=-COIN_DBL_MAX;
    columnUpper[iSequence]=COIN_DBL_MAX;
    cost[iSequence]=0.0;
    model->nonLinearCost()->setOne(iSequence,solution[iSequence],
                               columnLower[iSequence],
                               columnUpper[iSequence],0.0);
    model->setStatus(iSequence,ClpSimplex::basic);
    rhsOffset_[i+numberStaticRows_]=0.0;
  }
  if (currentNumberActiveSets!=numberActiveSets_||numberDeleted) {
    // deal with pivotVariable
    int * pivotVariable = model->pivotVariable();
    int sequence=firstDynamic_;
    int iRow=0;
    int base1 = firstDynamic_;
    int base2 = lastDynamic_;
    int base3 = numberColumns+numberStaticRows_;
    int numberRows = numberStaticRows_+currentNumberActiveSets;
    while (sequence<firstAvailable_) {
      int iPivot=pivotVariable[iRow];
      while (iPivot<base1||(iPivot>=base2&&iPivot<base3)) {
      iPivot= pivotVariable[++iRow];
      }
      pivotVariable[iRow++]=sequence;
      sequence++;
    }
    // move normal basic ones down
    int iPut=iRow;
    for (;iRow<numberRows;iRow++) {
      int iPivot=pivotVariable[iRow];
      if (iPivot<base1||(iPivot>=base2&&iPivot<base3)) {
      pivotVariable[iPut++]=iPivot;
      }
    }
    // and basic others
    for (i=0;i<numberActiveSets_;i++) {
      if (model->getRowStatus(i+numberStaticRows_)==ClpSimplex::basic) {
      pivotVariable[iPut++]=i+base3;
      }
    }
    for (i=numberActiveSets_;i<currentNumberActiveSets;i++) {
      pivotVariable[iPut++]=i+base3;
    }
    assert (iPut==numberRows);
  }
#ifdef CLP_DEBUG
#if 0
  printf("row for set 244 is %d, row status %d value %g ",toIndex_[244],status_[244],
       keyValue(244));
  int jj=startSet_[244];
  while (jj>=0) {
    printf("var %d status %d ",jj,dynamicStatus_[jj]);
    jj=next_[jj];
  }
  printf("\n");
#endif
#ifdef CLP_DEBUG
 {
   // debug coding to analyze set
   int i;
   int kSet=-6;
   if (kSet>=0) {
     int * back = new int [numberGubColumns_];
     for (i=0;i<numberGubColumns_;i++)
       back[i]=-1;
     for (i=firstDynamic_;i<firstAvailable_;i++)
       back[id_[i-firstDynamic_]]= i;
     int sequence = startSet_[kSet];
     if (toIndex_[kSet]<0){
       printf("Not in - Set %d - slack status %d, key %d\n",kSet,status_[kSet],keyVariable_[kSet]);
       while (sequence>=0) {
       printf("( %d status %d ) ",sequence,getDynamicStatus(sequence));
       sequence = next_[sequence];
       }
     } else {
       int iRow=numberStaticRows_+toIndex_[kSet];
       printf("In - Set %d - slack status %d, key %d offset %g slack %g\n",kSet,status_[kSet],keyVariable_[kSet],
            rhsOffset_[iRow],model->solutionRegion(0)[iRow]);
       while (sequence>=0) {
       int iBack = back[sequence];
       printf("( %d status %d value %g) ",sequence,getDynamicStatus(sequence),model->solutionRegion()[iBack]);
       sequence = next_[sequence];
       }
     }
     printf("\n");
     delete [] back;
   }
 }
#endif
  int n=numberActiveSets_;
  for (i=0;i<numberSets_;i++) {
    if (toIndex_[i]<0) {
      //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
      n++;
    }
  }
  assert (n==numberSets_);
#endif
  return 1;
}
void 
01529 ClpDynamicMatrix::times(double scalar,
                     const double * x, double * y) const
{
  if (model_->specialOptions()!=16) {
    ClpPackedMatrix::times(scalar,x,y);
  } else {
    int iRow;
    const double * element =  matrix_->getElements();
    const int * row = matrix_->getIndices();
    const CoinBigIndex * startColumn = matrix_->getVectorStarts();
    const int * length = matrix_->getVectorLengths();
    int * pivotVariable = model_->pivotVariable();
    for (iRow=0;iRow<numberStaticRows_+numberActiveSets_;iRow++) {
      y[iRow] -= scalar*rhsOffset_[iRow];
      int iColumn = pivotVariable[iRow];
      if (iColumn<lastDynamic_) {
      CoinBigIndex j;
      double value = scalar*x[iColumn];
      if (value) {
        for (j=startColumn[iColumn];
             j<startColumn[iColumn]+length[iColumn];j++) {
          int jRow=row[j];
          y[jRow] += value*element[j];
        }
      }
      }
    }
  }
}
// Modifies rhs offset
void 
01560 ClpDynamicMatrix::modifyOffset(int sequence, double amount)
{
  if (amount) {
    assert (rhsOffset_);
    CoinBigIndex j;
    for (j=startColumn_[sequence];j<startColumn_[sequence+1];j++) {
      int iRow = row_[j];
      rhsOffset_[iRow] += amount*element_[j];
    }
  }
}
// Gets key value when none in small
double 
01573 ClpDynamicMatrix::keyValue(int iSet) const
{
  double value=0.0;
  if (toIndex_[iSet]<0) {
    int key = keyVariable_[iSet];
    if (key<maximumGubColumns_) {
      if (getStatus(iSet)==ClpSimplex::atLowerBound)
      value=lowerSet_[iSet];
      else
      value=upperSet_[iSet];
      int numberKey=0;
      int j= startSet_[iSet];
      while (j>=0) {
      DynamicStatus status = getDynamicStatus(j);
      assert (status!=inSmall);
      if (status==soloKey) {
        numberKey++;
      } else if (status==atUpperBound) {
        value -= columnUpper_[j];
      } else if (columnLower_) {
        value -= columnLower_[j];
      }
      j = next_[j]; //onto next in set
      }
      assert (numberKey==1);
    } else {
      int j= startSet_[iSet];
      while (j>=0) {
      DynamicStatus status = getDynamicStatus(j);
      assert (status!=inSmall);
      assert (status!=soloKey);
      if (status==atUpperBound) {
        value += columnUpper_[j];
      } else if (columnLower_) {
        value += columnLower_[j];
      }
      j = next_[j]; //onto next in set
      }
    }
  }
  return value;
}
// Switches off dj checking each factorization (for BIG models)
void 
01617 ClpDynamicMatrix::switchOffCheck()
{
  noCheck_=0;
  infeasibilityWeight_=0.0;
}
/* Creates a variable.  This is called after partial pricing and may modify matrix.
   May update bestSequence.
*/
void 
01626 ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence)
{
  int numberRows = model->numberRows();
  int slackOffset = lastDynamic_+numberRows;
  int structuralOffset = slackOffset+numberSets_;
  int bestSequence2=savedBestSequence_-structuralOffset;
  if (bestSequence>=slackOffset) {
    double * columnLower = model->lowerRegion();
    double * columnUpper = model->upperRegion();
    double * solution = model->solutionRegion();
    double * reducedCost = model->djRegion();
    const double * duals = model->dualRowSolution();
    if (toIndex_[savedBestSet_]<0) {
      // need to put key into basis
      int newRow = numberActiveSets_+numberStaticRows_;
      model->dualRowSolution()[newRow]=savedBestGubDual_;
      double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
      toIndex_[savedBestSet_]=numberActiveSets_;
      fromIndex_[numberActiveSets_++]=savedBestSet_;
      int iSequence=lastDynamic_+newRow;
      // we need to get lower and upper correct
      double shift=0.0;
      int j= startSet_[savedBestSet_];
      while (j>=0) {
      if (getDynamicStatus(j)==atUpperBound) 
        shift += columnUpper_[j];
      else if (getDynamicStatus(j)==atLowerBound&&columnLower_)
        shift += columnLower_[j];
      j = next_[j]; //onto next in set
      }
      if (lowerSet_[savedBestSet_]>-1.0e20)
      columnLower[iSequence] = lowerSet_[savedBestSet_];
      else
      columnLower[iSequence]=-COIN_DBL_MAX;
      if (upperSet_[savedBestSet_]<1.0e20)
      columnUpper[iSequence] = upperSet_[savedBestSet_];
      else
      columnUpper[iSequence]=COIN_DBL_MAX;
#ifdef CLP_DEBUG
      if (model_->logLevel()==63) {
      printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
             bestSequence2,savedBestSet_,keyVariable_[savedBestSet_],
             columnLower[iSequence],columnUpper[iSequence],valueOfKey);
      int j= startSet_[savedBestSet_];
      while (j>=0) {
        if (getDynamicStatus(j)==atUpperBound) 
          printf("%d atup ",j);
        else if (getDynamicStatus(j)==atLowerBound)
          printf("%d atlo ",j);
        else if (getDynamicStatus(j)==soloKey)
          printf("%d solo ",j);
        else 
          abort();
        j = next_[j]; //onto next in set
      }
      printf("\n");
      }
#endif
      if (keyVariable_[savedBestSet_]<maximumGubColumns_) {
      // slack not key
      model_->pivotVariable()[newRow]=firstAvailable_;
      backToPivotRow_[firstAvailable_]=newRow;
      model->setStatus(iSequence,getStatus(savedBestSet_));
      model->djRegion()[iSequence] = savedBestGubDual_;
      if (getStatus(savedBestSet_)==ClpSimplex::atUpperBound)
        solution[iSequence] = columnUpper[iSequence];
      else
        solution[iSequence] = columnLower[iSequence];
      // create variable and pivot in
      int key=keyVariable_[savedBestSet_];
      setDynamicStatus(key,inSmall);
      double * element =  matrix_->getMutableElements();
      int * row = matrix_->getMutableIndices();
      CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
      int * length = matrix_->getMutableVectorLengths();
      CoinBigIndex numberElements = startColumn[firstAvailable_];
      int numberThis = startColumn_[key+1]-startColumn_[key]+1;
      if (numberElements+numberThis>numberElements_) {
        // need to redo
        numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
        matrix_->reserve(lastDynamic_,numberElements_);
        element =  matrix_->getMutableElements();
        row = matrix_->getMutableIndices();
        // these probably okay but be safe
        startColumn = matrix_->getMutableVectorStarts();
        length = matrix_->getMutableVectorLengths();
      }
      // already set startColumn[firstAvailable_]=numberElements;
      length[firstAvailable_]=numberThis;
      model->costRegion()[firstAvailable_]=cost_[key];
      CoinBigIndex base = startColumn_[key];
      for (int j=0;j<numberThis-1;j++) {
        row[numberElements]=row_[base+j];
        element[numberElements++]=element_[base+j];
      }
      row[numberElements]=newRow;
      element[numberElements++]=1.0;
      id_[firstAvailable_-firstDynamic_]=key;
      model->setObjectiveOffset(model->objectiveOffset()+cost_[key]*
                          valueOfKey);
      model->solutionRegion()[firstAvailable_]= valueOfKey;
      model->setStatus(firstAvailable_,ClpSimplex::basic);
      // ***** need to adjust effective rhs
      if (!columnLower_&&!columnUpper_) {
        columnLower[firstAvailable_] = 0.0;
        columnUpper[firstAvailable_] = COIN_DBL_MAX;
      }  else {
        if (columnLower_) 
          columnLower[firstAvailable_] = columnLower_[key];
        else
          columnLower[firstAvailable_] = 0.0;
        if (columnUpper_)
          columnUpper[firstAvailable_] = columnUpper_[key];
        else
          columnUpper[firstAvailable_] = COIN_DBL_MAX;
      }
      model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
                               columnLower[firstAvailable_],
                               columnUpper[firstAvailable_],cost_[key]);
      startColumn[firstAvailable_+1]=numberElements;
      reducedCost[firstAvailable_] = 0.0;
      modifyOffset(key,valueOfKey);
      rhsOffset_[newRow] = -shift; // sign?
#ifdef CLP_DEBUG
      model->rowArray(1)->checkClear();
#endif
      // now pivot in
      unpack(model,model->rowArray(1),firstAvailable_);
      model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(1));
      double alpha = model->rowArray(1)->denseVector()[newRow];
      int updateStatus = 
        model->factorization()->replaceColumn(model,
                                    model->rowArray(2),
                                    model->rowArray(1),
                                    newRow, alpha);
      model->rowArray(1)->clear();
      if (updateStatus) {
        if (updateStatus==3) {
          // out of memory
          // increase space if not many iterations
          if (model->factorization()->pivots()<
            0.5*model->factorization()->maximumPivots()&&
            model->factorization()->pivots()<400)
            model->factorization()->areaFactor(
                               model->factorization()->areaFactor() * 1.1);
        } else {
          printf("Bad returncode %d from replaceColumn\n",updateStatus);
        }
        bestSequence=-1;
        return;
      }
      // firstAvailable_ only finally updated if good pivot (in updatePivot)
      // otherwise it reverts to firstAvailableBefore_
      firstAvailable_++;
      } else {
      // slack key 
      model->setStatus(iSequence,ClpSimplex::basic);
      model->djRegion()[iSequence]=0.0;
      solution[iSequence]=shift;
      rhsOffset_[newRow] = -shift; // sign?
      }
      // correct slack
      model->costRegion()[iSequence] = 0.0;
      model->nonLinearCost()->setOne(iSequence,solution[iSequence],columnLower[iSequence],
                             columnUpper[iSequence],0.0);
    }
    if (savedBestSequence_>=structuralOffset) {
      // recompute dj and create
      double value=cost_[bestSequence2]-savedBestGubDual_;
      for (CoinBigIndex jBigIndex=startColumn_[bestSequence2];
         jBigIndex<startColumn_[bestSequence2+1];jBigIndex++) {
      int jRow=row_[jBigIndex];
      value -= duals[jRow]*element_[jBigIndex];
      }
      int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
      double * element =  matrix_->getMutableElements();
      int * row = matrix_->getMutableIndices();
      CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
      int * length = matrix_->getMutableVectorLengths();
      CoinBigIndex numberElements = startColumn[firstAvailable_];
      int numberThis = startColumn_[bestSequence2+1]-startColumn_[bestSequence2]+1;
      if (numberElements+numberThis>numberElements_) {
      // need to redo
      numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
      matrix_->reserve(lastDynamic_,numberElements_);
      element =  matrix_->getMutableElements();
      row = matrix_->getMutableIndices();
      // these probably okay but be safe
      startColumn = matrix_->getMutableVectorStarts();
      length = matrix_->getMutableVectorLengths();
      }
      // already set startColumn[firstAvailable_]=numberElements;
      length[firstAvailable_]=numberThis;
      model->costRegion()[firstAvailable_]=cost_[bestSequence2];
      CoinBigIndex base = startColumn_[bestSequence2];
      for (int j=0;j<numberThis-1;j++) {
      row[numberElements]=row_[base+j];
      element[numberElements++]=element_[base+j];
      }
      row[numberElements]=gubRow;
      element[numberElements++]=1.0;
      id_[firstAvailable_-firstDynamic_]=bestSequence2;
      //printf("best %d\n",bestSequence2);
      model->solutionRegion()[firstAvailable_]=0.0;
      if (!columnLower_&&!columnUpper_) {
      model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
      columnLower[firstAvailable_] = 0.0;
      columnUpper[firstAvailable_] = COIN_DBL_MAX;
      }  else {
      DynamicStatus status = getDynamicStatus(bestSequence2);
      if (columnLower_) 
        columnLower[firstAvailable_] = columnLower_[bestSequence2];
      else
        columnLower[firstAvailable_] = 0.0;
      if (columnUpper_)
        columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
      else
        columnUpper[firstAvailable_] = COIN_DBL_MAX;
      if (status==atLowerBound) {
        solution[firstAvailable_]=columnLower[firstAvailable_];
        model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
      } else {
        solution[firstAvailable_]=columnUpper[firstAvailable_];
        model->setStatus(firstAvailable_,ClpSimplex::atUpperBound);
      }
      }
      model->setObjectiveOffset(model->objectiveOffset()+cost_[bestSequence2]*
                        solution[firstAvailable_]);
      model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
                             columnLower[firstAvailable_],
                             columnUpper[firstAvailable_],cost_[bestSequence2]);
      bestSequence=firstAvailable_;
      // firstAvailable_ only updated if good pivot (in updatePivot)
      startColumn[firstAvailable_+1]=numberElements;
      //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_);
      reducedCost[bestSequence] = value;
    } else {
      // slack - row must just have been created
      assert (toIndex_[savedBestSet_]==numberActiveSets_-1);
      int newRow = numberStaticRows_+numberActiveSets_-1;
      bestSequence = lastDynamic_ + newRow;
      reducedCost[bestSequence] = savedBestGubDual_;
    }
  }
  // clear for next iteration
  savedBestSequence_=-1;
}
// Returns reduced cost of a variable
double 
01875 ClpDynamicMatrix::reducedCost(ClpSimplex * model,int sequence) const
{
  int numberRows = model->numberRows();
  int slackOffset = lastDynamic_+numberRows;
  if (sequence<slackOffset)
    return model->djRegion()[sequence];
  else
    return savedBestDj_;
}
// Does gub crash
void 
01886 ClpDynamicMatrix::gubCrash()
{
  // Do basis - cheapest or slack if feasible
  int longestSet=0;
  int iSet;
  for (iSet=0;iSet<numberSets_;iSet++) {
    int n=0;
    int j= startSet_[iSet];
    while (j>=0) {
      n++;
      j=next_[j];
    }
    longestSet = CoinMax(longestSet,n);
  }
  double * upper = new double[longestSet+1];
  double * cost = new double[longestSet+1];
  double * lower = new double[longestSet+1];
  double * solution = new double[longestSet+1];
  int * back = new int[longestSet+1];
  double tolerance = model_->primalTolerance();
  double objectiveOffset = 0.0;
  for (iSet=0;iSet<numberSets_;iSet++) {
    int iBasic=-1;
    double value=0.0;
    // find cheapest
    int numberInSet = 0;
    int j=startSet_[iSet];
    while (j>=0) {
      if (!columnLower_) 
      lower[numberInSet]=0.0;
      else
      lower[numberInSet]= columnLower_[j];
      if (!columnUpper_) 
      upper[numberInSet]=COIN_DBL_MAX;
      else
      upper[numberInSet]= columnUpper_[j];
      back[numberInSet++]=j;
      j = next_[j];
    }
    CoinFillN(solution,numberInSet,0.0);
    // and slack
    iBasic=numberInSet;
    solution[iBasic]=-value;
    lower[iBasic]=-upperSet_[iSet];
    upper[iBasic]=-lowerSet_[iSet];
    int kphase;
    if (value<lowerSet_[iSet]-tolerance||value>upperSet_[iSet]+tolerance) {
      // infeasible
      kphase=0;
      // remember bounds are flipped so opposite to natural
      if (value<lowerSet_[iSet]-tolerance)
      cost[iBasic]=1.0;
      else
      cost[iBasic]=-1.0;
      CoinZeroN(cost,numberInSet);
      double dualTolerance =model_->dualTolerance();
      for (int iphase =kphase;iphase<2;iphase++) {
      if (iphase) {
        cost[numberInSet]=0.0;
        for (int j=0;j<numberInSet;j++)
          cost[j]= cost_[back[j]];
      }
      // now do one row lp
      bool improve=true;
      while (improve) {
        improve=false;
        double dual = cost[iBasic];
        int chosen =-1;
        double best=dualTolerance;
        int way=0;
        for (int i=0;i<=numberInSet;i++) {
          double dj = cost[i]-dual;
          double improvement =0.0;
          double distance=0.0;
          if (iphase||i<numberInSet)
            assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
          if (dj>dualTolerance)
            improvement = dj*(solution[i]-lower[i]);
          else if (dj<-dualTolerance)
            improvement = dj*(solution[i]-upper[i]);
          if (improvement>best) {
            best=improvement;
            chosen=i;
            if (dj<0.0) {
            way = 1;
            distance = upper[i]-solution[i];
            } else {
            way = -1;
            distance = solution[i]-lower[i];
            }
          }
        }
        if (chosen>=0) {
          improve=true;
          // now see how far
          if (way>0) {
            // incoming increasing so basic decreasing
            // if phase 0 then go to nearest bound
            double distance=upper[chosen]-solution[chosen];
            double basicDistance;
            if (!iphase) {
            assert (iBasic==numberInSet);
            assert (solution[iBasic]>upper[iBasic]);
            basicDistance = solution[iBasic]-upper[iBasic];
            } else {
            basicDistance = solution[iBasic]-lower[iBasic];
            }
            // need extra coding for unbounded
            assert (CoinMin(distance,basicDistance)<1.0e20);
            if (distance>basicDistance) {
            // incoming becomes basic
            solution[chosen] += basicDistance;
            if (!iphase) 
              solution[iBasic]=upper[iBasic];
            else 
              solution[iBasic]=lower[iBasic];
            iBasic = chosen;
            } else {
            // flip
            solution[chosen]=upper[chosen];
            solution[iBasic] -= distance;
            }
          } else {
            // incoming decreasing so basic increasing
            // if phase 0 then go to nearest bound
            double distance=solution[chosen]-lower[chosen];
            double basicDistance;
            if (!iphase) {
            assert (iBasic==numberInSet);
            assert (solution[iBasic]<lower[iBasic]);
            basicDistance = lower[iBasic]-solution[iBasic];
            } else {
            basicDistance = upper[iBasic]-solution[iBasic];
            }
            // need extra coding for unbounded - for now just exit
            if (CoinMin(distance,basicDistance)>1.0e20) {
            printf("unbounded on set %d\n",iSet);
            iphase=1;
            iBasic=numberInSet;
            break;
            }
            if (distance>basicDistance) {
            // incoming becomes basic
            solution[chosen] -= basicDistance;
            if (!iphase) 
              solution[iBasic]=lower[iBasic];
            else 
              solution[iBasic]=upper[iBasic];
            iBasic = chosen;
            } else {
            // flip
            solution[chosen]=lower[chosen];
            solution[iBasic] += distance;
            }
          }
          if (!iphase) {
            if(iBasic<numberInSet)
            break; // feasible
            else if (solution[iBasic]>=lower[iBasic]&&
                   solution[iBasic]<=upper[iBasic])
            break; // feasible (on flip)
          }
        }
      }
      }
    }
    // do solution i.e. bounds
    if (columnLower_||columnUpper_) {
      for (int j=0;j<numberInSet;j++) {
      if (j!=iBasic) {
        objectiveOffset += solution[j]*cost[j];
        if (columnLower_&&columnUpper_) {
          if (fabs(solution[j]-columnLower_[back[j]])>
            fabs(solution[j]-columnUpper_[back[j]]))
            setDynamicStatus(back[j],atUpperBound);
        } else if (columnUpper_&&solution[j]>0.0) {
          setDynamicStatus(back[j],atUpperBound);
        } else {
          setDynamicStatus(back[j],atLowerBound);
          assert(!solution[j]);
        }
      }
      }
    }
    // convert iBasic back and do bounds
    if (iBasic==numberInSet) {
      // slack basic
      setStatus(iSet,ClpSimplex::basic);
      iBasic=iSet+maximumGubColumns_;
    } else {
      iBasic = back[iBasic];
      setDynamicStatus(iBasic,soloKey);
      // remember bounds flipped
      if (upper[numberInSet]==lower[numberInSet]) 
      setStatus(iSet,ClpSimplex::isFixed);
      else if (solution[numberInSet]==upper[numberInSet])
      setStatus(iSet,ClpSimplex::atLowerBound);
      else if (solution[numberInSet]==lower[numberInSet])
      setStatus(iSet,ClpSimplex::atUpperBound);
      else 
      abort();
    }
    keyVariable_[iSet]=iBasic;
  }
  model_->setObjectiveOffset(objectiveOffset_-objectiveOffset);
  delete [] lower;
  delete [] solution;
  delete [] upper;
  delete [] cost;
  delete [] back;
  // make sure matrix is in good shape
  matrix_->orderMatrix();
}
// Populates initial matrix from dynamic status
void 
02101 ClpDynamicMatrix::initialProblem()
{
  int iSet;
  double * element =  matrix_->getMutableElements();
  int * row = matrix_->getMutableIndices();
  CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
  int * length = matrix_->getMutableVectorLengths();
  double * cost = model_->objective();
  double * solution = model_->primalColumnSolution();
  double * columnLower = model_->columnLower();
  double * columnUpper = model_->columnUpper();
  double * rowSolution = model_->primalRowSolution();
  double * rowLower = model_->rowLower();
  double * rowUpper = model_->rowUpper();
  CoinBigIndex numberElements=startColumn[firstDynamic_];

  firstAvailable_ = firstDynamic_;
  numberActiveSets_=0;
  for (iSet=0;iSet<numberSets_;iSet++) {
    toIndex_[iSet]=-1;
    int numberActive=0;
    int whichKey=-1;
    if (getStatus(iSet)==ClpSimplex::basic)
      whichKey = maximumGubColumns_+iSet;
    else
      whichKey = -1;
    int j= startSet_[iSet];
    while (j>=0) {
      assert (getDynamicStatus(j)!=soloKey||whichKey==-1);
      if (getDynamicStatus(j)==inSmall)
      numberActive++;
      else if (getDynamicStatus(j)==soloKey)
      whichKey=j;
      j = next_[j]; //onto next in set
    }
    if (getStatus(iSet)==ClpSimplex::basic&&numberActive)
      numberActive++;
    if (numberActive) {
      assert (numberActive>1);
      int iRow = numberActiveSets_+numberStaticRows_;
      rowSolution[iRow]=0.0;
      double lowerValue;
      if (lowerSet_[iSet]>-1.0e20)
      lowerValue = lowerSet_[iSet];
      else
      lowerValue=-COIN_DBL_MAX;
      double upperValue;
      if (upperSet_[iSet]<1.0e20)
      upperValue = upperSet_[iSet];
      else
      upperValue=COIN_DBL_MAX;
      rowLower[iRow]=lowerValue;
      rowUpper[iRow]=upperValue;
      if (getStatus(iSet)==ClpSimplex::basic) {
      model_->setRowStatus(iRow,ClpSimplex::basic);
      rowSolution[iRow]=0.0;
      } else if (getStatus(iSet)==ClpSimplex::atLowerBound) {
      model_->setRowStatus(iRow,ClpSimplex::atLowerBound);
      rowSolution[iRow]=lowerValue;
      } else {
      model_->setRowStatus(iRow,ClpSimplex::atUpperBound);
      rowSolution[iRow]=upperValue;
      }
      j= startSet_[iSet];
      while (j>=0) {
      DynamicStatus status = getDynamicStatus(j);
      if (status==inSmall) {
        int numberThis = startColumn_[j+1]-startColumn_[j]+1;
        if (numberElements+numberThis>numberElements_) {
          // need to redo
          numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
          matrix_->reserve(lastDynamic_,numberElements_);
          element =  matrix_->getMutableElements();
          row = matrix_->getMutableIndices();
          // these probably okay but be safe
          startColumn = matrix_->getMutableVectorStarts();
          length = matrix_->getMutableVectorLengths();
        }
        length[firstAvailable_]=numberThis;
        cost[firstAvailable_]=cost_[j];
        CoinBigIndex base = startColumn_[j];
        for (int k=0;k<numberThis-1;k++) {
          row[numberElements]=row_[base+k];
          element[numberElements++]=element_[base+k];
        }
        row[numberElements]=iRow;
        element[numberElements++]=1.0;
        id_[firstAvailable_-firstDynamic_]=j;
        solution[firstAvailable_]=0.0;
        model_->setStatus(firstAvailable_,ClpSimplex::basic);
        if (!columnLower_&&!columnUpper_) {
          columnLower[firstAvailable_] = 0.0;
          columnUpper[firstAvailable_] = COIN_DBL_MAX;
        }  else {
          if (columnLower_) 
            columnLower[firstAvailable_] = columnLower_[j];
          else
            columnLower[firstAvailable_] = 0.0;
          if (columnUpper_)
            columnUpper[firstAvailable_] = columnUpper_[j];
          else
            columnUpper[firstAvailable_] = COIN_DBL_MAX;
          if (status==atLowerBound) {
            solution[firstAvailable_]=columnLower[firstAvailable_];
          } else {
            solution[firstAvailable_]=columnUpper[firstAvailable_];
          }
        }
        firstAvailable_++;
        startColumn[firstAvailable_]=numberElements;
      }
      j = next_[j]; //onto next in set
      }
      model_->setRowStatus(numberActiveSets_+numberStaticRows_,getStatus(iSet));
      toIndex_[iSet]=numberActiveSets_;
      fromIndex_[numberActiveSets_++]=iSet;
    }
    assert (toIndex_[iSet]>=0||whichKey>=0);
    keyVariable_[iSet]=whichKey;
  }
  return;
}
// Adds in a column to gub structure (called from descendant)
int 
02225 ClpDynamicMatrix::addColumn(int numberEntries,const int * row, const double * element,
                      double cost, double lower, double upper,int iSet,
                      DynamicStatus status)
{
  // check if already in
  int j=startSet_[iSet];
  while (j>=0) {
    if (startColumn_[j+1]-startColumn_[j]==numberEntries) {
      const int * row2 = row_+startColumn_[j];
      const double * element2 = element_ + startColumn_[j];
      bool same=true;
      for (int k=0;k<numberEntries;k++) {
      if (row[k]!=row2[k]||element[k]!=element2[k]) {
        same=false;
        break;
      }
      }
      if (same) {
      bool odd=false;
      if (cost!=cost_[j])
        odd =true;
      if (columnLower_&&lower!=columnLower_[j])
        odd=true;
      if (columnUpper_&&upper!=columnUpper_[j])
        odd=true;
      if (odd)
        printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
             cost,lower,upper,cost_[j],
             columnLower_ ? columnLower_[j] : 0.0,
             columnUpper_ ? columnUpper_[j] : 1.0e100);
      else
        return j;
      }
    }
    j = next_[j];
  }

  if (numberGubColumns_==maximumGubColumns_||
      startColumn_[numberGubColumns_]+numberEntries>maximumElements_) {
    CoinBigIndex j;
    int i;
    int put=0;
    int numberElements=0;
    CoinBigIndex start=0;
    // compress - leave ones at ub and basic
    int * which = new int [numberGubColumns_];
    for (i=0;i<numberGubColumns_;i++) {
      CoinBigIndex end=startColumn_[i+1];
      if (getDynamicStatus(i)!=atLowerBound) {
      // keep in
      for (j=start;j<end;j++) {
        row_[numberElements]=row_[j];
        element_[numberElements++]=element_[j];
      }
      startColumn_[put+1]=numberElements;
      cost_[put]=cost_[i];
      if (columnLower_)
        columnLower_[put]=columnLower_[i];
      if (columnUpper_)
        columnUpper_[put]=columnUpper_[i];
      dynamicStatus_[put]=dynamicStatus_[i];
      id_[put]=id_[i];
      which[i]=put;
      put++;
      } else {
      // out
      which[i]=-1;
      }
      start=end;
    }
    // now redo startSet_ and next_
    int * newNext = new int [maximumGubColumns_];
    for (int jSet=0;jSet<numberSets_;jSet++) {
      int sequence = startSet_[jSet];
      while (which[sequence]<0) {
      // out
      assert (next_[sequence]>=0);
      sequence=next_[sequence];
      }
      startSet_[jSet]=which[sequence];
      int last = which[sequence];
      while (next_[sequence]>=0) {
      sequence = next_[sequence];
      if(which[sequence]>=0) {
        // keep
        newNext[last]=which[sequence];
        last=which[sequence];
      }
      }
      newNext[last]=-jSet-1;
    }
    delete [] next_;
    next_=newNext;
    delete [] which;
    abort();
  }
  CoinBigIndex start = startColumn_[numberGubColumns_];
  CoinMemcpyN(row,numberEntries,row_+start);
  CoinMemcpyN(element,numberEntries,element_+start);
  startColumn_[numberGubColumns_+1]=start+numberEntries;
  cost_[numberGubColumns_]=cost;
  if (columnLower_) 
    columnLower_[numberGubColumns_]=lower;
  else
    assert (!lower);
  if (columnUpper_) 
    columnUpper_[numberGubColumns_]=upper;
  else
    assert (upper>1.0e20);
  setDynamicStatus(numberGubColumns_,status);
  // Do next_
  j=startSet_[iSet];
  startSet_[iSet]=numberGubColumns_;
  next_[numberGubColumns_]=j;
  numberGubColumns_++;
  return numberGubColumns_-1;
}
// Returns which set a variable is in
int 
02344 ClpDynamicMatrix::whichSet (int sequence) const
{
  while (next_[sequence]>=0)
    sequence = next_[sequence];
  int iSet = - next_[sequence]-1;
  return iSet;
}

Generated by  Doxygen 1.6.0   Back to index