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

ClpGubMatrix.cpp

/* $Id: ClpGubMatrix.cpp 1458 2009-11-05 12:34:07Z forrest $ */
// Copyright (C) 2002, 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 "ClpGubMatrix.hpp"
//#include "ClpGubDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
00029 ClpGubMatrix::ClpGubMatrix () 
  : ClpPackedMatrix(),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    infeasibilityWeight_(0.0),
    start_(NULL),
    end_(NULL),
    lower_(NULL),
    upper_(NULL),
    status_(NULL),
    saveStatus_(NULL),
    savedKeyVariable_(NULL),
    backward_(NULL),
    backToPivotRow_(NULL),
    changeCost_(NULL),
    keyVariable_(NULL),
    next_(NULL),
    toIndex_(NULL),
    fromIndex_(NULL),
    model_(NULL),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    noCheck_(-1),
    numberSets_(0),
    saveNumber_(0),
    possiblePivotKey_(0),
    gubSlackIn_(-1),
    firstGub_(0),
    lastGub_(0),
    gubType_(0)
{
  setType(16);
}

//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
00068 ClpGubMatrix::ClpGubMatrix (const ClpGubMatrix & rhs) 
: ClpPackedMatrix(rhs)
{  
  numberSets_ = rhs.numberSets_;
  saveNumber_ = rhs.saveNumber_;
  possiblePivotKey_ = rhs.possiblePivotKey_;
  gubSlackIn_ = rhs.gubSlackIn_;
  start_ = ClpCopyOfArray(rhs.start_,numberSets_);
  end_ = ClpCopyOfArray(rhs.end_,numberSets_);
  lower_ = ClpCopyOfArray(rhs.lower_,numberSets_);
  upper_ = ClpCopyOfArray(rhs.upper_,numberSets_);
  status_ = ClpCopyOfArray(rhs.status_,numberSets_);
  saveStatus_ = ClpCopyOfArray(rhs.saveStatus_,numberSets_);
  savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_,numberSets_);
  int numberColumns = getNumCols();
  backward_ = ClpCopyOfArray(rhs.backward_,numberColumns);
  backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns);
  changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_);
  fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1);
  keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
  // find longest set
  int * longest = new int[numberSets_];
  CoinZeroN(longest,numberSets_);
  int j;
  for (j=0;j<numberColumns;j++) {
    int iSet = backward_[j];
    if (iSet>=0)
      longest[iSet]++;
  }
  int length = 0;
  for (j=0;j<numberSets_;j++)
    length = CoinMax(length,longest[j]);
  next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length);
  toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
  sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
  infeasibilityWeight_=rhs.infeasibilityWeight_;
  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
  noCheck_ = rhs.noCheck_;
  firstGub_ = rhs.firstGub_;
  lastGub_ = rhs.lastGub_;
  gubType_ = rhs.gubType_;
  model_ = rhs.model_;
}

//-------------------------------------------------------------------
// assign matrix (for space reasons)
//-------------------------------------------------------------------
00119 ClpGubMatrix::ClpGubMatrix (CoinPackedMatrix * rhs) 
  : ClpPackedMatrix(rhs),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    infeasibilityWeight_(0.0),
    start_(NULL),
    end_(NULL),
    lower_(NULL),
    upper_(NULL),
    status_(NULL),
    saveStatus_(NULL),
    savedKeyVariable_(NULL),
    backward_(NULL),
    backToPivotRow_(NULL),
    changeCost_(NULL),
    keyVariable_(NULL),
    next_(NULL),
    toIndex_(NULL),
    fromIndex_(NULL),
    model_(NULL),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    noCheck_(-1),
    numberSets_(0),
    saveNumber_(0),
    possiblePivotKey_(0),
    gubSlackIn_(-1),
    firstGub_(0),
    lastGub_(0),
    gubType_(0)
{  
  setType(16);
}

/* This takes over ownership (for space reasons) and is the
   real constructor*/
00157 ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix * matrix, int numberSets,
                     const int * start, const int * end,
                     const double * lower, const double * upper,
                     const unsigned char * status)
  : ClpPackedMatrix(matrix->matrix()),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    saveNumber_(0),
    possiblePivotKey_(0),
    gubSlackIn_(-1)
{
  model_=NULL;
  numberSets_ = numberSets;
  start_ = ClpCopyOfArray(start,numberSets_);
  end_ = ClpCopyOfArray(end,numberSets_);
  lower_ = ClpCopyOfArray(lower,numberSets_);
  upper_ = ClpCopyOfArray(upper,numberSets_);
  // Check valid and ordered
  int last=-1;
  int numberColumns = matrix_->getNumCols();
  int numberRows = matrix_->getNumRows();
  backward_ = new int[numberColumns];
  backToPivotRow_ = new int[numberColumns];
  changeCost_ = new double [numberRows+numberSets_];
  keyVariable_ = new int[numberSets_];
  // signal to need new ordering
  next_ = NULL;
  for (int iColumn=0;iColumn<numberColumns;iColumn++) 
    backward_[iColumn]=-1;

  int iSet;
  for (iSet=0;iSet<numberSets_;iSet++) {
    // set key variable as slack
    keyVariable_[iSet]=iSet+numberColumns;
    if (start_[iSet]<0||start_[iSet]>=numberColumns)
      throw CoinError("Index out of range","constructor","ClpGubMatrix");
    if (end_[iSet]<0||end_[iSet]>numberColumns)
      throw CoinError("Index out of range","constructor","ClpGubMatrix");
    if (end_[iSet]<=start_[iSet])
      throw CoinError("Empty or negative set","constructor","ClpGubMatrix");
    if (start_[iSet]<last)
      throw CoinError("overlapping or non-monotonic sets","constructor","ClpGubMatrix");
    last=end_[iSet];
    int j;
    for (j=start_[iSet];j<end_[iSet];j++)
      backward_[j]=iSet;
  }
  // Find type of gub
  firstGub_=numberColumns+1;
  lastGub_=-1;
  int i;
  for (i=0;i<numberColumns;i++) {
    if (backward_[i]>=0) {
      firstGub_ = CoinMin(firstGub_,i);
      lastGub_ = CoinMax(lastGub_,i);
    }
  }
  gubType_=0;
  // adjust lastGub_
  if (lastGub_>0)
    lastGub_++;
  for (i=firstGub_;i<lastGub_;i++) {
    if (backward_[i]<0) {
      gubType_=1;
      printf("interior non gub %d\n",i);
      break;
    }
  }
  if (status) {
    status_ = ClpCopyOfArray(status,numberSets_);
  } 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);
    }
  }
  saveStatus_= new unsigned char [numberSets_];
  memset(saveStatus_,0,numberSets_);
  savedKeyVariable_= new int [numberSets_];
  memset(savedKeyVariable_,0,numberSets_*sizeof(int));
  noCheck_ = -1;
  infeasibilityWeight_=0.0;
}

00248 ClpGubMatrix::ClpGubMatrix (const CoinPackedMatrix & rhs) 
  : ClpPackedMatrix(rhs),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    infeasibilityWeight_(0.0),
    start_(NULL),
    end_(NULL),
    lower_(NULL),
    upper_(NULL),
    status_(NULL),
    saveStatus_(NULL),
    savedKeyVariable_(NULL),
    backward_(NULL),
    backToPivotRow_(NULL),
    changeCost_(NULL),
    keyVariable_(NULL),
    next_(NULL),
    toIndex_(NULL),
    fromIndex_(NULL),
    model_(NULL),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    noCheck_(-1),
    numberSets_(0),
    saveNumber_(0),
    possiblePivotKey_(0),
    gubSlackIn_(-1),
    firstGub_(0),
    lastGub_(0),
    gubType_(0)
{  
  setType(16);
  
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
00288 ClpGubMatrix::~ClpGubMatrix ()
{
  delete [] start_;
  delete [] end_;
  delete [] lower_;
  delete [] upper_;
  delete [] status_;
  delete [] saveStatus_;
  delete [] savedKeyVariable_;
  delete [] backward_;
  delete [] backToPivotRow_;
  delete [] changeCost_;
  delete [] keyVariable_;
  delete [] next_;
  delete [] toIndex_;
  delete [] fromIndex_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpGubMatrix &
ClpGubMatrix::operator=(const ClpGubMatrix& rhs)
{
  if (this != &rhs) {
    ClpPackedMatrix::operator=(rhs);
    delete [] start_;
    delete [] end_;
    delete [] lower_;
    delete [] upper_;
    delete [] status_;
    delete [] saveStatus_;
    delete [] savedKeyVariable_;
    delete [] backward_;
    delete [] backToPivotRow_;
    delete [] changeCost_;
    delete [] keyVariable_;
    delete [] next_;
    delete [] toIndex_;
    delete [] fromIndex_;
    numberSets_ = rhs.numberSets_;
    saveNumber_ = rhs.saveNumber_;
    possiblePivotKey_ = rhs.possiblePivotKey_;
    gubSlackIn_ = rhs.gubSlackIn_;
    start_ = ClpCopyOfArray(rhs.start_,numberSets_);
    end_ = ClpCopyOfArray(rhs.end_,numberSets_);
    lower_ = ClpCopyOfArray(rhs.lower_,numberSets_);
    upper_ = ClpCopyOfArray(rhs.upper_,numberSets_);
    status_ = ClpCopyOfArray(rhs.status_,numberSets_);
    saveStatus_ = ClpCopyOfArray(rhs.saveStatus_,numberSets_);
    savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_,numberSets_);
    int numberColumns = getNumCols();
    backward_ = ClpCopyOfArray(rhs.backward_,numberColumns);
    backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,numberColumns);
    changeCost_ = ClpCopyOfArray(rhs.changeCost_,getNumRows()+numberSets_);
    fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+numberSets_+1);
    keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
    // find longest set
    int * longest = new int[numberSets_];
    CoinZeroN(longest,numberSets_);
    int j;
    for (j=0;j<numberColumns;j++) {
      int iSet = backward_[j];
      if (iSet>=0)
      longest[iSet]++;
    }
    int length = 0;
    for (j=0;j<numberSets_;j++)
      length = CoinMax(length,longest[j]);
    next_ = ClpCopyOfArray(rhs.next_,numberColumns+numberSets_+2*length);
    toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
    sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
    sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
    sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
    sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
    infeasibilityWeight_=rhs.infeasibilityWeight_;
    numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
    numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
    noCheck_ = rhs.noCheck_;
    firstGub_ = rhs.firstGub_;
    lastGub_ = rhs.lastGub_;
    gubType_ = rhs.gubType_;
    model_ = rhs.model_;
  }
  return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
00377 ClpMatrixBase * ClpGubMatrix::clone() const
{
  return new ClpGubMatrix(*this);
}
/* Subset clone (without gaps).  Duplicates are allowed
   and order is as given */
ClpMatrixBase * 
00384 ClpGubMatrix::subsetClone (int numberRows, const int * whichRows,
                     int numberColumns, 
                     const int * whichColumns) const 
{
  return new ClpGubMatrix(*this, numberRows, whichRows,
                           numberColumns, whichColumns);
}
/* Returns a new matrix in reverse order without gaps
   Is allowed to return NULL if doesn't want to have row copy */
ClpMatrixBase * 
00394 ClpGubMatrix::reverseOrderedCopy() const 
{
  return NULL;
}
int 
00399 ClpGubMatrix::hiddenRows() const
{ 
  return numberSets_;
}
/* Subset constructor (without gaps).  Duplicates are allowed
   and order is as given */
00405 ClpGubMatrix::ClpGubMatrix (
                      const ClpGubMatrix & rhs,
                      int numberRows, const int * whichRows,
                      int numberColumns, const int * whichColumns)
  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns)
{
  // Assuming no gub rows deleted
  // We also assume all sets in same order
  // Get array with backward pointers
  int numberColumnsOld = rhs.matrix_->getNumCols();
  int * array = new int [ numberColumnsOld];
  int i;
  for (i=0;i<numberColumnsOld;i++)
    array[i]=-1;
  for (int iSet=0;iSet<numberSets_;iSet++) {
    for (int j=start_[iSet];j<end_[iSet];j++)
      array[j]=iSet;
  }
  numberSets_=-1;
  int lastSet=-1;
  bool inSet=false;
  for (i=0;i<numberColumns;i++) {
    int iColumn = whichColumns[i];
    int iSet=array[iColumn];
    if (iSet<0) {
      inSet=false;
    } else {
      if (!inSet) {
      // start of new set but check okay
      if (iSet<=lastSet)
        throw CoinError("overlapping or non-monotonic sets","subset constructor","ClpGubMatrix");
      lastSet = iSet;
      numberSets_++;
      start_[numberSets_]=i;
      end_[numberSets_]=i+1;
      lower_[numberSets_]=lower_[iSet];
      upper_[numberSets_]=upper_[iSet];
      inSet=true;
      } else {
      if (iSet<lastSet) {
        throw CoinError("overlapping or non-monotonic sets","subset constructor","ClpGubMatrix");
      } else if (iSet==lastSet) {
        end_[numberSets_]=i+1;
      } else {
        // new set
        lastSet = iSet;
        numberSets_++;
        start_[numberSets_]=i;
        end_[numberSets_]=i+1;
        lower_[numberSets_]=lower_[iSet];
        upper_[numberSets_]=upper_[iSet];
      }
      }
    }
  }
  delete [] array;
  numberSets_++; // adjust
  // Find type of gub
  firstGub_=numberColumns+1;
  lastGub_=-1;
  for (i=0;i<numberColumns;i++) {
    if (backward_[i]>=0) {
      firstGub_ = CoinMin(firstGub_,i);
      lastGub_ = CoinMax(lastGub_,i);
    }
  }
  if (lastGub_>0)
    lastGub_++;
  gubType_=0;
  for (i=firstGub_;i<lastGub_;i++) {
    if (backward_[i]<0) {
      gubType_=1;
      break;
    }
  }

  // Make sure key is feasible if only key in set
}
ClpGubMatrix::ClpGubMatrix (
                   const CoinPackedMatrix & rhs,
                   int numberRows, const int * whichRows,
                   int numberColumns, const int * whichColumns)
  : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns),
    sumDualInfeasibilities_(0.0),
    sumPrimalInfeasibilities_(0.0),
    sumOfRelaxedDualInfeasibilities_(0.0),
    sumOfRelaxedPrimalInfeasibilities_(0.0),
    start_(NULL),
    end_(NULL),
    lower_(NULL),
    upper_(NULL),
    backward_(NULL),
    backToPivotRow_(NULL),
    changeCost_(NULL),
    keyVariable_(NULL),
    next_(NULL),
    toIndex_(NULL),
    fromIndex_(NULL),
    numberDualInfeasibilities_(0),
    numberPrimalInfeasibilities_(0),
    numberSets_(0),
    saveNumber_(0),
    possiblePivotKey_(0),
    gubSlackIn_(-1),
    firstGub_(0),
    lastGub_(0),
    gubType_(0)
{
  setType(16);
}
/* Return <code>x * A + y</code> in <code>z</code>. 
      Squashes small elements and knows about ClpSimplex */
void 
00518 ClpGubMatrix::transposeTimes(const ClpSimplex * model, double scalar,
                        const CoinIndexedVector * rowArray,
                        CoinIndexedVector * y,
                        CoinIndexedVector * columnArray) const
{
  columnArray->clear();
  double * pi = rowArray->denseVector();
  int numberNonZero=0;
  int * index = columnArray->getIndices();
  double * array = columnArray->denseVector();
  int numberInRowArray = rowArray->getNumElements();
  // maybe I need one in OsiSimplex
  double zeroTolerance = model->zeroTolerance();
  int numberRows = model->numberRows();
  ClpPackedMatrix* rowCopy =
    dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
  bool packed = rowArray->packedMode();
  double factor = 0.3;
  // We may not want to do by row if there may be cache problems
  int numberColumns = model->numberColumns();
  // It would be nice to find L2 cache size - for moment 512K
  // Be slightly optimistic
  if (numberColumns*sizeof(double)>1000000) {
    if (numberRows*10<numberColumns)
      factor=0.1;
    else if (numberRows*4<numberColumns)
      factor=0.15;
    else if (numberRows*2<numberColumns)
      factor=0.2;
    //if (model->numberIterations()%50==0)
    //printf("%d nonzero\n",numberInRowArray);
  }
  // reduce for gub
  factor *= 0.5;
  assert (!y->getNumElements());
  if (numberInRowArray>factor*numberRows||!rowCopy) {
    // do by column
    int iColumn;
    // get matrix data pointers
    const int * row = matrix_->getIndices();
    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    const int * columnLength = matrix_->getVectorLengths(); 
    const double * elementByColumn = matrix_->getElements();
    const double * rowScale = model->rowScale();
    int numberColumns = model->numberColumns();
    int iSet = -1;
    double djMod=0.0;
    if (packed) {
      // need to expand pi into y
      assert(y->capacity()>=numberRows);
      double * piOld = pi;
      pi = y->denseVector();
      const int * whichRow = rowArray->getIndices();
      int i;
      if (!rowScale) {
      // modify pi so can collapse to one loop
      for (i=0;i<numberInRowArray;i++) {
        int iRow = whichRow[i];
        pi[iRow]=scalar*piOld[i];
      }
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
        if (backward_[iColumn]!=iSet) {
          // get pi on gub row
          iSet = backward_[iColumn];
          if (iSet>=0) {
            int iBasic = keyVariable_[iSet];
            if (iBasic<numberColumns) {
            // get dj without 
            assert (model->getStatus(iBasic)==ClpSimplex::basic);
            djMod=0.0;
            for (CoinBigIndex j=columnStart[iBasic];
                 j<columnStart[iBasic]+columnLength[iBasic];j++) {
              int jRow=row[j];
              djMod -= pi[jRow]*elementByColumn[j];
            }
            } else {
            djMod = 0.0;
            }
          } else {
            djMod = 0.0;
          }
        }
        double value = -djMod;
        CoinBigIndex j;
        for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
          int iRow = row[j];
          value += pi[iRow]*elementByColumn[j];
        }
        if (fabs(value)>zeroTolerance) {
          array[numberNonZero]=value;
          index[numberNonZero++]=iColumn;
        }
      }
      } else {
      // scaled
      // modify pi so can collapse to one loop
      for (i=0;i<numberInRowArray;i++) {
        int iRow = whichRow[i];
        pi[iRow]=scalar*piOld[i]*rowScale[iRow];
      }
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
        if (backward_[iColumn]!=iSet) {
          // get pi on gub row
          iSet = backward_[iColumn];
          if (iSet>=0) {
            int iBasic = keyVariable_[iSet];
            if (iBasic<numberColumns) {
            // get dj without 
            assert (model->getStatus(iBasic)==ClpSimplex::basic);
            djMod=0.0;
            // scaled
            for (CoinBigIndex j=columnStart[iBasic];
                 j<columnStart[iBasic]+columnLength[iBasic];j++) {
              int jRow=row[j];
              djMod -= pi[jRow]*elementByColumn[j]*rowScale[jRow];
            }
            } else {
            djMod = 0.0;
            }
          } else {
            djMod = 0.0;
          }
        }
        double value = -djMod;
        CoinBigIndex j;
        const double * columnScale = model->columnScale();
        for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
          int iRow = row[j];
          value += pi[iRow]*elementByColumn[j];
        }
        value *= columnScale[iColumn];
        if (fabs(value)>zeroTolerance) {
          array[numberNonZero]=value;
          index[numberNonZero++]=iColumn;
        }
      }
      }
      // zero out
      for (i=0;i<numberInRowArray;i++) {
      int iRow = whichRow[i];
      pi[iRow]=0.0;
      }
    } else {
      // code later
      assert (packed);
      if (!rowScale) {
      if (scalar==-1.0) {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j];
          }
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=-value;
          }
        }
      } else if (scalar==1.0) {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j];
          }
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=value;
          }
        }
      } else {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j];
          }
          value *= scalar;
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=value;
          }
        }
      }
      } else {
      // scaled
      if (scalar==-1.0) {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          const double * columnScale = model->columnScale();
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
          }
          value *= columnScale[iColumn];
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=-value;
          }
        }
      } else if (scalar==1.0) {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          const double * columnScale = model->columnScale();
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
          }
          value *= columnScale[iColumn];
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=value;
          }
        }
      } else {
        for (iColumn=0;iColumn<numberColumns;iColumn++) {
          double value = 0.0;
          CoinBigIndex j;
          const double * columnScale = model->columnScale();
          for (j=columnStart[iColumn];
             j<columnStart[iColumn]+columnLength[iColumn];j++) {
            int iRow = row[j];
            value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
          }
          value *= scalar*columnScale[iColumn];
          if (fabs(value)>zeroTolerance) {
            index[numberNonZero++]=iColumn;
            array[iColumn]=value;
          }
        }
      }
      }
    }
    columnArray->setNumElements(numberNonZero);
    y->setNumElements(0);
  } else {
    // do by row
    transposeTimesByRow(model, scalar, rowArray, y, columnArray);
  }
  if (packed)
    columnArray->setPackedMode(true);
  if (0) {
    columnArray->checkClean();
    int numberNonZero=columnArray->getNumElements();;
    int * index = columnArray->getIndices();
    double * array = columnArray->denseVector();
    int i;
    for (i=0;i<numberNonZero;i++) {
      int j=index[i];
      double value;
      if (packed)
      value=array[i];
      else
      value=array[j];
      printf("Ti %d %d %g\n",i,j,value);
    }
  }
}
/* Return <code>x * A + y</code> in <code>z</code>. 
      Squashes small elements and knows about ClpSimplex */
void 
00791 ClpGubMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
                        const CoinIndexedVector * rowArray,
                        CoinIndexedVector * y,
                        CoinIndexedVector * columnArray) const
{
  // Do packed part
  ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray);
  if (numberSets_) {
    /* what we need to do is do by row as normal but get list of sets touched
       and then update those ones */
    abort();
  }
}
/* Return <code>x *A in <code>z</code> but
   just for indices in y. */
void 
00807 ClpGubMatrix::subsetTransposeTimes(const ClpSimplex * model,
                        const CoinIndexedVector * rowArray,
                        const CoinIndexedVector * y,
                        CoinIndexedVector * columnArray) const
{
  columnArray->clear();
  double * pi = rowArray->denseVector();
  double * array = columnArray->denseVector();
  int jColumn;
  // get matrix data pointers
  const int * row = matrix_->getIndices();
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * columnLength = matrix_->getVectorLengths(); 
  const double * elementByColumn = matrix_->getElements();
  const double * rowScale = model->rowScale();
  int numberToDo = y->getNumElements();
  const int * which = y->getIndices();
  assert (!rowArray->packedMode());
  columnArray->setPacked();
  int numberTouched=0;
  if (!rowScale) {
    for (jColumn=0;jColumn<numberToDo;jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j;
      for (j=columnStart[iColumn];
         j<columnStart[iColumn]+columnLength[iColumn];j++) {
      int iRow = row[j];
      value += pi[iRow]*elementByColumn[j];
      }
      array[jColumn]=value;
      if (value) {
      int iSet = backward_[iColumn];
      if (iSet>=0) {
        int iBasic = keyVariable_[iSet];
        if (iBasic==iColumn) {
          toIndex_[iSet]=jColumn;
          fromIndex_[numberTouched++]=iSet;
        }
      }
      }
    }
  } else {
    // scaled
    for (jColumn=0;jColumn<numberToDo;jColumn++) {
      int iColumn = which[jColumn];
      double value = 0.0;
      CoinBigIndex j;
      const double * columnScale = model->columnScale();
      for (j=columnStart[iColumn];
         j<columnStart[iColumn]+columnLength[iColumn];j++) {
      int iRow = row[j];
      value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
      }
      value *= columnScale[iColumn];
      array[jColumn]=value;
      if (value) {
      int iSet = backward_[iColumn];
      if (iSet>=0) {
        int iBasic = keyVariable_[iSet];
        if (iBasic==iColumn) {
          toIndex_[iSet]=jColumn;
          fromIndex_[numberTouched++]=iSet;
        }
      }
      }
    }
  }
  // adjust djs
  for (jColumn=0;jColumn<numberToDo;jColumn++) {
    int iColumn = which[jColumn];
    int iSet = backward_[iColumn];
    if (iSet>=0) {
      int kColumn = toIndex_[iSet];
      if (kColumn>=0)
      array[jColumn] -= array[kColumn];
    }
  }
  // and clear basic
  for (int j=0;j<numberTouched;j++) {
    int iSet = fromIndex_[j];
    int kColumn = toIndex_[iSet];
    toIndex_[iSet]=-1;
    array[kColumn]=0.0;
  }
}
/// returns number of elements in column part of basis,
CoinBigIndex 
00895 ClpGubMatrix::countBasis(const int * whichColumn, 
                   int & numberColumnBasic)
{
  int i;
  int numberColumns = getNumCols();
  const int * columnLength = matrix_->getVectorLengths(); 
  int numberRows = getNumRows();
  int numberBasic=0;
  CoinBigIndex numberElements=0;
  int lastSet=-1;
  int key=-1;
  int keyLength=-1;
  double * work = new double[numberRows];
  CoinZeroN(work,numberRows);
  char * mark = new char[numberRows];
  CoinZeroN(mark,numberRows);
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * row = matrix_->getIndices();
  const double * elementByColumn = matrix_->getElements();
  //ClpGubDynamicMatrix* gubx =
  //dynamic_cast< ClpGubDynamicMatrix*>(this);
  //int * id = gubx->id();
  // just count 
  for (i=0;i<numberColumnBasic;i++) {
    int iColumn = whichColumn[i];
    int iSet = backward_[iColumn];
    int length = columnLength[iColumn];
    if (iSet<0||keyVariable_[iSet]>=numberColumns) {
      numberElements += length;
      numberBasic++;
      //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length);
    } else {
      // in gub set
      if (iColumn!=keyVariable_[iSet]) {
      numberBasic++;
      CoinBigIndex j;
      // not key 
      if (lastSet<iSet) {
        // erase work
        if (key>=0) {
          for (j=columnStart[key];j<columnStart[key]+keyLength;j++)
            work[row[j]]=0.0;
        }
        key=keyVariable_[iSet];
        lastSet=iSet;
        keyLength = columnLength[key];
        for (j=columnStart[key];j<columnStart[key]+keyLength;j++)
          work[row[j]]=elementByColumn[j];
      }
      int extra=keyLength;
      for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
        int iRow = row[j];
        double keyValue = work[iRow];
        double value=elementByColumn[j];
        if (!keyValue) {
          if (fabs(value)>1.0e-20)
            extra++;
        } else {
          value -= keyValue;
          if (fabs(value)<=1.0e-20)
            extra--;
        }
      }
      numberElements+=extra;
        //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra);
      }
    }
  }
  delete [] work;
  delete [] mark;
  // update number of column basic
  numberColumnBasic = numberBasic;
  return numberElements;
}
void
00970 ClpGubMatrix::fillBasis(ClpSimplex * model,
                   const int * whichColumn, 
                   int & numberColumnBasic,
                   int * indexRowU, int * start,
                   int * rowCount, int * columnCount,
                   CoinFactorizationDouble * elementU)
{
  int i;
  int numberColumns = getNumCols();
  const int * columnLength = matrix_->getVectorLengths(); 
  int numberRows = getNumRows();
  assert (next_ ||!elementU) ;
  CoinBigIndex numberElements=start[0];
  int lastSet=-1;
  int key=-1;
  int keyLength=-1;
  double * work = new double[numberRows];
  CoinZeroN(work,numberRows);
  char * mark = new char[numberRows];
  CoinZeroN(mark,numberRows);
  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
  const int * row = matrix_->getIndices();
  const double * elementByColumn = matrix_->getElements();
  const double * rowScale = model->rowScale();
  int numberBasic=0;
  if (0) {
    printf("%d basiccolumns\n",numberColumnBasic);
    int i;
    for (i=0;i<numberSets_;i++) {
      int k=keyVariable_[i];
      if (k<numberColumns) {
      printf("key %d on set %d, %d elements\n",k,i,columnStart[k+1]-columnStart[k]);
      for (int j=columnStart[k];j<columnStart[k+1];j++)
        printf("row %d el %g\n",row[j],elementByColumn[j]);
      } else {
      printf("slack key on set %d\n",i);
      }
    }
  }
  // fill
  if (!rowScale) {
    // no scaling
    for (i=0;i<numberColumnBasic;i++) {
      int iColumn = whichColumn[i];
      int iSet = backward_[iColumn];
      int length = columnLength[iColumn];
      if (0) {
      int k=iColumn;
      printf("column %d in set %d, %d elements\n",k,iSet,columnStart[k+1]-columnStart[k]);
      for (int j=columnStart[k];j<columnStart[k+1];j++)
        printf("row %d el %g\n",row[j],elementByColumn[j]);
      }
      CoinBigIndex j;
      if (iSet<0||keyVariable_[iSet]>=numberColumns) {
      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
        double value = elementByColumn[j];
        if (fabs(value)>1.0e-20) {
          int iRow = row[j];
          indexRowU[numberElements]=iRow;
          rowCount[iRow]++;
          elementU[numberElements++]=value;
        }
      }
      // end of column
      columnCount[numberBasic]=numberElements-start[numberBasic];
      numberBasic++;
      start[numberBasic]=numberElements;
      } else {
      // in gub set
      if (iColumn!=keyVariable_[iSet]) {
        // not key 
        if (lastSet!=iSet) {
          // erase work
          if (key>=0) {
            for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
            int iRow=row[j];
            work[iRow]=0.0;
            mark[iRow]=0;
            }
          }
          key=keyVariable_[iSet];
          lastSet=iSet;
          keyLength = columnLength[key];
          for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
            int iRow=row[j];
            work[iRow]=elementByColumn[j];
            mark[iRow]=1;
          }
        }
        for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
          int iRow = row[j];
          double value=elementByColumn[j];
          if (mark[iRow]) {
            mark[iRow]=0;
            double keyValue = work[iRow];
            value -= keyValue;
          }
          if (fabs(value)>1.0e-20) {
            indexRowU[numberElements]=iRow;
            rowCount[iRow]++;
            elementU[numberElements++]=value;
          }
        }
        for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
          int iRow = row[j];
          if (mark[iRow]) {
            double value = -work[iRow];
            if (fabs(value)>1.0e-20) {
            indexRowU[numberElements]=iRow;
            rowCount[iRow]++;
            elementU[numberElements++]=value;
            }
          } else {
            // just put back mark
            mark[iRow]=1;
          }
        }
        // end of column
        columnCount[numberBasic]=numberElements-start[numberBasic];
        numberBasic++;
        start[numberBasic]=numberElements;
      }
      }
    }
  } else {
    // scaling
    const double * columnScale = model->columnScale();
    for (i=0;i<numberColumnBasic;i++) {
      int iColumn = whichColumn[i];
      int iSet = backward_[iColumn];
      int length = columnLength[iColumn];
      CoinBigIndex j;
      if (iSet<0||keyVariable_[iSet]>=numberColumns) {
      double scale = columnScale[iColumn];
      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
        int iRow = row[j];
        double value = elementByColumn[j]*scale*rowScale[iRow];
        if (fabs(value)>1.0e-20) {
          indexRowU[numberElements]=iRow;
          rowCount[iRow]++;
          elementU[numberElements++]=value;
        }
      }
      // end of column
      columnCount[numberBasic]=numberElements-start[numberBasic];
      numberBasic++;
      start[numberBasic]=numberElements;
      } else {
      // in gub set
      if (iColumn!=keyVariable_[iSet]) {
        double scale = columnScale[iColumn];
        // not key 
        if (lastSet<iSet) {
          // erase work
          if (key>=0) {
            for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
            int iRow=row[j];
            work[iRow]=0.0;
            mark[iRow]=0;
            }
          }
          key=keyVariable_[iSet];
          lastSet=iSet;
          keyLength = columnLength[key];
          double scale = columnScale[key];
          for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
            int iRow=row[j];
            work[iRow]=elementByColumn[j]*scale*rowScale[iRow];
            mark[iRow]=1;
          }
        }
        for (j=columnStart[iColumn];j<columnStart[iColumn]+length;j++) {
          int iRow = row[j];
          double value=elementByColumn[j]*scale*rowScale[iRow];
          if (mark[iRow]) {
            mark[iRow]=0;
            double keyValue = work[iRow];
            value -= keyValue;
          }
          if (fabs(value)>1.0e-20) {
            indexRowU[numberElements]=iRow;
            rowCount[iRow]++;
            elementU[numberElements++]=value;
          }
        }
        for (j=columnStart[key];j<columnStart[key]+keyLength;j++) {
          int iRow = row[j];
          if (mark[iRow]) {
            double value = -work[iRow];
            if (fabs(value)>1.0e-20) {
            indexRowU[numberElements]=iRow;
            rowCount[iRow]++;
            elementU[numberElements++]=value;
            }
          } else {
            // just put back mark
            mark[iRow]=1;
          }
        }
        // end of column
        columnCount[numberBasic]=numberElements-start[numberBasic];
        numberBasic++;
        start[numberBasic]=numberElements;
      }
      }
    }
  }
  delete [] work;
  delete [] mark;
  // update number of column basic
  numberColumnBasic = numberBasic;
}
/* Unpacks a column into an CoinIndexedvector
 */
void 
01185 ClpGubMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
               int iColumn) const 
{
  assert (iColumn<model->numberColumns());
  // Do packed part
  ClpPackedMatrix::unpack(model,rowArray,iColumn);
  int iSet = backward_[iColumn];
  if (iSet>=0) {
    int iBasic = keyVariable_[iSet];
    if (iBasic <model->numberColumns()) {
      add(model,rowArray,iBasic,-1.0);
    }
  }
}
/* Unpacks a column into a CoinIndexedVector
** in packed format
Note that model is NOT const.  Bounds and objective could
be modified if doing column generation (just for this variable) */
void 
01204 ClpGubMatrix::unpackPacked(ClpSimplex * model,
                      CoinIndexedVector * rowArray,
                      int iColumn) const
{
  int numberColumns = model->numberColumns();
  if (iColumn<numberColumns) {
    // Do packed part
    ClpPackedMatrix::unpackPacked(model,rowArray,iColumn);
    int iSet = backward_[iColumn];
    if (iSet>=0) {
      // columns are in order
      int iBasic = keyVariable_[iSet];
      if (iBasic <numberColumns) {
      int number = rowArray->getNumElements();
      const double * rowScale = model->rowScale();
      const int * row = matrix_->getIndices();
      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
      const int * columnLength = matrix_->getVectorLengths(); 
      const double * elementByColumn = matrix_->getElements();
      double * array = rowArray->denseVector();
      int * index = rowArray->getIndices();
      CoinBigIndex i;
      int numberOld=number;
      int lastIndex=0;
      int next=index[lastIndex];
      if (!rowScale) {
        for (i=columnStart[iBasic];
             i<columnStart[iBasic]+columnLength[iBasic];i++) {
          int iRow = row[i];
          while (iRow>next) {
            lastIndex++;
            if (lastIndex==numberOld)
            next=matrix_->getNumRows();
            else
            next=index[lastIndex];
          }
          if (iRow<next) {
            array[number]=-elementByColumn[i];
            index[number++]=iRow;
          } else {
            assert (iRow==next);
            array[lastIndex] -= elementByColumn[i];
            if (!array[lastIndex])
            array[lastIndex]=1.0e-100;
          }
        }
      } else {
        // apply scaling
        double scale = model->columnScale()[iBasic];
        for (i=columnStart[iBasic];
             i<columnStart[iBasic]+columnLength[iBasic];i++) {
          int iRow = row[i];
          while (iRow>next) {
            lastIndex++;
            if (lastIndex==numberOld)
            next=matrix_->getNumRows();
            else
            next=index[lastIndex];
          }
          if (iRow<next) {
            array[number]=-elementByColumn[i]*scale*rowScale[iRow];
            index[number++]=iRow;
          } else {
            assert (iRow==next);
            array[lastIndex] -=elementByColumn[i]*scale*rowScale[iRow];
            if (!array[lastIndex])
            array[lastIndex]=1.0e-100;
          }
        }
      }
      rowArray->setNumElements(number);
      }
    }
  } else {
    // key slack entering
    int iBasic = keyVariable_[gubSlackIn_];
    assert (iBasic <numberColumns);
    int number = 0;
    const double * rowScale = model->rowScale();
    const int * row = matrix_->getIndices();
    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
    const int * columnLength = matrix_->getVectorLengths(); 
    const double * elementByColumn = matrix_->getElements();
    double * array = rowArray->denseVector();
    int * index = rowArray->getIndices();
    CoinBigIndex i;
    if (!rowScale) {
      for (i=columnStart[iBasic];
         i<columnStart[iBasic]+columnLength[iBasic];i++) {
      int iRow = row[i];
      array[number]=elementByColumn[i];
      index[number++]=iRow;
      }
    } else {
      // apply scaling
      double scale = model->columnScale()[iBasic];
      for (i=columnStart[iBasic];
         i<columnStart[iBasic]+columnLength[iBasic];i++) {
      int iRow = row[i];
      array[number]=elementByColumn[i]*scale*rowScale[iRow];
      index[number++]=iRow;
      }
    }
    rowArray->setNumElements(number);
    rowArray->setPacked();
  }
}
/* Adds multiple of a column into an CoinIndexedvector
      You can use quickAdd to add to vector */
void 
01314 ClpGubMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
               int iColumn, double multiplier) const 
{
  assert (iColumn<model->numberColumns());
  // Do packed part
  ClpPackedMatrix::add(model,rowArray,iColumn,multiplier);
  int iSet = backward_[iColumn];
  if (iSet>=0&&iColumn!=keyVariable_[iSet]) {
    ClpPackedMatrix::add(model,rowArray,keyVariable_[iSet],-multiplier);
  }
}
/* Adds multiple of a column into an array */
void 
01327 ClpGubMatrix::add(const ClpSimplex * model,double * array,
               int iColumn, double multiplier) const 
{
  assert (iColumn<model->numberColumns());
  // Do packed part
  ClpPackedMatrix::add(model,array,iColumn,multiplier);
  if (iColumn<model->numberColumns()) {
    int iSet = backward_[iColumn];
    if (iSet>=0&&iColumn!=keyVariable_[iSet]&&keyVariable_[iSet]<model->numberColumns()) {
      ClpPackedMatrix::add(model,array,keyVariable_[iSet],-multiplier);
    }
  }
}
// Partial pricing 
void 
01342 ClpGubMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
                        int & bestSequence, int & numberWanted)
{
  numberWanted=currentWanted_;
  if (numberSets_) {
    // Do packed part before gub
    int numberColumns = matrix_->getNumCols();
    double ratio = static_cast<double> (firstGub_)/
      static_cast<double> (numberColumns);
    ClpPackedMatrix::partialPricing(model,startFraction*ratio,
                            endFraction*ratio,bestSequence,numberWanted);
    if (numberWanted||minimumGoodReducedCosts_<-1) {
      // do gub
      const double * element =matrix_->getElements();
      const int * row = matrix_->getIndices();
      const CoinBigIndex * startColumn = matrix_->getVectorStarts();
      const int * length = matrix_->getVectorLengths();
      const double * rowScale = model->rowScale();
      const double * columnScale = model->columnScale();
      int iSequence;
      CoinBigIndex j;
      double tolerance=model->currentDualTolerance();
      double * reducedCost = model->djRegion();
      const double * duals = model->dualRowSolution();
      const double * cost = model->costRegion();
      double bestDj;
      int numberColumns = model->numberColumns();
      int numberRows = model->numberRows();
      if (bestSequence>=0)
      bestDj = fabs(this->reducedCost(model,bestSequence));
      else
      bestDj=tolerance;
      int sequenceOut = model->sequenceOut();
      int saveSequence = bestSequence;
      int startG = firstGub_+ static_cast<int> (startFraction* (lastGub_-firstGub_));
      int endG = firstGub_+ static_cast<int> (endFraction* (lastGub_-firstGub_));
      endG = CoinMin(lastGub_,endG+1);
      // If nothing found yet can go all the way to end
      int endAll = endG;
      if (bestSequence<0&&!startG)
      endAll = lastGub_;
      int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; 
      int minNeg = minimumGoodReducedCosts_==-1 ? 5 : minimumGoodReducedCosts_;
      int nSets=0;
      int iSet = -1;
      double djMod=0.0;
      double infeasibilityCost = model->infeasibilityCost();
      if (rowScale) {
      double bestDjMod=0.0;
      // scaled
      for (iSequence=startG;iSequence<endAll;iSequence++) {
        if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
          // give up
          numberWanted=0;
          break;
        } else if (iSequence==endG&&bestSequence>=0) {
          break;
        }
        if (backward_[iSequence]!=iSet) {
          // get pi on gub row
          iSet = backward_[iSequence];
          if (iSet>=0) {
            nSets++;
            int iBasic = keyVariable_[iSet];
            if (iBasic>=numberColumns) {
            djMod = - weight(iSet)*infeasibilityCost;
            } else {
            // get dj without 
            assert (model->getStatus(iBasic)==ClpSimplex::basic);
            djMod=0.0;
            // scaled
            for (j=startColumn[iBasic];
                 j<startColumn[iBasic]+length[iBasic];j++) {
              int jRow=row[j];
              djMod -= duals[jRow]*element[j]*rowScale[jRow];
            }
            // allow for scaling
            djMod +=  cost[iBasic]/columnScale[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 = numberRows+numberColumns+iSet;
                  bestDjMod = djMod;
                  } 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 = numberRows+numberColumns+iSet;
                  bestDjMod = djMod;
                  } else {
                  // just to make sure we don't exit before got something
                  numberWanted++;
                  abort();
                  }
                }
              }
            }
            }
          } else {
            // not in set
            djMod=0.0;
          }
        }
        if (iSequence!=sequenceOut) {
          double value;
          ClpSimplex::Status status = model->getStatus(iSequence);
          
          switch(status) {
            
          case ClpSimplex::basic:
          case ClpSimplex::isFixed:
            break;
          case ClpSimplex::isFree:
          case ClpSimplex::superBasic:
            value=-djMod;
            // scaled
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j]*rowScale[jRow];
            }
            value = fabs(cost[iSequence] +value*columnScale[iSequence]);
            if (value>FREE_ACCEPT*tolerance) {
            numberWanted--;
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          case ClpSimplex::atUpperBound:
            value=-djMod;
            // scaled
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j]*rowScale[jRow];
            }
            value = cost[iSequence] +value*columnScale[iSequence];
            if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          case ClpSimplex::atLowerBound:
            value=-djMod;
            // scaled
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j]*rowScale[jRow];
            }
            value = -(cost[iSequence] +value*columnScale[iSequence]);
            if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          }
        }
        if (!numberWanted)
          break;
      }
      if (bestSequence!=saveSequence) {
        if (bestSequence<numberRows+numberColumns) {
          // recompute dj
          double value=bestDjMod;
          // scaled
          for (j=startColumn[bestSequence];
             j<startColumn[bestSequence]+length[bestSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j]*rowScale[jRow];
          }
          reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
          gubSlackIn_=-1;
        } else {
          // slack - make last column
          gubSlackIn_= bestSequence - numberRows - numberColumns;
          bestSequence = numberColumns + 2*numberRows;
          reducedCost[bestSequence] = bestDjMod;
          model->setStatus(bestSequence,getStatus(gubSlackIn_));
          if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
            model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
          else
            model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
          model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
          model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
          model->costRegion()[bestSequence] = 0.0;
        }
        savedBestSequence_ = bestSequence;
        savedBestDj_ = reducedCost[savedBestSequence_];
      }
      } else {
      double bestDjMod=0.0;
      //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
      //     startG,endG,numberWanted);
      for (iSequence=startG;iSequence<endG;iSequence++) {
        if (numberWanted+minNeg<originalWanted_&&nSets>minSet) {
          // give up
          numberWanted=0;
          break;
        } else if (iSequence==endG&&bestSequence>=0) {
          break;
        }
        if (backward_[iSequence]!=iSet) {
          // get pi on gub row
          iSet = backward_[iSequence];
          if (iSet>=0) {
            nSets++;
            int iBasic = keyVariable_[iSet];
            if (iBasic>=numberColumns) {
            djMod = - weight(iSet)*infeasibilityCost;
            } else {
            // get dj without 
            assert (model->getStatus(iBasic)==ClpSimplex::basic);
            djMod=0.0;

            for (j=startColumn[iBasic];
                 j<startColumn[iBasic]+length[iBasic];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 = numberRows+numberColumns+iSet;
                  bestDjMod = djMod;
                  } 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 = numberRows+numberColumns+iSet;
                  bestDjMod = djMod;
                  } else {
                  // just to make sure we don't exit before got something
                  numberWanted++;
                  abort();
                  }
                }
              }
            }
            }
          } else {
            // not in set
            djMod=0.0;
          }
        }
        if (iSequence!=sequenceOut) {
          double value;
          ClpSimplex::Status status = model->getStatus(iSequence);
          
          switch(status) {
            
          case ClpSimplex::basic:
          case ClpSimplex::isFixed:
            break;
          case ClpSimplex::isFree:
          case ClpSimplex::superBasic:
            value=cost[iSequence]-djMod;
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j];
            }
            value = fabs(value);
            if (value>FREE_ACCEPT*tolerance) {
            numberWanted--;
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          case ClpSimplex::atUpperBound:
            value=cost[iSequence]-djMod;
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j];
            }
            if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          case ClpSimplex::atLowerBound:
            value=cost[iSequence]-djMod;
            for (j=startColumn[iSequence];
               j<startColumn[iSequence]+length[iSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j];
            }
            value = -value;
            if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
              // check flagged variable and correct dj
              if (!model->flagged(iSequence)) {
                bestDj=value;
                bestSequence = iSequence;
                bestDjMod = djMod;
              } else {
                // just to make sure we don't exit before got something
                numberWanted++;
              }
            }
            }
            break;
          }
        }
        if (!numberWanted)
          break;
      }
      if (bestSequence!=saveSequence) {
        if (bestSequence<numberRows+numberColumns) {
          // recompute dj
          double value=cost[bestSequence]-bestDjMod;
          for (j=startColumn[bestSequence];
             j<startColumn[bestSequence]+length[bestSequence];j++) {
            int jRow=row[j];
            value -= duals[jRow]*element[j];
          }
          //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
          reducedCost[bestSequence] = value;
          gubSlackIn_=-1;
        } else {
          // slack - make last column
          gubSlackIn_= bestSequence - numberRows - numberColumns;
          bestSequence = numberColumns + 2*numberRows;
          reducedCost[bestSequence] = bestDjMod;
          //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
          model->setStatus(bestSequence,getStatus(gubSlackIn_));
          if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
            model->solutionRegion()[bestSequence] = upper_[gubSlackIn_];
          else
            model->solutionRegion()[bestSequence] = lower_[gubSlackIn_];
          model->lowerRegion()[bestSequence] = lower_[gubSlackIn_];
          model->upperRegion()[bestSequence] = upper_[gubSlackIn_];
          model->costRegion()[bestSequence] = 0.0;
        }
      }
      }
      // See if may be finished
      if (startG==firstGub_&&bestSequence<0)
      infeasibilityWeight_=model_->infeasibilityCost();
      else if (bestSequence>=0) 
      infeasibilityWeight_=-1.0;
    }
    if (numberWanted) {
      // Do packed part after gub
      double offset = static_cast<double> (lastGub_)/
      static_cast<double> (numberColumns); 
      double ratio = static_cast<double> (numberColumns)/
      static_cast<double> (numberColumns)-offset;
      double start2 = offset + ratio*startFraction;
      double end2 = CoinMin(1.0,offset + ratio*endFraction+1.0e-6);
      ClpPackedMatrix::partialPricing(model,start2,end2,bestSequence,numberWanted);
    }
  } else {
    // no gub
    ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
  }
  if (bestSequence>=0)
    infeasibilityWeight_=-1.0; // not optimal
  currentWanted_=numberWanted;
}
/* expands an updated column to allow for extra rows which the main
   solver does not know about and returns number added. 
*/
int 
01795 ClpGubMatrix::extendUpdated(ClpSimplex * model,CoinIndexedVector * update, int mode)
{
  // I think we only need to bother about sets with two in basis or incoming set
  int number = update->getNumElements();
  double * array = update->denseVector();
  int * index = update->getIndices();
  int i;
  assert (!number||update->packedMode());
  int * pivotVariable = model->pivotVariable();
  int numberRows = model->numberRows();
  int numberColumns = model->numberColumns();
  int numberTotal = numberRows+numberColumns;
  int sequenceIn = model->sequenceIn();
  int returnCode=0;
  int iSetIn;
  if (sequenceIn<numberColumns) {
    iSetIn = backward_[sequenceIn];
    gubSlackIn_ = -1; // in case set
  } else if (sequenceIn<numberRows+numberColumns) {
    iSetIn = -1;
    gubSlackIn_ = -1; // in case set
  } else {
    iSetIn = gubSlackIn_;
  }
  double * lower = model->lowerRegion();
  double * upper = model->upperRegion();
  double * cost = model->costRegion();
  double * solution = model->solutionRegion();
  int number2=number;
  if (!mode) {
    double primalTolerance = model->primalTolerance();
    double infeasibilityCost = model->infeasibilityCost();
    // extend
    saveNumber_ = number;
    for (i=0;i<number;i++) {
      int iRow = index[i];
      int iPivot = pivotVariable[iRow];
      if (iPivot<numberColumns) {
      int iSet = backward_[iPivot];
      if (iSet>=0) {
        // two (or more) in set
        int iIndex =toIndex_[iSet];
        double otherValue=array[i];
        double value;
        if (iIndex<0) {
          toIndex_[iSet]=number2;
          int iNew = number2-number;
          fromIndex_[number2-number]=iSet;
          iIndex=number2;
          index[number2]=numberRows+iNew;
          // do key stuff
          int iKey = keyVariable_[iSet];
          if (iKey<numberColumns) {
            // Save current cost of key
            changeCost_[number2-number] = cost[iKey];
            if (iSet!=iSetIn)
            value = 0.0;
            else if (iSetIn!=gubSlackIn_)
            value = 1.0;
            else
            value =-1.0;
            pivotVariable[numberRows+iNew]=iKey;
            // Do I need to recompute?
            double sol;
            assert (getStatus(iSet)!=ClpSimplex::basic);
            if (getStatus(iSet)==ClpSimplex::atLowerBound)
            sol = lower_[iSet];
            else
            sol = upper_[iSet];
            if ((gubType_&8)!=0) {
            int iColumn =next_[iKey];
            // sum all non-key variables
            while(iColumn>=0) {
              sol -= solution[iColumn];
              iColumn=next_[iColumn];
            }
            } else {
            int stop = -(iKey+1);
            int iColumn =next_[iKey];
            // sum all non-key variables
            while(iColumn!=stop) {
              if (iColumn<0)
                iColumn = -iColumn-1;
              sol -= solution[iColumn];
              iColumn=next_[iColumn];
            }
            }
            solution[iKey]=sol;
            if (model->algorithm()>0)
            model->nonLinearCost()->setOne(iKey,sol);
            //assert (fabs(sol-solution[iKey])<1.0e-3);
          } else {
            // gub slack is basic
            // Save current cost of key
            changeCost_[number2-number]= -weight(iSet)*infeasibilityCost;
            otherValue = - otherValue; //allow for - sign on slack
            if (iSet!=iSetIn)
            value = 0.0;
            else
            value = -1.0;
            pivotVariable[numberRows+iNew]=iNew+numberTotal;
            model->djRegion()[iNew+numberTotal]=0.0;
            double sol=0.0;
            if ((gubType_&8)!=0) {
            int iColumn =next_[iKey];
            // sum all non-key variables
            while(iColumn>=0) {
              sol += solution[iColumn];
              iColumn=next_[iColumn];
            }
            } else {
            int stop = -(iKey+1);
            int iColumn =next_[iKey];
            // sum all non-key variables
            while(iColumn!=stop) {
              if (iColumn<0)
                iColumn = -iColumn-1;
              sol += solution[iColumn];
              iColumn=next_[iColumn];
            }
            }
            solution[iNew+numberTotal]=sol;
            // and do cost in nonLinearCost
            if (model->algorithm()>0)
            model->nonLinearCost()->setOne(iNew+numberTotal,sol,lower_[iSet],upper_[iSet]);
            if (sol>upper_[iSet]+primalTolerance) {
            setAbove(iSet);
            lower[iNew+numberTotal]=upper_[iSet];
            upper[iNew+numberTotal]=COIN_DBL_MAX;
            } else if (sol<lower_[iSet]-primalTolerance) {
            setBelow(iSet);
            lower[iNew+numberTotal]=-COIN_DBL_MAX;
            upper[iNew+numberTotal]=lower_[iSet];
            } else {
            setFeasible(iSet);
            lower[iNew+numberTotal]=lower_[iSet];
            upper[iNew+numberTotal]=upper_[iSet];
            }
            cost[iNew+numberTotal]=weight(iSet)*infeasibilityCost;
          }
          number2++;
        } else {
          value = array[iIndex];
          int iKey = keyVariable_[iSet];
          if (iKey>=numberColumns) 
            otherValue = - otherValue; //allow for - sign on slack
        }
        value -= otherValue;
        array[iIndex]=value;
      }
      }
    }
    if (iSetIn>=0&&toIndex_[iSetIn]<0) {
      // Do incoming
      update->setPacked(); // just in case no elements
      toIndex_[iSetIn]=number2;
      int iNew = number2-number;
      fromIndex_[number2-number]=iSetIn;
      // Save current cost of key
      double currentCost;
      int key=keyVariable_[iSetIn];
      if (key<numberColumns) 
      currentCost = cost[key];
      else
      currentCost = -weight(iSetIn)*infeasibilityCost;
      changeCost_[number2-number]=currentCost;
      index[number2]=numberRows+iNew;
      // do key stuff
      int iKey = keyVariable_[iSetIn];
      if (iKey<numberColumns) {
      if (gubSlackIn_<0)
        array[number2]=1.0;
      else
        array[number2]=-1.0;
      pivotVariable[numberRows+iNew]=iKey;
      // Do I need to recompute?
      double sol;
      assert (getStatus(iSetIn)!=ClpSimplex::basic);
      if (getStatus(iSetIn)==ClpSimplex::atLowerBound)
        sol = lower_[iSetIn];
      else
        sol = upper_[iSetIn];
      if ((gubType_&8)!=0) {
        int iColumn =next_[iKey];
        // sum all non-key variables
        while(iColumn>=0) {
          sol -= solution[iColumn];
          iColumn=next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(iKey+1);
        int iColumn =next_[iKey];
        // sum all non-key variables
        while(iColumn!=stop) {
          if (iColumn<0)
            iColumn = -iColumn-1;
          sol -= solution[iColumn];
          iColumn=next_[iColumn];
        }
      }
      solution[iKey]=sol;
      if (model->algorithm()>0)
        model->nonLinearCost()->setOne(iKey,sol);
      //assert (fabs(sol-solution[iKey])<1.0e-3);
      } else {
      // gub slack is basic
      array[number2]=-1.0;
      pivotVariable[numberRows+iNew]=iNew+numberTotal;
      model->djRegion()[iNew+numberTotal]=0.0;
      double sol=0.0;
      if ((gubType_&8)!=0) {
        int iColumn =next_[iKey];
        // sum all non-key variables
        while(iColumn>=0) {
          sol += solution[iColumn];
          iColumn=next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(iKey+1);
        int iColumn =next_[iKey];
        // sum all non-key variables
        while(iColumn!=stop) {
          if (iColumn<0)
            iColumn = -iColumn-1;
          sol += solution[iColumn];
          iColumn=next_[iColumn];
        }
      }
      solution[iNew+numberTotal]=sol;
      // and do cost in nonLinearCost
      if (model->algorithm()>0)
        model->nonLinearCost()->setOne(iNew+numberTotal,sol,lower_[iSetIn],upper_[iSetIn]);
      if (sol>upper_[iSetIn]+primalTolerance) {
        setAbove(iSetIn);
        lower[iNew+numberTotal]=upper_[iSetIn];
        upper[iNew+numberTotal]=COIN_DBL_MAX;
      } else if (sol<lower_[iSetIn]-primalTolerance) {
        setBelow(iSetIn);
        lower[iNew+numberTotal]=-COIN_DBL_MAX;
        upper[iNew+numberTotal]=lower_[iSetIn];
      } else {
        setFeasible(iSetIn);
        lower[iNew+numberTotal]=lower_[iSetIn];
        upper[iNew+numberTotal]=upper_[iSetIn];
      }
      cost[iNew+numberTotal]=weight(iSetIn)*infeasibilityCost;
      }
      number2++;
    }
    // mark end
    fromIndex_[number2-number]=-1;
    returnCode = number2-number;
    // make sure lower_ upper_ adjusted
    synchronize(model,9);
  } else {
    // take off?
    if (number>saveNumber_) {
      // clear
      double theta = model->theta();
      double * solution = model->solutionRegion();
      for (i=saveNumber_;i<number;i++) {
      int iRow = index[i];
      int iColumn = pivotVariable[iRow];
#ifdef CLP_DEBUG_PRINT
      printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n",
             iColumn,fromIndex_[i-saveNumber_],lower[iColumn],upper[iColumn],array[i],
             solution[iColumn],solution[iColumn]-model->theta()*array[i],model->theta());
#endif
      double value = array[i];
      array[i]=0.0;
      int iSet=fromIndex_[i-saveNumber_];
      toIndex_[iSet]=-1;
      if (iSet==iSetIn&&iColumn<numberColumns) {
        // update as may need value
        solution[iColumn] -= theta*value;
      }
      }
    }
#ifdef CLP_DEBUG
    for (i=0;i<numberSets_;i++)
      assert(toIndex_[i]==-1);
#endif
    number2= saveNumber_;
  }
  update->setNumElements(number2);
  return returnCode;
}
/*
     utility primal function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void 
02090 ClpGubMatrix::primalExpanded(ClpSimplex * model,int mode)
{
  int numberColumns = model->numberColumns();
  switch (mode) {
    // If key variable then slot in gub rhs so will get correct contribution
  case 0:
    {
      int i;
      double * solution = model->solutionRegion();
      ClpSimplex::Status iStatus;
      for (i=0;i<numberSets_;i++) {
      int iColumn = keyVariable_[i];
      if (iColumn<numberColumns) {
        // key is structural - where is slack
        iStatus = getStatus(i);
        assert (iStatus!=ClpSimplex::basic);
        if (iStatus==ClpSimplex::atLowerBound)
          solution[iColumn]=lower_[i];
        else
          solution[iColumn]=upper_[i];
      }
      }
    }
    break;
    // Compute values of key variables
  case 1:
    {
      int i;
      double * solution = model->solutionRegion();
      ClpSimplex::Status iStatus;
      //const int * columnLength = matrix_->getVectorLengths(); 
      //const CoinBigIndex * columnStart = matrix_->getVectorStarts();
      //const int * row = matrix_->getIndices();
      //const double * elementByColumn = matrix_->getElements();
      //int * pivotVariable = model->pivotVariable();
      sumPrimalInfeasibilities_=0.0;
      numberPrimalInfeasibilities_=0;
      double primalTolerance = model->primalTolerance();
      double relaxedTolerance=primalTolerance;
      // we can't really trust infeasibilities if there is primal error
      double error = CoinMin(1.0e-2,model->largestPrimalError());
      // allow tolerance at least slightly bigger than standard
      relaxedTolerance = relaxedTolerance +  error;
      // but we will be using difference
      relaxedTolerance -= primalTolerance;
      sumOfRelaxedPrimalInfeasibilities_ = 0.0;
      for (i=0;i<numberSets_;i++) { // Could just be over basics (esp if no bounds)
      int kColumn = keyVariable_[i];
      double value=0.0;
      if ((gubType_&8)!=0) {
        int iColumn =next_[kColumn];
        // sum all non-key variables
        while(iColumn>=0) {
          value+=solution[iColumn];
          iColumn=next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(kColumn+1);
        int iColumn =next_[kColumn];
        // sum all non-key variables
        while(iColumn!=stop) {
          if (iColumn<0)
            iColumn = -iColumn-1;
          value += solution[iColumn];
          iColumn=next_[iColumn];
        }
      }
      if (kColumn<numberColumns) {
        // make sure key is basic - so will be skipped in values pass
        model->setStatus(kColumn,ClpSimplex::basic);
        // feasibility will be done later
        assert (getStatus(i)!=ClpSimplex::basic);
        if (getStatus(i)==ClpSimplex::atUpperBound)
          solution[kColumn] = upper_[i]-value;
        else
          solution[kColumn] = lower_[i]-value;
        //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]);
      } else {
        // slack is key
        iStatus = getStatus(i);
        assert (iStatus==ClpSimplex::basic);
        double infeasibility=0.0;
        if (value>upper_[i]+primalTolerance) {
          infeasibility=value-upper_[i]-primalTolerance;
          setAbove(i);
        } else if (value<lower_[i]-primalTolerance) {
          infeasibility=lower_[i]-value-primalTolerance ;
          setBelow(i);
        } else {
          setFeasible(i);
        }
        //printf("Value of key slack for set %d is %g\n",i,value);
        if (infeasibility>0.0) {
          sumPrimalInfeasibilities_ += infeasibility;
          if (infeasibility>relaxedTolerance) 
            sumOfRelaxedPrimalInfeasibilities_ += infeasibility;
          numberPrimalInfeasibilities_ ++;
        }
      }
      }
    }
    break;
    // Report on infeasibilities of key variables
  case 2:
    {
      model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities()+
                               sumPrimalInfeasibilities_);
      model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities()+
                               numberPrimalInfeasibilities_);
      model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities()+
                               sumOfRelaxedPrimalInfeasibilities_);
    }
    break;
  }
}
/*
     utility dual function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
void 
02212 ClpGubMatrix::dualExpanded(ClpSimplex * model,
                      CoinIndexedVector * array,
                     double * /*other*/,int mode)
{
  switch (mode) {
    // modify costs before transposeUpdate
  case 0:
    {
      int i;
      double * cost = model->costRegion();
      ClpSimplex::Status iStatus;
      // not dual values yet
      //assert (!other);
      //double * work = array->denseVector();
      double infeasibilityCost = model->infeasibilityCost();
      int * pivotVariable = model->pivotVariable();
      int numberRows = model->numberRows();
      int numberColumns = model->numberColumns();
      for (i=0;i<numberRows;i++) {
      int iPivot = pivotVariable[i];
      if (iPivot<numberColumns) {
        int iSet = backward_[iPivot];
        if (iSet>=0) {
          int kColumn = keyVariable_[iSet];
          double costValue;
          if (kColumn<numberColumns) {
            // structural has cost
            costValue = cost[kColumn];
          } else {
            // slack is key
            iStatus = getStatus(iSet);
            assert (iStatus==ClpSimplex::basic);
            // negative as -1.0 for slack
            costValue=-weight(iSet)*infeasibilityCost;
          }
          array->add(i,-costValue); // was work[i]-costValue
        }
      }
      }
    }
    break;
    // create duals for key variables (without check on dual infeasible)
  case 1:
    {
      // If key slack then dual 0.0 (if feasible)
      // dj for key is zero so that defines dual on set
      int i;
      double * dj = model->djRegion();
      int numberColumns = model->numberColumns();
      double infeasibilityCost = model->infeasibilityCost();
      for (i=0;i<numberSets_;i++) {
      int kColumn = keyVariable_[i];
      if (kColumn<numberColumns) {
        // dj without set
        double value = dj[kColumn];
        // Now subtract out from all 
        dj[kColumn] =0.0;
        int iColumn =next_[kColumn];
        // modify all non-key variables
        while(iColumn>=0) {
          dj[iColumn]-=value;
          iColumn=next_[iColumn];
        }
      } else {
        // slack key - may not be feasible
        assert (getStatus(i)==ClpSimplex::basic);
        // negative as -1.0 for slack
        double value=-weight(i)*infeasibilityCost;
        if (value) {
          int iColumn =next_[kColumn];
          // modify all non-key variables basic
          while(iColumn>=0) {
            dj[iColumn]-=value;
            iColumn=next_[iColumn];
          }
        }
      }
      }
    }
    break;
    // as 1 but check slacks and compute djs
  case 2:
    {
      // 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;
      // make sure fromIndex will not confuse pricing
      fromIndex_[0]=-1;
      possiblePivotKey_=-1;
      // Create array
      int numberColumns = model->numberColumns();
      int * pivotVariable = model->pivotVariable();
      int numberRows = model->numberRows();
      for (i=0;i<numberRows;i++) {
      int iPivot = pivotVariable[i];
      if (iPivot<numberColumns)
        backToPivotRow_[iPivot]=i;
      }
      if (noCheck_>=0) {
      if (infeasibilityWeight_!=model->infeasibilityCost()) {
        // don't bother checking
        sumDualInfeasibilities_=100.0;
        numberDualInfeasibilities_=1;
        sumOfRelaxedDualInfeasibilities_ =100.0;
        return;
      }
      }
      double * dj = model->djRegion();
      double * dual = model->dualRowSolution();
      double * cost = model->costRegion();
      ClpSimplex::Status iStatus;
      const int * columnLength = matrix_->getVectorLengths(); 
      const CoinBigIndex * columnStart = matrix_->getVectorStarts();
      const int * row = matrix_->getIndices();
      const double * elementByColumn = matrix_->getElements();
      double infeasibilityCost = model->infeasibilityCost();
      sumDualInfeasibilities_=0.0;
      numberDualInfeasibilities_=0;
      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;
      sumOfRelaxedDualInfeasibilities_ = 0.0;
      for (i=0;i<numberSets_;i++) {
      int kColumn = keyVariable_[i];
      if (kColumn<numberColumns) {
        // dj without set
        double value = cost[kColumn];
        for (CoinBigIndex j=columnStart[kColumn];
             j<columnStart[kColumn]+columnLength[kColumn];j++) {
          int iRow = row[j];
          value -= dual[iRow]*elementByColumn[j];
        }
        // Now subtract out from all 
        dj[kColumn] -= value;
        int stop = -(kColumn+1);
        kColumn = next_[kColumn];
        while (kColumn!=stop) {
          if (kColumn<0)
            kColumn = -kColumn-1;
          double djValue = dj[kColumn]-value;
          dj[kColumn] = djValue;;
          double infeasibility=0.0;
          iStatus = model->getStatus(kColumn);
          if (iStatus==ClpSimplex::atLowerBound) {
            if (djValue<-dualTolerance) 
            infeasibility=-djValue-dualTolerance;
          } else if (iStatus==ClpSimplex::atUpperBound) {
            // at upper bound
            if (djValue>dualTolerance) 
            infeasibility=djValue-dualTolerance;
          }
          if (infeasibility>0.0) {
            sumDualInfeasibilities_ += infeasibility;
            if (infeasibility>relaxedTolerance) 
            sumOfRelaxedDualInfeasibilities_ += infeasibility;
            numberDualInfeasibilities_ ++;
          }
          kColumn = next_[kColumn];
        }
        // check slack
        iStatus = getStatus(i);
        assert (iStatus!=ClpSimplex::basic);
        double infeasibility=0.0;
        // dj of slack is -(-1.0)value
        if (iStatus==ClpSimplex::atLowerBound) {
          if (value<-dualTolerance) 
            infeasibility=-value-dualTolerance;
        } else if (iStatus==ClpSimplex::atUpperBound) {
          // at upper bound
          if (value>dualTolerance) 
            infeasibility=value-dualTolerance;
        }
        if (infeasibility>0.0) {
          sumDualInfeasibilities_ += infeasibility;
          if (infeasibility>relaxedTolerance) 
            sumOfRelaxedDualInfeasibilities_ += infeasibility;
          numberDualInfeasibilities_ ++;
        }
      } else {
        // slack key - may not be feasible
        assert (getStatus(i)==ClpSimplex::basic);
        // negative as -1.0 for slack
        double value=-weight(i)*infeasibilityCost;
        if (value) {
        // Now subtract out from all 
          int kColumn = i+numberColumns;
          int stop = -(kColumn+1);
          kColumn = next_[kColumn];
          while (kColumn!=stop) {
            if (kColumn<0)
            kColumn = -kColumn-1;
            double djValue = dj[kColumn]-value;
            dj[kColumn] = djValue;;
            double infeasibility=0.0;
            iStatus = model->getStatus(kColumn);
            if (iStatus==ClpSimplex::atLowerBound) {
            if (djValue<-dualTolerance) 
              infeasibility=-djValue-dualTolerance;
            } else if (iStatus==ClpSimplex::atUpperBound) {
            // at upper bound
            if (djValue>dualTolerance) 
              infeasibility=djValue-dualTolerance;
            }
            if (infeasibility>0.0) {
            sumDualInfeasibilities_ += infeasibility;
            if (infeasibility>relaxedTolerance) 
              sumOfRelaxedDualInfeasibilities_ += infeasibility;
            numberDualInfeasibilities_ ++;
            }
            kColumn = next_[kColumn];
          }
        }
      }
      }
      // and get statistics for column generation
      synchronize(model,4);
      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:
    {
      // First compute new costs etc for interesting gubs
      int iLook=0;
      int iSet=fromIndex_[0];
      double primalTolerance = model->primalTolerance();
      const double * cost = model->costRegion();
      double * solution = model->solutionRegion();
      double infeasibilityCost = model->infeasibilityCost();
      int numberColumns = model->numberColumns();
      int numberChanged=0;
      int * pivotVariable = model->pivotVariable();
      while (iSet>=0) {
      int key=keyVariable_[iSet];
      double value=0.0;
      // sum over all except key
      if ((gubType_&8)!=0) {
        int iColumn =next_[key];
        // sum all non-key variables
        while(iColumn>=0) {
          value += solution[iColumn];
          iColumn=next_[iColumn];
        }
      } else {
        // bounds exist - sum over all except key
        int stop = -(key+1);
        int iColumn =next_[key];
        // sum all non-key variables
        while(iColumn!=stop) {
          if (iColumn<0)
            iColumn = -iColumn-1;
          value += solution[iColumn];
          iColumn=next_[iColumn];
        }
      }
      double costChange;
      double oldCost = changeCost_[iLook];
      if (key<numberColumns) {
        assert (getStatus(iSet)!=ClpSimplex::basic);
        double sol;
        if (getStatus(iSet)==ClpSimplex::atUpperBound)
          sol = upper_[iSet]-value;
        else
          sol = lower_[iSet]-value;
        solution[key]=sol;
        // fix up cost
        model->nonLinearCost()->setOne(key,sol);
#ifdef CLP_DEBUG_PRINT
        printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n",key,iSet,sol,
             cost[key],oldCost);
#endif
        costChange = cost[key]-oldCost;
      } else {
        // slack is key
        if (value>upper_[iSet]+primalTolerance) {
          setAbove(iSet);
        } else if (value<lower_[iSet]-primalTolerance) {
          setBelow(iSet);
        } else {
          setFeasible(iSet);
        }
        // negative as -1.0 for slack
        costChange=-weight(iSet)*infeasibilityCost-oldCost;
#ifdef CLP_DEBUG_PRINT
        printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n",iSet,value,
             weight(iSet)*infeasibilityCost,oldCost);
#endif
      }
      if (costChange) {
        fromIndex_[numberChanged]=iSet;
        toIndex_[iSet]=numberChanged;
        changeCost_[numberChanged++]=costChange;
      }
      iSet = fromIndex_[++iLook];
      }
      if (numberChanged||possiblePivotKey_>=0) {
      // first do those in list already
      int number = array->getNumElements();
      array->setPacked();
      int i;
      double * work = array->denseVector();
      int * which = array->getIndices();
      for (i=0;i<number;i++) {
        int iRow = which[i];
        int iPivot = pivotVariable[iRow];
        if (iPivot<numberColumns) {
          int iSet = backward_[iPivot];
          if (iSet>=0&&toIndex_[iSet]>=0) {
            double newValue = work[i]+changeCost_[toIndex_[iSet]];
            if (!newValue)
            newValue =1.0e-100;
            work[i]=newValue;
            // mark as done
            backward_[iPivot]=-1;
          }
        }
        if (possiblePivotKey_==iRow) {
          double newValue = work[i]-model->dualIn();
          if (!newValue)
            newValue =1.0e-100;
          work[i]=newValue;
          possiblePivotKey_=-1;
        }
      }
      // now do rest and clean up
      for (i=0;i<numberChanged;i++) {
        int iSet = fromIndex_[i];
        int key=keyVariable_[iSet];
        int iColumn = next_[key];
        double change = changeCost_[i];
        while (iColumn>=0) {
          if (backward_[iColumn]>=0) {
            int iRow = backToPivotRow_[iColumn];
            assert (iRow>=0);
            work[number] = change;
            if (possiblePivotKey_==iRow) {
            double newValue = work[number]-model->dualIn();
            if (!newValue)
              newValue =1.0e-100;
            work[number]=newValue;
            possiblePivotKey_=-1;
            }
            which[number++]=iRow;
          } else {
            // reset
            backward_[iColumn]=iSet;
          }
          iColumn=next_[iColumn];
        }
        toIndex_[iSet]=-1;
      }
      if (possiblePivotKey_>=0) {
        work[number] = -model->dualIn();
        which[number++]=possiblePivotKey_;
        possiblePivotKey_=-1;
      }
      fromIndex_[0]=-1;
      array->setNumElements(number);
      }
    }
    break;
  }
}
// This is local to Gub to allow synchronization when status is good
int 
02594 ClpGubMatrix::synchronize(ClpSimplex *, int)
{
  return 0;
}
/*
     general utility function for dealing with dynamic constraints
     mode=n see ClpGubMatrix.hpp for definition
     Remember to update here when settled down
*/
int
02604 ClpGubMatrix::generalExpanded(ClpSimplex * model,int mode,int &number)
{
  int returnCode=0;
  int numberColumns = model->numberColumns();
  switch (mode) {
    // Fill in pivotVariable but not for key variables
  case 0:
    {
      if (!next_ ) {
      // do ordering
      assert (!rhsOffset_);
      // create and do gub crash
      useEffectiveRhs(model,false);
      }
      int i;
      int numberBasic=number;
      // 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) {
        int iSet = backward_[i];
        if (iSet<0||i!=keyVariable_[iSet])
          pivotVariable[numberBasic++]=i;
      }
      }
      number = numberBasic;
      if (model->numberIterations())
      assert (number==model->numberRows());
    }
    break;
    // Make all key variables basic
  case 1:
    {
      int i;
      for (i=0;i<numberSets_;i++) {
      int iColumn = keyVariable_[i];
      if (iColumn<numberColumns)
        model->setColumnStatus(iColumn,ClpSimplex::basic);
      }
    }
    break;
    // Do initial extra rows + maximum basic
  case 2:
    {
      returnCode= getNumRows()+1;
      number = model->numberRows()+numberSets_;
    }
    break;
    // Before normal replaceColumn
  case 3:
    {
      int sequenceIn = model->sequenceIn();
      int sequenceOut = model->sequenceOut();
      int numberColumns = model->numberColumns();
      int numberRows = model->numberRows();
      int pivotRow = model->pivotRow();
      if (gubSlackIn_>=0)
      assert (sequenceIn>numberRows+numberColumns);
      if (sequenceIn==sequenceOut) 
      return -1;
      int iSetIn=-1;
      int iSetOut=-1;
      if (sequenceOut<numberColumns) {
      iSetOut = backward_[sequenceOut];
      } else if (sequenceOut>=numberRows+numberColumns) {
      assert (pivotRow>=numberRows);
      int iExtra = pivotRow-numberRows;
      assert (iExtra>=0);
      if (iSetOut<0)
        iSetOut = fromIndex_[iExtra];
      else
        assert(iSetOut == fromIndex_[iExtra]);
      }
      if (sequenceIn<numberColumns) {
      iSetIn = backward_[sequenceIn];
      } else if (gubSlackIn_>=0) {
      iSetIn = gubSlackIn_;
      }
      possiblePivotKey_=-1;
      number=0; // say do ordinary
      int * pivotVariable = model->pivotVariable();
      if (pivotRow>=numberRows) {
      int iExtra = pivotRow-numberRows;
      //const int * length = matrix_->getVectorLengths();

      assert (sequenceOut>=numberRows+numberColumns||
            sequenceOut==keyVariable_[iSetOut]);
      int incomingColumn = sequenceIn; // to be used in updates
      if (iSetIn!=iSetOut) {
        // We need to find a possible pivot for incoming
        // look through rowArray_[1]
        int n = model->rowArray(1)->getNumElements();
        int * which = model->rowArray(1)->getIndices();
        double * array = model->rowArray(1)->denseVector();
        double bestAlpha = 1.0e-5;
        //int shortest=numberRows+1;
        for (int i=0;i<n;i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (iPivot<numberColumns&&backward_[iPivot]==iSetOut) {
            if (fabs(array[i])>fabs(bestAlpha)) {
            bestAlpha = array[i];
            possiblePivotKey_=iRow;
            }
          }
        }
        assert (possiblePivotKey_>=0); // could set returnCode=4
        number=1;
        if (sequenceIn>=numberRows+numberColumns) {
          number=3;
          // need swap as gub slack in and must become key
          // is this best way
          int key=keyVariable_[iSetIn];
          assert (key<numberColumns);
          // check other basic
          int iColumn = next_[key];
          // set new key to be used by unpack
          keyVariable_[iSetIn]=iSetIn+numberColumns;
          // change cost in changeCost
          {
            int iLook=0;
            int iSet=fromIndex_[0];
            while (iSet>=0) {
            if (iSet==iSetIn) {
              changeCost_[iLook]=0.0;
              break;
            }
            iSet = fromIndex_[++iLook];
            }
          }
          while (iColumn>=0) {
            if (iColumn!=sequenceOut) {
            // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
            printf("TTTTTry 5\n");
#endif
            int iRow = backToPivotRow_[iColumn];
            assert (iRow>=0);
            unpack(model,model->rowArray(3),iColumn);
            model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
            double alpha = model->rowArray(3)->denseVector()[iRow];
            //if (!alpha)
            //printf("zero alpha a\n");
            int updateStatus = model->factorization()->replaceColumn(model,
                                                       model->rowArray(2),
                                                       model->rowArray(3),
                                                       iRow, alpha);
            returnCode = CoinMax(updateStatus, returnCode);
            model->rowArray(3)->clear();
            if (returnCode)
              break;
            }
            iColumn=next_[iColumn];
          }
          if (!returnCode) {
            // now factorization looks as if key is out
            // pivot back in
#ifdef CLP_DEBUG_PRINT
            printf("TTTTTry 6\n");
#endif
            unpack(model,model->rowArray(3),key);
            model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
            pivotRow = possiblePivotKey_;
            double alpha = model->rowArray(3)->denseVector()[pivotRow];
            //if (!alpha)
            //printf("zero alpha b\n");
            int updateStatus = model->factorization()->replaceColumn(model,
                                                       model->rowArray(2),
                                                       model->rowArray(3),
                                                       pivotRow, alpha);
            returnCode = CoinMax(updateStatus, returnCode);
            model->rowArray(3)->clear();
          }
          // restore key
          keyVariable_[iSetIn]=key;
          // now alternate column can replace key on out 
          incomingColumn = pivotVariable[possiblePivotKey_];
        } else {
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTTry 4 %d\n",possiblePivotKey_);
#endif
          int updateStatus = model->factorization()->replaceColumn(model,
                                                     model->rowArray(2),
                                                     model->rowArray(1),
                                                     possiblePivotKey_, 
                                                     bestAlpha);
          returnCode = CoinMax(updateStatus, returnCode);
          incomingColumn = pivotVariable[possiblePivotKey_];
        }
        
        //returnCode=4; // need swap
      } else {
        // key swap 
        number=-1;
      }
      int key=keyVariable_[iSetOut];
      if (key<numberColumns)
        assert(key==sequenceOut);
      // check if any other basic
      int iColumn = next_[key];
      if (returnCode)
        iColumn = -1; // skip if error on previous
      // set new key to be used by unpack
      if (incomingColumn<numberColumns)
        keyVariable_[iSetOut]=incomingColumn;
      else
        keyVariable_[iSetOut]=iSetIn+numberColumns;
      double * cost = model->costRegion();
      if (possiblePivotKey_<0) {
        double dj = model->djRegion()[sequenceIn]-cost[sequenceIn];
        changeCost_[iExtra] = -dj;
#ifdef CLP_DEBUG_PRINT
        printf("modifying changeCost %d by %g - cost %g\n",iExtra, dj,cost[sequenceIn]);
#endif
      }
      while (iColumn>=0) {
        if (iColumn!=incomingColumn) {
          number=-2;
          // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTTry 1\n");
#endif
          int iRow = backToPivotRow_[iColumn];
          assert (iRow>=0&&iRow<numberRows);
          unpack(model,model->rowArray(3),iColumn);
          model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
          double * array = model->rowArray(3)->denseVector();
          double alpha = array[iRow];
          //if (!alpha)
          //printf("zero alpha d\n");
          int updateStatus = model->factorization()->replaceColumn(model,
                                                     model->rowArray(2),
                                                     model->rowArray(3),
                                                     iRow, alpha);
          returnCode = CoinMax(updateStatus, returnCode);
          model->rowArray(3)->clear();
          if (returnCode)
            break;
        }
        iColumn=next_[iColumn];
      }
      // restore key
      keyVariable_[iSetOut]=key;
      } else if (sequenceIn>=numberRows+numberColumns) {
      number=2;
      //returnCode=4;
      // need swap as gub slack in and must become key
      // is this best way
      int key=keyVariable_[iSetIn];
      assert (key<numberColumns);
      // check other basic
      int iColumn = next_[key];
      // set new key to be used by unpack
      keyVariable_[iSetIn]=iSetIn+numberColumns;
      // change cost in changeCost
      {
        int iLook=0;
        int iSet=fromIndex_[0];
        while (iSet>=0) {
          if (iSet==iSetIn) {
            changeCost_[iLook]=0.0;
            break;
          }
          iSet = fromIndex_[++iLook];
        }
      }
      while (iColumn>=0) {
        if (iColumn!=sequenceOut) {
          // need partial ftran and skip accuracy check in replaceColumn
#ifdef CLP_DEBUG_PRINT
          printf("TTTTTry 2\n");
#endif
          int iRow = backToPivotRow_[iColumn];
          assert (iRow>=0);
          unpack(model,model->rowArray(3),iColumn);
          model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
          double alpha = model->rowArray(3)->denseVector()[iRow];
          //if (!alpha)
          //printf("zero alpha e\n");
          int updateStatus = model->factorization()->replaceColumn(model,
                                                     model->rowArray(2),
                                                     model->rowArray(3),
                                                     iRow, alpha);
          returnCode = CoinMax(updateStatus, returnCode);
          model->rowArray(3)->clear();
          if (returnCode)
            break;
        }
        iColumn=next_[iColumn];
      }
      if (!returnCode) {
        // now factorization looks as if key is out
        // pivot back in
#ifdef CLP_DEBUG_PRINT
        printf("TTTTTry 3\n");
#endif
        unpack(model,model->rowArray(3),key);
        model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(3));
        double alpha = model->rowArray(3)->denseVector()[pivotRow];
        //if (!alpha)
        //printf("zero alpha f\n");
        int updateStatus = model->factorization()->replaceColumn(model,
                                                   model->rowArray(2),
                                                   model->rowArray(3),
                                                   pivotRow, alpha);
        returnCode = CoinMax(updateStatus, returnCode);
        model->rowArray(3)->clear();
      }
      // restore key
      keyVariable_[iSetIn]=key;
      } else {
      // normal - but might as well do here
      returnCode = model->factorization()->replaceColumn(model,
                                             model->rowArray(2),
                                             model->rowArray(1),
                                             model->pivotRow(),
                                             model->alpha());
      }
    }
#ifdef CLP_DEBUG_PRINT
    printf("Update type after %d - status %d - pivot row %d\n",
         number,returnCode,model->pivotRow());
#endif
    // see if column generation says time to re-factorize
    returnCode = CoinMax(returnCode,synchronize(model,5));
    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:
    {
      synchronize(model,0);
      CoinMemcpyN(status_,numberSets_,saveStatus_);
      CoinMemcpyN(keyVariable_,numberSets_,savedKeyVariable_);
    }
    break;
    // restore status
  case 6:
    {
      CoinMemcpyN(saveStatus_,numberSets_,status_);
      CoinMemcpyN(savedKeyVariable_,numberSets_,keyVariable_);
      // restore firstAvailable_
      synchronize(model,7);
      // redo next_
      int i;
      int * last = new int[numberSets_];
      for (i=0;i<numberSets_;i++) {
      int iKey=keyVariable_[i];
      assert(iKey>=numberColumns||backward_[iKey]==i);
      last[i]=iKey;
      // make sure basic
      //if (iKey<numberColumns)
      //model->setStatus(iKey,ClpSimplex::basic);
      }
      for (i=0;i<numberColumns;i++){
      int iSet = backward_[i];
      if (iSet>=0) {
        next_[last[iSet]]=i;
        last[iSet]=i;
      }
      }
      for (i=0;i<numberSets_;i++) {
      next_[last[i]]=-(keyVariable_[i]+1);
      redoSet(model,keyVariable_[i],keyVariable_[i],i);
      }
      delete [] last;
      // redo pivotVariable
      int * pivotVariable = model->pivotVariable();
      int iRow;
      int numberBasic=0;
      int numberRows = model->numberRows();
      for (iRow=0;iRow<numberRows;iRow++) {
      if (model->getRowStatus(iRow)==ClpSimplex::basic) {
        numberBasic++;
        pivotVariable[iRow]=iRow+numberColumns;
      } else {
        pivotVariable[iRow]=-1;
      }
      }
      i=0;
      int iColumn;
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (model->getStatus(iColumn)==ClpSimplex::basic) {
        int iSet = backward_[iColumn];
        if (iSet<0||keyVariable_[iSet]!=iColumn) {
          while (pivotVariable[i]>=0) {
            i++;
            assert (i<numberRows);
          }
          pivotVariable[i]=iColumn;
          backToPivotRow_[iColumn]=i;
          numberBasic++;
        }
      }
      }
      assert (numberBasic==numberRows);
      rhsOffset(model,true);
    }
    break;
    // flag a variable
  case 7:
    {
      assert (number==model->sequenceIn());
      synchronize(model,1);
      synchronize(model,8);
    }
    break;
    // unflag all variables
  case 8:
    {
      returnCode=synchronize(model,2);
    }
    break;
    // redo costs in primal
  case 9:
    {
      returnCode=synchronize(model,3);
    }
    break;
    // return 1 if there may be changing bounds on variable (column generation)
  case 10:
    {
      returnCode=synchronize(model,6);
    }
    break;
    // make sure set is clean
  case 11:
    {
      assert (number==model->sequenceIn());
      returnCode=synchronize(model,8);
    }
    break;
  default:
    break;
  }
  return returnCode;
}
// Sets up an effective RHS and does gub crash if needed
void 
03049 ClpGubMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
{
  // Do basis - cheapest or slack if feasible (unless cheapest set)
  int longestSet=0;
  int iSet;
  for (iSet=0;iSet<numberSets_;iSet++) 
    longestSet = CoinMax(longestSet,end_[iSet]-start_[iSet]);
    
  double * upper = new double[longestSet+1];
  double * cost = new double[longestSet+1];
  double * lower = new double[longestSet+1];
  double * solution = new double[longestSet+1];
  assert (!next_);
  int numberColumns = getNumCols();
  const int * columnLength = matrix_->getVectorLengths(); 
  const double * columnLower = model->lowerRegion();
  const double * columnUpper = model->upperRegion();
  double * columnSolution = model->solutionRegion();
  const double * objective = model->costRegion();
  int numberRows = getNumRows();
  toIndex_ = new int[numberSets_];
  for (iSet=0;iSet<numberSets_;iSet++) 
    toIndex_[iSet]=-1;
  fromIndex_ = new int [getNumRows()+1];
  double tolerance = model->primalTolerance();
  bool noNormalBounds=true;
  gubType_ &= ~8;
  bool gotBasis=false;
  for (iSet=0;iSet<numberSets_;iSet++) {
    if (keyVariable_[iSet]<numberColumns)
      gotBasis=true;
    CoinBigIndex j;
    CoinBigIndex iStart = start_[iSet];
    CoinBigIndex iEnd=end_[iSet];
    for (j=iStart;j<iEnd;j++) {
      if (columnLower[j]&&columnLower[j]>-1.0e20)
      noNormalBounds=false;
      if (columnUpper[j]&&columnUpper[j]<1.0e20)
      noNormalBounds=false;
    }
  }
  if (noNormalBounds)
    gubType_ |= 8;
  if (!gotBasis) {
    for (iSet=0;iSet<numberSets_;iSet++) {
      CoinBigIndex j;
      int numberBasic=0;
      int iBasic=-1;
      CoinBigIndex iStart = start_[iSet];
      CoinBigIndex iEnd=end_[iSet];
      // find one with smallest length
      int smallest=numberRows+1;
      double value=0.0;
      for (j=iStart;j<iEnd;j++) {
      if (model->getStatus(j)== ClpSimplex::basic) {
        if (columnLength[j]<smallest) {
          smallest=columnLength[j];
          iBasic=j;
        }
        numberBasic++;
      }
      value += columnSolution[j];
      }
      bool done=false;
      if (numberBasic>1||(numberBasic==1&&getStatus(iSet)==ClpSimplex::basic)) {
      if (getStatus(iSet)==ClpSimplex::basic) 
        iBasic = iSet+numberColumns;// slack key - use
      done=true;
      } else if (numberBasic==1) {
      // see if can be key
      double thisSolution = columnSolution[iBasic];
      if (thisSolution>columnUpper[iBasic]) {
        value -= thisSolution-columnUpper[iBasic];
        thisSolution = columnUpper[iBasic];
        columnSolution[iBasic]=thisSolution;
      }
      if (thisSolution<columnLower[iBasic]) {
        value -= thisSolution-columnLower[iBasic];
        thisSolution = columnLower[iBasic];
        columnSolution[iBasic]=thisSolution;
      }
      // try setting slack to a bound
      assert (upper_[iSet]<1.0e20||lower_[iSet]>-1.0e20);
      double cost1 = COIN_DBL_MAX;
      int whichBound=-1;
      if (upper_[iSet]<1.0e20) {
        // try slack at ub
        double newBasic = thisSolution +upper_[iSet]-value;
        if (newBasic>=columnLower[iBasic]-tolerance&&
            newBasic<=columnUpper[iBasic]+tolerance) {
          // can go
          whichBound=1;
          cost1 = newBasic*objective[iBasic];
          // But if exact then may be good solution
          if (fabs(upper_[iSet]-value)<tolerance)
            cost1=-COIN_DBL_MAX;
        }
      }
      if (lower_[iSet]>-1.0e20) {
        // try slack at lb
        double newBasic = thisSolution +lower_[iSet]-value;
        if (newBasic>=columnLower[iBasic]-tolerance&&
            newBasic<=columnUpper[iBasic]+tolerance) {
          // can go but is it cheaper
          double cost2 = newBasic*objective[iBasic];
          // But if exact then may be good solution
          if (fabs(lower_[iSet]-value)<tolerance)
            cost2=-COIN_DBL_MAX;
          if (cost2<cost1)
            whichBound=0;
        }
      }
      if (whichBound!=-1) {
        // key
        done=true;
        if (whichBound) {
          // slack to upper
          columnSolution[iBasic]=thisSolution + upper_[iSet]-value;
          setStatus(iSet,ClpSimplex::atUpperBound);
        } else {
          // slack to lower
          columnSolution[iBasic]=thisSolution + lower_[iSet]-value;
          setStatus(iSet,ClpSimplex::atLowerBound);
        }
      }
      }
      if (!done) {
      if (!cheapest) {
        // see if slack can be key
        if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
          done=true;
          setStatus(iSet,ClpSimplex::basic);
          iBasic=iSet+numberColumns;
        }
      }
      if (!done) {
        // set non basic if there was one
        if (iBasic>=0)
          model->setStatus(iBasic,ClpSimplex::atLowerBound);
        // find cheapest
        int numberInSet = iEnd-iStart;
        CoinMemcpyN(columnLower+iStart,numberInSet,lower);
        CoinMemcpyN(columnUpper+iStart,numberInSet,upper);
        CoinMemcpyN(columnSolution+iStart,numberInSet,solution);
        // and slack
        iBasic=numberInSet;
        solution[iBasic]=-value;
        lower[iBasic]=-upper_[iSet];
        upper[iBasic]=-lower_[iSet];
        int kphase;
        if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
          // feasible
          kphase=1;
          cost[iBasic]=0.0;
          CoinMemcpyN(objective+iStart,numberInSet,cost);
        } else {
          // infeasible
          kphase=0;
          // remember bounds are flipped so opposite to natural
          if (value<lower_[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;
            CoinMemcpyN(objective+iStart,numberInSet,cost);
          }
          // 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)
            }
            }
          }
        }
        // convert iBasic back and do bounds
        if (iBasic==numberInSet) {
          // slack basic
          setStatus(iSet,ClpSimplex::basic);
          iBasic=iSet+numberColumns;
        } else {
          iBasic += start_[iSet];
          model->setStatus(iBasic,ClpSimplex::basic);
          // 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();
        }
        for (j=iStart;j<iEnd;j++) {
          if (model->getStatus(j)!=ClpSimplex::basic) {
            int inSet=j-iStart;
            columnSolution[j]=solution[inSet];
            if (upper[inSet]==lower[inSet]) 
            model->setStatus(j,ClpSimplex::isFixed);
            else if (solution[inSet]==upper[inSet])
            model->setStatus(j,ClpSimplex::atUpperBound);
            else if (solution[inSet]==lower[inSet])
            model->setStatus(j,ClpSimplex::atLowerBound);
          }
        }
      }
      } 
      keyVariable_[iSet]=iBasic;
    }
  }
  delete [] lower;
  delete [] solution;
  delete [] upper;
  delete [] cost;
  // make sure matrix is in good shape
  matrix_->orderMatrix();
  // create effective rhs
  delete [] rhsOffset_;
  rhsOffset_ = new double[numberRows];
  delete [] next_;
  next_ = new int[numberColumns+numberSets_+2*longestSet];
  char * mark = new char[numberColumns];
  memset(mark,0,numberColumns);
  for (int iColumn=0;iColumn<numberColumns;iColumn++) 
    next_[iColumn]=COIN_INT_MAX;
  int i;
  int * keys = new int[numberSets_];
  for (i=0;i<numberSets_;i++) 
    keys[i]=COIN_INT_MAX;
  // set up chains
  for (i=0;i<numberColumns;i++){
    if (model->getStatus(i)==ClpSimplex::basic) 
      mark[i]=1;
    int iSet = backward_[i];
    if (iSet>=0) {
      int iNext = keys[iSet];
      next_[i]=iNext;
      keys[iSet]=i;
    }
  }
  for (i=0;i<numberSets_;i++) {
    int j;
    if (getStatus(i)!=ClpSimplex::basic) {
      // make sure fixed if it is
      if (upper_[i]==lower_[i])
      setStatus(i,ClpSimplex::isFixed);
      // slack not key - choose one with smallest length
      int smallest=numberRows+1;
      int key=-1;
      j = keys[i];
      if (j!=COIN_INT_MAX) {
      while (1) {
        if (mark[j]&&columnLength[j]<smallest&&!gotBasis) {
          key=j;
          smallest=columnLength[j];
        }
        if (next_[j]!=COIN_INT_MAX) {
          j = next_[j];
        } else {
          // correct end
          next_[j]=-(keys[i]+1);
          break;
        }
      }
      } else {
      next_[i+numberColumns] = -(numberColumns+i+1);
      }  
      if (gotBasis)
      key =keyVariable_[i];
      if (key>=0) {
      keyVariable_[i]=key;
      } else {
      // nothing basic - make slack key
      //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
      // fudge to avoid const problem
      status_[i]=1;
      }
    } else {
      // slack key
      keyVariable_[i]=numberColumns+i;
      int j;
      double sol=0.0;
      j = keys[i];
      if (j!=COIN_INT_MAX) {
      while (1) {
        sol += columnSolution[j];
        if (next_[j]!=COIN_INT_MAX) {
          j = next_[j];
        } else {
          // correct end
          next_[j]=-(keys[i]+1);
          break;
        }
      }
      } else {
      next_[i+numberColumns] = -(numberColumns+i+1);
      }
      if (sol>upper_[i]+tolerance) {
      setAbove(i);
      } else if (sol<lower_[i]-tolerance) {
      setBelow(i);
      } else {
      setFeasible(i);
      }
    }
    // Create next_
    int key = keyVariable_[i];
    redoSet(model,key,keys[i],i);
  }
  delete [] keys;
  delete [] mark;
  rhsOffset(model,true);
}
// redoes next_ for a set.  
void 
03463 ClpGubMatrix::redoSet(ClpSimplex * model, int newKey, int oldKey, int iSet)
{
  int numberColumns = model->numberColumns();
  int * save = next_+numberColumns+numberSets_;
  int number=0;
  int stop = -(oldKey+1);
  int j=next_[oldKey];
  while (j!=stop) {
    if (j<0)
      j = -j-1;
    if (j!=newKey)
      save[number++]=j;
    j=next_[j];
  }
  // and add oldkey
  if (newKey!=oldKey)
    save[number++]=oldKey;
  // now do basic
  int lastMarker = -(newKey+1);
  keyVariable_[iSet]=newKey;
  next_[newKey]=lastMarker;
  int last = newKey;
  for ( j=0;j<number;j++) {
    int iColumn=save[j];
    if (iColumn<numberColumns) {
      if (model->getStatus(iColumn)==ClpSimplex::basic) {
      next_[last]=iColumn;
      next_[iColumn]=lastMarker;
      last = iColumn;
      }
    }
  }
  // now add in non-basic
  for ( j=0;j<number;j++) {
    int iColumn=save[j];
    if (iColumn<numberColumns) {
      if (model->getStatus(iColumn)!=ClpSimplex::basic) {
      next_[last]=-(iColumn+1);
      next_[iColumn]=lastMarker;
      last = iColumn;
      }
    }
  }

}
/* 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 * 
03512 ClpGubMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,bool 
#ifdef CLP_DEBUG
check
#endif
)
{
  //forceRefresh=true;
  if (rhsOffset_) {
#ifdef CLP_DEBUG
    if (check) {
      // no need - but check anyway
      // zero out basic
      int numberRows = model->numberRows();
      int numberColumns = model->numberColumns();
      double * solution = new double [numberColumns];
      double * rhs = new double[numberRows];
      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
      CoinZeroN(rhs,numberRows);
      int iRow;
      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
        solution[iColumn]=0.0;
      }
      for (int iSet=0;iSet<numberSets_;iSet++) {
      int iColumn = keyVariable_[iSet];
      if (iColumn<numberColumns) 
        solution[iColumn]=0.0;
      }
      times(-1.0,solution,rhs);
      delete [] solution;
      const double * columnSolution = model->solutionRegion();
      // and now subtract out non basic
      ClpSimplex::Status iStatus;
      for (int iSet=0;iSet<numberSets_;iSet++) {
      int iColumn = keyVariable_[iSet];
      if (iColumn<numberColumns) {
        double b=0.0;
        // key is structural - where is slack
        iStatus = getStatus(iSet);
        assert (iStatus!=ClpSimplex::basic);
        if (iStatus==ClpSimplex::atLowerBound)
          b=lower_[iSet];
        else
          b=upper_[iSet];
        // subtract out others at bounds
        if ((gubType_&8)==0) {
          int stop = -(iColumn+1);
          int jColumn =next_[iColumn];
          // sum all non-basic variables - first skip basic
          while(jColumn>=0) 
            jColumn=next_[jColumn];
          while(jColumn!=stop) {
            assert (jColumn<0);
            jColumn = -jColumn-1;
            b -= columnSolution[jColumn];
            jColumn=next_[jColumn];
          }
        }
        // subtract out
        ClpPackedMatrix::add(model,rhs,iColumn,-b);
      }
      }
      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]);
      }
      delete [] rhs;
    }
#endif
    if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
                   lastRefresh_+refreshFrequency_)) {
      // zero out basic
      int numberRows = model->numberRows();
      int numberColumns = model->numberColumns();
      double * solution = new double [numberColumns];
      CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
      CoinZeroN(rhsOffset_,numberRows);
      for (int iColumn=0;iColumn<numberColumns;iColumn++) {
      if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
        solution[iColumn]=0.0;
      }
      int iSet;
      for ( iSet=0;iSet<numberSets_;iSet++) {
      int iColumn = keyVariable_[iSet];
      if (iColumn<numberColumns) 
        solution[iColumn]=0.0;
      }
      times(-1.0,solution,rhsOffset_);
      delete [] solution;
      lastRefresh_ = model->numberIterations();
      const double * columnSolution = model->solutionRegion();
      // and now subtract out non basic
      ClpSimplex::Status iStatus;
      for ( iSet=0;iSet<numberSets_;iSet++) {
      int iColumn = keyVariable_[iSet];
      if (iColumn<numberColumns) {
        double b=0.0;
        // key is structural - where is slack
        iStatus = getStatus(iSet);
        assert (iStatus!=ClpSimplex::basic);
        if (iStatus==ClpSimplex::atLowerBound)
          b=lower_[iSet];
        else
          b=upper_[iSet];
        // subtract out others at bounds
        if ((gubType_&8)==0) {
          int stop = -(iColumn+1);
          int jColumn =next_[iColumn];
          // sum all non-basic variables - first skip basic
          while(jColumn>=0) 
            jColumn=next_[jColumn];
          while(jColumn!=stop) {
            assert (jColumn<0);
            jColumn = -jColumn-1;
            b -= columnSolution[jColumn];
            jColumn=next_[jColumn];
          }
        }
        // subtract out
        if (b)
          ClpPackedMatrix::add(model,rhsOffset_,iColumn,-b);
      }
      }
    }
  }
  return rhsOffset_;
}
/* 
   update information for a pivot (and effective rhs)
*/
int 
03643 ClpGubMatrix::updatePivot(ClpSimplex * model,double oldInValue, double /*oldOutValue*/)
{
  int sequenceIn = model->sequenceIn();
  int sequenceOut = model->sequenceOut();
  double * solution = model->solutionRegion();
  int numberColumns = model->numberColumns();
  int numberRows = model->numberRows();
  int pivotRow = model->pivotRow();
  int iSetIn;
  // Correct sequence in
  trueSequenceIn_=sequenceIn;
  if (sequenceIn<numberColumns) {
    iSetIn = backward_[sequenceIn];
  } else if (sequenceIn<numberColumns+numberRows) {
    iSetIn = -1;
  } else {
    iSetIn = gubSlackIn_;
    trueSequenceIn_ = numberColumns+numberRows+iSetIn;
  }
  int iSetOut=-1;
  trueSequenceOut_=sequenceOut;
  if (sequenceOut<numberColumns) {
    iSetOut = backward_[sequenceOut];
  } else if (sequenceOut>=numberRows+numberColumns) {
    assert (pivotRow>=numberRows);
    int iExtra = pivotRow-numberRows;
    assert (iExtra>=0);
    if (iSetOut<0)
      iSetOut = fromIndex_[iExtra];
    else
      assert(iSetOut == fromIndex_[iExtra]);
    trueSequenceOut_ = numberColumns+numberRows+iSetOut;
  }
  if (rhsOffset_) {
    // update effective rhs
    if (sequenceIn==sequenceOut) {
      assert (sequenceIn<numberRows+numberColumns); // should be easy to deal with
      if (sequenceIn<numberColumns)
      add(model,rhsOffset_,sequenceIn,oldInValue-solution[sequenceIn]);
    } else {
      if (sequenceIn<numberColumns) {
      // we need to test if WILL be key
      ClpPackedMatrix::add(model,rhsOffset_,sequenceIn,oldInValue);
      if (iSetIn>=0)  {
        // old contribution to rhsOffset_
        int key = keyVariable_[iSetIn];
        if (key<numberColumns) {
          double oldB=0.0;
          ClpSimplex::Status iStatus = getStatus(iSetIn);
          if (iStatus==ClpSimplex::atLowerBound)
            oldB=lower_[iSetIn];
          else
            oldB=upper_[iSetIn];
          // subtract out others at bounds
          if ((gubType_&8)==0) {
            int stop = -(key+1);
            int iColumn =next_[key];
            // skip basic
            while (iColumn>=0)
            iColumn = next_[iColumn];
            // sum all non-key variables
            while(iColumn!=stop) {
            assert (iColumn<0);
            iColumn = -iColumn-1;
            if (iColumn == sequenceIn) 
              oldB -= oldInValue;
            else if ( iColumn != sequenceOut )
              oldB -= solution[iColumn];
            iColumn=next_[iColumn];
            }
          }
          if (oldB)
            ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
        }
      }
      } else if (sequenceIn<numberRows+numberColumns) {
      //rhsOffset_[sequenceIn-numberColumns] -= oldInValue;
      } else {
#ifdef CLP_DEBUG_PRINT
      printf("** in is key slack %d\n",sequenceIn);
#endif
      // old contribution to rhsOffset_
      int key = keyVariable_[iSetIn];
      if (key<numberColumns) {
        double oldB=0.0;
        ClpSimplex::Status iStatus = getStatus(iSetIn);
        if (iStatus==ClpSimplex::atLowerBound)
          oldB=lower_[iSetIn];
        else
          oldB=upper_[iSetIn];
        // subtract out others at bounds
        if ((gubType_&8)==0) {
          int stop = -(key+1);
          int iColumn =next_[key];
          // skip basic
          while (iColumn>=0)
            iColumn = next_[iColumn];
          // sum all non-key variables
          while(iColumn!=stop) {
            assert (iColumn<0);
            iColumn = -iColumn-1;
            if ( iColumn != sequenceOut )
            oldB -= solution[iColumn];
            iColumn=next_[iColumn];
          }
        }
        if (oldB)
          ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
      }
      }
      if (sequenceOut<numberColumns) {
      ClpPackedMatrix::add(model,rhsOffset_,sequenceOut,-solution[sequenceOut]);
      if (iSetOut>=0) {
        // old contribution to rhsOffset_
        int key = keyVariable_[iSetOut];
        if (key<numberColumns&&iSetIn!=iSetOut) {
          double oldB=0.0;
          ClpSimplex::Status iStatus = getStatus(iSetOut);
          if (iStatus==ClpSimplex::atLowerBound)
            oldB=lower_[iSetOut];
          else
            oldB=upper_[iSetOut];
          // subtract out others at bounds
          if ((gubType_&8)==0) {
            int stop = -(key+1);
            int iColumn =next_[key];
            // skip basic
            while (iColumn>=0)
            iColumn = next_[iColumn];
            // sum all non-key variables
            while(iColumn!=stop) {
            assert (iColumn<0);
            iColumn = -iColumn-1;
            if (iColumn == sequenceIn) 
              oldB -= oldInValue;
            else if ( iColumn != sequenceOut )
              oldB -= solution[iColumn];
            iColumn=next_[iColumn];
            }
          }
          if (oldB)
            ClpPackedMatrix::add(model,rhsOffset_,key,oldB);
        }
      }
      } else if (sequenceOut<numberRows+numberColumns) {
      //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut];
      } else {
#ifdef CLP_DEBUG_PRINT
      printf("** out is key slack %d\n",sequenceOut);
#endif
      assert (pivotRow>=numberRows);
      }
    }
  }
  int * pivotVariable = model->pivotVariable();
  // may need to deal with key
  // Also need coding to mark/allow key slack entering
  if (pivotRow>=numberRows) {
    assert (sequenceOut>=numberRows+numberColumns||sequenceOut==keyVariable_[iSetOut]);
#ifdef CLP_DEBUG_PRINT
    if (sequenceIn>=numberRows+numberColumns)
      printf("key slack %d in, set out %d\n",gubSlackIn_,iSetOut);
    printf("** danger - key out for set %d in %d (set %d)\n",iSetOut,sequenceIn,
         iSetIn);
#endif   
    // if slack out mark correctly
    if (sequenceOut>=numberRows+numberColumns) {
      double value=model->valueOut();
      if (value==upper_[iSetOut]) {
      setStatus(iSetOut,ClpSimplex::atUpperBound);
      } else if (value==lower_[iSetOut]) {
      setStatus(iSetOut,ClpSimplex::atLowerBound);
      } else {
      if (fabs(value-upper_[iSetOut])<
          fabs(value-lower_[iSetOut])) {
        setStatus(iSetOut,ClpSimplex::atUpperBound);
      } else {
        setStatus(iSetOut,ClpSimplex::atLowerBound);
      }
      }
      if (upper_[iSetOut]==lower_[iSetOut])
      setStatus(iSetOut,ClpSimplex::isFixed);
      setFeasible(iSetOut);
    }
    if (iSetOut==iSetIn) {
      // key swap
      int key;
      if (sequenceIn>=numberRows+numberColumns) {
      key = numberColumns+iSetIn;
      setStatus(iSetIn,ClpSimplex::basic);
      } else {
      key = sequenceIn;
      }
      redoSet(model,key,keyVariable_[iSetIn],iSetIn);
    } else {
      // key was chosen
      assert (possiblePivotKey_>=0&&possiblePivotKey_<numberRows);
      int key=pivotVariable[possiblePivotKey_];
      // and set incoming here
      if (sequenceIn>=numberRows+numberColumns) {
      // slack in - so use old key
      sequenceIn = keyVariable_[iSetIn];
      model->setStatus(sequenceIn,ClpSimplex::basic);
      setStatus(iSetIn,ClpSimplex::basic);
      redoSet(model,iSetIn+numberColumns,keyVariable_[iSetIn],iSetIn);
      }
      //? do not do if iSetIn<0 ? as will be done later
      pivotVariable[possiblePivotKey_]=sequenceIn;
      if (sequenceIn<numberColumns)
      backToPivotRow_[sequenceIn]=possiblePivotKey_;
      redoSet(model,key,keyVariable_[iSetOut],iSetOut);
    }
  } else {
    if (sequenceOut<numberColumns) {
      if (iSetIn>=0&&iSetOut==iSetIn) {
      // key not out - only problem is if slack in
      int key;
      if (sequenceIn>=numberRows+numberColumns) {
        key = numberColumns+iSetIn;
        setStatus(iSetIn,ClpSimplex::basic);
        assert (pivotRow<numberRows);
        // must swap with current key
        int key=keyVariable_[iSetIn];
        model->setStatus(key,ClpSimplex::basic);
        pivotVariable[pivotRow]=key;
        backToPivotRow_[key]=pivotRow;
      } else {
        key = keyVariable_[iSetIn];
      }
      redoSet(model,key,keyVariable_[iSetIn],iSetIn);
      } else if (iSetOut>=0) {
      // just redo set
      int key=keyVariable_[iSetOut];;
      redoSet(model,key,keyVariable_[iSetOut],iSetOut);
      }
    }
  }
  if (iSetIn>=0&&iSetIn!=iSetOut) {
    int key=keyVariable_[iSetIn];
    if (sequenceIn == numberColumns+2*numberRows) {
      // key slack in
      assert (pivotRow<numberRows);
      // must swap with current key
      model->setStatus(key,ClpSimplex::basic);
      pivotVariable[pivotRow]=key;
      backToPivotRow_[key]=pivotRow;
      setStatus(iSetIn,ClpSimplex::basic);
      key = iSetIn+numberColumns;
    }
    // redo set to allow for new one
    redoSet(model,key,keyVariable_[iSetIn],iSetIn);
  }
  // update pivot 
  if (sequenceIn<numberColumns) {
    if (pivotRow<numberRows) {
      backToPivotRow_[sequenceIn]=pivotRow;
    } else {
      if (possiblePivotKey_>=0) {
      assert (possiblePivotKey_<numberRows);
      backToPivotRow_[sequenceIn]=possiblePivotKey_;
      pivotVariable[possiblePivotKey_]=sequenceIn;
      }
    }
  } else if (sequenceIn>=numberRows+numberColumns) {
    // key in - something should have been done before
    int key = keyVariable_[iSetIn];
    assert (key==numberColumns+iSetIn);
    //pivotVariable[pivotRow]=key;
    //backToPivotRow_[key]=pivotRow;
    //model->setStatus(key,ClpSimplex::basic);
    //key=numberColumns+iSetIn;
    setStatus(iSetIn,ClpSimplex::basic);
    redoSet(model,key,keyVariable_[iSetIn],iSetIn);
  }
#ifdef CLP_DEBUG
  {
    char * xx = new char[numberColumns+numberRows];
    memset(xx,0,numberRows+numberColumns);
    for (int i=0;i<numberRows;i++) {
      int iPivot = pivotVariable[i];
      assert (iPivot<numberRows+numberColumns);
      assert (!xx[iPivot]);
      xx[iPivot]=1;
      if (iPivot<numberColumns) {
      int iBack = backToPivotRow_[iPivot];
      assert (i==iBack);
      }
    }
    delete [] xx;
  }
#endif
  if (rhsOffset_) {
    // update effective rhs
    if (sequenceIn!=sequenceOut) {
      if (sequenceIn<numberColumns) {
      if (iSetIn>=0) {
        // new contribution to rhsOffset_
        int key = keyVariable_[iSetIn];
        if (key<numberColumns) {
          double newB=0.0;
          ClpSimplex::Status iStatus = getStatus(iSetIn);
          if (iStatus==ClpSimplex::atLowerBound)
            newB=lower_[iSetIn];
          else
            newB=upper_[iSetIn];
          // subtract out others at bounds
          if ((gubType_&8)==0) {
            int stop = -(key+1);
            int iColumn =next_[key];
            // skip basic
            while (iColumn>=0)
            iColumn = next_[iColumn];
            // sum all non-key variables
            while(iColumn!=stop) {
            assert (iColumn<0);
            iColumn = -iColumn-1;
            newB -= solution[iColumn];
            iColumn=next_[iColumn];
            }
          }
          if (newB)
            ClpPackedMatrix::add(model,rhsOffset_,key,-newB);
        }
      }
      }
      if (iSetOut>=0) {
      // new contribution to rhsOffset_
      int key = keyVariable_[iSetOut];
      if (key<numberColumns&&iSetIn!=iSetOut) {
        double newB=0.0;
        ClpSimplex::Status iStatus = getStatus(iSetOut);
        if (iStatus==ClpSimplex::atLowerBound)
          newB=lower_[iSetOut];
        else
          newB=upper_[iSetOut];
        // subtract out others at bounds
        if ((gubType_&8)==0) {
          int stop = -(key+1);
          int iColumn =next_[key];
          // skip basic
          while (iColumn>=0)
            iColumn = next_[iColumn];
          // sum all non-key variables
          while(iColumn!=stop) {
            assert (iColumn<0);
            iColumn = -iColumn-1;
            newB -= solution[iColumn];
            iColumn=next_[iColumn];
          }
        }
        if (newB)
          ClpPackedMatrix::add(model,rhsOffset_,key,-newB);
      }
      }
    }
  }
#ifdef CLP_DEBUG
  // debug
  {
    int i;
    char * xxxx = new char[numberColumns];
    memset(xxxx,0,numberColumns);
    for (i=0;i<numberRows;i++) {
      int iPivot = pivotVariable[i];
      assert (model->getStatus(iPivot)==ClpSimplex::basic);
      if (iPivot<numberColumns&&backward_[iPivot]>=0)
      xxxx[iPivot]=1;
    }
    double primalTolerance = model->primalTolerance();
    for (i=0;i<numberSets_;i++) {
      int key=keyVariable_[i];
      double value=0.0;
      // sum over all except key
      int iColumn =next_[key];
      // sum all non-key variables
      int k=0;
      int stop = -(key+1);
      while (iColumn!=stop) {
      if (iColumn<0)
        iColumn = -iColumn-1;
      value += solution[iColumn];
      k++;
      assert (k<100);
      assert (backward_[iColumn]==i);
      iColumn=next_[iColumn];
      }
      iColumn = next_[key];
      if (key<numberColumns) {
      // feasibility will be done later
      assert (getStatus(i)!=ClpSimplex::basic);
      double sol;
      if (getStatus(i)==ClpSimplex::atUpperBound)
        sol = upper_[i]-value;
      else
        sol = lower_[i]-value;
      //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol,
      //     cost[key]);
      //if (fabs(sol-solution[key])>1.0e-3)
      //printf("** stored value was %g\n",solution[key]);
      } else {
      // slack is key
      double infeasibility=0.0;
      if (value>upper_[i]+primalTolerance) {
        infeasibility=value-upper_[i]-primalTolerance;
        //setAbove(i);
      } else if (value<lower_[i]-primalTolerance) {
        infeasibility=lower_[i]-value-primalTolerance ;
        //setBelow(i);
      } else {
        //setFeasible(i);
      }
      //printf("xx Value of key slack for set %d is %g\n",i,value);
      }
      while (iColumn>=0) {
      assert (xxxx[iColumn]);
      xxxx[iColumn]=0;
      iColumn=next_[iColumn];
      }
    }
    for (i=0;i<numberColumns;i++) {
      if (i<numberColumns&&backward_[i]>=0) {
      assert (!xxxx[i]||i==keyVariable_[backward_[i]]);
      }
    }
    delete [] xxxx;
  }
#endif
  return 0;
}
// Switches off dj checking each factorization (for BIG models)
void 
04074 ClpGubMatrix::switchOffCheck()
{
  noCheck_=0;
  infeasibilityWeight_=0.0;
}
// Correct sequence in and out to give true value
void 
04081 ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/,int & sequenceIn, int & sequenceOut)
{
  if (sequenceIn!=-999) {
    sequenceIn = trueSequenceIn_;
    sequenceOut = trueSequenceOut_;
  }
}

Generated by  Doxygen 1.6.0   Back to index