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

ClpPrimalColumnSteepest.cpp

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

#include "CoinPragma.hpp"

#include "ClpSimplex.hpp"
#include "ClpPrimalColumnSteepest.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpFactorization.hpp"
#include "ClpNonLinearCost.hpp"
#include "ClpMessage.hpp"
#include "CoinHelperFunctions.hpp"
#include <stdio.h>
//#define CLP_DEBUG
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################

//-------------------------------------------------------------------
// Default Constructor 
//-------------------------------------------------------------------
00023 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode) 
  : ClpPrimalColumnPivot(),
    devex_(0.0),
    weights_(NULL),
    infeasible_(NULL),
    alternateWeights_(NULL),
    savedWeights_(NULL),
    reference_(NULL),
    state_(-1),
    mode_(mode),
    persistence_(normal),
    numberSwitched_(0),
    pivotSequence_(-1),
    savedPivotSequence_(-1),
    savedSequenceOut_(-1),
    sizeFactorization_(0)
{
  type_=2+64*mode;
}
//-------------------------------------------------------------------
// Copy constructor 
//-------------------------------------------------------------------
00045 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs) 
: ClpPrimalColumnPivot(rhs)
{  
  state_=rhs.state_;
  mode_ = rhs.mode_;
  persistence_ = rhs.persistence_;
  numberSwitched_ = rhs.numberSwitched_;
  model_ = rhs.model_;
  pivotSequence_ = rhs.pivotSequence_;
  savedPivotSequence_ = rhs.savedPivotSequence_;
  savedSequenceOut_ = rhs.savedSequenceOut_;
  sizeFactorization_ = rhs.sizeFactorization_;
  devex_ = rhs.devex_;
  if ((model_&&model_->whatsChanged()&1)!=0) {
    if (rhs.infeasible_) {
      infeasible_= new CoinIndexedVector(rhs.infeasible_);
    } else {
      infeasible_=NULL;
    }
    reference_=NULL;
    if (rhs.weights_) {
      assert(model_);
      int number = model_->numberRows()+model_->numberColumns();
      assert (number == rhs.model_->numberRows()+rhs.model_->numberColumns());
      weights_= new double[number];
      CoinMemcpyN(rhs.weights_,number,weights_);
      savedWeights_= new double[number];
      CoinMemcpyN(rhs.savedWeights_,number,savedWeights_);
      if (mode_!=1) {
        reference_ = CoinCopyOfArray(rhs.reference_,(number+31)>>5);
      }
    } else {
      weights_=NULL;
      savedWeights_=NULL;
    }
    if (rhs.alternateWeights_) {
      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    } else {
      alternateWeights_=NULL;
    }
  } else {
    infeasible_=NULL;
    reference_=NULL;
    weights_=NULL;
    savedWeights_=NULL;
    alternateWeights_=NULL;
  }
}

//-------------------------------------------------------------------
// Destructor 
//-------------------------------------------------------------------
00097 ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
{
  delete [] weights_;
  delete infeasible_;
  delete alternateWeights_;
  delete [] savedWeights_;
  delete [] reference_;
}

//----------------------------------------------------------------
// Assignment operator 
//-------------------------------------------------------------------
ClpPrimalColumnSteepest &
00110 ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
{
  if (this != &rhs) {
    ClpPrimalColumnPivot::operator=(rhs);
    state_=rhs.state_;
    mode_ = rhs.mode_;
    persistence_ = rhs.persistence_;
    numberSwitched_ = rhs.numberSwitched_;
    model_ = rhs.model_;
    pivotSequence_ = rhs.pivotSequence_;
    savedPivotSequence_ = rhs.savedPivotSequence_;
    savedSequenceOut_ = rhs.savedSequenceOut_;
    sizeFactorization_ = rhs.sizeFactorization_;
    devex_ = rhs.devex_;
    delete [] weights_;
    delete [] reference_;
    reference_=NULL;
    delete infeasible_;
    delete alternateWeights_;
    delete [] savedWeights_;
    savedWeights_ = NULL;
    if (rhs.infeasible_!=NULL) {
      infeasible_= new CoinIndexedVector(rhs.infeasible_);
    } else {
      infeasible_=NULL;
    }
    if (rhs.weights_!=NULL) {
      assert(model_);
      int number = model_->numberRows()+model_->numberColumns();
      assert (number == rhs.model_->numberRows()+rhs.model_->numberColumns());
      weights_= new double[number];
      CoinMemcpyN(rhs.weights_,number,weights_);
      savedWeights_= new double[number];
      CoinMemcpyN(rhs.savedWeights_,number,savedWeights_);
      if (mode_!=1) {
      reference_ = CoinCopyOfArray(rhs.reference_,(number+31)>>5);
      }
    } else {
      weights_=NULL;
    }
    if (rhs.alternateWeights_!=NULL) {
      alternateWeights_= new CoinIndexedVector(rhs.alternateWeights_);
    } else {
      alternateWeights_=NULL;
    }
  }
  return *this;
}
// These have to match ClpPackedMatrix version
#define TRY_NORM 1.0e-4
#define ADD_ONE 1.0
// Returns pivot column, -1 if none
/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
      that is just +-weight where a feasibility changed.  It also has 
      reduced cost from last iteration in pivot row*/
int 
00166 ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
                            CoinIndexedVector * spareRow1,
                            CoinIndexedVector * spareRow2,
                            CoinIndexedVector * spareColumn1,
                            CoinIndexedVector * spareColumn2)
{
  assert(model_);
  if (model_->nonLinearCost()->lookBothWays()||model_->algorithm()==2) {
    // Do old way
    updates->expand();
    return pivotColumnOldMethod(updates,spareRow1,spareRow2,
                        spareColumn1,spareColumn2);
  }
  int number=0;
  int * index;
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  int anyUpdates;
  double * infeas = infeasible_->denseVector();

  // Local copy of mode so can decide what to do
  int switchType;
  if (mode_==4) 
    switchType = 5-numberSwitched_;
  else if (mode_>=10)
    switchType=3;
  else
    switchType = mode_;
  /* switchType - 
     0 - all exact devex
     1 - all steepest
     2 - some exact devex
     3 - auto some exact devex
     4 - devex
     5 - dantzig
     10 - can go to mini-sprint
  */
  // Look at gub
#if 1
  model_->clpMatrix()->dualExpanded(model_,updates,NULL,4);
#else
  updates->clear();
  model_->computeDuals(NULL);
#endif
  if (updates->getNumElements()>1) {
    // would have to have two goes for devex, three for steepest
    anyUpdates=2;
  } else if (updates->getNumElements()) {
    if (updates->getIndices()[0]==pivotRow&&fabs(updates->denseVector()[0])>1.0e-6) {
      // reasonable size
      anyUpdates=1;
      //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
      //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
    } else {
      // too small
      anyUpdates=2;
    }
  } else if (pivotSequence_>=0){
    // just after re-factorization
    anyUpdates=-1;
  } else {
    // sub flip - nothing to do
    anyUpdates=0;
  }
  int sequenceOut = model_->sequenceOut();
  if (switchType==5) {
    // If known matrix then we will do partial pricing
    if (model_->clpMatrix()->canDoPartialPricing()) {
      pivotSequence_=-1;
      pivotRow=-1;
      // See if to switch
      int numberRows = model_->numberRows();
      int numberWanted=10;
      int numberColumns = model_->numberColumns();
      int numberHiddenRows = model_->clpMatrix()->hiddenRows();
      double ratio = static_cast<double> (sizeFactorization_+numberHiddenRows)/
      static_cast<double> (numberRows + 2*numberHiddenRows);
      // Number of dual infeasibilities at last invert
      int numberDual = model_->numberDualInfeasibilities();
      int numberLook = CoinMin(numberDual,numberColumns/10);
      if (ratio<1.0) {
      numberWanted = 100;
      numberLook /= 20;
      numberWanted = CoinMax(numberWanted,numberLook);
      } else if (ratio<3.0) {
      numberWanted = 500;
      numberLook /= 15;
      numberWanted = CoinMax(numberWanted,numberLook);
      } else if (ratio<4.0||mode_==5) {
      numberWanted = 1000;
      numberLook /= 10;
      numberWanted = CoinMax(numberWanted,numberLook);
      } else if (mode_!=5) {
      switchType=4;
      // initialize
      numberSwitched_++;
      // Make sure will re-do
      delete [] weights_;
      weights_=NULL;
      model_->computeDuals(NULL);
      saveWeights(model_,4);
      anyUpdates=0;
      printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
      }
      if (switchType==5) {
      numberLook *= 5; // needs tuning for gub
        if (model_->numberIterations()%1000==0&&model_->logLevel()>1) {
        printf("numels %d ratio %g wanted %d look %d\n",
             sizeFactorization_,ratio,numberWanted,numberLook);
        }
      // Update duals and row djs
      // Do partial pricing
      return partialPricing(updates,spareRow2,
                        numberWanted,numberLook);
      }
    }
  }
  if (switchType==5) {
    if (anyUpdates>0) {
      justDjs(updates,spareRow2,
      spareColumn1,spareColumn2);
    }
  } else if (anyUpdates==1) {
    if (switchType<4) {
      // exact etc when can use dj 
      djsAndSteepest(updates,spareRow2,
      spareColumn1,spareColumn2);
    } else {
      // devex etc when can use dj 
      djsAndDevex(updates,spareRow2,
      spareColumn1,spareColumn2);
    }
  } else if (anyUpdates==-1) {
    if (switchType<4) {
      // exact etc when djs okay 
      justSteepest(updates,spareRow2,
      spareColumn1,spareColumn2);
    } else {
      // devex etc when djs okay 
      justDevex(updates,spareRow2,
      spareColumn1,spareColumn2);
    }
  } else if (anyUpdates==2) {
    if (switchType<4) {
      // exact etc when have to use pivot
      djsAndSteepest2(updates,spareRow2,
      spareColumn1,spareColumn2);
    } else {
      // devex etc when have to use pivot
      djsAndDevex2(updates,spareRow2,
      spareColumn1,spareColumn2);
    }
  } 
#ifdef CLP_DEBUG
  alternateWeights_->checkClear();
#endif
  // make sure outgoing from last iteration okay
  if (sequenceOut>=0) {
    ClpSimplex::Status status = model_->getStatus(sequenceOut);
    double value = model_->reducedCost(sequenceOut);
    
    switch(status) {
      
    case ClpSimplex::basic:
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) { 
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
      break;
    case ClpSimplex::atUpperBound:
      if (value>tolerance) {
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
      break;
    case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
    }
  }
  // update of duals finished - now do pricing
  // See what sort of pricing
  int numberWanted=10;
  number = infeasible_->getNumElements();
  int numberColumns = model_->numberColumns();
  if (switchType==5) {
    pivotSequence_=-1;
    pivotRow=-1;
    // See if to switch
    int numberRows = model_->numberRows();
    // ratio is done on number of columns here
    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
    double ratio = static_cast<double> (sizeFactorization_)/static_cast<double> (numberRows);
    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
    if (ratio<1.0) {
      numberWanted = CoinMax(100,number/200);
    } else if (ratio<2.0-1.0) {
      numberWanted = CoinMax(500,number/40);
    } else if (ratio<4.0-3.0||mode_==5) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/30);
    } else if (mode_!=5) {
      switchType=4;
      // initialize
      numberSwitched_++;
      // Make sure will re-do
      delete [] weights_;
      weights_=NULL;
      saveWeights(model_,4);
      printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
    }
    //if (model_->numberIterations()%1000==0)
    //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
  }
  int numberRows = model_->numberRows();
  // ratio is done on number of rows here
  double ratio = static_cast<double> (sizeFactorization_)/static_cast<double> (numberRows);
  if(switchType==4) {
    // Still in devex mode
    // Go to steepest if lot of iterations?
    if (ratio<5.0) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/20);
    } else if (ratio<7.0) {
      numberWanted = CoinMax(2000,number/5);
      numberWanted = CoinMax(numberWanted,numberColumns/10);
    } else {
      // we can zero out
      updates->clear();
      spareColumn1->clear();
      switchType=3;
      // initialize
      pivotSequence_=-1;
      pivotRow=-1;
      numberSwitched_++;
      // Make sure will re-do
      delete [] weights_;
      weights_=NULL;
      saveWeights(model_,4);
      printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
      updates->clear();
    }
    if (model_->numberIterations()%1000==0)
    printf("numels %d ratio %g wanted %d type x\n",sizeFactorization_,ratio,numberWanted);
  } 
  if (switchType<4) {
    if (switchType<2 ) {
      numberWanted = number+1;
    } else if (switchType==2) {
      numberWanted = CoinMax(2000,number/8);
    } else {
      if (ratio<1.0) {
      numberWanted = CoinMax(2000,number/20);
      } else if (ratio<5.0) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/40);
      } else if (ratio<10.0) {
      numberWanted = CoinMax(2000,number/8);
      numberWanted = CoinMax(numberWanted,numberColumns/20);
      } else {
      ratio = number * (ratio/80.0);
      if (ratio>number) {
        numberWanted=number+1;
      } else {
        numberWanted = CoinMax(2000,static_cast<int> (ratio));
        numberWanted = CoinMax(numberWanted,numberColumns/10);
      }
      }
    }
    //if (model_->numberIterations()%1000==0)
    //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
    //switchType);
  }


  double bestDj = 1.0e-30;
  int bestSequence=-1;

  int i,iSequence;
  index = infeasible_->getIndices();
  number = infeasible_->getNumElements();
  // Re-sort infeasible every 100 pivots
  // Not a good idea
  if (0&&model_->factorization()->pivots()>0&&
      (model_->factorization()->pivots()%100)==0) {
    int nLook = model_->numberRows()+numberColumns;
    number=0;
    for (i=0;i<nLook;i++) {
      if (infeas[i]) { 
      if (fabs(infeas[i])>COIN_INDEXED_TINY_ELEMENT) 
        index[number++]=i;
      else
        infeas[i]=0.0;
      }
    }
    infeasible_->setNumElements(number);
  }
  if(model_->numberIterations()<model_->lastBadIteration()+200&&
     model_->factorization()->pivots()>10) {
    // we can't really trust infeasibilities if there is dual error
    double checkTolerance = 1.0e-8;
    if (model_->largestDualError()>checkTolerance)
      tolerance *= model_->largestDualError()/checkTolerance;
    // But cap
    tolerance = CoinMin(1000.0,tolerance);
  }
#ifdef CLP_DEBUG
  if (model_->numberDualInfeasibilities()==1) 
    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
         model_->largestDualError(),model_,model_->messageHandler(),
         number);
#endif
  // stop last one coming immediately
  double saveOutInfeasibility=0.0;
  if (sequenceOut>=0) {
    saveOutInfeasibility = infeas[sequenceOut];
    infeas[sequenceOut]=0.0;
  }
  if (model_->factorization()->pivots()&&model_->numberPrimalInfeasibilities())
    tolerance = CoinMax(tolerance,1.0e-10*model_->infeasibilityCost());
  tolerance *= tolerance; // as we are using squares

  int iPass;
  // Setup two passes
  int start[4];
  start[1]=number;
  start[2]=0;
  double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
  start[0]=static_cast<int> (dstart);
  start[3]=start[0];
  //double largestWeight=0.0;
  //double smallestWeight=1.0e100;
  for (iPass=0;iPass<2;iPass++) {
    int end = start[2*iPass+1];
    if (switchType<5) {
      for (i=start[2*iPass];i<end;i++) {
      iSequence = index[i];
      double value = infeas[iSequence];
      double weight = weights_[iSequence];
      if (value>tolerance) {
        //weight=1.0;
        if (value>bestDj*weight) {
          // check flagged variable and correct dj
          if (!model_->flagged(iSequence)) {
            bestDj=value/weight;
            bestSequence = iSequence;
          } else {
            // just to make sure we don't exit before got something
            numberWanted++;
          }
        }
        numberWanted--;
      }
      if (!numberWanted)
        break;
      }
    } else {
      // Dantzig
      for (i=start[2*iPass];i<end;i++) {
      iSequence = index[i];
      double value = infeas[iSequence];
      if (value>tolerance) {
        if (value>bestDj) {
          // check flagged variable and correct dj
          if (!model_->flagged(iSequence)) {
            bestDj=value;
            bestSequence = iSequence;
          } else {
            // just to make sure we don't exit before got something
            numberWanted++;
          }
        }
        numberWanted--;
      }
      if (!numberWanted)
        break;
      }
    }
    if (!numberWanted)
      break;
  }
  model_->clpMatrix()->setSavedBestSequence(bestSequence);
  if (bestSequence>=0)
    model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
  if (sequenceOut>=0) {
    infeas[sequenceOut]=saveOutInfeasibility;
  }
  /*if (model_->numberIterations()%100==0)
    printf("%d best %g\n",bestSequence,bestDj);*/
  
#ifndef NDEBUG
  if (bestSequence>=0) {
    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
      assert(model_->reducedCost(bestSequence)<0.0);
    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound) {
      assert(model_->reducedCost(bestSequence)>0.0);
    }
  }
#endif
  return bestSequence;
}
// Just update djs
void 
00597 ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int iSection,j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  double * infeas = infeasible_->denseVector();
  //updates->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  
  // put row of tableau in rowArray and columnArray (packed mode)
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  // normal
  for (iSection=0;iSection<2;iSection++) {
    
    reducedCost=model_->djRegion(iSection);
    int addSequence;
    
    if (!iSection) {
      number = updates->getNumElements();
      index = updates->getIndices();
      updateBy = updates->denseVector();
      addSequence = model_->numberColumns();
    } else {
      number = spareColumn1->getNumElements();
      index = spareColumn1->getIndices();
      updateBy = spareColumn1->denseVector();
      addSequence = 0;
    }
    
    for (j=0;j<number;j++) {
      int iSequence = index[j];
      double value = reducedCost[iSequence];
      value -= updateBy[j];
      updateBy[j]=0.0;
      reducedCost[iSequence] = value;
      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
      
      switch(status) {

      case ClpSimplex::basic:
      infeasible_->zero(iSequence+addSequence);
      case ClpSimplex::isFixed:
      break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) {
        // we are going to bias towards free (but only if reasonable)
        value *= FREE_BIAS;
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atUpperBound:
      if (value>tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      }
    }
  }
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
  if (pivotRow>=0) {
    // make sure infeasibility on incoming is 0.0
    int sequenceIn = model_->sequenceIn();
    infeasible_->zero(sequenceIn);
  }
}
// Update djs, weights for Devex
void 
00701 ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  // for weights update we use pivotSequence
  // unset in case sub flip
  assert (pivotSequence_>=0);
  assert (model_->pivotVariable()[pivotSequence_]==model_->sequenceIn());
  pivotSequence_=-1;
  double * infeas = infeasible_->denseVector();
  //updates->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  // and we can see if reference
  double referenceIn=0.0;
  int sequenceIn = model_->sequenceIn();
  if (mode_!=1&&reference(sequenceIn))
    referenceIn=1.0;
  // save outgoing weight round update
  double outgoingWeight=0.0;
  int sequenceOut = model_->sequenceOut();
  if (sequenceOut>=0)
    outgoingWeight=weights_[sequenceOut];
    
  double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
  // put row of tableau in rowArray and columnArray (packed mode)
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  // update weights
  double * weight;
  int numberColumns = model_->numberColumns();
  // rows
  reducedCost=model_->djRegion(0);
  int addSequence=model_->numberColumns();;
    
  number = updates->getNumElements();
  index = updates->getIndices();
  updateBy = updates->denseVector();
  weight = weights_+numberColumns;
  // Devex
  for (j=0;j<number;j++) {
    double thisWeight;
    double pivot;
    double value3;
    int iSequence = index[j];
    double value = reducedCost[iSequence];
    double value2 = updateBy[j];
    updateBy[j]=0.0;
    value -= value2;
    reducedCost[iSequence] = value;
    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    
    switch(status) {
      
    case ClpSimplex::basic:
      infeasible_->zero(iSequence+addSequence);
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (fabs(value)>FREE_ACCEPT*tolerance) {
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
      break;
    case ClpSimplex::atUpperBound:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (value>tolerance) {
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
      break;
    case ClpSimplex::atLowerBound:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (value<-tolerance) {
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
    }
  }
  
  // columns
  weight = weights_;
  
  scaleFactor = -scaleFactor;
  reducedCost=model_->djRegion(1);
  number = spareColumn1->getNumElements();
  index = spareColumn1->getIndices();
  updateBy = spareColumn1->denseVector();

  // Devex
  
  for (j=0;j<number;j++) {
    double thisWeight;
    double pivot;
    double value3;
    int iSequence = index[j];
    double value = reducedCost[iSequence];
    double value2 = updateBy[j];
    value -= value2;
    updateBy[j]=0.0;
    reducedCost[iSequence] = value;
    ClpSimplex::Status status = model_->getStatus(iSequence);
    
    switch(status) {
      
    case ClpSimplex::basic:
      infeasible_->zero(iSequence);
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (fabs(value)>FREE_ACCEPT*tolerance) {
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
      break;
    case ClpSimplex::atUpperBound:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (value>tolerance) {
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
      break;
    case ClpSimplex::atLowerBound:
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      value3 = pivot * pivot*devex_;
      if (reference(iSequence))
      value3 += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value3);
      if (value<-tolerance) {
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
    }
  }
  // restore outgoing weight
  if (sequenceOut>=0)
    weights_[sequenceOut]=outgoingWeight;
  // make sure infeasibility on incoming is 0.0
  infeasible_->zero(sequenceIn);
  spareRow2->setNumElements(0);
  //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
  // check for accuracy
  int iCheck=892;
  //printf("weight for iCheck is %g\n",weights_[iCheck]);
  int numberRows = model_->numberRows();
  //int numberColumns = model_->numberColumns();
  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
      !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
  }
#endif
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
}
// Update djs, weights for Steepest
void 
00940 ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  // for weights update we use pivotSequence
  // unset in case sub flip
  assert (pivotSequence_>=0);
  assert (model_->pivotVariable()[pivotSequence_]==model_->sequenceIn());
  double * infeas = infeasible_->denseVector();
  double scaleFactor = 1.0/updates->denseVector()[0]; // as formula is with 1.0
  assert (updates->getIndices()[0]==pivotSequence_);
  pivotSequence_=-1;
  //updates->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  //alternateWeights_->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,
                                     alternateWeights_);
  // and we can see if reference
  int sequenceIn = model_->sequenceIn();
  double referenceIn;
  if (mode_!=1) {
    if(reference(sequenceIn))
      referenceIn=1.0;
    else
      referenceIn=0.0;
  } else {
    referenceIn=-1.0;
  }
  // save outgoing weight round update
  double outgoingWeight=0.0;
  int sequenceOut = model_->sequenceOut();
  if (sequenceOut>=0)
    outgoingWeight=weights_[sequenceOut];
  // update row weights here so we can scale alternateWeights_
  // update weights
  double * weight;
  double * other = alternateWeights_->denseVector();
  int numberColumns = model_->numberColumns();
  // rows
  reducedCost=model_->djRegion(0);
  int addSequence=model_->numberColumns();;
    
  number = updates->getNumElements();
  index = updates->getIndices();
  updateBy = updates->denseVector();
  weight = weights_+numberColumns;
  
  for (j=0;j<number;j++) {
    double thisWeight;
    double pivot;
    double modification;
    double pivotSquared;
    int iSequence = index[j];
    double value2 = updateBy[j];
    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    double value;
    
    switch(status) {
      
    case ClpSimplex::basic:
      infeasible_->zero(iSequence+addSequence);
      reducedCost[iSequence] = 0.0;
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      value = reducedCost[iSequence] - value2;
      modification = other[iSequence];
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      pivotSquared = pivot * pivot;
      
      thisWeight += pivotSquared * devex_ + pivot * modification;
      reducedCost[iSequence] = value;
      if (thisWeight<TRY_NORM) {
      if (mode_==1) {
        // steepest
        thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
      } else {
        // exact
        thisWeight = referenceIn*pivotSquared;
        if (reference(iSequence+numberColumns))
          thisWeight += 1.0;
        thisWeight = CoinMax(thisWeight,TRY_NORM);
      }
      }
      weight[iSequence] = thisWeight;
      if (fabs(value)>FREE_ACCEPT*tolerance) {
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
      break;
    case ClpSimplex::atUpperBound:
      value = reducedCost[iSequence] - value2;
      modification = other[iSequence];
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      pivotSquared = pivot * pivot;
      
      thisWeight += pivotSquared * devex_ + pivot * modification;
      reducedCost[iSequence] = value;
      if (thisWeight<TRY_NORM) {
      if (mode_==1) {
        // steepest
        thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
      } else {
        // exact
        thisWeight = referenceIn*pivotSquared;
        if (reference(iSequence+numberColumns))
          thisWeight += 1.0;
        thisWeight = CoinMax(thisWeight,TRY_NORM);
      }
      }
      weight[iSequence] = thisWeight;
      if (value>tolerance) {
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
      break;
    case ClpSimplex::atLowerBound:
      value = reducedCost[iSequence] - value2;
      modification = other[iSequence];
      thisWeight = weight[iSequence];
      // row has -1 
      pivot = value2*scaleFactor;
      pivotSquared = pivot * pivot;
      
      thisWeight += pivotSquared * devex_ + pivot * modification;
      reducedCost[iSequence] = value;
      if (thisWeight<TRY_NORM) {
      if (mode_==1) {
        // steepest
        thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
      } else {
        // exact
        thisWeight = referenceIn*pivotSquared;
        if (reference(iSequence+numberColumns))
          thisWeight += 1.0;
        thisWeight = CoinMax(thisWeight,TRY_NORM);
      }
      }
      weight[iSequence] = thisWeight;
      if (value<-tolerance) {
      // store square in list
      if (infeas[iSequence+addSequence])
        infeas[iSequence+addSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
      infeasible_->zero(iSequence+addSequence);
      }
    }
  }
  // put row of tableau in rowArray and columnArray (packed)
  // get subset which have nonzero tableau elements
  transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,
                  -scaleFactor);
  // zero updateBy
  CoinZeroN(updateBy,number);
  alternateWeights_->clear();
  // columns
  assert (scaleFactor);
  reducedCost=model_->djRegion(1);
  number = spareColumn1->getNumElements();
  index = spareColumn1->getIndices();
  updateBy = spareColumn1->denseVector();

  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double value = reducedCost[iSequence];
    double value2 = updateBy[j];
    updateBy[j]=0.0;
    value -= value2;
    reducedCost[iSequence] = value;
    ClpSimplex::Status status = model_->getStatus(iSequence);
    
    switch(status) {
      
    case ClpSimplex::basic:
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) {
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
      break;
    case ClpSimplex::atUpperBound:
      if (value>tolerance) {
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
      break;
    case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
      // store square in list
      if (infeas[iSequence])
        infeas[iSequence] = value*value; // already there
      else
        infeasible_->quickAdd(iSequence,value*value);
      } else {
      infeasible_->zero(iSequence);
      }
    }
  }
  // restore outgoing weight
  if (sequenceOut>=0)
    weights_[sequenceOut]=outgoingWeight;
  // make sure infeasibility on incoming is 0.0
  infeasible_->zero(sequenceIn);
  spareColumn2->setNumElements(0);
  //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
  // check for accuracy
  int iCheck=892;
  //printf("weight for iCheck is %g\n",weights_[iCheck]);
  int numberRows = model_->numberRows();
  //int numberColumns = model_->numberColumns();
  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
      !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
  }
#endif
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
}
// Update djs, weights for Devex
void 
01208 ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int iSection,j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  // dj could be very small (or even zero - take care)
  double dj = model_->dualIn();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  double * infeas = infeasible_->denseVector();
  //updates->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  
  // put row of tableau in rowArray and columnArray
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  // normal
  for (iSection=0;iSection<2;iSection++) {
    
    reducedCost=model_->djRegion(iSection);
    int addSequence;
    
    if (!iSection) {
      number = updates->getNumElements();
      index = updates->getIndices();
      updateBy = updates->denseVector();
      addSequence = model_->numberColumns();
    } else {
      number = spareColumn1->getNumElements();
      index = spareColumn1->getIndices();
      updateBy = spareColumn1->denseVector();
      addSequence = 0;
    }
    
    for (j=0;j<number;j++) {
      int iSequence = index[j];
      double value = reducedCost[iSequence];
      value -= updateBy[j];
      updateBy[j]=0.0;
      reducedCost[iSequence] = value;
      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
      
      switch(status) {

      case ClpSimplex::basic:
      infeasible_->zero(iSequence+addSequence);
      case ClpSimplex::isFixed:
      break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) {
        // we are going to bias towards free (but only if reasonable)
        value *= FREE_BIAS;
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atUpperBound:
      if (value>tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      }
    }
  }
  // They are empty
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
  // make sure infeasibility on incoming is 0.0
  int sequenceIn = model_->sequenceIn();
  infeasible_->zero(sequenceIn);
  // for weights update we use pivotSequence
  if (pivotSequence_>=0) {
    pivotRow = pivotSequence_;
    // unset in case sub flip
    pivotSequence_=-1;
    // make sure infeasibility on incoming is 0.0
    const int * pivotVariable = model_->pivotVariable();
    sequenceIn = pivotVariable[pivotRow];
    infeasible_->zero(sequenceIn);
    // and we can see if reference
    double referenceIn=0.0;
    if (mode_!=1&&reference(sequenceIn))
      referenceIn=1.0;
    // save outgoing weight round update
    double outgoingWeight=0.0;
    int sequenceOut = model_->sequenceOut();
    if (sequenceOut>=0)
      outgoingWeight=weights_[sequenceOut];
    // update weights
    updates->setNumElements(0);
    spareColumn1->setNumElements(0);
    // might as well set dj to 1
    dj = 1.0;
    updates->insert(pivotRow,-dj);
    model_->factorization()->updateColumnTranspose(spareRow2,updates);
    // put row of tableau in rowArray and columnArray
    model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
    double * weight;
    int numberColumns = model_->numberColumns();
    // rows
    number = updates->getNumElements();
    index = updates->getIndices();
    updateBy = updates->denseVector();
    weight = weights_+numberColumns;
    
    assert (devex_>0.0);
    for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = - updateBy[iSequence];
      updateBy[iSequence]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
      value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
    }
    
    // columns
    weight = weights_;
    
    number = spareColumn1->getNumElements();
    index = spareColumn1->getIndices();
    updateBy = spareColumn1->denseVector();
    for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = updateBy[iSequence];
      updateBy[iSequence]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence))
      value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
    }
    // restore outgoing weight
    if (sequenceOut>=0)
      weights_[sequenceOut]=outgoingWeight;
    spareColumn2->setNumElements(0);
    //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
    // check for accuracy
    int iCheck=892;
    //printf("weight for iCheck is %g\n",weights_[iCheck]);
    int numberRows = model_->numberRows();
    //int numberColumns = model_->numberColumns();
    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
    }
#endif
    updates->setNumElements(0);
    spareColumn1->setNumElements(0);
  }
}
// Update djs, weights for Steepest
void 
01399 ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int iSection,j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  // dj could be very small (or even zero - take care)
  double dj = model_->dualIn();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  double * infeas = infeasible_->denseVector();
  //updates->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  
  // put row of tableau in rowArray and columnArray
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  // normal
  for (iSection=0;iSection<2;iSection++) {
    
    reducedCost=model_->djRegion(iSection);
    int addSequence;
    
    if (!iSection) {
      number = updates->getNumElements();
      index = updates->getIndices();
      updateBy = updates->denseVector();
      addSequence = model_->numberColumns();
    } else {
      number = spareColumn1->getNumElements();
      index = spareColumn1->getIndices();
      updateBy = spareColumn1->denseVector();
      addSequence = 0;
    }
    
    for (j=0;j<number;j++) {
      int iSequence = index[j];
      double value = reducedCost[iSequence];
      value -= updateBy[j];
      updateBy[j]=0.0;
      reducedCost[iSequence] = value;
      ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
      
      switch(status) {

      case ClpSimplex::basic:
      infeasible_->zero(iSequence+addSequence);
      case ClpSimplex::isFixed:
      break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) {
        // we are going to bias towards free (but only if reasonable)
        value *= FREE_BIAS;
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atUpperBound:
      if (value>tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      break;
      case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
        // store square in list
        if (infeas[iSequence+addSequence])
          infeas[iSequence+addSequence] = value*value; // already there
        else
          infeasible_->quickAdd(iSequence+addSequence,value*value);
      } else {
        infeasible_->zero(iSequence+addSequence);
      }
      }
    }
  }
  // we can zero out as will have to get pivot row
  // ***** move
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
  if (pivotRow>=0) {
    // make sure infeasibility on incoming is 0.0
    int sequenceIn = model_->sequenceIn();
    infeasible_->zero(sequenceIn);
  }
  // for weights update we use pivotSequence
  pivotRow = pivotSequence_;
  // unset in case sub flip
  pivotSequence_=-1;
  if (pivotRow>=0) {
    // make sure infeasibility on incoming is 0.0
    const int * pivotVariable = model_->pivotVariable();
    int sequenceIn = pivotVariable[pivotRow];
    assert (sequenceIn==model_->sequenceIn());
    infeasible_->zero(sequenceIn);
    // and we can see if reference
    double referenceIn;
    if (mode_!=1) {
      if(reference(sequenceIn))
        referenceIn=1.0;
      else
        referenceIn=0.0;
    } else {
      referenceIn=-1.0;
    }
    // save outgoing weight round update
    double outgoingWeight=0.0;
    int sequenceOut = model_->sequenceOut();
    if (sequenceOut>=0)
      outgoingWeight=weights_[sequenceOut];
    // update weights
    updates->setNumElements(0);
    spareColumn1->setNumElements(0);
    // might as well set dj to 1
    dj = -1.0;
    updates->createPacked(1,&pivotRow,&dj);
    model_->factorization()->updateColumnTranspose(spareRow2,updates);
    bool needSubset =  (mode_<4||numberSwitched_>1||mode_>=10);

    double * weight;
    double * other = alternateWeights_->denseVector();
    int numberColumns = model_->numberColumns();
    // rows
    number = updates->getNumElements();
    index = updates->getIndices();
    updateBy = updates->denseVector();
    weight = weights_+numberColumns;
    if (needSubset) {
      // now update weight update array
      model_->factorization()->updateColumnTranspose(spareRow2,
                                         alternateWeights_);
      // do alternateWeights_ here so can scale
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      assert (iSequence>=0&&iSequence<model_->numberRows());
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = - updateBy[j];
      double modification = other[iSequence];
      double pivotSquared = pivot * pivot;

      thisWeight += pivotSquared * devex_ + pivot * modification;
      if (thisWeight<TRY_NORM) {
        if (mode_==1) {
          // steepest
          thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
        } else {
          // exact
          thisWeight = referenceIn*pivotSquared;
          if (reference(iSequence+numberColumns))
            thisWeight += 1.0;
          thisWeight = CoinMax(thisWeight,TRY_NORM);
        }
      }
      weight[iSequence] = thisWeight;
      }
      transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,0.0);
    } else {
      // put row of tableau in rowArray and columnArray
      model_->clpMatrix()->transposeTimes(model_,-1.0,
                                          updates,spareColumn2,spareColumn1);
    }
    
    if (needSubset) {
      CoinZeroN(updateBy,number);
    } else if (mode_==4) {
      // Devex
      assert (devex_>0.0);
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = -updateBy[j];
      updateBy[j]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
        value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
      }
    }
    
    // columns
    weight = weights_;
    
    number = spareColumn1->getNumElements();
    index = spareColumn1->getIndices();
    updateBy = spareColumn1->denseVector();
    if (needSubset) {
      // Exact - already done
    } else if (mode_==4) {
      // Devex
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = updateBy[j];
      updateBy[j]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence))
        value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
      }
    }
    // restore outgoing weight
    if (sequenceOut>=0)
      weights_[sequenceOut]=outgoingWeight;
    alternateWeights_->clear();
    spareColumn2->setNumElements(0);
    //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
    // check for accuracy
    int iCheck=892;
    //printf("weight for iCheck is %g\n",weights_[iCheck]);
    int numberRows = model_->numberRows();
    //int numberColumns = model_->numberColumns();
    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
    }
#endif
  }
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
}
// Updates two arrays for steepest
void 
01646 ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
                                         const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
                                         CoinIndexedVector * spare,
                                         double scaleFactor)
{
  // see if reference
  int sequenceIn = model_->sequenceIn();
  double referenceIn;
  if (mode_!=1) {
    if(reference(sequenceIn))
      referenceIn=1.0;
    else
      referenceIn=0.0;
  } else {
    referenceIn=-1.0;
  }
  if (model_->clpMatrix()->canCombine(model_,pi1)) {
    // put row of tableau in rowArray and columnArray
    model_->clpMatrix()->transposeTimes2(model_,pi1,dj1,pi2,spare,referenceIn, devex_,
                                         reference_,
                                         weights_,scaleFactor);
  } else {
    // put row of tableau in rowArray and columnArray
    model_->clpMatrix()->transposeTimes(model_,-1.0,
                              pi1,dj2,dj1);
    // get subset which have nonzero tableau elements
    model_->clpMatrix()->subsetTransposeTimes(model_,pi2,dj1,dj2);
    bool killDjs = (scaleFactor==0.0);
    if (!scaleFactor)
      scaleFactor=1.0;
    // columns
    double * weight = weights_;
    
    int number = dj1->getNumElements();
    const int * index = dj1->getIndices();
    double * updateBy = dj1->denseVector();
    double * updateBy2 = dj2->denseVector();
    
    for (int j=0;j<number;j++) {
      double thisWeight;
      double pivot;
      double pivotSquared;
      int iSequence = index[j];
      double value2 = updateBy[j];
      if (killDjs)
        updateBy[j]=0.0;
      double modification=updateBy2[j];
      updateBy2[j]=0.0;
      ClpSimplex::Status status = model_->getStatus(iSequence);
      
      if (status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
        thisWeight = weight[iSequence];
        pivot = value2*scaleFactor;
        pivotSquared = pivot * pivot;
        
        thisWeight += pivotSquared * devex_ + pivot * modification;
        if (thisWeight<TRY_NORM) {
          if (referenceIn<0.0) {
            // steepest
            thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
          } else {
            // exact
            thisWeight = referenceIn*pivotSquared;
            if (reference(iSequence))
              thisWeight += 1.0;
            thisWeight = CoinMax(thisWeight,TRY_NORM);
          }
        }
        weight[iSequence] = thisWeight;
      }
    }
  }
  dj2->setNumElements(0);
}
// Update weights for Devex
void 
01722 ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int j;
  int number=0;
  int * index;
  double * updateBy;
  // dj could be very small (or even zero - take care)
  double dj = model_->dualIn();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  // for weights update we use pivotSequence
  pivotRow = pivotSequence_;
  assert (pivotRow>=0);
  // make sure infeasibility on incoming is 0.0
  const int * pivotVariable = model_->pivotVariable();
  int sequenceIn = pivotVariable[pivotRow];
  infeasible_->zero(sequenceIn);
  // and we can see if reference
  double referenceIn=0.0;
  if (mode_!=1&&reference(sequenceIn))
    referenceIn=1.0;
  // save outgoing weight round update
  double outgoingWeight=0.0;
  int sequenceOut = model_->sequenceOut();
  if (sequenceOut>=0)
    outgoingWeight=weights_[sequenceOut];
  assert (!updates->getNumElements());
  assert (!spareColumn1->getNumElements());
  // unset in case sub flip
  pivotSequence_=-1;
  // might as well set dj to 1
  dj = -1.0;
  updates->createPacked(1,&pivotRow,&dj);
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  // put row of tableau in rowArray and columnArray
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  double * weight;
  int numberColumns = model_->numberColumns();
  // rows
  number = updates->getNumElements();
  index = updates->getIndices();
  updateBy = updates->denseVector();
  weight = weights_+numberColumns;
  
  // Devex
  assert (devex_>0.0);
  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double thisWeight = weight[iSequence];
    // row has -1 
    double pivot = - updateBy[j];
    updateBy[j]=0.0;
    double value = pivot * pivot*devex_;
    if (reference(iSequence+numberColumns))
      value += 1.0;
    weight[iSequence] = CoinMax(0.99*thisWeight,value);
  }
  
  // columns
  weight = weights_;
  
  number = spareColumn1->getNumElements();
  index = spareColumn1->getIndices();
  updateBy = spareColumn1->denseVector();
  // Devex
  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double thisWeight = weight[iSequence];
    // row has -1 
    double pivot = updateBy[j];
    updateBy[j]=0.0;
    double value = pivot * pivot*devex_;
    if (reference(iSequence))
      value += 1.0;
    weight[iSequence] = CoinMax(0.99*thisWeight,value);
  }
  // restore outgoing weight
  if (sequenceOut>=0)
    weights_[sequenceOut]=outgoingWeight;
  spareColumn2->setNumElements(0);
  //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
  // check for accuracy
  int iCheck=892;
  //printf("weight for iCheck is %g\n",weights_[iCheck]);
  int numberRows = model_->numberRows();
  //int numberColumns = model_->numberColumns();
  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
      !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
  }
#endif
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
}
// Update weights for Steepest
void 
01829 ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
             CoinIndexedVector * spareRow2,
             CoinIndexedVector * spareColumn1,
             CoinIndexedVector * spareColumn2)
{
  int j;
  int number=0;
  int * index;
  double * updateBy;
  // dj could be very small (or even zero - take care)
  double dj = model_->dualIn();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  // for weights update we use pivotSequence
  pivotRow = pivotSequence_;
  // unset in case sub flip
  pivotSequence_=-1;
  assert (pivotRow>=0);
  // make sure infeasibility on incoming is 0.0
  const int * pivotVariable = model_->pivotVariable();
  int sequenceIn = pivotVariable[pivotRow];
  infeasible_->zero(sequenceIn);
  // and we can see if reference
  double referenceIn=0.0;
  if (mode_!=1&&reference(sequenceIn))
    referenceIn=1.0;
  // save outgoing weight round update
  double outgoingWeight=0.0;
  int sequenceOut = model_->sequenceOut();
  if (sequenceOut>=0)
    outgoingWeight=weights_[sequenceOut];
  assert (!updates->getNumElements());
  assert (!spareColumn1->getNumElements());
  // update weights
  //updates->setNumElements(0);
  //spareColumn1->setNumElements(0);
  // might as well set dj to 1
  dj = -1.0;
  updates->createPacked(1,&pivotRow,&dj);
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  // put row of tableau in rowArray and columnArray
  model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
  double * weight;
  double * other = alternateWeights_->denseVector();
  int numberColumns = model_->numberColumns();
  // rows
  number = updates->getNumElements();
  index = updates->getIndices();
  updateBy = updates->denseVector();
  weight = weights_+numberColumns;
  
  // Exact
  // now update weight update array
  //alternateWeights_->scanAndPack();
  model_->factorization()->updateColumnTranspose(spareRow2,
                                     alternateWeights_);
  // get subset which have nonzero tableau elements
  model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
                                  spareColumn1,
                                  spareColumn2);
  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double thisWeight = weight[iSequence];
    // row has -1 
    double pivot = -updateBy[j];
    updateBy[j]=0.0;
    double modification = other[iSequence];
    double pivotSquared = pivot * pivot;
    
    thisWeight += pivotSquared * devex_ + pivot * modification;
    if (thisWeight<TRY_NORM) {
      if (mode_==1) {
      // steepest
      thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
      } else {
      // exact
      thisWeight = referenceIn*pivotSquared;
      if (reference(iSequence+numberColumns))
        thisWeight += 1.0;
      thisWeight = CoinMax(thisWeight,TRY_NORM);
      }
    }
    weight[iSequence] = thisWeight;
  }
  
  // columns
  weight = weights_;
  
  number = spareColumn1->getNumElements();
  index = spareColumn1->getIndices();
  updateBy = spareColumn1->denseVector();
  // Exact
  double * updateBy2 = spareColumn2->denseVector();
  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double thisWeight = weight[iSequence];
    double pivot = updateBy[j];
    updateBy[j]=0.0;
    double modification = updateBy2[j];
    updateBy2[j]=0.0;
    double pivotSquared = pivot * pivot;
    
    thisWeight += pivotSquared * devex_ + pivot * modification;
    if (thisWeight<TRY_NORM) {
      if (mode_==1) {
      // steepest
      thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
      } else {
      // exact
      thisWeight = referenceIn*pivotSquared;
      if (reference(iSequence))
        thisWeight += 1.0;
      thisWeight = CoinMax(thisWeight,TRY_NORM);
      }
    }
    weight[iSequence] = thisWeight;
  }
  // restore outgoing weight
  if (sequenceOut>=0)
    weights_[sequenceOut]=outgoingWeight;
  alternateWeights_->clear();
  spareColumn2->setNumElements(0);
  //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
  // check for accuracy
  int iCheck=892;
  //printf("weight for iCheck is %g\n",weights_[iCheck]);
  int numberRows = model_->numberRows();
  //int numberColumns = model_->numberColumns();
  for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
    if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
      !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
  }
#endif
  updates->setNumElements(0);
  spareColumn1->setNumElements(0);
}
// Returns pivot column, -1 if none
int 
01975 ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
                            CoinIndexedVector * ,
                            CoinIndexedVector * spareRow2,
                            CoinIndexedVector * spareColumn1,
                            CoinIndexedVector * spareColumn2)
{
  assert(model_);
  int iSection,j;
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  // dj could be very small (or even zero - take care)
  double dj = model_->dualIn();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  int pivotRow = model_->pivotRow();
  int anyUpdates;
  double * infeas = infeasible_->denseVector();

  // Local copy of mode so can decide what to do
  int switchType;
  if (mode_==4) 
    switchType = 5-numberSwitched_;
  else if (mode_>=10)
    switchType=3;
  else
    switchType = mode_;
  /* switchType - 
     0 - all exact devex
     1 - all steepest
     2 - some exact devex
     3 - auto some exact devex
     4 - devex
     5 - dantzig
  */
  if (updates->getNumElements()) {
    // would have to have two goes for devex, three for steepest
    anyUpdates=2;
    // add in pivot contribution
    if (pivotRow>=0) 
      updates->add(pivotRow,-dj);
  } else if (pivotRow>=0) {
    if (fabs(dj)>1.0e-15) {
      // some dj
      updates->insert(pivotRow,-dj);
      if (fabs(dj)>1.0e-6) {
      // reasonable size
      anyUpdates=1;
      } else {
      // too small
      anyUpdates=2;
      }
    } else {
      // zero dj
      anyUpdates=-1;
    }
  } else if (pivotSequence_>=0){
    // just after re-factorization
    anyUpdates=-1;
  } else {
    // sub flip - nothing to do
    anyUpdates=0;
  }

  if (anyUpdates>0) {
    model_->factorization()->updateColumnTranspose(spareRow2,updates);
    
    // put row of tableau in rowArray and columnArray
    model_->clpMatrix()->transposeTimes(model_,-1.0,
                              updates,spareColumn2,spareColumn1);
    // normal
    for (iSection=0;iSection<2;iSection++) {
      
      reducedCost=model_->djRegion(iSection);
      int addSequence;
      
      if (!iSection) {
      number = updates->getNumElements();
      index = updates->getIndices();
      updateBy = updates->denseVector();
      addSequence = model_->numberColumns();
      } else {
      number = spareColumn1->getNumElements();
      index = spareColumn1->getIndices();
      updateBy = spareColumn1->denseVector();
      addSequence = 0;
      }
      if (!model_->nonLinearCost()->lookBothWays()) {

      for (j=0;j<number;j++) {
        int iSequence = index[j];
        double value = reducedCost[iSequence];
        value -= updateBy[iSequence];
        reducedCost[iSequence] = value;
        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
        
        switch(status) {
          
        case ClpSimplex::basic:
          infeasible_->zero(iSequence+addSequence);
        case ClpSimplex::isFixed:
          break;
        case ClpSimplex::isFree:
        case ClpSimplex::superBasic:
          if (fabs(value)>FREE_ACCEPT*tolerance) {
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atUpperBound:
          if (value>tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atLowerBound:
          if (value<-tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
        }
      }
      } else {
      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
      // We can go up OR down
      for (j=0;j<number;j++) {
        int iSequence = index[j];
        double value = reducedCost[iSequence];
        value -= updateBy[iSequence];
        reducedCost[iSequence] = value;
        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
        
        switch(status) {
          
        case ClpSimplex::basic:
          infeasible_->zero(iSequence+addSequence);
        case ClpSimplex::isFixed:
          break;
        case ClpSimplex::isFree:
        case ClpSimplex::superBasic:
          if (fabs(value)>FREE_ACCEPT*tolerance) {
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atUpperBound:
          if (value>tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            // look other way - change up should be negative
            value -= nonLinear->changeUpInCost(iSequence+addSequence);
            if (value>-tolerance) {
            infeasible_->zero(iSequence+addSequence);
            } else {
            // store square in list
            if (infeas[iSequence+addSequence])
              infeas[iSequence+addSequence] = value*value; // already there
            else
              infeasible_->quickAdd(iSequence+addSequence,value*value);
            }
          }
          break;
        case ClpSimplex::atLowerBound:
          if (value<-tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            // look other way - change down should be positive
            value -= nonLinear->changeDownInCost(iSequence+addSequence);
            if (value<tolerance) {
            infeasible_->zero(iSequence+addSequence);
            } else {
            // store square in list
            if (infeas[iSequence+addSequence])
              infeas[iSequence+addSequence] = value*value; // already there
            else
              infeasible_->quickAdd(iSequence+addSequence,value*value);
            }
          }
        }
      }
      }
    }
    if (anyUpdates==2) {
      // we can zero out as will have to get pivot row
      updates->clear();
      spareColumn1->clear();
    }
    if (pivotRow>=0) {
      // make sure infeasibility on incoming is 0.0
      int sequenceIn = model_->sequenceIn();
      infeasible_->zero(sequenceIn);
    }
  }
  // make sure outgoing from last iteration okay
  int sequenceOut = model_->sequenceOut();
  if (sequenceOut>=0) {
    ClpSimplex::Status status = model_->getStatus(sequenceOut);
    double value = model_->reducedCost(sequenceOut);
    
    switch(status) {
      
    case ClpSimplex::basic:
    case ClpSimplex::isFixed:
      break;
    case ClpSimplex::isFree:
    case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) { 
      // we are going to bias towards free (but only if reasonable)
      value *= FREE_BIAS;
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
      break;
    case ClpSimplex::atUpperBound:
      if (value>tolerance) {
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
      break;
    case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
      // store square in list
      if (infeas[sequenceOut])
        infeas[sequenceOut] = value*value; // already there
      else
        infeasible_->quickAdd(sequenceOut,value*value);
      } else {
      infeasible_->zero(sequenceOut);
      }
    }
  }

  // If in quadratic re-compute all
  if (model_->algorithm()==2) {
    for (iSection=0;iSection<2;iSection++) {
      
      reducedCost=model_->djRegion(iSection);
      int addSequence;
      int iSequence;
      
      if (!iSection) {
      number = model_->numberRows();
      addSequence = model_->numberColumns();
      } else {
      number = model_->numberColumns();
      addSequence = 0;
      }
      
      if (!model_->nonLinearCost()->lookBothWays()) {
      for (iSequence=0;iSequence<number;iSequence++) {
        double value = reducedCost[iSequence];
        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
        
        switch(status) {
          
        case ClpSimplex::basic:
          infeasible_->zero(iSequence+addSequence);
        case ClpSimplex::isFixed:
          break;
        case ClpSimplex::isFree:
        case ClpSimplex::superBasic:
          if (fabs(value)>tolerance) {
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atUpperBound:
          if (value>tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atLowerBound:
          if (value<-tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
        }
      }
      } else {
      // we can go both ways
      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
      for (iSequence=0;iSequence<number;iSequence++) {
        double value = reducedCost[iSequence];
        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
        
        switch(status) {
          
        case ClpSimplex::basic:
          infeasible_->zero(iSequence+addSequence);
        case ClpSimplex::isFixed:
          break;
        case ClpSimplex::isFree:
        case ClpSimplex::superBasic:
          if (fabs(value)>tolerance) {
            // we are going to bias towards free (but only if reasonable)
            value *= FREE_BIAS;
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            infeasible_->zero(iSequence+addSequence);
          }
          break;
        case ClpSimplex::atUpperBound:
          if (value>tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            // look other way - change up should be negative
            value -= nonLinear->changeUpInCost(iSequence+addSequence);
            if (value>-tolerance) {
            infeasible_->zero(iSequence+addSequence);
            } else {
            // store square in list
            if (infeas[iSequence+addSequence])
              infeas[iSequence+addSequence] = value*value; // already there
            else
              infeasible_->quickAdd(iSequence+addSequence,value*value);
            }
          }
          break;
        case ClpSimplex::atLowerBound:
          if (value<-tolerance) {
            // store square in list
            if (infeas[iSequence+addSequence])
            infeas[iSequence+addSequence] = value*value; // already there
            else
            infeasible_->quickAdd(iSequence+addSequence,value*value);
          } else {
            // look other way - change down should be positive
            value -= nonLinear->changeDownInCost(iSequence+addSequence);
            if (value<tolerance) {
            infeasible_->zero(iSequence+addSequence);
            } else {
            // store square in list
            if (infeas[iSequence+addSequence])
              infeas[iSequence+addSequence] = value*value; // already there
            else
              infeasible_->quickAdd(iSequence+addSequence,value*value);
            }
          }
        }
      }
      }
    }
  }
  // See what sort of pricing
  int numberWanted=10;
  number = infeasible_->getNumElements();
  int numberColumns = model_->numberColumns();
  if (switchType==5) {
    // we can zero out
    updates->clear();
    spareColumn1->clear();
    pivotSequence_=-1;
    pivotRow=-1;
    // See if to switch
    int numberRows = model_->numberRows();
    // ratio is done on number of columns here
    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
    double ratio = static_cast<double> (sizeFactorization_)/static_cast<double> (numberRows);
    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
    if (ratio<0.1) {
      numberWanted = CoinMax(100,number/200);
    } else if (ratio<0.3) {
      numberWanted = CoinMax(500,number/40);
    } else if (ratio<0.5||mode_==5) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/30);
    } else if (mode_!=5) {
      switchType=4;
      // initialize
      numberSwitched_++;
      // Make sure will re-do
      delete [] weights_;
      weights_=NULL;
      saveWeights(model_,4);
      printf("switching to devex %d nel ratio %g\n",sizeFactorization_,ratio);
      updates->clear();
    }
    if (model_->numberIterations()%1000==0)
      printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
  }
  if(switchType==4) {
    // Still in devex mode
    int numberRows = model_->numberRows();
    // ratio is done on number of rows here
    double ratio = static_cast<double> (sizeFactorization_)/static_cast<double> (numberRows);
    // Go to steepest if lot of iterations?
    if (ratio<1.0) {
      numberWanted = CoinMax(2000,number/20);
    } else if (ratio<5.0) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/20);
    } else {
      // we can zero out
      updates->clear();
      spareColumn1->clear();
      switchType=3;
      // initialize
      pivotSequence_=-1;
      pivotRow=-1;
      numberSwitched_++;
      // Make sure will re-do
      delete [] weights_;
      weights_=NULL;
      saveWeights(model_,4);
      printf("switching to exact %d nel ratio %g\n",sizeFactorization_,ratio);
      updates->clear();
    }
    if (model_->numberIterations()%1000==0)
      printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
  } 
  if (switchType<4) {
    if (switchType<2 ) {
      numberWanted = number+1;
    } else if (switchType==2) {
      numberWanted = CoinMax(2000,number/8);
    } else {
      double ratio = static_cast<double> (sizeFactorization_)/static_cast<double> (model_->numberRows());
      if (ratio<1.0) {
      numberWanted = CoinMax(2000,number/20);
      } else if (ratio<5.0) {
      numberWanted = CoinMax(2000,number/10);
      numberWanted = CoinMax(numberWanted,numberColumns/20);
      } else if (ratio<10.0) {
      numberWanted = CoinMax(2000,number/8);
      numberWanted = CoinMax(numberWanted,numberColumns/20);
      } else {
      ratio = number * (ratio/80.0);
      if (ratio>number) {
        numberWanted=number+1;
      } else {
        numberWanted = CoinMax(2000,static_cast<int> (ratio));
        numberWanted = CoinMax(numberWanted,numberColumns/10);
      }
      }
    }
  }
  // for weights update we use pivotSequence
  pivotRow = pivotSequence_;
  // unset in case sub flip
  pivotSequence_=-1;
  if (pivotRow>=0) {
    // make sure infeasibility on incoming is 0.0
    const int * pivotVariable = model_->pivotVariable();
    int sequenceIn = pivotVariable[pivotRow];
    infeasible_->zero(sequenceIn);
    // and we can see if reference
    double referenceIn=0.0;
    if (switchType!=1&&reference(sequenceIn))
      referenceIn=1.0;
    // save outgoing weight round update
    double outgoingWeight=0.0;
    if (sequenceOut>=0)
      outgoingWeight=weights_[sequenceOut];
    // update weights
    if (anyUpdates!=1) {
      updates->setNumElements(0);
      spareColumn1->setNumElements(0);
      // might as well set dj to 1
      dj = 1.0;
      updates->insert(pivotRow,-dj);
      model_->factorization()->updateColumnTranspose(spareRow2,updates);
      // put row of tableau in rowArray and columnArray
      model_->clpMatrix()->transposeTimes(model_,-1.0,
                                updates,spareColumn2,spareColumn1);
    }
    double * weight;
    double * other = alternateWeights_->denseVector();
    int numberColumns = model_->numberColumns();
    double scaleFactor = -1.0/dj; // as formula is with 1.0
    // rows
    number = updates->getNumElements();
    index = updates->getIndices();
    updateBy = updates->denseVector();
    weight = weights_+numberColumns;
      
    if (switchType<4) {
      // Exact
      // now update weight update array
      model_->factorization()->updateColumnTranspose(spareRow2,
                                         alternateWeights_);
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = updateBy[iSequence]*scaleFactor;
      updateBy[iSequence]=0.0;
      double modification = other[iSequence];
      double pivotSquared = pivot * pivot;

      thisWeight += pivotSquared * devex_ + pivot * modification;
      if (thisWeight<TRY_NORM) {
        if (switchType==1) {
          // steepest
          thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
        } else {
          // exact
          thisWeight = referenceIn*pivotSquared;
          if (reference(iSequence+numberColumns))
            thisWeight += 1.0;
          thisWeight = CoinMax(thisWeight,TRY_NORM);
        }
      }
      weight[iSequence] = thisWeight;
      }
    } else if (switchType==4) {
      // Devex
      assert (devex_>0.0);
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = updateBy[iSequence]*scaleFactor;
      updateBy[iSequence]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence+numberColumns))
        value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
      }
    }
    
    // columns
    weight = weights_;
    
    scaleFactor = -scaleFactor;
    
    number = spareColumn1->getNumElements();
    index = spareColumn1->getIndices();
    updateBy = spareColumn1->denseVector();
    if (switchType<4) {
      // Exact
      // get subset which have nonzero tableau elements
      model_->clpMatrix()->subsetTransposeTimes(model_,alternateWeights_,
                                    spareColumn1,
                                    spareColumn2);
      double * updateBy2 = spareColumn2->denseVector();
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      double pivot = updateBy[iSequence]*scaleFactor;
      updateBy[iSequence]=0.0;
      double modification = updateBy2[j];
      updateBy2[j]=0.0;
      double pivotSquared = pivot * pivot;

      thisWeight += pivotSquared * devex_ + pivot * modification;
      if (thisWeight<TRY_NORM) {
        if (switchType==1) {
          // steepest
          thisWeight = CoinMax(TRY_NORM,ADD_ONE+pivotSquared);
        } else {
          // exact
          thisWeight = referenceIn*pivotSquared;
          if (reference(iSequence))
            thisWeight += 1.0;
          thisWeight = CoinMax(thisWeight,TRY_NORM);
        }
      }
      weight[iSequence] = thisWeight;
      }
    } else if (switchType==4) {
      // Devex
      for (j=0;j<number;j++) {
      int iSequence = index[j];
      double thisWeight = weight[iSequence];
      // row has -1 
      double pivot = updateBy[iSequence]*scaleFactor;
      updateBy[iSequence]=0.0;
      double value = pivot * pivot*devex_;
      if (reference(iSequence))
        value += 1.0;
      weight[iSequence] = CoinMax(0.99*thisWeight,value);
      }
    }
    // restore outgoing weight
    if (sequenceOut>=0)
      weights_[sequenceOut]=outgoingWeight;
    alternateWeights_->clear();
    spareColumn2->setNumElements(0);
    //#define SOME_DEBUG_1
#ifdef SOME_DEBUG_1
    // check for accuracy
    int iCheck=892;
    //printf("weight for iCheck is %g\n",weights_[iCheck]);
    int numberRows = model_->numberRows();
    //int numberColumns = model_->numberColumns();
    for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) {
      if (model_->getStatus(iCheck)!=ClpSimplex::basic&&
        !model_->getStatus(iCheck)!=ClpSimplex::isFixed)
      checkAccuracy(iCheck,1.0e-1,updates,spareRow2);
    }
#endif
    updates->setNumElements(0);
    spareColumn1->setNumElements(0);
  }

  // update of duals finished - now do pricing


  double bestDj = 1.0e-30;
  int bestSequence=-1;

  int i,iSequence;
  index = infeasible_->getIndices();
  number = infeasible_->getNumElements();
  if(model_->numberIterations()<model_->lastBadIteration()+200) {
    // we can't really trust infeasibilities if there is dual error
    double checkTolerance = 1.0e-8;
    if (!model_->factorization()->pivots())
      checkTolerance = 1.0e-6;
    if (model_->largestDualError()>checkTolerance)
      tolerance *= model_->largestDualError()/checkTolerance;
    // But cap
    tolerance = CoinMin(1000.0,tolerance);
  }
#ifdef CLP_DEBUG
  if (model_->numberDualInfeasibilities()==1) 
    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
         model_->largestDualError(),model_,model_->messageHandler(),
         number);
#endif
  // stop last one coming immediately
  double saveOutInfeasibility=0.0;
  if (sequenceOut>=0) {
    saveOutInfeasibility = infeas[sequenceOut];
    infeas[sequenceOut]=0.0;
  }
  tolerance *= tolerance; // as we are using squares

  int iPass;
  // Setup two passes
  int start[4];
  start[1]=number;
  start[2]=0;
  double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
  start[0]=static_cast<int> (dstart);
  start[3]=start[0];
  //double largestWeight=0.0;
  //double smallestWeight=1.0e100;
  for (iPass=0;iPass<2;iPass++) {
    int end = start[2*iPass+1];
    if (switchType<5) {
      for (i=start[2*iPass];i<end;i++) {
      iSequence = index[i];
      double value = infeas[iSequence];
      if (value>tolerance) {
        double weight = weights_[iSequence];
        //weight=1.0;
        if (value>bestDj*weight) {
          // check flagged variable and correct dj
          if (!model_->flagged(iSequence)) {
            bestDj=value/weight;
            bestSequence = iSequence;
          } else {
            // just to make sure we don't exit before got something
            numberWanted++;
          }
        }
      }
      numberWanted--;
      if (!numberWanted)
        break;
      }
    } else {
      // Dantzig
      for (i=start[2*iPass];i<end;i++) {
      iSequence = index[i];
      double value = infeas[iSequence];
      if (value>tolerance) {
        if (value>bestDj) {
          // check flagged variable and correct dj
          if (!model_->flagged(iSequence)) {
            bestDj=value;
            bestSequence = iSequence;
          } else {
            // just to make sure we don't exit before got something
            numberWanted++;
          }
        }
      }
      numberWanted--;
      if (!numberWanted)
        break;
      }
    }
    if (!numberWanted)
      break;
  }
  if (sequenceOut>=0) {
    infeas[sequenceOut]=saveOutInfeasibility;
  }
  /*if (model_->numberIterations()%100==0)
    printf("%d best %g\n",bestSequence,bestDj);*/
  reducedCost=model_->djRegion();
  model_->clpMatrix()->setSavedBestSequence(bestSequence);
  if (bestSequence>=0)
    model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
  
#ifdef CLP_DEBUG
  if (bestSequence>=0) {
    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
      assert(model_->reducedCost(bestSequence)<0.0);
    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
      assert(model_->reducedCost(bestSequence)>0.0);
  }
#endif
  return bestSequence;
}
// Called when maximum pivots changes
void 
02757 ClpPrimalColumnSteepest::maximumPivotsChanged()
{
  if (alternateWeights_&&
      alternateWeights_->capacity()!=model_->numberRows()+
      model_->factorization()->maximumPivots()) {
    delete alternateWeights_;
    alternateWeights_ = new CoinIndexedVector();
    // enough space so can use it for factorization
    alternateWeights_->reserve(model_->numberRows()+
                           model_->factorization()->maximumPivots());
  }
}
/* 
   1) before factorization
   2) after factorization
   3) just redo infeasibilities
   4) restore weights
   5) at end of values pass (so need initialization)
*/
void 
02777 ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
{
  model_ = model;
  if (mode_==4||mode_==5) {
    if (mode==1&&!weights_) 
      numberSwitched_=0; // Reset
  }
  // alternateWeights_ is defined as indexed but is treated oddly
  // at times
  int numberRows = model_->numberRows();
  int numberColumns = model_->numberColumns();
  const int * pivotVariable = model_->pivotVariable();
  bool doInfeasibilities=true;
  if (mode==1) {
    if(weights_) {
      // Check if size has changed
      if (infeasible_->capacity()==numberRows+numberColumns&&
      alternateWeights_->capacity()==numberRows+
        model_->factorization()->maximumPivots()) {
      //alternateWeights_->clear();
      if (pivotSequence_>=0&&pivotSequence_<numberRows) {
        // save pivot order
        CoinMemcpyN(pivotVariable,
             numberRows,alternateWeights_->getIndices());
        // change from pivot row number to sequence number
        pivotSequence_=pivotVariable[pivotSequence_];
      } else {
        pivotSequence_ = -1;
      }
      state_=1;
      } else {
      // size has changed - clear everything
      delete [] weights_;
      weights_=NULL;
      delete infeasible_;
      infeasible_=NULL;
      delete alternateWeights_;
      alternateWeights_=NULL;
      delete [] savedWeights_;
      savedWeights_=NULL;
      delete [] reference_;
      reference_=NULL;
      state_=-1;
      pivotSequence_=-1;
      }
    }
  } else if (mode==2||mode==4||mode==5) {
    // restore
    if (!weights_||state_==-1||mode==5) {
      // Partial is only allowed with certain types of matrix
      if ((mode_!=4&&mode_!=5)||numberSwitched_||!model_->clpMatrix()->canDoPartialPricing()) {
      // initialize weights
      delete [] weights_;
      delete alternateWeights_;
      weights_ = new double[numberRows+numberColumns];
      alternateWeights_ = new CoinIndexedVector();
      // enough space so can use it for factorization
      alternateWeights_->reserve(numberRows+
                           model_->factorization()->maximumPivots());
      initializeWeights();
      // create saved weights 
      delete [] savedWeights_;
      savedWeights_ = CoinCopyOfArray(weights_,numberRows+numberColumns);
      // just do initialization
      mode=3;
      } else {
      // Partial pricing
      // use region as somewhere to save non-fixed slacks
      // set up infeasibilities
      if (!infeasible_) {
        infeasible_ = new CoinIndexedVector();
        infeasible_->reserve(numberColumns+numberRows);
      }
      infeasible_->clear();
      int number = model_->numberRows() + model_->numberColumns();
      int iSequence;
      int numberLook=0;
      int * which = infeasible_->getIndices();
      for (iSequence=model_->numberColumns();iSequence<number;iSequence++) {
        ClpSimplex::Status status = model_->getStatus(iSequence);
        if (status!=ClpSimplex::isFixed)
          which[numberLook++]=iSequence;
      }
      infeasible_->setNumElements(numberLook);
      doInfeasibilities=false;
      }
      savedPivotSequence_=-2;
      savedSequenceOut_ = -2;
      
    } else {
      if (mode!=4) {
      // save
      CoinMemcpyN(weights_,(numberRows+numberColumns),savedWeights_);
      savedPivotSequence_= pivotSequence_;
      savedSequenceOut_ = model_->sequenceOut();
      } else {
      // restore
      CoinMemcpyN(savedWeights_,(numberRows+numberColumns),weights_);
      // was - but I think should not be
      //pivotSequence_= savedPivotSequence_;
      //model_->setSequenceOut(savedSequenceOut_); 
      pivotSequence_= -1;
      model_->setSequenceOut(-1);
      // indices are wrong so clear by hand
      //alternateWeights_->clear();
      CoinZeroN(alternateWeights_->denseVector(),
              alternateWeights_->capacity());
      alternateWeights_->setNumElements(0);
      }
    }
    state_=0;
    // set up infeasibilities
    if (!infeasible_) {
      infeasible_ = new CoinIndexedVector();
      infeasible_->reserve(numberColumns+numberRows);
    }
  }
  if (mode>=2&&mode!=5) {
    if (mode!=3) {
      if (pivotSequence_>=0) {
      // restore pivot row
      int iRow;
      // permute alternateWeights
      double * temp = model_->rowArray(3)->denseVector();;
      double * work = alternateWeights_->denseVector();
      int * savePivotOrder = model_->rowArray(3)->getIndices();
      int * oldPivotOrder = alternateWeights_->getIndices();
      for (iRow=0;iRow<numberRows;iRow++) {
        int iPivot=oldPivotOrder[iRow];
        temp[iPivot]=work[iRow];
          savePivotOrder[iRow]=iPivot;
      }
      int number=0;
      int found=-1;
      int * which = oldPivotOrder;
      // find pivot row and re-create alternateWeights
      for (iRow=0;iRow<numberRows;iRow++) {
        int iPivot=pivotVariable[iRow];
        if (iPivot==pivotSequence_) 
          found=iRow;
        work[iRow]=temp[iPivot];
        if (work[iRow])
          which[number++]=iRow;
      }
      alternateWeights_->setNumElements(number);
#ifdef CLP_DEBUG
      // Can happen but I should clean up
      assert(found>=0);
#endif
      pivotSequence_ = found;
      for (iRow=0;iRow<numberRows;iRow++) {
        int iPivot=savePivotOrder[iRow];
        temp[iPivot]=0.0;
        }
      } else {
      // Just clean up
      if (alternateWeights_)
        alternateWeights_->clear();
      }
    }
    // Save size of factorization
    if (!model->factorization()->pivots())
      sizeFactorization_ = model_->factorization()->numberElements();
    if(!doInfeasibilities)
      return; // don't disturb infeasibilities
    infeasible_->clear();
    double tolerance=model_->currentDualTolerance();
    int number = model_->numberRows() + model_->numberColumns();
    int iSequence;
   
    double * reducedCost = model_->djRegion();
      
    if (!model_->nonLinearCost()->lookBothWays()) {
      for (iSequence=0;iSequence<number;iSequence++) {
      double value = reducedCost[iSequence];
      ClpSimplex::Status status = model_->getStatus(iSequence);

      switch(status) {
        
      case ClpSimplex::basic:
      case ClpSimplex::isFixed:
        break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
        if (fabs(value)>FREE_ACCEPT*tolerance) { 
          // we are going to bias towards free (but only if reasonable)
          value *= FREE_BIAS;
          // store square in list
          infeasible_->quickAdd(iSequence,value*value);
        }
        break;
      case ClpSimplex::atUpperBound:
        if (value>tolerance) {
          infeasible_->quickAdd(iSequence,value*value);
        }
        break;
      case ClpSimplex::atLowerBound:
        if (value<-tolerance) {
          infeasible_->quickAdd(iSequence,value*value);
        }
      }
      }
    } else {
      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
      // can go both ways
      for (iSequence=0;iSequence<number;iSequence++) {
      double value = reducedCost[iSequence];
      ClpSimplex::Status status = model_->getStatus(iSequence);

      switch(status) {
        
      case ClpSimplex::basic:
      case ClpSimplex::isFixed:
        break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
        if (fabs(value)>FREE_ACCEPT*tolerance) { 
          // we are going to bias towards free (but only if reasonable)
          value *= FREE_BIAS;
          // store square in list
          infeasible_->quickAdd(iSequence,value*value);
        }
        break;
      case ClpSimplex::atUpperBound:
        if (value>tolerance) {
          infeasible_->quickAdd(iSequence,value*value);
        } else {
          // look other way - change up should be negative
          value -= nonLinear->changeUpInCost(iSequence);
          if (value<-tolerance) {
            // store square in list
            infeasible_->quickAdd(iSequence,value*value);
          }
        }
        break;
      case ClpSimplex::atLowerBound:
        if (value<-tolerance) {
          infeasible_->quickAdd(iSequence,value*value);
        } else {
          // look other way - change down should be positive
          value -= nonLinear->changeDownInCost(iSequence);
          if (value>tolerance) {
            // store square in list
            infeasible_->quickAdd(iSequence,value*value);
          }
        }
      }
      }
    }
  }
}
// Gets rid of last update
void 
03030 ClpPrimalColumnSteepest::unrollWeights()
{
  if ((mode_==4||mode_==5)&&!numberSwitched_)
    return;
  double * saved = alternateWeights_->denseVector();
  int number = alternateWeights_->getNumElements();
  int * which = alternateWeights_->getIndices();
  int i;
  for (i=0;i<number;i++) {
    int iRow = which[i];
    weights_[iRow]=saved[iRow];
    saved[iRow]=0.0;
  }
  alternateWeights_->setNumElements(0);
}
  
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
03049 ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
{
  if (CopyData) {
    return new ClpPrimalColumnSteepest(*this);
  } else {
    return new ClpPrimalColumnSteepest();
  }
}
void
03058 ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
{
  // Local copy of mode so can decide what to do
  int switchType = mode_;
  if (mode_==4&&numberSwitched_)
    switchType=3;
  else if (mode_==4||mode_==5)
    return;
  int number=input->getNumElements();
  int * which = input->getIndices();
  double * work = input->denseVector();
  int newNumber = 0;
  int * newWhich = alternateWeights_->getIndices();
  double * newWork = alternateWeights_->denseVector();
  int i;
  int sequenceIn = model_->sequenceIn();
  int sequenceOut = model_->sequenceOut();
  const int * pivotVariable = model_->pivotVariable();

  int pivotRow = model_->pivotRow();
  pivotSequence_ = pivotRow;

  devex_ =0.0;
  // Can't create alternateWeights_ as packed as needed unpacked
  if (!input->packedMode()) {
    if (pivotRow>=0) {
      if (switchType==1) {
      for (i=0;i<number;i++) {
        int iRow = which[i];
        devex_ += work[iRow]*work[iRow];
        newWork[iRow] = -2.0*work[iRow];
      }
      newWork[pivotRow] = -2.0*CoinMax(devex_,0.0);
      devex_+=ADD_ONE;
      weights_[sequenceOut]=1.0+ADD_ONE;
      CoinMemcpyN(which,number,newWhich);
      alternateWeights_->setNumElements(number);
      } else {
      if ((mode_!=4&&mode_!=5)||numberSwitched_>1) {
        for (i=0;i<number;i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (reference(iPivot)) {
            devex_ += work[iRow]*work[iRow];
            newWork[iRow] = -2.0*work[iRow];
            newWhich[newNumber++]=iRow;
          }
        }
        if (!newWork[pivotRow]&&devex_>0.0)
          newWhich[newNumber++]=pivotRow; // add if not already in
        newWork[pivotRow] = -2.0*CoinMax(devex_,0.0);
      } else {
        for (i=0;i<number;i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (reference(iPivot)) 
            devex_ += work[iRow]*work[iRow];
        }
      }
      if (reference(sequenceIn)) {
        devex_+=1.0;
      } else {
      }
      if (reference(sequenceOut)) {
        weights_[sequenceOut]=1.0+1.0;
      } else {
        weights_[sequenceOut]=1.0;
      }
      alternateWeights_->setNumElements(newNumber);
      }
    } else {
      if (switchType==1) {
      for (i=0;i<number;i++) {
        int iRow = which[i];
        devex_ += work[iRow]*work[iRow];
      }
      devex_ += ADD_ONE;
      } else {
      for (i=0;i<number;i++) {
        int iRow = which[i];
        int iPivot = pivotVariable[iRow];
        if (reference(iPivot)) {
          devex_ += work[iRow]*work[iRow];
        }
      }
      if (reference(sequenceIn)) 
        devex_+=1.0;
      }
    }
  } else {
    // packed input
    if (pivotRow>=0) {
      if (switchType==1) {
      for (i=0;i<number;i++) {
        int iRow = which[i];
        devex_ += work[i]*work[i];
        newWork[iRow] = -2.0*work[i];
      }
      newWork[pivotRow] = -2.0*CoinMax(devex_,0.0);
      devex_+=ADD_ONE;
      weights_[sequenceOut]=1.0+ADD_ONE;
      CoinMemcpyN(which,number,newWhich);
      alternateWeights_->setNumElements(number);
      } else {
      if ((mode_!=4&&mode_!=5)||numberSwitched_>1) {
        for (i=0;i<number;i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (reference(iPivot)) {
            devex_ += work[i]*work[i];
            newWork[iRow] = -2.0*work[i];
            newWhich[newNumber++]=iRow;
          }
        }
        if (!newWork[pivotRow]&&devex_>0.0)
          newWhich[newNumber++]=pivotRow; // add if not already in
        newWork[pivotRow] = -2.0*CoinMax(devex_,0.0);
      } else {
        for (i=0;i<number;i++) {
          int iRow = which[i];
          int iPivot = pivotVariable[iRow];
          if (reference(iPivot)) 
            devex_ += work[i]*work[i];
        }
      }
      if (reference(sequenceIn)) {
        devex_+=1.0;
      } else {
      }
      if (reference(sequenceOut)) {
        weights_[sequenceOut]=1.0+1.0;
      } else {
        weights_[sequenceOut]=1.0;
      }
      alternateWeights_->setNumElements(newNumber);
      }
    } else {
      if (switchType==1) {
      for (i=0;i<number;i++) {
        devex_ += work[i]*work[i];
      }
      devex_ += ADD_ONE;
      } else {
      for (i=0;i<number;i++) {
        int iRow = which[i];
        int iPivot = pivotVariable[iRow];
        if (reference(iPivot)) {
          devex_ += work[i]*work[i];
        }
      }
      if (reference(sequenceIn)) 
        devex_+=1.0;
      }
    }
  }
  double oldDevex = weights_[sequenceIn];
#ifdef CLP_DEBUG
  if ((model_->messageHandler()->logLevel()&32))
    printf("old weight %g, new %g\n",oldDevex,devex_);
#endif
  double check = CoinMax(devex_,oldDevex)+0.1;
  weights_[sequenceIn] = devex_;
  double testValue=0.1;
  if (mode_==4&&numberSwitched_==1)
    testValue = 0.5;
  if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
#ifdef CLP_DEBUG
    if ((model_->messageHandler()->logLevel()&48)==16)
      printf("old weight %g, new %g\n",oldDevex,devex_);
#endif
    //printf("old weight %g, new %g\n",oldDevex,devex_);
    testValue=0.99;
    if (mode_==1) 
      testValue=1.01e1; // make unlikely to do if steepest
    else if (mode_==4&&numberSwitched_==1)
      testValue=0.9;
    double difference = fabs(devex_-oldDevex);
    if ( difference > testValue * check ) {
      // need to redo
      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
                              *model_->messagesPointer())
                                <<oldDevex<<devex_
                                <<CoinMessageEol;
      initializeWeights();
    }
  }
  if (pivotRow>=0) {
    // set outgoing weight here
    weights_[model_->sequenceOut()]=devex_/(model_->alpha()*model_->alpha());
  }
}
// Checks accuracy - just for debug
void 
03251 ClpPrimalColumnSteepest::checkAccuracy(int sequence,
                               double relativeTolerance,
                               CoinIndexedVector * rowArray1,
                               CoinIndexedVector * rowArray2)
{
  if ((mode_==4||mode_==5)&&!numberSwitched_)
    return;
  model_->unpack(rowArray1,sequence);
  model_->factorization()->updateColumn(rowArray2,rowArray1);
  int number=rowArray1->getNumElements();
  int * which = rowArray1->getIndices();
  double * work = rowArray1->denseVector();
  const int * pivotVariable = model_->pivotVariable();
  
  double devex =0.0;
  int i;

  if (mode_==1) {
    for (i=0;i<number;i++) {
      int iRow = which[i];
      devex += work[iRow]*work[iRow];
      work[iRow]=0.0;
    }
    devex += ADD_ONE;
  } else {
    for (i=0;i<number;i++) {
      int iRow = which[i];
      int iPivot = pivotVariable[iRow];
      if (reference(iPivot)) {
      devex += work[iRow]*work[iRow];
      }
      work[iRow]=0.0;
    }
    if (reference(sequence)) 
      devex+=1.0;
  }

  double oldDevex = weights_[sequence];
  double check = CoinMax(devex,oldDevex);;
  if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
    printf("check %d old weight %g, new %g\n",sequence,oldDevex,devex);
    // update so won't print again
    weights_[sequence]=devex;
  }
  rowArray1->setNumElements(0);
}

// Initialize weights
void 
03300 ClpPrimalColumnSteepest::initializeWeights()
{
  int numberRows = model_->numberRows();
  int numberColumns = model_->numberColumns();
  int number = numberRows + numberColumns;
  int iSequence;
  if (mode_!=1) {
    // initialize to 1.0 
    // and set reference framework
    if (!reference_) {
      int nWords = (number+31)>>5;
      reference_ = new unsigned int[nWords];
      CoinZeroN(reference_,nWords);
    }
    
    for (iSequence=0;iSequence<number;iSequence++) {
      weights_[iSequence]=1.0;
      if (model_->getStatus(iSequence)==ClpSimplex::basic) {
      setReference(iSequence,false);
      } else {
      setReference(iSequence,true);
      }
    }
  } else {
    CoinIndexedVector * temp = new CoinIndexedVector();
    temp->reserve(numberRows+
              model_->factorization()->maximumPivots());
    double * array = alternateWeights_->denseVector();
    int * which = alternateWeights_->getIndices();
      
    for (iSequence=0;iSequence<number;iSequence++) {
      weights_[iSequence]=1.0+ADD_ONE;
      if (model_->getStatus(iSequence)!=ClpSimplex::basic&&
        model_->getStatus(iSequence)!=ClpSimplex::isFixed) {
      model_->unpack(alternateWeights_,iSequence);
      double value=ADD_ONE;
      model_->factorization()->updateColumn(temp,alternateWeights_);
      int number = alternateWeights_->getNumElements();     int j;
      for (j=0;j<number;j++) {
        int iRow=which[j];
        value+=array[iRow]*array[iRow];
        array[iRow]=0.0;
      }
      alternateWeights_->setNumElements(0);
      weights_[iSequence] = value;
      }
    }
    delete temp;
  }
}
// Gets rid of all arrays
void 
03352 ClpPrimalColumnSteepest::clearArrays()
{
  if (persistence_==normal) {
    delete [] weights_;
    weights_=NULL;
    delete infeasible_;
    infeasible_ = NULL;
    delete alternateWeights_;
    alternateWeights_ = NULL;
    delete [] savedWeights_;
    savedWeights_ = NULL;
    delete [] reference_;
    reference_ = NULL;
  }
  pivotSequence_=-1;
  state_ = -1;
  savedPivotSequence_ = -1;
  savedSequenceOut_ = -1;  
  devex_ = 0.0;
}
// Returns true if would not find any column
bool 
03374 ClpPrimalColumnSteepest::looksOptimal() const
{
  if (looksOptimal_)
    return true; // user overrode
  //**** THIS MUST MATCH the action coding above
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  if(model_->numberIterations()<model_->lastBadIteration()+200) {
    // we can't really trust infeasibilities if there is dual error
    double checkTolerance = 1.0e-8;
    if (!model_->factorization()->pivots())
      checkTolerance = 1.0e-6;
    if (model_->largestDualError()>checkTolerance)
      tolerance *= model_->largestDualError()/checkTolerance;
    // But cap
    tolerance = CoinMin(1000.0,tolerance);
  }
  int number = model_->numberRows() + model_->numberColumns();
  int iSequence;
  
  double * reducedCost = model_->djRegion();
  int numberInfeasible=0;
  if (!model_->nonLinearCost()->lookBothWays()) {
    for (iSequence=0;iSequence<number;iSequence++) {
      double value = reducedCost[iSequence];
      ClpSimplex::Status status = model_->getStatus(iSequence);
      
      switch(status) {

      case ClpSimplex::basic:
      case ClpSimplex::isFixed:
      break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) 
        numberInfeasible++;
      break;
      case ClpSimplex::atUpperBound:
      if (value>tolerance) 
        numberInfeasible++;
      break;
      case ClpSimplex::atLowerBound:
      if (value<-tolerance) 
        numberInfeasible++;
      }
    }
  } else {
    ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
    // can go both ways
    for (iSequence=0;iSequence<number;iSequence++) {
      double value = reducedCost[iSequence];
      ClpSimplex::Status status = model_->getStatus(iSequence);
      
      switch(status) {

      case ClpSimplex::basic:
      case ClpSimplex::isFixed:
      break;
      case ClpSimplex::isFree:
      case ClpSimplex::superBasic:
      if (fabs(value)>FREE_ACCEPT*tolerance) 
        numberInfeasible++;
      break;
      case ClpSimplex::atUpperBound:
      if (value>tolerance) {
        numberInfeasible++;
      } else {
        // look other way - change up should be negative
        value -= nonLinear->changeUpInCost(iSequence);
        if (value<-tolerance) 
          numberInfeasible++;
      }
      break;
      case ClpSimplex::atLowerBound:
      if (value<-tolerance) {
        numberInfeasible++;
      } else {
        // look other way - change down should be positive
        value -= nonLinear->changeDownInCost(iSequence);
        if (value>tolerance) 
          numberInfeasible++;
      }
      }
    }
  }
  return numberInfeasible==0;
}
/* Returns number of extra columns for sprint algorithm - 0 means off.
   Also number of iterations before recompute
*/
int 
03469 ClpPrimalColumnSteepest::numberSprintColumns(int & numberIterations) const
{
  numberIterations=0;
  int numberAdd=0;
  if (!numberSwitched_&&mode_>=10) {
    numberIterations = CoinMin(2000,model_->numberRows()/5);
    numberIterations = CoinMax(numberIterations,model_->factorizationFrequency());
    numberIterations = CoinMax(numberIterations,500);
    if (mode_==10) {
      numberAdd=CoinMax(300,model_->numberColumns()/10);
      numberAdd=CoinMax(numberAdd,model_->numberRows()/5);
      // fake all
      //numberAdd=1000000;
      numberAdd = CoinMin(numberAdd,model_->numberColumns());
    } else {
      abort();
    }
  }
  return numberAdd;
}
// Switch off sprint idea
void 
03491 ClpPrimalColumnSteepest::switchOffSprint()
{
  numberSwitched_=10;
}
// Update djs doing partial pricing (dantzig)
int 
03497 ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
                              CoinIndexedVector * spareRow2,
                              int numberWanted,
                              int numberLook)
{
  int number=0;
  int * index;
  double * updateBy;
  double * reducedCost;
  double saveTolerance = model_->currentDualTolerance();
  double tolerance=model_->currentDualTolerance();
  // we can't really trust infeasibilities if there is dual error
  // this coding has to mimic coding in checkDualSolution
  double error = CoinMin(1.0e-2,model_->largestDualError());
  // allow tolerance at least slightly bigger than standard
  tolerance = tolerance +  error;
  if(model_->numberIterations()<model_->lastBadIteration()+200) {
    // we can't really trust infeasibilities if there is dual error
    double checkTolerance = 1.0e-8;
    if (!model_->factorization()->pivots())
      checkTolerance = 1.0e-6;
    if (model_->largestDualError()>checkTolerance)
      tolerance *= model_->largestDualError()/checkTolerance;
    // But cap
    tolerance = CoinMin(1000.0,tolerance);
  }
  if (model_->factorization()->pivots()&&model_->numberPrimalInfeasibilities())
    tolerance = CoinMax(tolerance,1.0e-10*model_->infeasibilityCost());
  // So partial pricing can use
  model_->setCurrentDualTolerance(tolerance);
  model_->factorization()->updateColumnTranspose(spareRow2,updates);
  int numberColumns = model_->numberColumns();

  // Rows
  reducedCost=model_->djRegion(0);
  int addSequence;
    
  number = updates->getNumElements();
  index = updates->getIndices();
  updateBy = updates->denseVector();
  addSequence = numberColumns;
  int j;  
  double * duals = model_->dualRowSolution();
  for (j=0;j<number;j++) {
    int iSequence = index[j];
    double value = duals[iSequence];
    value -= updateBy[j];
    updateBy[j]=0.0;
    duals[iSequence] = value;
  }
  //#define CLP_DEBUG
#ifdef CLP_DEBUG
  // check duals
  {
    int numberRows = model_->numberRows();
    //work space
    CoinIndexedVector arrayVector;
    arrayVector.reserve(numberRows+1000);
    CoinIndexedVector workSpace;
    workSpace.reserve(numberRows+1000);
    
    
    int iRow;
    double * array = arrayVector.denseVector();
    int * index = arrayVector.getIndices();
    int number=0;
    int * pivotVariable = model_->pivotVariable();
    double * cost = model_->costRegion();
    for (iRow=0;iRow<numberRows;iRow++) {
      int iPivot=pivotVariable[iRow];
      double value = cost[iPivot];
      if (value) {
      array[iRow]=value;
      index[number++]=iRow;
      }
    }
    arrayVector.setNumElements(number);
    // Extended duals before "updateTranspose"
    model_->clpMatrix()->dualExpanded(model_,&arrayVector,NULL,0);
    
    // Btran basic costs
    model_->factorization()->updateColumnTranspose(&workSpace,&arrayVector);
    
    // now look at dual solution
    for (iRow=0;iRow<numberRows;iRow++) {
      // slack
      double value = array[iRow];
      if (fabs(duals[iRow]-value)>1.0e-3)
      printf("bad row %d old dual %g new %g\n",iRow,duals[iRow],value);
      //duals[iRow]=value;
    }
  }
#endif
#undef CLP_DEBUG
  double bestDj = tolerance;
  int bestSequence=-1;

  const double * cost = model_->costRegion(1);

  model_->clpMatrix()->setOriginalWanted(numberWanted);
  model_->clpMatrix()->setCurrentWanted(numberWanted);
  int iPassR=0,iPassC=0;
  // Setup two passes
  // This biases towards picking row variables
  // This probably should be fixed
  int startR[4];
  const int * which = infeasible_->getIndices();
  int nSlacks = infeasible_->getNumElements();
  startR[1]=nSlacks;
  startR[2]=0;
  double randomR = model_->randomNumberGenerator()->randomDouble();
  double dstart = static_cast<double> (nSlacks) * randomR;
  startR[0]=static_cast<int> (dstart);
  startR[3]=startR[0];
  double startC[4];
  startC[1]=1.0;
  startC[2]=0;
  double randomC = model_->randomNumberGenerator()->randomDouble();
  startC[0]=randomC;
  startC[3]=randomC;
  reducedCost = model_->djRegion(1);
  int sequenceOut = model_->sequenceOut();
  double * duals2 = duals-numberColumns;
  int chunk = CoinMin(1024,(numberColumns+nSlacks)/32);
  if (model_->numberIterations()%1000==0&&model_->logLevel()>1) {
    printf("%d wanted, nSlacks %d, chunk %d\n",numberWanted,nSlacks,chunk);
    int i;
    for (i=0;i<4;i++)
      printf("start R %d C %g ",startR[i],startC[i]);
    printf("\n");
  }
  chunk = CoinMax(chunk,256);
  bool finishedR=false,finishedC=false;
  bool doingR = randomR>randomC;
  //doingR=false;
  int saveNumberWanted = numberWanted;
  while (!finishedR||!finishedC) {
    if (finishedR)
      doingR=false;
    if (doingR) {
      int saveSequence = bestSequence;
      int start = startR[iPassR];
      int end = CoinMin(startR[iPassR+1],start+chunk/2);
      int jSequence;
      for (jSequence=start;jSequence<end;jSequence++) {
      int iSequence=which[jSequence];
      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 = fabs(cost[iSequence]+duals2[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;
            } else {
              // just to make sure we don't exit before got something
              numberWanted++;
            }
            }
          }
          break;
        case ClpSimplex::atUpperBound:
          value = cost[iSequence]+duals2[iSequence];
          if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
            // check flagged variable and correct dj
            if (!model_->flagged(iSequence)) {
              bestDj=value;
              bestSequence = iSequence;
            } else {
              // just to make sure we don't exit before got something
              numberWanted++;
            }
            }
          }
          break;
        case ClpSimplex::atLowerBound:
          value = -(cost[iSequence]+duals2[iSequence]);
          if (value>tolerance) {
            numberWanted--;
            if (value>bestDj) {
            // check flagged variable and correct dj
            if (!model_->flagged(iSequence)) {
              bestDj=value;
              bestSequence = iSequence;
            } else {
              // just to make sure we don't exit before got something
              numberWanted++;
            }
            }
          }
          break;
        }
      }
      if (!numberWanted)
        break;
      }
      numberLook -= (end-start);
      if (numberLook<0&&(10*(saveNumberWanted-numberWanted)>saveNumberWanted))
      numberWanted=0; // give up
      if (saveSequence!=bestSequence) {
      // dj
      reducedCost[bestSequence] = cost[bestSequence] +duals[bestSequence-numberColumns];
      bestDj=fabs(reducedCost[bestSequence]);
      model_->clpMatrix()->setSavedBestSequence(bestSequence);
      model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
      }
      model_->clpMatrix()->setCurrentWanted(numberWanted);
      if (!numberWanted) 
      break;
      doingR=false;
      // update start
      startR[iPassR]=jSequence;
      if (jSequence>=startR[iPassR+1]) {
      if (iPassR)
        finishedR=true;
      else
        iPassR=2;
      }
    }
    if (finishedC)
      doingR=true;
    if (!doingR) {
      int saveSequence = bestSequence;
      // Columns
      double start = startC[iPassC];
      // If we put this idea back then each function needs to update endFraction **
#if 0
      double dchunk = (static_cast<double> chunk)/(static_cast<double> numberColumns);
      double end = CoinMin(startC[iPassC+1],start+dchunk);;
#else
      double end=startC[iPassC+1]; // force end
#endif
      model_->clpMatrix()->partialPricing(model_,start,end,bestSequence,numberWanted);
      numberWanted=model_->clpMatrix()->currentWanted();
      numberLook -= static_cast<int> ((end-start)*numberColumns);
      if (numberLook<0&&(10*(saveNumberWanted-numberWanted)>saveNumberWanted))
      numberWanted=0; // give up
      if (saveSequence!=bestSequence) {
      // dj
      bestDj=fabs(model_->clpMatrix()->reducedCost(model_,bestSequence));
      }
      if (!numberWanted)
      break;
      doingR=true;
      // update start
      startC[iPassC]=end;
      if (end>=startC[iPassC+1]-1.0e-8) {
      if (iPassC)
        finishedC=true;
      else
        iPassC=2;
      }
    }
  }
  updates->setNumElements(0);

  // Restore tolerance
  model_->setCurrentDualTolerance(saveTolerance);
  // Now create variable if column generation
  model_->clpMatrix()->createVariable(model_,bestSequence);
#ifndef NDEBUG
  if (bestSequence>=0) {
    if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
      assert(model_->reducedCost(bestSequence)<0.0);
    if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
      assert(model_->reducedCost(bestSequence)>0.0);
  }
#endif
  return bestSequence;
}

Generated by  Doxygen 1.6.0   Back to index