Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/AbcPrimalColumnSteepest.cpp

    r2078 r2385  
    2424// Default Constructor
    2525//-------------------------------------------------------------------
    26 AbcPrimalColumnSteepest::AbcPrimalColumnSteepest (int mode)
    27      : AbcPrimalColumnPivot(),
    28        devex_(0.0),
    29        weights_(NULL),
    30        infeasible_(NULL),
    31        alternateWeights_(NULL),
    32        savedWeights_(NULL),
    33        reference_(NULL),
    34        state_(-1),
    35        mode_(mode),
    36        persistence_(normal),
    37        numberSwitched_(0),
    38        pivotSequence_(-1),
    39        savedPivotSequence_(-1),
    40        savedSequenceOut_(-1),
    41       sizeFactorization_(0)
     26AbcPrimalColumnSteepest::AbcPrimalColumnSteepest(int mode)
     27  : AbcPrimalColumnPivot()
     28  , devex_(0.0)
     29  , weights_(NULL)
     30  , infeasible_(NULL)
     31  , alternateWeights_(NULL)
     32  , savedWeights_(NULL)
     33  , reference_(NULL)
     34  , state_(-1)
     35  , mode_(mode)
     36  , persistence_(normal)
     37  , numberSwitched_(0)
     38  , pivotSequence_(-1)
     39  , savedPivotSequence_(-1)
     40  , savedSequenceOut_(-1)
     41  , sizeFactorization_(0)
    4242{
    43      type_ = 2 + 64 * mode;
     43  type_ = 2 + 64 * mode;
    4444}
    4545//-------------------------------------------------------------------
    4646// Copy constructor
    4747//-------------------------------------------------------------------
    48 AbcPrimalColumnSteepest::AbcPrimalColumnSteepest (const AbcPrimalColumnSteepest & rhs)
    49      : AbcPrimalColumnPivot(rhs)
     48AbcPrimalColumnSteepest::AbcPrimalColumnSteepest(const AbcPrimalColumnSteepest &rhs)
     49  : AbcPrimalColumnPivot(rhs)
    5050{
    51      state_ = rhs.state_;
    52      mode_ = rhs.mode_;
    53      persistence_ = rhs.persistence_;
    54      numberSwitched_ = rhs.numberSwitched_;
    55      model_ = rhs.model_;
    56      pivotSequence_ = rhs.pivotSequence_;
    57      savedPivotSequence_ = rhs.savedPivotSequence_;
    58      savedSequenceOut_ = rhs.savedSequenceOut_;
    59      sizeFactorization_ = rhs.sizeFactorization_;
    60      devex_ = rhs.devex_;
    61      if ((model_ && model_->whatsChanged() & 1) != 0) {
    62           if (rhs.infeasible_) {
    63                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    64           } else {
    65                infeasible_ = NULL;
    66           }
    67           reference_ = NULL;
    68           if (rhs.weights_) {
    69                assert(model_);
    70                int number = model_->numberRows() + model_->numberColumns();
    71                assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
    72                weights_ = new double[number];
    73                CoinMemcpyN(rhs.weights_, number, weights_);
    74                savedWeights_ = new double[number];
    75                CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
    76                if (mode_ != 1) {
    77                     reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
    78                }
    79           } else {
    80                weights_ = NULL;
    81                savedWeights_ = NULL;
    82           }
    83           if (rhs.alternateWeights_) {
     51  state_ = rhs.state_;
     52  mode_ = rhs.mode_;
     53  persistence_ = rhs.persistence_;
     54  numberSwitched_ = rhs.numberSwitched_;
     55  model_ = rhs.model_;
     56  pivotSequence_ = rhs.pivotSequence_;
     57  savedPivotSequence_ = rhs.savedPivotSequence_;
     58  savedSequenceOut_ = rhs.savedSequenceOut_;
     59  sizeFactorization_ = rhs.sizeFactorization_;
     60  devex_ = rhs.devex_;
     61  if ((model_ && model_->whatsChanged() & 1) != 0) {
     62    if (rhs.infeasible_) {
     63      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     64    } else {
     65      infeasible_ = NULL;
     66    }
     67    reference_ = NULL;
     68    if (rhs.weights_) {
     69      assert(model_);
     70      int number = model_->numberRows() + model_->numberColumns();
     71      assert(number == rhs.model_->numberRows() + rhs.model_->numberColumns());
     72      weights_ = new double[number];
     73      CoinMemcpyN(rhs.weights_, number, weights_);
     74      savedWeights_ = new double[number];
     75      CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
     76      if (mode_ != 1) {
     77        reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
     78      }
     79    } else {
     80      weights_ = NULL;
     81      savedWeights_ = NULL;
     82    }
     83    if (rhs.alternateWeights_) {
    8484#ifndef ABC_USE_COIN_FACTORIZATIONz
    85                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    86                // zero
    87                CoinZeroN(alternateWeights_->getIndices(),
    88                         alternateWeights_->capacity());
     85      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     86      // zero
     87      CoinZeroN(alternateWeights_->getIndices(),
     88        alternateWeights_->capacity());
    8989#else
    9090#endif
    91           } else {
    92                alternateWeights_ = NULL;
    93           }
    94      } else {
    95           infeasible_ = NULL;
    96           reference_ = NULL;
    97           weights_ = NULL;
    98           savedWeights_ = NULL;
    99           alternateWeights_ = NULL;
    100      }
     91    } else {
     92      alternateWeights_ = NULL;
     93    }
     94  } else {
     95    infeasible_ = NULL;
     96    reference_ = NULL;
     97    weights_ = NULL;
     98    savedWeights_ = NULL;
     99    alternateWeights_ = NULL;
     100  }
    101101}
    102102
     
    104104// Destructor
    105105//-------------------------------------------------------------------
    106 AbcPrimalColumnSteepest::~AbcPrimalColumnSteepest ()
     106AbcPrimalColumnSteepest::~AbcPrimalColumnSteepest()
    107107{
    108      delete [] weights_;
    109      delete infeasible_;
    110      delete alternateWeights_;
    111      delete [] savedWeights_;
    112      delete [] reference_;
     108  delete[] weights_;
     109  delete infeasible_;
     110  delete alternateWeights_;
     111  delete[] savedWeights_;
     112  delete[] reference_;
    113113}
    114114
     
    117117//-------------------------------------------------------------------
    118118AbcPrimalColumnSteepest &
    119 AbcPrimalColumnSteepest::operator=(const AbcPrimalColumnSteepest& rhs)
     119AbcPrimalColumnSteepest::operator=(const AbcPrimalColumnSteepest &rhs)
    120120{
    121      if (this != &rhs) {
    122           AbcPrimalColumnPivot::operator=(rhs);
    123           state_ = rhs.state_;
    124           mode_ = rhs.mode_;
    125           persistence_ = rhs.persistence_;
    126           numberSwitched_ = rhs.numberSwitched_;
    127           model_ = rhs.model_;
    128           pivotSequence_ = rhs.pivotSequence_;
    129           savedPivotSequence_ = rhs.savedPivotSequence_;
    130           savedSequenceOut_ = rhs.savedSequenceOut_;
    131           sizeFactorization_ = rhs.sizeFactorization_;
    132           devex_ = rhs.devex_;
    133           delete [] weights_;
    134           delete [] reference_;
    135           reference_ = NULL;
    136           delete infeasible_;
    137           delete alternateWeights_;
    138           delete [] savedWeights_;
    139           savedWeights_ = NULL;
    140           if (rhs.infeasible_ != NULL) {
    141                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    142           } else {
    143                infeasible_ = NULL;
    144           }
    145           if (rhs.weights_ != NULL) {
    146                assert(model_);
    147                int number = model_->numberRows() + model_->numberColumns();
    148                assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
    149                weights_ = new double[number];
    150                CoinMemcpyN(rhs.weights_, number, weights_);
    151                savedWeights_ = new double[number];
    152                CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
    153                if (mode_ != 1) {
    154                     reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
    155                }
    156           } else {
    157                weights_ = NULL;
    158           }
    159           if (rhs.alternateWeights_ != NULL) {
     121  if (this != &rhs) {
     122    AbcPrimalColumnPivot::operator=(rhs);
     123    state_ = rhs.state_;
     124    mode_ = rhs.mode_;
     125    persistence_ = rhs.persistence_;
     126    numberSwitched_ = rhs.numberSwitched_;
     127    model_ = rhs.model_;
     128    pivotSequence_ = rhs.pivotSequence_;
     129    savedPivotSequence_ = rhs.savedPivotSequence_;
     130    savedSequenceOut_ = rhs.savedSequenceOut_;
     131    sizeFactorization_ = rhs.sizeFactorization_;
     132    devex_ = rhs.devex_;
     133    delete[] weights_;
     134    delete[] reference_;
     135    reference_ = NULL;
     136    delete infeasible_;
     137    delete alternateWeights_;
     138    delete[] savedWeights_;
     139    savedWeights_ = NULL;
     140    if (rhs.infeasible_ != NULL) {
     141      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     142    } else {
     143      infeasible_ = NULL;
     144    }
     145    if (rhs.weights_ != NULL) {
     146      assert(model_);
     147      int number = model_->numberRows() + model_->numberColumns();
     148      assert(number == rhs.model_->numberRows() + rhs.model_->numberColumns());
     149      weights_ = new double[number];
     150      CoinMemcpyN(rhs.weights_, number, weights_);
     151      savedWeights_ = new double[number];
     152      CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
     153      if (mode_ != 1) {
     154        reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
     155      }
     156    } else {
     157      weights_ = NULL;
     158    }
     159    if (rhs.alternateWeights_ != NULL) {
    160160#ifndef ABC_USE_COIN_FACTORIZATIONz
    161                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    162                // zero
    163                CoinZeroN(alternateWeights_->getIndices(),
    164                         alternateWeights_->capacity());
     161      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     162      // zero
     163      CoinZeroN(alternateWeights_->getIndices(),
     164        alternateWeights_->capacity());
    165165#else
    166166#endif
    167           } else {
    168                alternateWeights_ = NULL;
    169           }
    170      }
    171      return *this;
     167    } else {
     168      alternateWeights_ = NULL;
     169    }
     170  }
     171  return *this;
    172172}
    173173// These have to match ClpPackedMatrix version
     
    178178        that is just +-weight where a feasibility changed.  It also has
    179179        reduced cost from last iteration in pivot row*/
    180 int
    181 AbcPrimalColumnSteepest::pivotColumn(CoinPartitionedVector * updates,
    182                                      CoinPartitionedVector * spareRow2,
    183                                      CoinPartitionedVector * spareColumn1)
     180int AbcPrimalColumnSteepest::pivotColumn(CoinPartitionedVector *updates,
     181  CoinPartitionedVector *spareRow2,
     182  CoinPartitionedVector *spareColumn1)
    184183{
    185      assert(model_);
    186      int number = 0;
    187      int * index;
    188      double tolerance = model_->currentDualTolerance();
    189      // we can't really trust infeasibilities if there is dual error
    190      // this coding has to mimic coding in checkDualSolution
    191      double error = CoinMin(1.0e-2, model_->largestDualError());
    192      // allow tolerance at least slightly bigger than standard
    193      tolerance = tolerance + error;
    194      int pivotRow = model_->pivotRow();
    195      int anyUpdates;
    196      double * infeas = infeasible_->denseVector();
    197 
    198      // Local copy of mode so can decide what to do
    199      int switchType;
    200      if (mode_ == 4)
    201           switchType = 5 - numberSwitched_;
    202      else if (mode_ >= 10)
    203           switchType = 3;
    204      else
    205           switchType = mode_;
    206      /* switchType -
     184  assert(model_);
     185  int number = 0;
     186  int *index;
     187  double tolerance = model_->currentDualTolerance();
     188  // we can't really trust infeasibilities if there is dual error
     189  // this coding has to mimic coding in checkDualSolution
     190  double error = CoinMin(1.0e-2, model_->largestDualError());
     191  // allow tolerance at least slightly bigger than standard
     192  tolerance = tolerance + error;
     193  int pivotRow = model_->pivotRow();
     194  int anyUpdates;
     195  double *infeas = infeasible_->denseVector();
     196
     197  // Local copy of mode so can decide what to do
     198  int switchType;
     199  if (mode_ == 4)
     200    switchType = 5 - numberSwitched_;
     201  else if (mode_ >= 10)
     202    switchType = 3;
     203  else
     204    switchType = mode_;
     205  /* switchType -
    207206        0 - all exact devex
    208207        1 - all steepest
     
    213212        10 - can go to mini-sprint
    214213     */
    215      if (updates->getNumElements() > 1) {
    216           // would have to have two goes for devex, three for steepest
    217           anyUpdates = 2;
    218      } else if (updates->getNumElements()) {
    219           if (updates->getIndices()[0] == pivotRow
    220               && fabs(updates->denseVector()[pivotRow]) > 1.0e-6) {
    221             //&& fabs(updates->denseVector()[pivotRow]) < 1.0e6) {
    222                // reasonable size
    223                anyUpdates = 1;
    224                //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
    225                //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
    226           } else {
    227                // too small
    228                anyUpdates = 2;
     214  if (updates->getNumElements() > 1) {
     215    // would have to have two goes for devex, three for steepest
     216    anyUpdates = 2;
     217  } else if (updates->getNumElements()) {
     218    if (updates->getIndices()[0] == pivotRow
     219      && fabs(updates->denseVector()[pivotRow]) > 1.0e-6) {
     220      //&& fabs(updates->denseVector()[pivotRow]) < 1.0e6) {
     221      // reasonable size
     222      anyUpdates = 1;
     223      //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
     224      //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
     225    } else {
     226      // too small
     227      anyUpdates = 2;
     228    }
     229  } else if (pivotSequence_ >= 0) {
     230    // just after re-factorization
     231    anyUpdates = -1;
     232  } else {
     233    // sub flip - nothing to do
     234    anyUpdates = 0;
     235  }
     236  int sequenceOut = model_->sequenceOut();
     237  if (switchType == 5) {
     238    pivotSequence_ = -1;
     239    pivotRow = -1;
     240    // See if to switch
     241    int numberRows = model_->numberRows();
     242    int numberWanted = 10;
     243    int numberColumns = model_->numberColumns();
     244    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     245    // Number of dual infeasibilities at last invert
     246    int numberDual = model_->numberDualInfeasibilities();
     247    int numberLook = CoinMin(numberDual, numberColumns / 10);
     248    if (ratio < 1.0) {
     249      numberWanted = 100;
     250      numberLook /= 20;
     251      numberWanted = CoinMax(numberWanted, numberLook);
     252    } else if (ratio < 3.0) {
     253      numberWanted = 500;
     254      numberLook /= 15;
     255      numberWanted = CoinMax(numberWanted, numberLook);
     256    } else if (ratio < 4.0 || mode_ == 5) {
     257      numberWanted = 1000;
     258      numberLook /= 10;
     259      numberWanted = CoinMax(numberWanted, numberLook);
     260    } else if (mode_ != 5) {
     261      switchType = 4;
     262      // initialize
     263      numberSwitched_++;
     264      // Make sure will re-do
     265      delete[] weights_;
     266      weights_ = NULL;
     267      //work space
     268      int whichArray[2];
     269      for (int i = 0; i < 2; i++)
     270        whichArray[i] = model_->getAvailableArray();
     271      CoinIndexedVector *array1 = model_->usefulArray(whichArray[0]);
     272      CoinIndexedVector *array2 = model_->usefulArray(whichArray[1]);
     273      model_->computeDuals(NULL, array1, array2);
     274      for (int i = 0; i < 2; i++)
     275        model_->setAvailableArray(whichArray[i]);
     276      saveWeights(model_, 4);
     277      anyUpdates = 0;
     278      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
     279    }
     280    if (switchType == 5) {
     281      numberLook *= 5; // needs tuning for gub
     282      if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
     283        COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
     284          sizeFactorization_, ratio, numberWanted, numberLook));
     285      }
     286      // Update duals and row djs
     287      // Do partial pricing
     288      return partialPricing(updates, numberWanted, numberLook);
     289    }
     290  }
     291  if (switchType == 5) {
     292    if (anyUpdates > 0) {
     293      justDjs(updates, spareColumn1);
     294    }
     295  } else if (anyUpdates == 1) {
     296    if (switchType < 4) {
     297      // exact etc when can use dj
     298      doSteepestWork(updates, spareRow2, spareColumn1, 2);
     299    } else {
     300      // devex etc when can use dj
     301      djsAndDevex(updates, spareRow2,
     302        spareColumn1);
     303    }
     304  } else if (anyUpdates == -1) {
     305    if (switchType < 4) {
     306      // exact etc when djs okay
     307      //if (model_->maximumIterations()!=1000)
     308      doSteepestWork(updates, spareRow2, spareColumn1, 1);
     309    } else {
     310      // devex etc when djs okay
     311      justDevex(updates,
     312        spareColumn1);
     313    }
     314  } else if (anyUpdates == 2) {
     315    if (switchType < 4) {
     316      // exact etc when have to use pivot
     317      doSteepestWork(updates, spareRow2, spareColumn1, 3);
     318    } else {
     319      // devex etc when have to use pivot
     320      djsAndDevex2(updates,
     321        spareColumn1);
     322    }
     323  }
     324#ifdef CLP_DEBUG
     325  alternateWeights_->checkClear();
     326#endif
     327  // make sure outgoing from last iteration okay
     328  if (sequenceOut >= 0) {
     329    AbcSimplex::Status status = model_->getInternalStatus(sequenceOut);
     330    double value = model_->reducedCost(sequenceOut);
     331
     332    switch (status) {
     333
     334    case AbcSimplex::basic:
     335    case AbcSimplex::isFixed:
     336      break;
     337    case AbcSimplex::isFree:
     338    case AbcSimplex::superBasic:
     339      if (fabs(value) > FREE_ACCEPT * tolerance) {
     340        // we are going to bias towards free (but only if reasonable)
     341        value *= FREE_BIAS;
     342        // store square in list
     343        if (infeas[sequenceOut])
     344          infeas[sequenceOut] = value * value; // already there
     345        else
     346          infeasible_->quickAdd(sequenceOut, value * value);
     347      } else {
     348        infeasible_->zero(sequenceOut);
     349      }
     350      break;
     351    case AbcSimplex::atUpperBound:
     352      if (value > tolerance) {
     353        // store square in list
     354        if (infeas[sequenceOut])
     355          infeas[sequenceOut] = value * value; // already there
     356        else
     357          infeasible_->quickAdd(sequenceOut, value * value);
     358      } else {
     359        infeasible_->zero(sequenceOut);
     360      }
     361      break;
     362    case AbcSimplex::atLowerBound:
     363      if (value < -tolerance) {
     364        // store square in list
     365        if (infeas[sequenceOut])
     366          infeas[sequenceOut] = value * value; // already there
     367        else
     368          infeasible_->quickAdd(sequenceOut, value * value);
     369      } else {
     370        infeasible_->zero(sequenceOut);
     371      }
     372    }
     373  }
     374  // update of duals finished - now do pricing
     375  // See what sort of pricing
     376  int numberWanted = 10;
     377  number = infeasible_->getNumElements();
     378  int numberColumns = model_->numberColumns();
     379  if (switchType == 5) {
     380    pivotSequence_ = -1;
     381    pivotRow = -1;
     382    // See if to switch
     383    int numberRows = model_->numberRows();
     384    // ratio is done on number of columns here
     385    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
     386    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     387    if (ratio < 1.0) {
     388      numberWanted = CoinMax(100, number / 200);
     389    } else if (ratio < 2.0 - 1.0) {
     390      numberWanted = CoinMax(500, number / 40);
     391    } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
     392      numberWanted = CoinMax(2000, number / 10);
     393      numberWanted = CoinMax(numberWanted, numberColumns / 30);
     394    } else if (mode_ != 5) {
     395      switchType = 4;
     396      // initialize
     397      numberSwitched_++;
     398      // Make sure will re-do
     399      delete[] weights_;
     400      weights_ = NULL;
     401      saveWeights(model_, 4);
     402      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
     403    }
     404    //if (model_->numberIterations()%1000==0)
     405    //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
     406  }
     407  int numberRows = model_->numberRows();
     408  // ratio is done on number of rows here
     409  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     410  if (switchType == 4) {
     411    // Still in devex mode
     412    // Go to steepest if lot of iterations?
     413    if (ratio < 5.0) {
     414      numberWanted = CoinMax(2000, number / 10);
     415      numberWanted = CoinMax(numberWanted, numberColumns / 20);
     416    } else if (ratio < 7.0) {
     417      numberWanted = CoinMax(2000, number / 5);
     418      numberWanted = CoinMax(numberWanted, numberColumns / 10);
     419    } else {
     420      // we can zero out
     421      updates->clear();
     422      spareColumn1->clear();
     423      switchType = 3;
     424      // initialize
     425      pivotSequence_ = -1;
     426      pivotRow = -1;
     427      numberSwitched_++;
     428      // Make sure will re-do
     429      delete[] weights_;
     430      weights_ = NULL;
     431      saveWeights(model_, 4);
     432      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
     433      updates->clear();
     434    }
     435    if (model_->numberIterations() % 1000 == 0)
     436      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
     437  }
     438  if (switchType < 4) {
     439    if (switchType < 2) {
     440      numberWanted = number + 1;
     441    } else if (switchType == 2) {
     442      numberWanted = CoinMax(2000, number / 8);
     443    } else {
     444      if (ratio < 1.0) {
     445        numberWanted = CoinMax(2000, number / 20);
     446      } else if (ratio < 5.0) {
     447        numberWanted = CoinMax(2000, number / 10);
     448        numberWanted = CoinMax(numberWanted, numberColumns / 40);
     449      } else if (ratio < 10.0) {
     450        numberWanted = CoinMax(2000, number / 8);
     451        numberWanted = CoinMax(numberWanted, numberColumns / 20);
     452      } else {
     453        ratio = number * (ratio / 80.0);
     454        if (ratio > number) {
     455          numberWanted = number + 1;
     456        } else {
     457          numberWanted = CoinMax(2000, static_cast< int >(ratio));
     458          numberWanted = CoinMax(numberWanted, numberColumns / 10);
     459        }
     460      }
     461    }
     462    //if (model_->numberIterations()%1000==0)
     463    //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
     464    //switchType);
     465  }
     466
     467  double bestDj = 1.0e-30;
     468  int bestSequence = -1;
     469
     470  int i, iSequence;
     471  index = infeasible_->getIndices();
     472  number = infeasible_->getNumElements();
     473  if (model_->numberIterations() < model_->lastBadIteration() + 200 && model_->factorization()->pivots() > 10) {
     474    // we can't really trust infeasibilities if there is dual error
     475    double checkTolerance = 1.0e-8;
     476    if (model_->largestDualError() > checkTolerance)
     477      tolerance *= model_->largestDualError() / checkTolerance;
     478    // But cap
     479    tolerance = CoinMin(1000.0, tolerance);
     480  }
     481#ifdef CLP_DEBUG
     482  if (model_->numberDualInfeasibilities() == 1)
     483    printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
     484      model_->largestDualError(), model_, model_->messageHandler(),
     485      number);
     486#endif
     487  // stop last one coming immediately
     488  double saveOutInfeasibility = 0.0;
     489  if (sequenceOut >= 0) {
     490    saveOutInfeasibility = infeas[sequenceOut];
     491    infeas[sequenceOut] = 0.0;
     492  }
     493  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
     494    tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
     495  tolerance *= tolerance; // as we are using squares
     496
     497  int iPass;
     498  // Setup two passes
     499  int start[4];
     500  start[1] = number;
     501  start[2] = 0;
     502  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
     503  start[0] = static_cast< int >(dstart);
     504  start[3] = start[0];
     505  //double largestWeight=0.0;
     506  //double smallestWeight=1.0e100;
     507  for (iPass = 0; iPass < 2; iPass++) {
     508    int end = start[2 * iPass + 1];
     509    if (switchType < 5) {
     510      for (i = start[2 * iPass]; i < end; i++) {
     511        iSequence = index[i];
     512        double value = infeas[iSequence];
     513        double weight = weights_[iSequence];
     514        if (value > tolerance) {
     515          //weight=1.0;
     516          if (value > bestDj * weight) {
     517            // check flagged variable and correct dj
     518            if (!model_->flagged(iSequence)) {
     519              bestDj = value / weight;
     520              bestSequence = iSequence;
     521            } else {
     522              // just to make sure we don't exit before got something
     523              numberWanted++;
     524            }
    229525          }
    230      } else if (pivotSequence_ >= 0) {
    231           // just after re-factorization
    232           anyUpdates = -1;
    233      } else {
    234           // sub flip - nothing to do
    235           anyUpdates = 0;
    236      }
    237      int sequenceOut = model_->sequenceOut();
    238      if (switchType == 5) {
    239                pivotSequence_ = -1;
    240                pivotRow = -1;
    241                // See if to switch
    242                int numberRows = model_->numberRows();
    243                int numberWanted = 10;
    244                int numberColumns = model_->numberColumns();
    245                double ratio = static_cast<double> (sizeFactorization_) /
    246                               static_cast<double> (numberRows);
    247                // Number of dual infeasibilities at last invert
    248                int numberDual = model_->numberDualInfeasibilities();
    249                int numberLook = CoinMin(numberDual, numberColumns / 10);
    250                if (ratio < 1.0) {
    251                     numberWanted = 100;
    252                     numberLook /= 20;
    253                     numberWanted = CoinMax(numberWanted, numberLook);
    254                } else if (ratio < 3.0) {
    255                     numberWanted = 500;
    256                     numberLook /= 15;
    257                     numberWanted = CoinMax(numberWanted, numberLook);
    258                } else if (ratio < 4.0 || mode_ == 5) {
    259                     numberWanted = 1000;
    260                     numberLook /= 10;
    261                     numberWanted = CoinMax(numberWanted, numberLook);
    262                } else if (mode_ != 5) {
    263                     switchType = 4;
    264                     // initialize
    265                     numberSwitched_++;
    266                     // Make sure will re-do
    267                     delete [] weights_;
    268                     weights_ = NULL;
    269                     //work space
    270                     int whichArray[2];
    271                     for (int i=0;i<2;i++)
    272                       whichArray[i]=model_->getAvailableArray();
    273                     CoinIndexedVector * array1 = model_->usefulArray(whichArray[0]);
    274                     CoinIndexedVector * array2 = model_->usefulArray(whichArray[1]);
    275                     model_->computeDuals(NULL,array1,array2);
    276                     for (int i=0;i<2;i++)
    277                       model_->setAvailableArray(whichArray[i]);
    278                     saveWeights(model_, 4);
    279                     anyUpdates = 0;
    280                     COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
    281                }
    282                if (switchType == 5) {
    283                     numberLook *= 5; // needs tuning for gub
    284                     if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
    285                          COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
    286                                                   sizeFactorization_, ratio, numberWanted, numberLook));
    287                     }
    288                     // Update duals and row djs
    289                     // Do partial pricing
    290                     return partialPricing(updates, numberWanted, numberLook);
    291                }
    292      }
    293      if (switchType == 5) {
    294           if (anyUpdates > 0) {
    295                justDjs(updates,spareColumn1);
     526          numberWanted--;
     527        }
     528        if (!numberWanted)
     529          break;
     530      }
     531    } else {
     532      // Dantzig
     533      for (i = start[2 * iPass]; i < end; i++) {
     534        iSequence = index[i];
     535        double value = infeas[iSequence];
     536        if (value > tolerance) {
     537          if (value > bestDj) {
     538            // check flagged variable and correct dj
     539            if (!model_->flagged(iSequence)) {
     540              bestDj = value;
     541              bestSequence = iSequence;
     542            } else {
     543              // just to make sure we don't exit before got something
     544              numberWanted++;
     545            }
    296546          }
    297      } else if (anyUpdates == 1) {
    298           if (switchType < 4) {
    299                // exact etc when can use dj
    300             doSteepestWork(updates,spareRow2,spareColumn1,2);
    301           } else {
    302                // devex etc when can use dj
    303                djsAndDevex(updates, spareRow2,
    304                            spareColumn1);
    305           }
    306      } else if (anyUpdates == -1) {
    307           if (switchType < 4) {
    308                // exact etc when djs okay
    309             //if (model_->maximumIterations()!=1000)
    310             doSteepestWork(updates,spareRow2,spareColumn1,1);
    311           } else {
    312                // devex etc when djs okay
    313                justDevex(updates,
    314                          spareColumn1);
    315           }
    316      } else if (anyUpdates == 2) {
    317           if (switchType < 4) {
    318                // exact etc when have to use pivot
    319             doSteepestWork(updates,spareRow2,spareColumn1,3);
    320           } else {
    321                // devex etc when have to use pivot
    322                djsAndDevex2(updates,
    323                             spareColumn1);
    324           }
    325      }
    326 #ifdef CLP_DEBUG
    327      alternateWeights_->checkClear();
    328 #endif
    329      // make sure outgoing from last iteration okay
    330      if (sequenceOut >= 0) {
    331           AbcSimplex::Status status = model_->getInternalStatus(sequenceOut);
    332           double value = model_->reducedCost(sequenceOut);
    333 
    334           switch(status) {
    335 
    336           case AbcSimplex::basic:
    337           case AbcSimplex::isFixed:
    338                break;
    339           case AbcSimplex::isFree:
    340           case AbcSimplex::superBasic:
    341                if (fabs(value) > FREE_ACCEPT * tolerance) {
    342                     // we are going to bias towards free (but only if reasonable)
    343                     value *= FREE_BIAS;
    344                     // store square in list
    345                     if (infeas[sequenceOut])
    346                          infeas[sequenceOut] = value * value; // already there
    347                     else
    348                          infeasible_->quickAdd(sequenceOut, value * value);
    349                } else {
    350                     infeasible_->zero(sequenceOut);
    351                }
    352                break;
    353           case AbcSimplex::atUpperBound:
    354                if (value > tolerance) {
    355                     // store square in list
    356                     if (infeas[sequenceOut])
    357                          infeas[sequenceOut] = value * value; // already there
    358                     else
    359                          infeasible_->quickAdd(sequenceOut, value * value);
    360                } else {
    361                     infeasible_->zero(sequenceOut);
    362                }
    363                break;
    364           case AbcSimplex::atLowerBound:
    365                if (value < -tolerance) {
    366                     // store square in list
    367                     if (infeas[sequenceOut])
    368                          infeas[sequenceOut] = value * value; // already there
    369                     else
    370                          infeasible_->quickAdd(sequenceOut, value * value);
    371                } else {
    372                     infeasible_->zero(sequenceOut);
    373                }
    374           }
    375      }
    376      // update of duals finished - now do pricing
    377      // See what sort of pricing
    378      int numberWanted = 10;
    379      number = infeasible_->getNumElements();
    380      int numberColumns = model_->numberColumns();
    381      if (switchType == 5) {
    382           pivotSequence_ = -1;
    383           pivotRow = -1;
    384           // See if to switch
    385           int numberRows = model_->numberRows();
    386           // ratio is done on number of columns here
    387           //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
    388           double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    389           if (ratio < 1.0) {
    390                numberWanted = CoinMax(100, number / 200);
    391           } else if (ratio < 2.0 - 1.0) {
    392                numberWanted = CoinMax(500, number / 40);
    393           } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
    394                numberWanted = CoinMax(2000, number / 10);
    395                numberWanted = CoinMax(numberWanted, numberColumns / 30);
    396           } else if (mode_ != 5) {
    397                switchType = 4;
    398                // initialize
    399                numberSwitched_++;
    400                // Make sure will re-do
    401                delete [] weights_;
    402                weights_ = NULL;
    403                saveWeights(model_, 4);
    404                COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
    405           }
    406           //if (model_->numberIterations()%1000==0)
    407           //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
    408      }
    409      int numberRows = model_->numberRows();
    410      // ratio is done on number of rows here
    411      double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    412      if(switchType == 4) {
    413           // Still in devex mode
    414           // Go to steepest if lot of iterations?
    415           if (ratio < 5.0) {
    416                numberWanted = CoinMax(2000, number / 10);
    417                numberWanted = CoinMax(numberWanted, numberColumns / 20);
    418           } else if (ratio < 7.0) {
    419                numberWanted = CoinMax(2000, number / 5);
    420                numberWanted = CoinMax(numberWanted, numberColumns / 10);
    421           } else {
    422                // we can zero out
    423                updates->clear();
    424                spareColumn1->clear();
    425                switchType = 3;
    426                // initialize
    427                pivotSequence_ = -1;
    428                pivotRow = -1;
    429                numberSwitched_++;
    430                // Make sure will re-do
    431                delete [] weights_;
    432                weights_ = NULL;
    433                saveWeights(model_, 4);
    434                COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
    435                updates->clear();
    436           }
    437           if (model_->numberIterations() % 1000 == 0)
    438             COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
    439      }
    440      if (switchType < 4) {
    441           if (switchType < 2 ) {
    442                numberWanted = number + 1;
    443           } else if (switchType == 2) {
    444                numberWanted = CoinMax(2000, number / 8);
    445           } else {
    446                if (ratio < 1.0) {
    447                     numberWanted = CoinMax(2000, number / 20);
    448                } else if (ratio < 5.0) {
    449                     numberWanted = CoinMax(2000, number / 10);
    450                     numberWanted = CoinMax(numberWanted, numberColumns / 40);
    451                } else if (ratio < 10.0) {
    452                     numberWanted = CoinMax(2000, number / 8);
    453                     numberWanted = CoinMax(numberWanted, numberColumns / 20);
    454                } else {
    455                     ratio = number * (ratio / 80.0);
    456                     if (ratio > number) {
    457                          numberWanted = number + 1;
    458                     } else {
    459                          numberWanted = CoinMax(2000, static_cast<int> (ratio));
    460                          numberWanted = CoinMax(numberWanted, numberColumns / 10);
    461                     }
    462                }
    463           }
    464           //if (model_->numberIterations()%1000==0)
    465           //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
    466           //switchType);
    467      }
    468 
    469 
    470      double bestDj = 1.0e-30;
    471      int bestSequence = -1;
    472 
    473      int i, iSequence;
    474      index = infeasible_->getIndices();
    475      number = infeasible_->getNumElements();
    476      if(model_->numberIterations() < model_->lastBadIteration() + 200 &&
    477                model_->factorization()->pivots() > 10) {
    478           // we can't really trust infeasibilities if there is dual error
    479           double checkTolerance = 1.0e-8;
    480           if (model_->largestDualError() > checkTolerance)
    481                tolerance *= model_->largestDualError() / checkTolerance;
    482           // But cap
    483           tolerance = CoinMin(1000.0, tolerance);
    484      }
    485 #ifdef CLP_DEBUG
    486      if (model_->numberDualInfeasibilities() == 1)
    487           printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
    488                  model_->largestDualError(), model_, model_->messageHandler(),
    489                  number);
    490 #endif
    491      // stop last one coming immediately
    492      double saveOutInfeasibility = 0.0;
    493      if (sequenceOut >= 0) {
    494           saveOutInfeasibility = infeas[sequenceOut];
    495           infeas[sequenceOut] = 0.0;
    496      }
    497      if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    498           tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    499      tolerance *= tolerance; // as we are using squares
    500 
    501      int iPass;
    502      // Setup two passes
    503      int start[4];
    504      start[1] = number;
    505      start[2] = 0;
    506      double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
    507      start[0] = static_cast<int> (dstart);
    508      start[3] = start[0];
    509      //double largestWeight=0.0;
    510      //double smallestWeight=1.0e100;
    511      for (iPass = 0; iPass < 2; iPass++) {
    512           int end = start[2*iPass+1];
    513           if (switchType < 5) {
    514                for (i = start[2*iPass]; i < end; i++) {
    515                     iSequence = index[i];
    516                     double value = infeas[iSequence];
    517                     double weight = weights_[iSequence];
    518                     if (value > tolerance) {
    519                          //weight=1.0;
    520                          if (value > bestDj * weight) {
    521                               // check flagged variable and correct dj
    522                               if (!model_->flagged(iSequence)) {
    523                                    bestDj = value / weight;
    524                                    bestSequence = iSequence;
    525                               } else {
    526                                    // just to make sure we don't exit before got something
    527                                    numberWanted++;
    528                               }
    529                          }
    530                          numberWanted--;
    531                     }
    532                     if (!numberWanted)
    533                          break;
    534                }
    535           } else {
    536                // Dantzig
    537                for (i = start[2*iPass]; i < end; i++) {
    538                     iSequence = index[i];
    539                     double value = infeas[iSequence];
    540                     if (value > tolerance) {
    541                          if (value > bestDj) {
    542                               // check flagged variable and correct dj
    543                               if (!model_->flagged(iSequence)) {
    544                                    bestDj = value;
    545                                    bestSequence = iSequence;
    546                               } else {
    547                                    // just to make sure we don't exit before got something
    548                                    numberWanted++;
    549                               }
    550                          }
    551                          numberWanted--;
    552                     }
    553                     if (!numberWanted)
    554                          break;
    555                }
    556           }
    557           if (!numberWanted)
    558                break;
    559      }
    560      if (sequenceOut >= 0) {
    561           infeas[sequenceOut] = saveOutInfeasibility;
    562      }
    563      /*if (model_->numberIterations()%100==0)
     547          numberWanted--;
     548        }
     549        if (!numberWanted)
     550          break;
     551      }
     552    }
     553    if (!numberWanted)
     554      break;
     555  }
     556  if (sequenceOut >= 0) {
     557    infeas[sequenceOut] = saveOutInfeasibility;
     558  }
     559  /*if (model_->numberIterations()%100==0)
    564560       printf("%d best %g\n",bestSequence,bestDj);*/
    565561
    566562#ifndef NDEBUG
    567      if (bestSequence >= 0) {
    568           if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
    569                assert(model_->reducedCost(bestSequence) < 0.0);
    570           if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound) {
    571                assert(model_->reducedCost(bestSequence) > 0.0);
    572           }
    573      }
     563  if (bestSequence >= 0) {
     564    if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
     565      assert(model_->reducedCost(bestSequence) < 0.0);
     566    if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound) {
     567      assert(model_->reducedCost(bestSequence) > 0.0);
     568    }
     569  }
    574570#endif
    575571#if 1
    576      if (model_->logLevel()==127) {
    577        double * reducedCost=model_->djRegion();
    578        for (int i=0;i<numberRows;i++)
    579          printf("row %d weight %g infeas %g dj %g\n",i,weights_[i],infeas[i]>1.0e-13 ? infeas[i]: 0.0,reducedCost[i]);
    580        for (int i=numberRows;i<numberRows+numberColumns;i++)
    581          printf("column %d weight %g infeas %g dj %g\n",i-numberRows,weights_[i],infeas[i]>1.0e-13 ? infeas[i]: 0.0,reducedCost[i]);
    582      }
     572  if (model_->logLevel() == 127) {
     573    double *reducedCost = model_->djRegion();
     574    for (int i = 0; i < numberRows; i++)
     575      printf("row %d weight %g infeas %g dj %g\n", i, weights_[i], infeas[i] > 1.0e-13 ? infeas[i] : 0.0, reducedCost[i]);
     576    for (int i = numberRows; i < numberRows + numberColumns; i++)
     577      printf("column %d weight %g infeas %g dj %g\n", i - numberRows, weights_[i], infeas[i] > 1.0e-13 ? infeas[i] : 0.0, reducedCost[i]);
     578  }
    583579#endif
    584580#if SOME_DEBUG_1
    585      if (switchType<4) {
    586        // check for accuracy
    587        int iCheck = 892;
    588        //printf("weight for iCheck is %g\n",weights_[iCheck]);
    589        int numberRows = model_->numberRows();
    590        int numberColumns = model_->numberColumns();
    591        for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    592          if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
    593              model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
    594            checkAccuracy(iCheck, 1.0e-1, updates);
    595        }
    596      }
     581  if (switchType < 4) {
     582    // check for accuracy
     583    int iCheck = 892;
     584    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     585    int numberRows = model_->numberRows();
     586    int numberColumns = model_->numberColumns();
     587    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     588      if (model_->getInternalStatus(iCheck) != AbcSimplex::basic && model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
     589        checkAccuracy(iCheck, 1.0e-1, updates);
     590    }
     591  }
    597592#endif
    598593#if 0
     
    603598     }
    604599#endif
    605      return bestSequence;
     600  return bestSequence;
    606601}
    607602/* Does steepest work
     
    612607   3 - both using extra array
    613608*/
    614 int
    615 AbcPrimalColumnSteepest::doSteepestWork(CoinPartitionedVector * updates,
    616                                         CoinPartitionedVector * spareRow2,
    617                                         CoinPartitionedVector * spareColumn1,
    618                                         int type)
     609int AbcPrimalColumnSteepest::doSteepestWork(CoinPartitionedVector *updates,
     610  CoinPartitionedVector *spareRow2,
     611  CoinPartitionedVector *spareColumn1,
     612  int type)
    619613{
    620   int pivotRow = model_->pivotRow(); 
    621   if (type==1) {
    622     pivotRow=pivotSequence_;
    623     if (pivotRow<0)
    624       type=-1;
    625   } else if (pivotSequence_<0) {
    626     assert (model_->sequenceIn()==model_->sequenceOut());
    627     type=0;
     614  int pivotRow = model_->pivotRow();
     615  if (type == 1) {
     616    pivotRow = pivotSequence_;
     617    if (pivotRow < 0)
     618      type = -1;
     619  } else if (pivotSequence_ < 0) {
     620    assert(model_->sequenceIn() == model_->sequenceOut());
     621    type = 0;
    628622  }
    629623  // see if reference
    630   int sequenceIn = pivotRow>=0 ? model_->sequenceIn() : -1;
     624  int sequenceIn = pivotRow >= 0 ? model_->sequenceIn() : -1;
    631625  // save outgoing weight round update
    632626  double outgoingWeight = 0.0;
     
    636630  double referenceIn;
    637631  if (mode_ != 1) {
    638     if(sequenceIn>=0&&reference(sequenceIn))
     632    if (sequenceIn >= 0 && reference(sequenceIn))
    639633      referenceIn = 1.0;
    640634    else
     
    643637    referenceIn = -1.0;
    644638  }
    645   double * infeasibilities=infeasible_->denseVector();
    646   if (sequenceIn>=0)
    647     infeasibilities[sequenceIn]=0.0;
    648   double * array=updates->denseVector();
     639  double *infeasibilities = infeasible_->denseVector();
     640  if (sequenceIn >= 0)
     641    infeasibilities[sequenceIn] = 0.0;
     642  double *array = updates->denseVector();
    649643  double scaleFactor;
    650   double devex=devex_;
    651   if (type==0) {
     644  double devex = devex_;
     645  if (type == 0) {
    652646    // just djs - to keep clean swap
    653     scaleFactor=1.0;
    654     CoinPartitionedVector * temp = updates;
    655     updates=spareRow2;
    656     spareRow2=temp;
    657     assert (!alternateWeights_->getNumElements());
    658     devex=0.0;
    659   } else if (type==1) {
     647    scaleFactor = 1.0;
     648    CoinPartitionedVector *temp = updates;
     649    updates = spareRow2;
     650    spareRow2 = temp;
     651    assert(!alternateWeights_->getNumElements());
     652    devex = 0.0;
     653  } else if (type == 1) {
    660654    // just steepest
    661655    updates->clear();
    662     scaleFactor=COIN_DBL_MAX;
     656    scaleFactor = COIN_DBL_MAX;
    663657    // might as well set dj to 1
    664658    double dj = -1.0;
    665     assert (pivotRow>=0);
     659    assert(pivotRow >= 0);
    666660    spareRow2->createUnpacked(1, &pivotRow, &dj);
    667   } else if (type==2) {
     661  } else if (type == 2) {
    668662    // using scaleFactor - swap
    669     assert (pivotRow>=0);
    670     scaleFactor=-array[pivotRow];
    671     array[pivotRow]=-1.0;
    672     CoinPartitionedVector * temp = updates;
    673     updates=spareRow2;
    674     spareRow2=temp;
    675   } else if (type==3) {
     663    assert(pivotRow >= 0);
     664    scaleFactor = -array[pivotRow];
     665    array[pivotRow] = -1.0;
     666    CoinPartitionedVector *temp = updates;
     667    updates = spareRow2;
     668    spareRow2 = temp;
     669  } else if (type == 3) {
    676670    // need two arrays
    677     scaleFactor=0.0;
     671    scaleFactor = 0.0;
    678672    // might as well set dj to 1
    679673    double dj = -1.0;
    680     assert (pivotRow>=0);
     674    assert(pivotRow >= 0);
    681675    spareRow2->createUnpacked(1, &pivotRow, &dj);
    682676  }
    683   if (type>=0) {
     677  if (type >= 0) {
    684678    // parallelize this
    685     if (type==0) {
     679    if (type == 0) {
    686680      model_->factorization()->updateColumnTranspose(*spareRow2);
    687     } else if (type<3) {
    688       cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2,0);
    689       model_->factorization()->updateColumnTransposeCpu(*alternateWeights_,1);
     681    } else if (type < 3) {
     682      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2, 0);
     683      model_->factorization()->updateColumnTransposeCpu(*alternateWeights_, 1);
    690684      cilk_sync;
    691685    } else {
    692       cilk_spawn model_->factorization()->updateColumnTransposeCpu(*updates,0);
    693       cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2,1);
    694       model_->factorization()->updateColumnTransposeCpu(*alternateWeights_,2);
     686      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*updates, 0);
     687      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2, 1);
     688      model_->factorization()->updateColumnTransposeCpu(*alternateWeights_, 2);
    695689      cilk_sync;
    696690    }
    697691    model_->abcMatrix()->primalColumnDouble(*spareRow2,
    698                                             *updates,
    699                                             *alternateWeights_,
    700                                             *spareColumn1,
    701                                             *infeasible_,referenceIn,devex,
    702                                             reference_,weights_,scaleFactor);
    703   }
    704   pivotSequence_=-1;
     692      *updates,
     693      *alternateWeights_,
     694      *spareColumn1,
     695      *infeasible_, referenceIn, devex,
     696      reference_, weights_, scaleFactor);
     697  }
     698  pivotSequence_ = -1;
    705699  // later do pricing here
    706700  // later move pricing into abcmatrix
     
    714708}
    715709// Just update djs
    716 void
    717 AbcPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
    718                                  CoinIndexedVector * spareColumn1)
     710void AbcPrimalColumnSteepest::justDjs(CoinIndexedVector *updates,
     711  CoinIndexedVector *spareColumn1)
    719712{
    720      int iSection, j;
    721      int number = 0;
    722      double multiplier;
    723      int * index;
    724      double * updateBy;
    725      double * reducedCost;
    726      double tolerance = model_->currentDualTolerance();
    727      // we can't really trust infeasibilities if there is dual error
    728      // this coding has to mimic coding in checkDualSolution
    729      double error = CoinMin(1.0e-2, model_->largestDualError());
    730      // allow tolerance at least slightly bigger than standard
    731      tolerance = tolerance +  error;
    732      int pivotRow = model_->pivotRow();
    733      double * infeas = infeasible_->denseVector();
    734      //updates->scanAndPack();
    735      model_->factorization()->updateColumnTranspose(*updates);
    736 
    737      // put row of tableau in rowArray and columnArray
    738      model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
    739      // normal
    740      reducedCost = model_->djRegion();
    741      for (iSection = 0; iSection < 2; iSection++) {
    742 
    743           int addSequence;
    744 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    745           double slack_multiplier;
    746 #endif
    747 
    748           if (!iSection) {
    749                number = updates->getNumElements();
    750                index = updates->getIndices();
    751                updateBy = updates->denseVector();
    752                addSequence = 0;
    753                multiplier=-1;
    754 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    755                slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
    756 #endif
    757           } else {
    758                number = spareColumn1->getNumElements();
    759                index = spareColumn1->getIndices();
    760                updateBy = spareColumn1->denseVector();
    761                addSequence = model_->maximumAbcNumberRows();
    762                multiplier=1;
    763 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    764                slack_multiplier = 1.0;
    765 #endif
    766           }
    767 
    768           for (j = 0; j < number; j++) {
    769                int iSequence = index[j];
    770                double tableauValue=updateBy[iSequence];
    771                updateBy[iSequence] = 0.0;
    772                iSequence += addSequence;
    773                double value = reducedCost[iSequence];
    774                value -= multiplier*tableauValue;
    775                reducedCost[iSequence] = value;
    776                AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    777 
    778                switch(status) {
    779 
    780                case AbcSimplex::basic:
    781                     infeasible_->zero(iSequence);
    782                case AbcSimplex::isFixed:
    783                     break;
    784                case AbcSimplex::isFree:
    785                case AbcSimplex::superBasic:
    786                     if (fabs(value) > FREE_ACCEPT * tolerance) {
    787                          // we are going to bias towards free (but only if reasonable)
    788                          value *= FREE_BIAS;
    789                          // store square in list
    790                          if (infeas[iSequence])
    791                               infeas[iSequence] = value * value; // already there
    792                          else
    793                               infeasible_->quickAdd(iSequence , value * value);
    794                     } else {
    795                          infeasible_->zero(iSequence);
    796                     }
    797                     break;
    798                case AbcSimplex::atUpperBound:
    799                     if (value > tolerance) {
     713  int iSection, j;
     714  int number = 0;
     715  double multiplier;
     716  int *index;
     717  double *updateBy;
     718  double *reducedCost;
     719  double tolerance = model_->currentDualTolerance();
     720  // we can't really trust infeasibilities if there is dual error
     721  // this coding has to mimic coding in checkDualSolution
     722  double error = CoinMin(1.0e-2, model_->largestDualError());
     723  // allow tolerance at least slightly bigger than standard
     724  tolerance = tolerance + error;
     725  int pivotRow = model_->pivotRow();
     726  double *infeas = infeasible_->denseVector();
     727  //updates->scanAndPack();
     728  model_->factorization()->updateColumnTranspose(*updates);
     729
     730  // put row of tableau in rowArray and columnArray
     731  model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
     732  // normal
     733  reducedCost = model_->djRegion();
     734  for (iSection = 0; iSection < 2; iSection++) {
     735
     736    int addSequence;
    800737#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    801                         value *= value*slack_multiplier;
     738    double slack_multiplier;
     739#endif
     740
     741    if (!iSection) {
     742      number = updates->getNumElements();
     743      index = updates->getIndices();
     744      updateBy = updates->denseVector();
     745      addSequence = 0;
     746      multiplier = -1;
     747#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     748      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
     749#endif
     750    } else {
     751      number = spareColumn1->getNumElements();
     752      index = spareColumn1->getIndices();
     753      updateBy = spareColumn1->denseVector();
     754      addSequence = model_->maximumAbcNumberRows();
     755      multiplier = 1;
     756#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     757      slack_multiplier = 1.0;
     758#endif
     759    }
     760
     761    for (j = 0; j < number; j++) {
     762      int iSequence = index[j];
     763      double tableauValue = updateBy[iSequence];
     764      updateBy[iSequence] = 0.0;
     765      iSequence += addSequence;
     766      double value = reducedCost[iSequence];
     767      value -= multiplier * tableauValue;
     768      reducedCost[iSequence] = value;
     769      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     770
     771      switch (status) {
     772
     773      case AbcSimplex::basic:
     774        infeasible_->zero(iSequence);
     775      case AbcSimplex::isFixed:
     776        break;
     777      case AbcSimplex::isFree:
     778      case AbcSimplex::superBasic:
     779        if (fabs(value) > FREE_ACCEPT * tolerance) {
     780          // we are going to bias towards free (but only if reasonable)
     781          value *= FREE_BIAS;
     782          // store square in list
     783          if (infeas[iSequence])
     784            infeas[iSequence] = value * value; // already there
     785          else
     786            infeasible_->quickAdd(iSequence, value * value);
     787        } else {
     788          infeasible_->zero(iSequence);
     789        }
     790        break;
     791      case AbcSimplex::atUpperBound:
     792        if (value > tolerance) {
     793#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     794          value *= value * slack_multiplier;
    802795#else
    803                         value *= value;
    804 #endif
    805                          // store square in list
    806                          if (infeas[iSequence])
    807                               infeas[iSequence] = value; // already there
    808                          else
    809                               infeasible_->quickAdd(iSequence, value);
    810                     } else {
    811                          infeasible_->zero(iSequence);
    812                     }
    813                     break;
    814                case AbcSimplex::atLowerBound:
    815                     if (value < -tolerance) {
     796          value *= value;
     797#endif
     798          // store square in list
     799          if (infeas[iSequence])
     800            infeas[iSequence] = value; // already there
     801          else
     802            infeasible_->quickAdd(iSequence, value);
     803        } else {
     804          infeasible_->zero(iSequence);
     805        }
     806        break;
     807      case AbcSimplex::atLowerBound:
     808        if (value < -tolerance) {
    816809#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    817                         value *= value*slack_multiplier;
     810          value *= value * slack_multiplier;
    818811#else
    819                         value *= value;
    820 #endif
    821                          // store square in list
    822                          if (infeas[iSequence])
    823                               infeas[iSequence] = value; // already there
    824                          else
    825                               infeasible_->quickAdd(iSequence, value);
    826                     } else {
    827                          infeasible_->zero(iSequence);
    828                     }
    829                }
    830           }
    831      }
    832      updates->setNumElements(0);
    833      spareColumn1->setNumElements(0);
    834      if (pivotRow >= 0) {
    835           // make sure infeasibility on incoming is 0.0
    836           int sequenceIn = model_->sequenceIn();
    837           infeasible_->zero(sequenceIn);
    838      }
     812          value *= value;
     813#endif
     814          // store square in list
     815          if (infeas[iSequence])
     816            infeas[iSequence] = value; // already there
     817          else
     818            infeasible_->quickAdd(iSequence, value);
     819        } else {
     820          infeasible_->zero(iSequence);
     821        }
     822      }
     823    }
     824  }
     825  updates->setNumElements(0);
     826  spareColumn1->setNumElements(0);
     827  if (pivotRow >= 0) {
     828    // make sure infeasibility on incoming is 0.0
     829    int sequenceIn = model_->sequenceIn();
     830    infeasible_->zero(sequenceIn);
     831  }
    839832}
    840833// Update djs, weights for Devex
    841 void
    842 AbcPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
    843                                      CoinIndexedVector * spareRow2,
    844                                      CoinIndexedVector * spareColumn1)
     834void AbcPrimalColumnSteepest::djsAndDevex(CoinIndexedVector *updates,
     835  CoinIndexedVector *spareRow2,
     836  CoinIndexedVector *spareColumn1)
    845837{
    846      int j;
    847      int number = 0;
    848      int * index;
    849      double * updateBy;
    850      double * reducedCost;
    851      double tolerance = model_->currentDualTolerance();
    852      // we can't really trust infeasibilities if there is dual error
    853      // this coding has to mimic coding in checkDualSolution
    854      double error = CoinMin(1.0e-2, model_->largestDualError());
    855      // allow tolerance at least slightly bigger than standard
    856      tolerance = tolerance + error;
    857      // for weights update we use pivotSequence
    858      // unset in case sub flip
    859      assert (pivotSequence_ >= 0);
    860      assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
    861      double scaleFactor = 1.0 / updates->denseVector()[pivotSequence_]; // as formula is with 1.0
    862      pivotSequence_ = -1;
    863      double * infeas = infeasible_->denseVector();
    864      //updates->scanAndPack();
    865      model_->factorization()->updateColumnTranspose(*updates);
    866      // and we can see if reference
    867      //double referenceIn = 0.0;
    868      int sequenceIn = model_->sequenceIn();
    869      //if (mode_ != 1 && reference(sequenceIn))
    870      //   referenceIn = 1.0;
    871      // save outgoing weight round update
    872      double outgoingWeight = 0.0;
    873      int sequenceOut = model_->sequenceOut();
    874      if (sequenceOut >= 0)
    875           outgoingWeight = weights_[sequenceOut];
    876 
    877      // put row of tableau in rowArray and columnArray
    878      model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
    879      // update weights
    880      double * weight;
    881      // rows
    882      reducedCost = model_->djRegion();
    883 
    884      number = updates->getNumElements();
    885      index = updates->getIndices();
    886      updateBy = updates->denseVector();
    887      weight = weights_;
    888      // Devex
    889      for (j = 0; j < number; j++) {
    890           double thisWeight;
    891           double pivot;
    892           double value3;
    893           int iSequence = index[j];
    894           double value = reducedCost[iSequence];
    895           double value2 = updateBy[iSequence];
    896           updateBy[iSequence] = 0.0;
    897           value -= -value2;
    898           reducedCost[iSequence] = value;
    899           AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    900 
    901           switch(status) {
    902 
    903           case AbcSimplex::basic:
    904                infeasible_->zero(iSequence);
    905           case AbcSimplex::isFixed:
    906                break;
    907           case AbcSimplex::isFree:
    908           case AbcSimplex::superBasic:
    909                thisWeight = weight[iSequence];
    910                // row has -1
    911                pivot = value2 * scaleFactor;
    912                value3 = pivot * pivot * devex_;
    913                if (reference(iSequence))
    914                     value3 += 1.0;
    915                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    916                if (fabs(value) > FREE_ACCEPT * tolerance) {
    917                     // we are going to bias towards free (but only if reasonable)
    918                     value *= FREE_BIAS;
    919                     // store square in list
    920                     if (infeas[iSequence])
    921                          infeas[iSequence] = value * value; // already there
    922                     else
    923                          infeasible_->quickAdd(iSequence , value * value);
    924                } else {
    925                     infeasible_->zero(iSequence);
    926                }
    927                break;
    928           case AbcSimplex::atUpperBound:
    929                thisWeight = weight[iSequence];
    930                // row has -1
    931                pivot = value2 * scaleFactor;
    932                value3 = pivot * pivot * devex_;
    933                if (reference(iSequence))
    934                     value3 += 1.0;
    935                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    936                if (value > tolerance) {
    937                     // store square in list
     838  int j;
     839  int number = 0;
     840  int *index;
     841  double *updateBy;
     842  double *reducedCost;
     843  double tolerance = model_->currentDualTolerance();
     844  // we can't really trust infeasibilities if there is dual error
     845  // this coding has to mimic coding in checkDualSolution
     846  double error = CoinMin(1.0e-2, model_->largestDualError());
     847  // allow tolerance at least slightly bigger than standard
     848  tolerance = tolerance + error;
     849  // for weights update we use pivotSequence
     850  // unset in case sub flip
     851  assert(pivotSequence_ >= 0);
     852  assert(model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
     853  double scaleFactor = 1.0 / updates->denseVector()[pivotSequence_]; // as formula is with 1.0
     854  pivotSequence_ = -1;
     855  double *infeas = infeasible_->denseVector();
     856  //updates->scanAndPack();
     857  model_->factorization()->updateColumnTranspose(*updates);
     858  // and we can see if reference
     859  //double referenceIn = 0.0;
     860  int sequenceIn = model_->sequenceIn();
     861  //if (mode_ != 1 && reference(sequenceIn))
     862  //   referenceIn = 1.0;
     863  // save outgoing weight round update
     864  double outgoingWeight = 0.0;
     865  int sequenceOut = model_->sequenceOut();
     866  if (sequenceOut >= 0)
     867    outgoingWeight = weights_[sequenceOut];
     868
     869  // put row of tableau in rowArray and columnArray
     870  model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
     871  // update weights
     872  double *weight;
     873  // rows
     874  reducedCost = model_->djRegion();
     875
     876  number = updates->getNumElements();
     877  index = updates->getIndices();
     878  updateBy = updates->denseVector();
     879  weight = weights_;
     880  // Devex
     881  for (j = 0; j < number; j++) {
     882    double thisWeight;
     883    double pivot;
     884    double value3;
     885    int iSequence = index[j];
     886    double value = reducedCost[iSequence];
     887    double value2 = updateBy[iSequence];
     888    updateBy[iSequence] = 0.0;
     889    value -= -value2;
     890    reducedCost[iSequence] = value;
     891    AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     892
     893    switch (status) {
     894
     895    case AbcSimplex::basic:
     896      infeasible_->zero(iSequence);
     897    case AbcSimplex::isFixed:
     898      break;
     899    case AbcSimplex::isFree:
     900    case AbcSimplex::superBasic:
     901      thisWeight = weight[iSequence];
     902      // row has -1
     903      pivot = value2 * scaleFactor;
     904      value3 = pivot * pivot * devex_;
     905      if (reference(iSequence))
     906        value3 += 1.0;
     907      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     908      if (fabs(value) > FREE_ACCEPT * tolerance) {
     909        // we are going to bias towards free (but only if reasonable)
     910        value *= FREE_BIAS;
     911        // store square in list
     912        if (infeas[iSequence])
     913          infeas[iSequence] = value * value; // already there
     914        else
     915          infeasible_->quickAdd(iSequence, value * value);
     916      } else {
     917        infeasible_->zero(iSequence);
     918      }
     919      break;
     920    case AbcSimplex::atUpperBound:
     921      thisWeight = weight[iSequence];
     922      // row has -1
     923      pivot = value2 * scaleFactor;
     924      value3 = pivot * pivot * devex_;
     925      if (reference(iSequence))
     926        value3 += 1.0;
     927      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     928      if (value > tolerance) {
     929        // store square in list
    938930#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    939                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
     931        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
    940932#else
    941                     value *= value;
    942 #endif
    943                     if (infeas[iSequence])
    944                          infeas[iSequence] = value; // already there
    945                     else
    946                          infeasible_->quickAdd(iSequence , value);
    947                } else {
    948                     infeasible_->zero(iSequence);
    949                }
    950                break;
    951           case AbcSimplex::atLowerBound:
    952                thisWeight = weight[iSequence];
    953                // row has -1
    954                pivot = value2 * scaleFactor;
    955                value3 = pivot * pivot * devex_;
    956                if (reference(iSequence))
    957                     value3 += 1.0;
    958                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    959                if (value < -tolerance) {
    960                     // store square in list
     933        value *= value;
     934#endif
     935        if (infeas[iSequence])
     936          infeas[iSequence] = value; // already there
     937        else
     938          infeasible_->quickAdd(iSequence, value);
     939      } else {
     940        infeasible_->zero(iSequence);
     941      }
     942      break;
     943    case AbcSimplex::atLowerBound:
     944      thisWeight = weight[iSequence];
     945      // row has -1
     946      pivot = value2 * scaleFactor;
     947      value3 = pivot * pivot * devex_;
     948      if (reference(iSequence))
     949        value3 += 1.0;
     950      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     951      if (value < -tolerance) {
     952        // store square in list
    961953#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    962                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
     954        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
    963955#else
    964                     value *= value;
    965 #endif
    966                     if (infeas[iSequence])
    967                          infeas[iSequence] = value; // already there
    968                     else
    969                          infeasible_->quickAdd(iSequence , value);
    970                } else {
    971                     infeasible_->zero(iSequence);
    972                }
    973           }
    974      }
    975 
    976      // columns
    977      int addSequence = model_->maximumAbcNumberRows();
    978 
    979      scaleFactor = -scaleFactor;
    980      number = spareColumn1->getNumElements();
    981      index = spareColumn1->getIndices();
    982      updateBy = spareColumn1->denseVector();
    983 
    984      // Devex
    985 
    986      for (j = 0; j < number; j++) {
    987           double thisWeight;
    988           double pivot;
    989           double value3;
    990           int iSequence = index[j];
    991           double value2=updateBy[iSequence];
    992           updateBy[iSequence] = 0.0;
    993           iSequence += addSequence;
    994           double value = reducedCost[iSequence];
    995           value -= value2;
    996           reducedCost[iSequence] = value;
    997           AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    998 
    999           switch(status) {
    1000 
    1001           case AbcSimplex::basic:
    1002                infeasible_->zero(iSequence);
    1003           case AbcSimplex::isFixed:
    1004                break;
    1005           case AbcSimplex::isFree:
    1006           case AbcSimplex::superBasic:
    1007                thisWeight = weight[iSequence];
    1008                // row has -1
    1009                pivot = value2 * scaleFactor;
    1010                value3 = pivot * pivot * devex_;
    1011                if (reference(iSequence))
    1012                     value3 += 1.0;
    1013                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1014                if (fabs(value) > FREE_ACCEPT * tolerance) {
    1015                     // we are going to bias towards free (but only if reasonable)
    1016                     value *= FREE_BIAS;
    1017                     // store square in list
    1018                     if (infeas[iSequence])
    1019                          infeas[iSequence] = value * value; // already there
    1020                     else
    1021                          infeasible_->quickAdd(iSequence, value * value);
    1022                } else {
    1023                     infeasible_->zero(iSequence);
    1024                }
    1025                break;
    1026           case AbcSimplex::atUpperBound:
    1027                thisWeight = weight[iSequence];
    1028                // row has -1
    1029                pivot = value2 * scaleFactor;
    1030                value3 = pivot * pivot * devex_;
    1031                if (reference(iSequence))
    1032                     value3 += 1.0;
    1033                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1034                if (value > tolerance) {
    1035                     // store square in list
    1036                     if (infeas[iSequence])
    1037                          infeas[iSequence] = value * value; // already there
    1038                     else
    1039                          infeasible_->quickAdd(iSequence, value * value);
    1040                } else {
    1041                     infeasible_->zero(iSequence);
    1042                }
    1043                break;
    1044           case AbcSimplex::atLowerBound:
    1045                thisWeight = weight[iSequence];
    1046                // row has -1
    1047                pivot = value2 * scaleFactor;
    1048                value3 = pivot * pivot * devex_;
    1049                if (reference(iSequence))
    1050                     value3 += 1.0;
    1051                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1052                if (value < -tolerance) {
    1053                     // store square in list
    1054                     if (infeas[iSequence])
    1055                          infeas[iSequence] = value * value; // already there
    1056                     else
    1057                          infeasible_->quickAdd(iSequence, value * value);
    1058                } else {
    1059                     infeasible_->zero(iSequence);
    1060                }
    1061           }
    1062      }
    1063      // restore outgoing weight
    1064      if (sequenceOut >= 0)
    1065           weights_[sequenceOut] = outgoingWeight;
    1066      // make sure infeasibility on incoming is 0.0
    1067      infeasible_->zero(sequenceIn);
    1068      spareRow2->setNumElements(0);
    1069      //#define SOME_DEBUG_1
     956        value *= value;
     957#endif
     958        if (infeas[iSequence])
     959          infeas[iSequence] = value; // already there
     960        else
     961          infeasible_->quickAdd(iSequence, value);
     962      } else {
     963        infeasible_->zero(iSequence);
     964      }
     965    }
     966  }
     967
     968  // columns
     969  int addSequence = model_->maximumAbcNumberRows();
     970
     971  scaleFactor = -scaleFactor;
     972  number = spareColumn1->getNumElements();
     973  index = spareColumn1->getIndices();
     974  updateBy = spareColumn1->denseVector();
     975
     976  // Devex
     977
     978  for (j = 0; j < number; j++) {
     979    double thisWeight;
     980    double pivot;
     981    double value3;
     982    int iSequence = index[j];
     983    double value2 = updateBy[iSequence];
     984    updateBy[iSequence] = 0.0;
     985    iSequence += addSequence;
     986    double value = reducedCost[iSequence];
     987    value -= value2;
     988    reducedCost[iSequence] = value;
     989    AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     990
     991    switch (status) {
     992
     993    case AbcSimplex::basic:
     994      infeasible_->zero(iSequence);
     995    case AbcSimplex::isFixed:
     996      break;
     997    case AbcSimplex::isFree:
     998    case AbcSimplex::superBasic:
     999      thisWeight = weight[iSequence];
     1000      // row has -1
     1001      pivot = value2 * scaleFactor;
     1002      value3 = pivot * pivot * devex_;
     1003      if (reference(iSequence))
     1004        value3 += 1.0;
     1005      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1006      if (fabs(value) > FREE_ACCEPT * tolerance) {
     1007        // we are going to bias towards free (but only if reasonable)
     1008        value *= FREE_BIAS;
     1009        // store square in list
     1010        if (infeas[iSequence])
     1011          infeas[iSequence] = value * value; // already there
     1012        else
     1013          infeasible_->quickAdd(iSequence, value * value);
     1014      } else {
     1015        infeasible_->zero(iSequence);
     1016      }
     1017      break;
     1018    case AbcSimplex::atUpperBound:
     1019      thisWeight = weight[iSequence];
     1020      // row has -1
     1021      pivot = value2 * scaleFactor;
     1022      value3 = pivot * pivot * devex_;
     1023      if (reference(iSequence))
     1024        value3 += 1.0;
     1025      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1026      if (value > tolerance) {
     1027        // store square in list
     1028        if (infeas[iSequence])
     1029          infeas[iSequence] = value * value; // already there
     1030        else
     1031          infeasible_->quickAdd(iSequence, value * value);
     1032      } else {
     1033        infeasible_->zero(iSequence);
     1034      }
     1035      break;
     1036    case AbcSimplex::atLowerBound:
     1037      thisWeight = weight[iSequence];
     1038      // row has -1
     1039      pivot = value2 * scaleFactor;
     1040      value3 = pivot * pivot * devex_;
     1041      if (reference(iSequence))
     1042        value3 += 1.0;
     1043      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1044      if (value < -tolerance) {
     1045        // store square in list
     1046        if (infeas[iSequence])
     1047          infeas[iSequence] = value * value; // already there
     1048        else
     1049          infeasible_->quickAdd(iSequence, value * value);
     1050      } else {
     1051        infeasible_->zero(iSequence);
     1052      }
     1053    }
     1054  }
     1055  // restore outgoing weight
     1056  if (sequenceOut >= 0)
     1057    weights_[sequenceOut] = outgoingWeight;
     1058  // make sure infeasibility on incoming is 0.0
     1059  infeasible_->zero(sequenceIn);
     1060  spareRow2->setNumElements(0);
     1061  //#define SOME_DEBUG_1
    10701062#ifdef SOME_DEBUG_1
    1071      // check for accuracy
    1072      int iCheck = 892;
    1073      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1074      int numberRows = model_->numberRows();
    1075      //int numberColumns = model_->numberColumns();
    1076      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1077           if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
    1078                     !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
    1079                checkAccuracy(iCheck, 1.0e-1, updates);
    1080      }
    1081 #endif
    1082      updates->setNumElements(0);
    1083      spareColumn1->setNumElements(0);
     1063  // check for accuracy
     1064  int iCheck = 892;
     1065  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1066  int numberRows = model_->numberRows();
     1067  //int numberColumns = model_->numberColumns();
     1068  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1069    if (model_->getInternalStatus(iCheck) != AbcSimplex::basic && !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
     1070      checkAccuracy(iCheck, 1.0e-1, updates);
     1071  }
     1072#endif
     1073  updates->setNumElements(0);
     1074  spareColumn1->setNumElements(0);
    10841075}
    10851076// Update djs, weights for Devex
    1086 void
    1087 AbcPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
    1088                                       CoinIndexedVector * spareColumn1)
     1077void AbcPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector *updates,
     1078  CoinIndexedVector *spareColumn1)
    10891079{
    1090      int iSection, j;
    1091      int number = 0;
    1092      double multiplier;
    1093      int * index;
    1094      double * updateBy;
    1095      double * reducedCost;
    1096      // dj could be very small (or even zero - take care)
    1097      double dj = model_->dualIn();
    1098      double tolerance = model_->currentDualTolerance();
    1099      // we can't really trust infeasibilities if there is dual error
    1100      // this coding has to mimic coding in checkDualSolution
    1101      double error = CoinMin(1.0e-2, model_->largestDualError());
    1102      // allow tolerance at least slightly bigger than standard
    1103      tolerance = tolerance +  error;
    1104      int pivotRow = model_->pivotRow();
    1105      double * infeas = infeasible_->denseVector();
    1106      //updates->scanAndPack();
    1107      model_->factorization()->updateColumnTranspose(*updates);
    1108 
    1109      // put row of tableau in rowArray and columnArray
    1110      model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
    1111      // normal
    1112      reducedCost = model_->djRegion();
    1113      for (iSection = 0; iSection < 2; iSection++) {
    1114 
    1115           int addSequence;
    1116 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1117           double slack_multiplier;
    1118 #endif
    1119 
    1120           if (!iSection) {
    1121                number = updates->getNumElements();
    1122                index = updates->getIndices();
    1123                updateBy = updates->denseVector();
    1124                addSequence = 0;
    1125                multiplier=-1;
    1126 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1127                slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
    1128 #endif
    1129           } else {
    1130                number = spareColumn1->getNumElements();
    1131                index = spareColumn1->getIndices();
    1132                updateBy = spareColumn1->denseVector();
    1133                addSequence = model_->maximumAbcNumberRows();
    1134                multiplier=1;
    1135 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1136                slack_multiplier = 1.0;
    1137 #endif
    1138           }
    1139 
    1140           for (j = 0; j < number; j++) {
    1141                int iSequence = index[j];
    1142                double tableauValue=updateBy[iSequence];
    1143                updateBy[iSequence] = 0.0;
    1144                iSequence += addSequence;
    1145                double value = reducedCost[iSequence];
    1146                value -= multiplier*tableauValue;
    1147                AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    1148 
    1149                switch(status) {
    1150 
    1151                case AbcSimplex::basic:
    1152                     infeasible_->zero(iSequence);
    1153                case AbcSimplex::isFixed:
    1154                     break;
    1155                case AbcSimplex::isFree:
    1156                case AbcSimplex::superBasic:
    1157                  reducedCost[iSequence] = value;
    1158                     if (fabs(value) > FREE_ACCEPT * tolerance) {
    1159                          // we are going to bias towards free (but only if reasonable)
    1160                          value *= FREE_BIAS;
    1161                          // store square in list
    1162                          if (infeas[iSequence])
    1163                               infeas[iSequence] = value * value; // already there
    1164                          else
    1165                               infeasible_->quickAdd(iSequence , value * value);
    1166                     } else {
    1167                          infeasible_->zero(iSequence );
    1168                     }
    1169                     break;
    1170                case AbcSimplex::atUpperBound:
    1171                  reducedCost[iSequence] = value;
    1172                     if (value > tolerance) {
     1080  int iSection, j;
     1081  int number = 0;
     1082  double multiplier;
     1083  int *index;
     1084  double *updateBy;
     1085  double *reducedCost;
     1086  // dj could be very small (or even zero - take care)
     1087  double dj = model_->dualIn();
     1088  double tolerance = model_->currentDualTolerance();
     1089  // we can't really trust infeasibilities if there is dual error
     1090  // this coding has to mimic coding in checkDualSolution
     1091  double error = CoinMin(1.0e-2, model_->largestDualError());
     1092  // allow tolerance at least slightly bigger than standard
     1093  tolerance = tolerance + error;
     1094  int pivotRow = model_->pivotRow();
     1095  double *infeas = infeasible_->denseVector();
     1096  //updates->scanAndPack();
     1097  model_->factorization()->updateColumnTranspose(*updates);
     1098
     1099  // put row of tableau in rowArray and columnArray
     1100  model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
     1101  // normal
     1102  reducedCost = model_->djRegion();
     1103  for (iSection = 0; iSection < 2; iSection++) {
     1104
     1105    int addSequence;
    11731106#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1174                         value *= value*slack_multiplier;
     1107    double slack_multiplier;
     1108#endif
     1109
     1110    if (!iSection) {
     1111      number = updates->getNumElements();
     1112      index = updates->getIndices();
     1113      updateBy = updates->denseVector();
     1114      addSequence = 0;
     1115      multiplier = -1;
     1116#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1117      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
     1118#endif
     1119    } else {
     1120      number = spareColumn1->getNumElements();
     1121      index = spareColumn1->getIndices();
     1122      updateBy = spareColumn1->denseVector();
     1123      addSequence = model_->maximumAbcNumberRows();
     1124      multiplier = 1;
     1125#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1126      slack_multiplier = 1.0;
     1127#endif
     1128    }
     1129
     1130    for (j = 0; j < number; j++) {
     1131      int iSequence = index[j];
     1132      double tableauValue = updateBy[iSequence];
     1133      updateBy[iSequence] = 0.0;
     1134      iSequence += addSequence;
     1135      double value = reducedCost[iSequence];
     1136      value -= multiplier * tableauValue;
     1137      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     1138
     1139      switch (status) {
     1140
     1141      case AbcSimplex::basic:
     1142        infeasible_->zero(iSequence);
     1143      case AbcSimplex::isFixed:
     1144        break;
     1145      case AbcSimplex::isFree:
     1146      case AbcSimplex::superBasic:
     1147        reducedCost[iSequence] = value;
     1148        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1149          // we are going to bias towards free (but only if reasonable)
     1150          value *= FREE_BIAS;
     1151          // store square in list
     1152          if (infeas[iSequence])
     1153            infeas[iSequence] = value * value; // already there
     1154          else
     1155            infeasible_->quickAdd(iSequence, value * value);
     1156        } else {
     1157          infeasible_->zero(iSequence);
     1158        }
     1159        break;
     1160      case AbcSimplex::atUpperBound:
     1161        reducedCost[iSequence] = value;
     1162        if (value > tolerance) {
     1163#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1164          value *= value * slack_multiplier;
    11751165#else
    1176                         value *= value;
    1177 #endif
    1178                          // store square in list
    1179                          if (infeas[iSequence])
    1180                               infeas[iSequence] = value; // already there
    1181                          else
    1182                               infeasible_->quickAdd(iSequence, value);
    1183                     } else {
    1184                          infeasible_->zero(iSequence);
    1185                     }
    1186                     break;
    1187                case AbcSimplex::atLowerBound:
    1188                 reducedCost[iSequence] = value;
    1189                     if (value < -tolerance) {
     1166          value *= value;
     1167#endif
     1168          // store square in list
     1169          if (infeas[iSequence])
     1170            infeas[iSequence] = value; // already there
     1171          else
     1172            infeasible_->quickAdd(iSequence, value);
     1173        } else {
     1174          infeasible_->zero(iSequence);
     1175        }
     1176        break;
     1177      case AbcSimplex::atLowerBound:
     1178        reducedCost[iSequence] = value;
     1179        if (value < -tolerance) {
    11901180#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1191                         value *= value*slack_multiplier;
     1181          value *= value * slack_multiplier;
    11921182#else
    1193                         value *= value;
    1194 #endif
    1195                          // store square in list
    1196                          if (infeas[iSequence])
    1197                               infeas[iSequence] = value; // already there
    1198                          else
    1199                               infeasible_->quickAdd(iSequence, value);
    1200                     } else {
    1201                          infeasible_->zero(iSequence);
    1202                     }
    1203                }
    1204           }
    1205      }
    1206      // They are empty
    1207      updates->setNumElements(0);
    1208      spareColumn1->setNumElements(0);
    1209      // make sure infeasibility on incoming is 0.0
    1210      int sequenceIn = model_->sequenceIn();
    1211      infeasible_->zero(sequenceIn);
    1212      // for weights update we use pivotSequence
    1213      if (pivotSequence_ >= 0) {
    1214           pivotRow = pivotSequence_;
    1215           // unset in case sub flip
    1216           pivotSequence_ = -1;
    1217           // make sure infeasibility on incoming is 0.0
    1218           const int * pivotVariable = model_->pivotVariable();
    1219           sequenceIn = pivotVariable[pivotRow];
    1220           infeasible_->zero(sequenceIn);
    1221           // and we can see if reference
    1222           //double referenceIn = 0.0;
    1223           //if (mode_ != 1 && reference(sequenceIn))
    1224           //   referenceIn = 1.0;
    1225           // save outgoing weight round update
    1226           double outgoingWeight = 0.0;
    1227           int sequenceOut = model_->sequenceOut();
    1228           if (sequenceOut >= 0)
    1229                outgoingWeight = weights_[sequenceOut];
    1230           // update weights
    1231           updates->setNumElements(0);
    1232           spareColumn1->setNumElements(0);
    1233           // might as well set dj to 1
    1234           dj = 1.0;
    1235           updates->insert(pivotRow, -dj);
    1236           model_->factorization()->updateColumnTranspose(*updates);
    1237           // put row of tableau in rowArray and columnArray
    1238           model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
    1239           double * weight;
    1240           int numberColumns = model_->numberColumns();
    1241           // rows
    1242           number = updates->getNumElements();
    1243           index = updates->getIndices();
    1244           updateBy = updates->denseVector();
    1245           weight = weights_;
    1246 
    1247           assert (devex_ > 0.0);
    1248           for (j = 0; j < number; j++) {
    1249                int iSequence = index[j];
    1250                double thisWeight = weight[iSequence];
    1251                // row has -1
    1252                double pivot = - updateBy[iSequence];
    1253                updateBy[iSequence] = 0.0;
    1254                double value = pivot * pivot * devex_;
    1255                if (reference(iSequence + numberColumns))
    1256                     value += 1.0;
    1257                weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1258           }
    1259 
    1260           // columns
    1261           number = spareColumn1->getNumElements();
    1262           index = spareColumn1->getIndices();
    1263           updateBy = spareColumn1->denseVector();
    1264           int addSequence = model_->maximumAbcNumberRows();
    1265           for (j = 0; j < number; j++) {
    1266                int iSequence = index[j];
    1267                double pivot=updateBy[iSequence];
    1268                updateBy[iSequence] = 0.0;
    1269                iSequence += addSequence;
    1270                double thisWeight = weight[iSequence];
    1271                // row has -1
    1272                double value = pivot * pivot * devex_;
    1273                if (reference(iSequence))
    1274                     value += 1.0;
    1275                weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1276           }
    1277           // restore outgoing weight
    1278           if (sequenceOut >= 0)
    1279                weights_[sequenceOut] = outgoingWeight;
    1280           //#define SOME_DEBUG_1
     1183          value *= value;
     1184#endif
     1185          // store square in list
     1186          if (infeas[iSequence])
     1187            infeas[iSequence] = value; // already there
     1188          else
     1189            infeasible_->quickAdd(iSequence, value);
     1190        } else {
     1191          infeasible_->zero(iSequence);
     1192        }
     1193      }
     1194    }
     1195  }
     1196  // They are empty
     1197  updates->setNumElements(0);
     1198  spareColumn1->setNumElements(0);
     1199  // make sure infeasibility on incoming is 0.0
     1200  int sequenceIn = model_->sequenceIn();
     1201  infeasible_->zero(sequenceIn);
     1202  // for weights update we use pivotSequence
     1203  if (pivotSequence_ >= 0) {
     1204    pivotRow = pivotSequence_;
     1205    // unset in case sub flip
     1206    pivotSequence_ = -1;
     1207    // make sure infeasibility on incoming is 0.0
     1208    const int *pivotVariable = model_->pivotVariable();
     1209    sequenceIn = pivotVariable[pivotRow];
     1210    infeasible_->zero(sequenceIn);
     1211    // and we can see if reference
     1212    //double referenceIn = 0.0;
     1213    //if (mode_ != 1 && reference(sequenceIn))
     1214    //   referenceIn = 1.0;
     1215    // save outgoing weight round update
     1216    double outgoingWeight = 0.0;
     1217    int sequenceOut = model_->sequenceOut();
     1218    if (sequenceOut >= 0)
     1219      outgoingWeight = weights_[sequenceOut];
     1220    // update weights
     1221    updates->setNumElements(0);
     1222    spareColumn1->setNumElements(0);
     1223    // might as well set dj to 1
     1224    dj = 1.0;
     1225    updates->insert(pivotRow, -dj);
     1226    model_->factorization()->updateColumnTranspose(*updates);
     1227    // put row of tableau in rowArray and columnArray
     1228    model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
     1229    double *weight;
     1230    int numberColumns = model_->numberColumns();
     1231    // rows
     1232    number = updates->getNumElements();
     1233    index = updates->getIndices();
     1234    updateBy = updates->denseVector();
     1235    weight = weights_;
     1236
     1237    assert(devex_ > 0.0);
     1238    for (j = 0; j < number; j++) {
     1239      int iSequence = index[j];
     1240      double thisWeight = weight[iSequence];
     1241      // row has -1
     1242      double pivot = -updateBy[iSequence];
     1243      updateBy[iSequence] = 0.0;
     1244      double value = pivot * pivot * devex_;
     1245      if (reference(iSequence + numberColumns))
     1246        value += 1.0;
     1247      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1248    }
     1249
     1250    // columns
     1251    number = spareColumn1->getNumElements();
     1252    index = spareColumn1->getIndices();
     1253    updateBy = spareColumn1->denseVector();
     1254    int addSequence = model_->maximumAbcNumberRows();
     1255    for (j = 0; j < number; j++) {
     1256      int iSequence = index[j];
     1257      double pivot = updateBy[iSequence];
     1258      updateBy[iSequence] = 0.0;
     1259      iSequence += addSequence;
     1260      double thisWeight = weight[iSequence];
     1261      // row has -1
     1262      double value = pivot * pivot * devex_;
     1263      if (reference(iSequence))
     1264        value += 1.0;
     1265      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1266    }
     1267    // restore outgoing weight
     1268    if (sequenceOut >= 0)
     1269      weights_[sequenceOut] = outgoingWeight;
     1270      //#define SOME_DEBUG_1
    12811271#ifdef SOME_DEBUG_1
    1282           // check for accuracy
    1283           int iCheck = 892;
    1284           //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1285           int numberRows = model_->numberRows();
    1286           //int numberColumns = model_->numberColumns();
    1287           for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1288                if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
    1289                          !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
    1290                     checkAccuracy(iCheck, 1.0e-1, updates);
    1291           }
    1292 #endif
    1293           updates->setNumElements(0);
    1294           spareColumn1->setNumElements(0);
    1295      }
     1272    // check for accuracy
     1273    int iCheck = 892;
     1274    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1275    int numberRows = model_->numberRows();
     1276    //int numberColumns = model_->numberColumns();
     1277    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1278      if (model_->getInternalStatus(iCheck) != AbcSimplex::basic && !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
     1279        checkAccuracy(iCheck, 1.0e-1, updates);
     1280    }
     1281#endif
     1282    updates->setNumElements(0);
     1283    spareColumn1->setNumElements(0);
     1284  }
    12961285}
    12971286// Update weights for Devex
    1298 void
    1299 AbcPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
    1300                                    CoinIndexedVector * spareColumn1)
     1287void AbcPrimalColumnSteepest::justDevex(CoinIndexedVector *updates,
     1288  CoinIndexedVector *spareColumn1)
    13011289{
    1302      int j;
    1303      int number = 0;
    1304      int * index;
    1305      double * updateBy;
    1306      // dj could be very small (or even zero - take care)
    1307      double dj = model_->dualIn();
    1308      double tolerance = model_->currentDualTolerance();
    1309      // we can't really trust infeasibilities if there is dual error
    1310      // this coding has to mimic coding in checkDualSolution
    1311      double error = CoinMin(1.0e-2, model_->largestDualError());
    1312      // allow tolerance at least slightly bigger than standard
    1313      tolerance = tolerance + error;
    1314      int pivotRow = model_->pivotRow();
    1315      // for weights update we use pivotSequence
    1316      pivotRow = pivotSequence_;
    1317      assert (pivotRow >= 0);
    1318      // make sure infeasibility on incoming is 0.0
    1319      const int * pivotVariable = model_->pivotVariable();
    1320      int sequenceIn = pivotVariable[pivotRow];
    1321      infeasible_->zero(sequenceIn);
    1322      // and we can see if reference
    1323      //double referenceIn = 0.0;
    1324      //if (mode_ != 1 && reference(sequenceIn))
    1325      //   referenceIn = 1.0;
    1326      // save outgoing weight round update
    1327      double outgoingWeight = 0.0;
    1328      int sequenceOut = model_->sequenceOut();
    1329      if (sequenceOut >= 0)
    1330           outgoingWeight = weights_[sequenceOut];
    1331      assert (!updates->getNumElements());
    1332      assert (!spareColumn1->getNumElements());
    1333      // unset in case sub flip
    1334      pivotSequence_ = -1;
    1335      // might as well set dj to 1
    1336      dj = -1.0;
    1337      updates->createUnpacked(1, &pivotRow, &dj);
    1338      model_->factorization()->updateColumnTranspose(*updates);
    1339      // put row of tableau in rowArray and columnArray
    1340      model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
    1341      double * weight;
    1342      // rows
    1343      number = updates->getNumElements();
    1344      index = updates->getIndices();
    1345      updateBy = updates->denseVector();
    1346      weight = weights_;
    1347 
    1348      // Devex
    1349      assert (devex_ > 0.0);
    1350      for (j = 0; j < number; j++) {
    1351           int iSequence = index[j];
    1352           double thisWeight = weight[iSequence];
    1353           // row has -1
    1354           double pivot = - updateBy[iSequence];
    1355           updateBy[iSequence] = 0.0;
    1356           double value = pivot * pivot * devex_;
    1357           if (reference(iSequence))
    1358                value += 1.0;
    1359           weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1360      }
    1361 
    1362      // columns
    1363      int addSequence = model_->maximumAbcNumberRows();
    1364      number = spareColumn1->getNumElements();
    1365      index = spareColumn1->getIndices();
    1366      updateBy = spareColumn1->denseVector();
    1367      // Devex
    1368      for (j = 0; j < number; j++) {
    1369           int iSequence = index[j];
    1370           double pivot=updateBy[iSequence];
    1371           updateBy[iSequence] = 0.0;
    1372           iSequence += addSequence;
    1373           double thisWeight = weight[iSequence];
    1374           // row has -1
    1375           double value = pivot * pivot * devex_;
    1376           if (reference(iSequence))
    1377                value += 1.0;
    1378           weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1379      }
    1380      // restore outgoing weight
    1381      if (sequenceOut >= 0)
    1382           weights_[sequenceOut] = outgoingWeight;
    1383      //#define SOME_DEBUG_1
     1290  int j;
     1291  int number = 0;
     1292  int *index;
     1293  double *updateBy;
     1294  // dj could be very small (or even zero - take care)
     1295  double dj = model_->dualIn();
     1296  double tolerance = model_->currentDualTolerance();
     1297  // we can't really trust infeasibilities if there is dual error
     1298  // this coding has to mimic coding in checkDualSolution
     1299  double error = CoinMin(1.0e-2, model_->largestDualError());
     1300  // allow tolerance at least slightly bigger than standard
     1301  tolerance = tolerance + error;
     1302  int pivotRow = model_->pivotRow();
     1303  // for weights update we use pivotSequence
     1304  pivotRow = pivotSequence_;
     1305  assert(pivotRow >= 0);
     1306  // make sure infeasibility on incoming is 0.0
     1307  const int *pivotVariable = model_->pivotVariable();
     1308  int sequenceIn = pivotVariable[pivotRow];
     1309  infeasible_->zero(sequenceIn);
     1310  // and we can see if reference
     1311  //double referenceIn = 0.0;
     1312  //if (mode_ != 1 && reference(sequenceIn))
     1313  //   referenceIn = 1.0;
     1314  // save outgoing weight round update
     1315  double outgoingWeight = 0.0;
     1316  int sequenceOut = model_->sequenceOut();
     1317  if (sequenceOut >= 0)
     1318    outgoingWeight = weights_[sequenceOut];
     1319  assert(!updates->getNumElements());
     1320  assert(!spareColumn1->getNumElements());
     1321  // unset in case sub flip
     1322  pivotSequence_ = -1;
     1323  // might as well set dj to 1
     1324  dj = -1.0;
     1325  updates->createUnpacked(1, &pivotRow, &dj);
     1326  model_->factorization()->updateColumnTranspose(*updates);
     1327  // put row of tableau in rowArray and columnArray
     1328  model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
     1329  double *weight;
     1330  // rows
     1331  number = updates->getNumElements();
     1332  index = updates->getIndices();
     1333  updateBy = updates->denseVector();
     1334  weight = weights_;
     1335
     1336  // Devex
     1337  assert(devex_ > 0.0);
     1338  for (j = 0; j < number; j++) {
     1339    int iSequence = index[j];
     1340    double thisWeight = weight[iSequence];
     1341    // row has -1
     1342    double pivot = -updateBy[iSequence];
     1343    updateBy[iSequence] = 0.0;
     1344    double value = pivot * pivot * devex_;
     1345    if (reference(iSequence))
     1346      value += 1.0;
     1347    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1348  }
     1349
     1350  // columns
     1351  int addSequence = model_->maximumAbcNumberRows();
     1352  number = spareColumn1->getNumElements();
     1353  index = spareColumn1->getIndices();
     1354  updateBy = spareColumn1->denseVector();
     1355  // Devex
     1356  for (j = 0; j < number; j++) {
     1357    int iSequence = index[j];
     1358    double pivot = updateBy[iSequence];
     1359    updateBy[iSequence] = 0.0;
     1360    iSequence += addSequence;
     1361    double thisWeight = weight[iSequence];
     1362    // row has -1
     1363    double value = pivot * pivot * devex_;
     1364    if (reference(iSequence))
     1365      value += 1.0;
     1366    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1367  }
     1368  // restore outgoing weight
     1369  if (sequenceOut >= 0)
     1370    weights_[sequenceOut] = outgoingWeight;
     1371    //#define SOME_DEBUG_1
    13841372#ifdef SOME_DEBUG_1
    1385      // check for accuracy
    1386      int iCheck = 892;
    1387      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1388      int numberRows = model_->numberRows();
    1389      //int numberColumns = model_->numberColumns();
    1390      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1391           if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
    1392                     !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
    1393                checkAccuracy(iCheck, 1.0e-1, updates);
    1394      }
    1395 #endif
    1396      updates->setNumElements(0);
    1397      spareColumn1->setNumElements(0);
     1373  // check for accuracy
     1374  int iCheck = 892;
     1375  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1376  int numberRows = model_->numberRows();
     1377  //int numberColumns = model_->numberColumns();
     1378  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1379    if (model_->getInternalStatus(iCheck) != AbcSimplex::basic && !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
     1380      checkAccuracy(iCheck, 1.0e-1, updates);
     1381  }
     1382#endif
     1383  updates->setNumElements(0);
     1384  spareColumn1->setNumElements(0);
    13981385}
    13991386// Called when maximum pivots changes
    1400 void
    1401 AbcPrimalColumnSteepest::maximumPivotsChanged()
     1387void AbcPrimalColumnSteepest::maximumPivotsChanged()
    14021388{
    1403      if (alternateWeights_ &&
    1404                alternateWeights_->capacity() != model_->numberRows() +
    1405                model_->factorization()->maximumPivots()) {
    1406           delete alternateWeights_;
    1407           alternateWeights_ = new CoinIndexedVector();
    1408           // enough space so can use it for factorization
     1389  if (alternateWeights_ && alternateWeights_->capacity() != model_->numberRows() + model_->factorization()->maximumPivots()) {
     1390    delete alternateWeights_;
     1391    alternateWeights_ = new CoinIndexedVector();
     1392    // enough space so can use it for factorization
    14091393#ifndef ABC_USE_COIN_FACTORIZATION
    1410           // enough for ordered
    1411           alternateWeights_->reserve(2*model_->numberRows() +
    1412                                      model_->factorization()->maximumPivots());
     1394    // enough for ordered
     1395    alternateWeights_->reserve(2 * model_->numberRows() + model_->factorization()->maximumPivots());
    14131396#else
    1414           int n=(2*model_->numberRows()+
    1415              model_->factorization()->maximumPivots()+7)&~3;
    1416           alternateWeights_->reserve(n);
    1417           // zero
    1418           CoinZeroN(alternateWeights_->getIndices(),
    1419                     alternateWeights_->capacity());
    1420 #endif
    1421      }
     1397    int n = (2 * model_->numberRows() + model_->factorization()->maximumPivots() + 7) & ~3;
     1398    alternateWeights_->reserve(n);
     1399    // zero
     1400    CoinZeroN(alternateWeights_->getIndices(),
     1401      alternateWeights_->capacity());
     1402#endif
     1403  }
    14221404}
    14231405/*
     
    14281410   5) at end of values pass (so need initialization)
    14291411*/
    1430 void
    1431 AbcPrimalColumnSteepest::saveWeights(AbcSimplex * model, int mode)
     1412void AbcPrimalColumnSteepest::saveWeights(AbcSimplex *model, int mode)
    14321413{
    14331414  model_ = model;
     
    14401421  int numberRows = model_->numberRows();
    14411422  int numberColumns = model_->numberColumns();
    1442   int maximumRows=model_->maximumAbcNumberRows();
    1443   const int * pivotVariable = model_->pivotVariable();
     1423  int maximumRows = model_->maximumAbcNumberRows();
     1424  const int *pivotVariable = model_->pivotVariable();
    14441425  bool doInfeasibilities = true;
    14451426  if (mode == 1) {
    1446     if(weights_) {
     1427    if (weights_) {
    14471428      // Check if size has changed
    1448       if (infeasible_->capacity() == numberRows + numberColumns &&
    1449           alternateWeights_->capacity() == 2*numberRows +
    1450           model_->factorization()->maximumPivots()) {
    1451         //alternateWeights_->clear();
    1452         if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
    1453           // save pivot order
    1454           CoinMemcpyN(pivotVariable,
    1455                       numberRows, alternateWeights_->getIndices());
    1456           // change from pivot row number to sequence number
    1457           pivotSequence_ = pivotVariable[pivotSequence_];
    1458         } else {
    1459           pivotSequence_ = -1;
    1460         }
    1461         state_ = 1;
     1429      if (infeasible_->capacity() == numberRows + numberColumns && alternateWeights_->capacity() == 2 * numberRows + model_->factorization()->maximumPivots()) {
     1430        //alternateWeights_->clear();
     1431        if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
     1432          // save pivot order
     1433          CoinMemcpyN(pivotVariable,
     1434            numberRows, alternateWeights_->getIndices());
     1435          // change from pivot row number to sequence number
     1436          pivotSequence_ = pivotVariable[pivotSequence_];
     1437        } else {
     1438          pivotSequence_ = -1;
     1439        }
     1440        state_ = 1;
    14621441      } else {
    1463         // size has changed - clear everything
    1464         delete [] weights_;
    1465         weights_ = NULL;
    1466         delete infeasible_;
    1467         infeasible_ = NULL;
    1468         delete alternateWeights_;
    1469         alternateWeights_ = NULL;
    1470         delete [] savedWeights_;
    1471         savedWeights_ = NULL;
    1472         delete [] reference_;
    1473         reference_ = NULL;
    1474         state_ = -1;
    1475         pivotSequence_ = -1;
     1442        // size has changed - clear everything
     1443        delete[] weights_;
     1444        weights_ = NULL;
     1445        delete infeasible_;
     1446        infeasible_ = NULL;
     1447        delete alternateWeights_;
     1448        alternateWeights_ = NULL;
     1449        delete[] savedWeights_;
     1450        savedWeights_ = NULL;
     1451        delete[] reference_;
     1452        reference_ = NULL;
     1453        state_ = -1;
     1454        pivotSequence_ = -1;
    14761455      }
    14771456    }
     
    14801459    if (!weights_ || state_ == -1 || mode == 5) {
    14811460      // Partial is only allowed with certain types of matrix
    1482       if ((mode_ != 4 && mode_ != 5) || numberSwitched_ ) {
    1483         // initialize weights
    1484         delete [] weights_;
    1485         delete alternateWeights_;
    1486         weights_ = new double[numberRows+numberColumns];
    1487         alternateWeights_ = new CoinIndexedVector();
    1488         // enough space so can use it for factorization
     1461      if ((mode_ != 4 && mode_ != 5) || numberSwitched_) {
     1462        // initialize weights
     1463        delete[] weights_;
     1464        delete alternateWeights_;
     1465        weights_ = new double[numberRows + numberColumns];
     1466        alternateWeights_ = new CoinIndexedVector();
     1467        // enough space so can use it for factorization
    14891468#ifndef ABC_USE_COIN_FACTORIZATION
    1490         alternateWeights_->reserve(2*numberRows +
    1491                                    model_->factorization()->maximumPivots());
     1469        alternateWeights_->reserve(2 * numberRows + model_->factorization()->maximumPivots());
    14921470#else
    1493         int n=(2*model_->numberRows()+
    1494            model_->factorization()->maximumPivots()+7)&~3;
    1495         alternateWeights_->reserve(n);
    1496         // zero
    1497         CoinZeroN(alternateWeights_->getIndices(),
    1498                   alternateWeights_->capacity());
    1499 #endif
    1500         initializeWeights();
    1501         // create saved weights
    1502         delete [] savedWeights_;
    1503         savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
    1504         // just do initialization
    1505         mode = 3;
     1471        int n = (2 * model_->numberRows() + model_->factorization()->maximumPivots() + 7) & ~3;
     1472        alternateWeights_->reserve(n);
     1473        // zero
     1474        CoinZeroN(alternateWeights_->getIndices(),
     1475          alternateWeights_->capacity());
     1476#endif
     1477        initializeWeights();
     1478        // create saved weights
     1479        delete[] savedWeights_;
     1480        savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
     1481        // just do initialization
     1482        mode = 3;
    15061483      } else {
    1507         // Partial pricing
    1508         // use region as somewhere to save non-fixed slacks
    1509         // set up infeasibilities
    1510         if (!infeasible_) {
    1511           infeasible_ = new CoinIndexedVector();
    1512           infeasible_->reserve(numberColumns + numberRows);
    1513         }
    1514         infeasible_->clear();
    1515         int number = model_->numberRows();
    1516         int iSequence;
    1517         int numberLook = 0;
    1518         int * which = infeasible_->getIndices();
    1519         for (iSequence = 0; iSequence < number; iSequence++) {
    1520           AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    1521           if (status != AbcSimplex::isFixed)
    1522             which[numberLook++] = iSequence;
    1523         }
    1524         infeasible_->setNumElements(numberLook);
    1525         doInfeasibilities = false;
     1484        // Partial pricing
     1485        // use region as somewhere to save non-fixed slacks
     1486        // set up infeasibilities
     1487        if (!infeasible_) {
     1488          infeasible_ = new CoinIndexedVector();
     1489          infeasible_->reserve(numberColumns + numberRows);
     1490        }
     1491        infeasible_->clear();
     1492        int number = model_->numberRows();
     1493        int iSequence;
     1494        int numberLook = 0;
     1495        int *which = infeasible_->getIndices();
     1496        for (iSequence = 0; iSequence < number; iSequence++) {
     1497          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     1498          if (status != AbcSimplex::isFixed)
     1499            which[numberLook++] = iSequence;
     1500        }
     1501        infeasible_->setNumElements(numberLook);
     1502        doInfeasibilities = false;
    15261503      }
    15271504      savedPivotSequence_ = -2;
     
    15291506    } else {
    15301507      if (mode != 4) {
    1531         // save
    1532         CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
    1533         savedPivotSequence_ = pivotSequence_;
    1534         savedSequenceOut_ = model_->sequenceOut();
     1508        // save
     1509        CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
     1510        savedPivotSequence_ = pivotSequence_;
     1511        savedSequenceOut_ = model_->sequenceOut();
    15351512      } else {
    1536         // restore
    1537         CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
    1538         // was - but I think should not be
    1539         //pivotSequence_= savedPivotSequence_;
    1540         //model_->setSequenceOut(savedSequenceOut_);
    1541         pivotSequence_ = -1;
    1542         model_->setSequenceOut(-1);
    1543         // indices are wrong so clear by hand
    1544         //alternateWeights_->clear();
    1545         CoinZeroN(alternateWeights_->denseVector(),
    1546                   alternateWeights_->capacity());
    1547         alternateWeights_->setNumElements(0);
     1513        // restore
     1514        CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
     1515        // was - but I think should not be
     1516        //pivotSequence_= savedPivotSequence_;
     1517        //model_->setSequenceOut(savedSequenceOut_);
     1518        pivotSequence_ = -1;
     1519        model_->setSequenceOut(-1);
     1520        // indices are wrong so clear by hand
     1521        //alternateWeights_->clear();
     1522        CoinZeroN(alternateWeights_->denseVector(),
     1523          alternateWeights_->capacity());
     1524        alternateWeights_->setNumElements(0);
    15481525      }
    15491526    }
     
    15581535    if (mode != 3) {
    15591536      if (pivotSequence_ >= 0) {
    1560         // restore pivot row
    1561         int iRow;
    1562         // permute alternateWeights
    1563         int iVector=model_->getAvailableArray();
    1564         double * temp = model_->usefulArray(iVector)->denseVector();;
    1565         double * work = alternateWeights_->denseVector();
    1566         int * savePivotOrder = model_->usefulArray(iVector)->getIndices();
    1567         int * oldPivotOrder = alternateWeights_->getIndices();
    1568         for (iRow = 0; iRow < numberRows; iRow++) {
    1569           int iPivot = oldPivotOrder[iRow];
    1570           temp[iPivot] = work[iRow];
    1571           savePivotOrder[iRow] = iPivot;
    1572         }
    1573         int number = 0;
    1574         int found = -1;
    1575         int * which = oldPivotOrder;
    1576         // find pivot row and re-create alternateWeights
    1577         for (iRow = 0; iRow < numberRows; iRow++) {
    1578           int iPivot = pivotVariable[iRow];
    1579           if (iPivot == pivotSequence_)
    1580             found = iRow;
    1581           work[iRow] = temp[iPivot];
    1582           if (work[iRow])
    1583             which[number++] = iRow;
    1584         }
    1585         alternateWeights_->setNumElements(number);
     1537        // restore pivot row
     1538        int iRow;
     1539        // permute alternateWeights
     1540        int iVector = model_->getAvailableArray();
     1541        double *temp = model_->usefulArray(iVector)->denseVector();
     1542        ;
     1543        double *work = alternateWeights_->denseVector();
     1544        int *savePivotOrder = model_->usefulArray(iVector)->getIndices();
     1545        int *oldPivotOrder = alternateWeights_->getIndices();
     1546        for (iRow = 0; iRow < numberRows; iRow++) {
     1547          int iPivot = oldPivotOrder[iRow];
     1548          temp[iPivot] = work[iRow];
     1549          savePivotOrder[iRow] = iPivot;
     1550        }
     1551        int number = 0;
     1552        int found = -1;
     1553        int *which = oldPivotOrder;
     1554        // find pivot row and re-create alternateWeights
     1555        for (iRow = 0; iRow < numberRows; iRow++) {
     1556          int iPivot = pivotVariable[iRow];
     1557          if (iPivot == pivotSequence_)
     1558            found = iRow;
     1559          work[iRow] = temp[iPivot];
     1560          if (work[iRow])
     1561            which[number++] = iRow;
     1562        }
     1563        alternateWeights_->setNumElements(number);
    15861564#ifdef CLP_DEBUG
    1587         // Can happen but I should clean up
    1588         assert(found >= 0);
    1589 #endif
    1590         pivotSequence_ = found;
    1591         for (iRow = 0; iRow < numberRows; iRow++) {
    1592           int iPivot = savePivotOrder[iRow];
    1593           temp[iPivot] = 0.0;
    1594         }
    1595         model_->setAvailableArray(iVector);
     1565        // Can happen but I should clean up
     1566        assert(found >= 0);
     1567#endif
     1568        pivotSequence_ = found;
     1569        for (iRow = 0; iRow < numberRows; iRow++) {
     1570          int iPivot = savePivotOrder[iRow];
     1571          temp[iPivot] = 0.0;
     1572        }
     1573        model_->setAvailableArray(iVector);
    15961574      } else {
    1597         // Just clean up
    1598         if (alternateWeights_)
    1599           alternateWeights_->clear();
     1575        // Just clean up
     1576        if (alternateWeights_)
     1577          alternateWeights_->clear();
    16001578      }
    16011579    }
     
    16031581    if (!model->factorization()->pivots())
    16041582      sizeFactorization_ = model_->factorization()->numberElements();
    1605     if(!doInfeasibilities)
     1583    if (!doInfeasibilities)
    16061584      return; // don't disturb infeasibilities
    16071585    infeasible_->clear();
     
    16091587    int number = model_->numberRows() + model_->numberColumns();
    16101588    int iSequence;
    1611    
    1612     double * reducedCost = model_->djRegion();
    1613 #ifndef CLP_PRIMAL_SLACK_MULTIPLIER 
     1589
     1590    double *reducedCost = model_->djRegion();
     1591#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
    16141592    for (iSequence = 0; iSequence < number; iSequence++) {
    16151593      double value = reducedCost[iSequence];
    16161594      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    1617      
    1618       switch(status) {
    1619        
     1595
     1596      switch (status) {
     1597
    16201598      case AbcSimplex::basic:
    16211599      case AbcSimplex::isFixed:
    1622         break;
     1600        break;
    16231601      case AbcSimplex::isFree:
    16241602      case AbcSimplex::superBasic:
    1625         if (fabs(value) > FREE_ACCEPT * tolerance) {
    1626           // we are going to bias towards free (but only if reasonable)
    1627           value *= FREE_BIAS;
    1628           // store square in list
    1629           infeasible_->quickAdd(iSequence, value * value);
    1630         }
    1631         break;
     1603        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1604          // we are going to bias towards free (but only if reasonable)
     1605          value *= FREE_BIAS;
     1606          // store square in list
     1607          infeasible_->quickAdd(iSequence, value * value);
     1608        }
     1609        break;
    16321610      case AbcSimplex::atUpperBound:
    1633         if (value > tolerance) {
    1634           infeasible_->quickAdd(iSequence, value * value);
    1635         }
    1636         break;
     1611        if (value > tolerance) {
     1612          infeasible_->quickAdd(iSequence, value * value);
     1613        }
     1614        break;
    16371615      case AbcSimplex::atLowerBound:
    1638         if (value < -tolerance) {
    1639           infeasible_->quickAdd(iSequence, value * value);
    1640         }
     1616        if (value < -tolerance) {
     1617          infeasible_->quickAdd(iSequence, value * value);
     1618        }
    16411619      }
    16421620    }
     
    16461624      double value = reducedCost[iSequence];
    16471625      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    1648      
    1649       switch(status) {
    1650        
     1626
     1627      switch (status) {
     1628
    16511629      case AbcSimplex::basic:
    16521630      case AbcSimplex::isFixed:
    1653         break;
     1631        break;
    16541632      case AbcSimplex::isFree:
    16551633      case AbcSimplex::superBasic:
    1656         if (fabs(value) > FREE_ACCEPT * tolerance) {
    1657           // we are going to bias towards free (but only if reasonable)
    1658           value *= FREE_BIAS;
    1659           // store square in list
    1660           infeasible_->quickAdd(iSequence, value * value);
    1661         }
    1662         break;
     1634        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1635          // we are going to bias towards free (but only if reasonable)
     1636          value *= FREE_BIAS;
     1637          // store square in list
     1638          infeasible_->quickAdd(iSequence, value * value);
     1639        }
     1640        break;
    16631641      case AbcSimplex::atUpperBound:
    1664         if (value > tolerance) {
    1665           infeasible_->quickAdd(iSequence, value * value);
    1666         }
    1667         break;
     1642        if (value > tolerance) {
     1643          infeasible_->quickAdd(iSequence, value * value);
     1644        }
     1645        break;
    16681646      case AbcSimplex::atLowerBound:
    1669         if (value < -tolerance) {
    1670           infeasible_->quickAdd(iSequence, value * value);
    1671         }
     1647        if (value < -tolerance) {
     1648          infeasible_->quickAdd(iSequence, value * value);
     1649        }
    16721650      }
    16731651    }
    16741652    // Rows
    1675     for (iSequence=0 ; iSequence < numberRows; iSequence++) {
     1653    for (iSequence = 0; iSequence < numberRows; iSequence++) {
    16761654      double value = reducedCost[iSequence];
    16771655      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    1678      
    1679       switch(status) {
    1680        
     1656
     1657      switch (status) {
     1658
    16811659      case AbcSimplex::basic:
    16821660      case AbcSimplex::isFixed:
    1683         break;
     1661        break;
    16841662      case AbcSimplex::isFree:
    16851663      case AbcSimplex::superBasic:
    1686         if (fabs(value) > FREE_ACCEPT * tolerance) {
    1687           // we are going to bias towards free (but only if reasonable)
    1688           value *= FREE_BIAS;
    1689           // store square in list
    1690           infeasible_->quickAdd(iSequence, value * value);
    1691         }
    1692         break;
     1664        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1665          // we are going to bias towards free (but only if reasonable)
     1666          value *= FREE_BIAS;
     1667          // store square in list
     1668          infeasible_->quickAdd(iSequence, value * value);
     1669        }
     1670        break;
    16931671      case AbcSimplex::atUpperBound:
    1694         if (value > tolerance) {
    1695           infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
    1696         }
    1697         break;
     1672        if (value > tolerance) {
     1673          infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
     1674        }
     1675        break;
    16981676      case AbcSimplex::atLowerBound:
    1699         if (value < -tolerance) {
    1700           infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
    1701         }
    1702       }
    1703     }
    1704 #endif
    1705   } 
     1677        if (value < -tolerance) {
     1678          infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
     1679        }
     1680      }
     1681    }
     1682#endif
     1683  }
    17061684}
    17071685// Gets rid of last update
    1708 void
    1709 AbcPrimalColumnSteepest::unrollWeights()
     1686void AbcPrimalColumnSteepest::unrollWeights()
    17101687{
    1711      if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
    1712           return;
    1713      double * saved = alternateWeights_->denseVector();
    1714      int number = alternateWeights_->getNumElements();
    1715      int * which = alternateWeights_->getIndices();
    1716      int i;
    1717      for (i = 0; i < number; i++) {
    1718           int iRow = which[i];
    1719           weights_[iRow] = saved[iRow];
    1720           saved[iRow] = 0.0;
    1721      }
    1722      alternateWeights_->setNumElements(0);
     1688  if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
     1689    return;
     1690  double *saved = alternateWeights_->denseVector();
     1691  int number = alternateWeights_->getNumElements();
     1692  int *which = alternateWeights_->getIndices();
     1693  int i;
     1694  for (i = 0; i < number; i++) {
     1695    int iRow = which[i];
     1696    weights_[iRow] = saved[iRow];
     1697    saved[iRow] = 0.0;
     1698  }
     1699  alternateWeights_->setNumElements(0);
    17231700}
    17241701
     
    17261703// Clone
    17271704//-------------------------------------------------------------------
    1728 AbcPrimalColumnPivot * AbcPrimalColumnSteepest::clone(bool CopyData) const
     1705AbcPrimalColumnPivot *AbcPrimalColumnSteepest::clone(bool CopyData) const
    17291706{
    1730      if (CopyData) {
    1731           return new AbcPrimalColumnSteepest(*this);
    1732      } else {
    1733           return new AbcPrimalColumnSteepest();
    1734      }
     1707  if (CopyData) {
     1708    return new AbcPrimalColumnSteepest(*this);
     1709  } else {
     1710    return new AbcPrimalColumnSteepest();
     1711  }
    17351712}
    1736 void
    1737 AbcPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
     1713void AbcPrimalColumnSteepest::updateWeights(CoinIndexedVector *input)
    17381714{
    1739      // Local copy of mode so can decide what to do
    1740      int switchType = mode_;
    1741      if (mode_ == 4 && numberSwitched_)
    1742           switchType = 3;
    1743      else if (mode_ == 4 || mode_ == 5)
    1744           return;
    1745      int number = input->getNumElements();
    1746      int * which = input->getIndices();
    1747      double * work = input->denseVector();
    1748      int newNumber = 0;
    1749      int * newWhich = alternateWeights_->getIndices();
    1750      double * newWork = alternateWeights_->denseVector();
    1751      int i;
    1752      int sequenceIn = model_->sequenceIn();
    1753      int sequenceOut = model_->sequenceOut();
    1754      const int * pivotVariable = model_->pivotVariable();
    1755 
    1756      int pivotRow = model_->pivotRow();
    1757      pivotSequence_ = pivotRow;
    1758 
    1759      devex_ = 0.0;
    1760           if (pivotRow >= 0) {
    1761                if (switchType == 1) {
    1762                     for (i = 0; i < number; i++) {
    1763                          int iRow = which[i];
    1764                          devex_ += work[iRow] * work[iRow];
    1765                          newWork[iRow] = -2.0 * work[iRow];
    1766                     }
    1767                     newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
    1768                     devex_ += ADD_ONE;
    1769                     weights_[sequenceOut] = 1.0 + ADD_ONE;
    1770                     CoinMemcpyN(which, number, newWhich);
    1771                     alternateWeights_->setNumElements(number);
    1772                } else {
    1773                     if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
    1774                          for (i = 0; i < number; i++) {
    1775                               int iRow = which[i];
    1776                               int iPivot = pivotVariable[iRow];
    1777                               if (reference(iPivot)) {
    1778                                    devex_ += work[iRow] * work[iRow];
    1779                                    newWork[iRow] = -2.0 * work[iRow];
    1780                                    newWhich[newNumber++] = iRow;
    1781                               }
    1782                          }
    1783                          if (!newWork[pivotRow] && devex_ > 0.0)
    1784                               newWhich[newNumber++] = pivotRow; // add if not already in
    1785                          newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
    1786                     } else {
    1787                          for (i = 0; i < number; i++) {
    1788                               int iRow = which[i];
    1789                               int iPivot = pivotVariable[iRow];
    1790                               if (reference(iPivot))
    1791                                    devex_ += work[iRow] * work[iRow];
    1792                          }
    1793                     }
    1794                     if (reference(sequenceIn)) {
    1795                          devex_ += 1.0;
    1796                     } else {
    1797                     }
    1798                     if (reference(sequenceOut)) {
    1799                          weights_[sequenceOut] = 1.0 + 1.0;
    1800                     } else {
    1801                          weights_[sequenceOut] = 1.0;
    1802                     }
    1803                     alternateWeights_->setNumElements(newNumber);
    1804                }
    1805           } else {
    1806                if (switchType == 1) {
    1807                     for (i = 0; i < number; i++) {
    1808                          int iRow = which[i];
    1809                          devex_ += work[iRow] * work[iRow];
    1810                     }
    1811                     devex_ += ADD_ONE;
    1812                } else {
    1813                     for (i = 0; i < number; i++) {
    1814                          int iRow = which[i];
    1815                          int iPivot = pivotVariable[iRow];
    1816                          if (reference(iPivot)) {
    1817                               devex_ += work[iRow] * work[iRow];
    1818                          }
    1819                     }
    1820                     if (reference(sequenceIn))
    1821                          devex_ += 1.0;
    1822                }
     1715  // Local copy of mode so can decide what to do
     1716  int switchType = mode_;
     1717  if (mode_ == 4 && numberSwitched_)
     1718    switchType = 3;
     1719  else if (mode_ == 4 || mode_ == 5)
     1720    return;
     1721  int number = input->getNumElements();
     1722  int *which = input->getIndices();
     1723  double *work = input->denseVector();
     1724  int newNumber = 0;
     1725  int *newWhich = alternateWeights_->getIndices();
     1726  double *newWork = alternateWeights_->denseVector();
     1727  int i;
     1728  int sequenceIn = model_->sequenceIn();
     1729  int sequenceOut = model_->sequenceOut();
     1730  const int *pivotVariable = model_->pivotVariable();
     1731
     1732  int pivotRow = model_->pivotRow();
     1733  pivotSequence_ = pivotRow;
     1734
     1735  devex_ = 0.0;
     1736  if (pivotRow >= 0) {
     1737    if (switchType == 1) {
     1738      for (i = 0; i < number; i++) {
     1739        int iRow = which[i];
     1740        devex_ += work[iRow] * work[iRow];
     1741        newWork[iRow] = -2.0 * work[iRow];
     1742      }
     1743      newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     1744      devex_ += ADD_ONE;
     1745      weights_[sequenceOut] = 1.0 + ADD_ONE;
     1746      CoinMemcpyN(which, number, newWhich);
     1747      alternateWeights_->setNumElements(number);
     1748    } else {
     1749      if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
     1750        for (i = 0; i < number; i++) {
     1751          int iRow = which[i];
     1752          int iPivot = pivotVariable[iRow];
     1753          if (reference(iPivot)) {
     1754            devex_ += work[iRow] * work[iRow];
     1755            newWork[iRow] = -2.0 * work[iRow];
     1756            newWhich[newNumber++] = iRow;
    18231757          }
    1824      double oldDevex = weights_[sequenceIn];
     1758        }
     1759        if (!newWork[pivotRow] && devex_ > 0.0)
     1760          newWhich[newNumber++] = pivotRow; // add if not already in
     1761        newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
     1762      } else {
     1763        for (i = 0; i < number; i++) {
     1764          int iRow = which[i];
     1765          int iPivot = pivotVariable[iRow];
     1766          if (reference(iPivot))
     1767            devex_ += work[iRow] * work[iRow];
     1768        }
     1769      }
     1770      if (reference(sequenceIn)) {
     1771        devex_ += 1.0;
     1772      } else {
     1773      }
     1774      if (reference(sequenceOut)) {
     1775        weights_[sequenceOut] = 1.0 + 1.0;
     1776      } else {
     1777        weights_[sequenceOut] = 1.0;
     1778      }
     1779      alternateWeights_->setNumElements(newNumber);
     1780    }
     1781  } else {
     1782    if (switchType == 1) {
     1783      for (i = 0; i < number; i++) {
     1784        int iRow = which[i];
     1785        devex_ += work[iRow] * work[iRow];
     1786      }
     1787      devex_ += ADD_ONE;
     1788    } else {
     1789      for (i = 0; i < number; i++) {
     1790        int iRow = which[i];
     1791        int iPivot = pivotVariable[iRow];
     1792        if (reference(iPivot)) {
     1793          devex_ += work[iRow] * work[iRow];
     1794        }
     1795      }
     1796      if (reference(sequenceIn))
     1797        devex_ += 1.0;
     1798    }
     1799  }
     1800  double oldDevex = weights_[sequenceIn];
    18251801#ifdef CLP_DEBUG
    1826      if ((model_->messageHandler()->logLevel() & 32))
    1827           printf("old weight %g, new %g\n", oldDevex, devex_);
    1828 #endif
    1829      double check = CoinMax(devex_, oldDevex) + 0.1;
    1830      weights_[sequenceIn] = devex_;
    1831      double testValue = 0.1;
    1832      if (mode_ == 4 && numberSwitched_ == 1)
    1833           testValue = 0.5;
    1834      if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
     1802  if ((model_->messageHandler()->logLevel() & 32))
     1803    printf("old weight %g, new %g\n", oldDevex, devex_);
     1804#endif
     1805  double check = CoinMax(devex_, oldDevex) + 0.1;
     1806  weights_[sequenceIn] = devex_;
     1807  double testValue = 0.1;
     1808  if (mode_ == 4 && numberSwitched_ == 1)
     1809    testValue = 0.5;
     1810  if (fabs(devex_ - oldDevex) > testValue * check) {
    18351811#ifdef CLP_DEBUG
    1836           if ((model_->messageHandler()->logLevel() & 48) == 16)
    1837                printf("old weight %g, new %g\n", oldDevex, devex_);
    1838 #endif
    1839           //printf("old weight %g, new %g\n",oldDevex,devex_);
    1840           testValue = 0.99;
    1841           if (mode_ == 1)
    1842                testValue = 1.01e1; // make unlikely to do if steepest
    1843           else if (mode_ == 4 && numberSwitched_ == 1)
    1844                testValue = 0.9;
    1845           double difference = fabs(devex_ - oldDevex);
    1846           if ( difference > testValue * check ) {
    1847                // need to redo
    1848                model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
    1849                                                  *model_->messagesPointer())
    1850                          << oldDevex << devex_
    1851                          << CoinMessageEol;
    1852                initializeWeights();
    1853           }
    1854      }
    1855      if (pivotRow >= 0) {
    1856           // set outgoing weight here
    1857           weights_[model_->sequenceOut()] = devex_ / (model_->alpha() * model_->alpha());
    1858      }
     1812    if ((model_->messageHandler()->logLevel() & 48) == 16)
     1813      printf("old weight %g, new %g\n", oldDevex, devex_);
     1814#endif
     1815    //printf("old weight %g, new %g\n",oldDevex,devex_);
     1816    testValue = 0.99;
     1817    if (mode_ == 1)
     1818      testValue = 1.01e1; // make unlikely to do if steepest
     1819    else if (mode_ == 4 && numberSwitched_ == 1)
     1820      testValue = 0.9;
     1821    double difference = fabs(devex_ - oldDevex);
     1822    if (difference > testValue * check) {
     1823      // need to redo
     1824      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
     1825        *model_->messagesPointer())
     1826        << oldDevex << devex_
     1827        << CoinMessageEol;
     1828      initializeWeights();
     1829    }
     1830  }
     1831  if (pivotRow >= 0) {
     1832    // set outgoing weight here
     1833    weights_[model_->sequenceOut()] = devex_ / (model_->alpha() * model_->alpha());
     1834  }
    18591835}
    18601836// Checks accuracy - just for debug
    1861 void
    1862 AbcPrimalColumnSteepest::checkAccuracy(int sequence,
    1863                                        double relativeTolerance,
    1864                                        CoinIndexedVector * rowArray1)
     1837void AbcPrimalColumnSteepest::checkAccuracy(int sequence,
     1838  double relativeTolerance,
     1839  CoinIndexedVector *rowArray1)
    18651840{
    1866      if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
    1867           return;
    1868      model_->unpack(*rowArray1, sequence);
    1869      model_->factorization()->updateColumn(*rowArray1);
    1870      int number = rowArray1->getNumElements();
    1871      int * which = rowArray1->getIndices();
    1872      double * work = rowArray1->denseVector();
    1873      const int * pivotVariable = model_->pivotVariable();
    1874 
    1875      double devex = 0.0;
    1876      int i;
    1877 
    1878      if (mode_ == 1) {
    1879           for (i = 0; i < number; i++) {
    1880                int iRow = which[i];
    1881                devex += work[iRow] * work[iRow];
    1882                work[iRow] = 0.0;
    1883           }
    1884           devex += ADD_ONE;
    1885      } else {
    1886           for (i = 0; i < number; i++) {
    1887                int iRow = which[i];
    1888                int iPivot = pivotVariable[iRow];
    1889                if (reference(iPivot)) {
    1890                     devex += work[iRow] * work[iRow];
    1891                }
    1892                work[iRow] = 0.0;
    1893           }
    1894           if (reference(sequence))
    1895                devex += 1.0;
    1896      }
    1897 
    1898      double oldDevex = weights_[sequence];
    1899      double check = CoinMax(devex, oldDevex);;
    1900      if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
    1901        COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
    1902           // update so won't print again
    1903           weights_[sequence] = devex;
    1904      }
    1905      rowArray1->setNumElements(0);
     1841  if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
     1842    return;
     1843  model_->unpack(*rowArray1, sequence);
     1844  model_->factorization()->updateColumn(*rowArray1);
     1845  int number = rowArray1->getNumElements();
     1846  int *which = rowArray1->getIndices();
     1847  double *work = rowArray1->denseVector();
     1848  const int *pivotVariable = model_->pivotVariable();
     1849
     1850  double devex = 0.0;
     1851  int i;
     1852
     1853  if (mode_ == 1) {
     1854    for (i = 0; i < number; i++) {
     1855      int iRow = which[i];
     1856      devex += work[iRow] * work[iRow];
     1857      work[iRow] = 0.0;
     1858    }
     1859    devex += ADD_ONE;
     1860  } else {
     1861    for (i = 0; i < number; i++) {
     1862      int iRow = which[i];
     1863      int iPivot = pivotVariable[iRow];
     1864      if (reference(iPivot)) {
     1865        devex += work[iRow] * work[iRow];
     1866      }
     1867      work[iRow] = 0.0;
     1868    }
     1869    if (reference(sequence))
     1870      devex += 1.0;
     1871  }
     1872
     1873  double oldDevex = weights_[sequence];
     1874  double check = CoinMax(devex, oldDevex);
     1875  ;
     1876  if (fabs(devex - oldDevex) > relativeTolerance * check) {
     1877    COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
     1878    // update so won't print again
     1879    weights_[sequence] = devex;
     1880  }
     1881  rowArray1->setNumElements(0);
    19061882}
    19071883
    19081884// Initialize weights
    1909 void
    1910 AbcPrimalColumnSteepest::initializeWeights()
     1885void AbcPrimalColumnSteepest::initializeWeights()
    19111886{
    1912      int numberRows = model_->numberRows();
    1913      int numberColumns = model_->numberColumns();
    1914      int number = numberRows + numberColumns;
    1915      int iSequence;
    1916      if (mode_ != 1) {
    1917           // initialize to 1.0
    1918           // and set reference framework
    1919           if (!reference_) {
    1920                int nWords = (number + 31) >> 5;
    1921                reference_ = new unsigned int[nWords];
    1922                CoinZeroN(reference_, nWords);
    1923           }
    1924 
    1925           for (iSequence = 0; iSequence < number; iSequence++) {
    1926                weights_[iSequence] = 1.0;
    1927                if (model_->getInternalStatus(iSequence) == AbcSimplex::basic) {
    1928                     setReference(iSequence, false);
    1929                } else {
    1930                     setReference(iSequence, true);
    1931                }
    1932           }
    1933      } else {
    1934           CoinIndexedVector * temp = new CoinIndexedVector();
     1887  int numberRows = model_->numberRows();
     1888  int numberColumns = model_->numberColumns();
     1889  int number = numberRows + numberColumns;
     1890  int iSequence;
     1891  if (mode_ != 1) {
     1892    // initialize to 1.0
     1893    // and set reference framework
     1894    if (!reference_) {
     1895      int nWords = (number + 31) >> 5;
     1896      reference_ = new unsigned int[nWords];
     1897      CoinZeroN(reference_, nWords);
     1898    }
     1899
     1900    for (iSequence = 0; iSequence < number; iSequence++) {
     1901      weights_[iSequence] = 1.0;
     1902      if (model_->getInternalStatus(iSequence) == AbcSimplex::basic) {
     1903        setReference(iSequence, false);
     1904      } else {
     1905        setReference(iSequence, true);
     1906      }
     1907    }
     1908  } else {
     1909    CoinIndexedVector *temp = new CoinIndexedVector();
    19351910#ifndef ABC_USE_COIN_FACTORIZATIONz
    1936           temp->reserve(numberRows +
    1937                         model_->factorization()->maximumPivots());
     1911    temp->reserve(numberRows + model_->factorization()->maximumPivots());
    19381912#else
    19391913#endif
    1940           double * array = alternateWeights_->denseVector();
    1941           int * which = alternateWeights_->getIndices();
    1942 
    1943           for (iSequence = 0; iSequence < number; iSequence++) {
    1944                weights_[iSequence] = 1.0 + ADD_ONE;
    1945                if (model_->getInternalStatus(iSequence) != AbcSimplex::basic &&
    1946                          model_->getInternalStatus(iSequence) != AbcSimplex::isFixed) {
    1947                     model_->unpack(*alternateWeights_, iSequence);
    1948                     double value = ADD_ONE;
    1949                     model_->factorization()->updateColumn(*alternateWeights_);
    1950                     int number = alternateWeights_->getNumElements();
    1951                     int j;
    1952                     for (j = 0; j < number; j++) {
    1953                          int iRow = which[j];
    1954                          value += array[iRow] * array[iRow];
    1955                          array[iRow] = 0.0;
    1956                     }
    1957                     alternateWeights_->setNumElements(0);
    1958                     weights_[iSequence] = value;
    1959                }
    1960           }
    1961           delete temp;
    1962      }
     1914    double *array = alternateWeights_->denseVector();
     1915    int *which = alternateWeights_->getIndices();
     1916
     1917    for (iSequence = 0; iSequence < number; iSequence++) {
     1918      weights_[iSequence] = 1.0 + ADD_ONE;
     1919      if (model_->getInternalStatus(iSequence) != AbcSimplex::basic && model_->getInternalStatus(iSequence) != AbcSimplex::isFixed) {
     1920        model_->unpack(*alternateWeights_, iSequence);
     1921        double value = ADD_ONE;
     1922        model_->factorization()->updateColumn(*alternateWeights_);
     1923        int number = alternateWeights_->getNumElements();
     1924        int j;
     1925        for (j = 0; j < number; j++) {
     1926          int iRow = which[j];
     1927          value += array[iRow] * array[iRow];
     1928          array[iRow] = 0.0;
     1929        }
     1930        alternateWeights_->setNumElements(0);
     1931        weights_[iSequence] = value;
     1932      }
     1933    }
     1934    delete temp;
     1935  }
    19631936}
    19641937// Gets rid of all arrays
    1965 void
    1966 AbcPrimalColumnSteepest::clearArrays()
     1938void AbcPrimalColumnSteepest::clearArrays()
    19671939{
    1968      if (persistence_ == normal) {
    1969           delete [] weights_;
    1970           weights_ = NULL;
    1971           delete infeasible_;
    1972           infeasible_ = NULL;
    1973           delete alternateWeights_;
    1974           alternateWeights_ = NULL;
    1975           delete [] savedWeights_;
    1976           savedWeights_ = NULL;
    1977           delete [] reference_;
    1978           reference_ = NULL;
    1979      }
    1980      pivotSequence_ = -1;
    1981      state_ = -1;
    1982      savedPivotSequence_ = -1;
    1983      savedSequenceOut_ = -1;
    1984      devex_ = 0.0;
     1940  if (persistence_ == normal) {
     1941    delete[] weights_;
     1942    weights_ = NULL;
     1943    delete infeasible_;
     1944    infeasible_ = NULL;
     1945    delete alternateWeights_;
     1946    alternateWeights_ = NULL;
     1947    delete[] savedWeights_;
     1948    savedWeights_ = NULL;
     1949    delete[] reference_;
     1950    reference_ = NULL;
     1951  }
     1952  pivotSequence_ = -1;
     1953  state_ = -1;
     1954  savedPivotSequence_ = -1;
     1955  savedSequenceOut_ = -1;
     1956  devex_ = 0.0;
    19851957}
    19861958// Returns true if would not find any column
    1987 bool
    1988 AbcPrimalColumnSteepest::looksOptimal() const
     1959bool AbcPrimalColumnSteepest::looksOptimal() const
    19891960{
    1990      if (looksOptimal_)
    1991           return true; // user overrode
    1992      //**** THIS MUST MATCH the action coding above
    1993      double tolerance = model_->currentDualTolerance();
    1994      // we can't really trust infeasibilities if there is dual error
    1995      // this coding has to mimic coding in checkDualSolution
    1996      double error = CoinMin(1.0e-2, model_->largestDualError());
    1997      // allow tolerance at least slightly bigger than standard
    1998      tolerance = tolerance + error;
    1999      if(model_->numberIterations() < model_->lastBadIteration() + 200) {
    2000           // we can't really trust infeasibilities if there is dual error
    2001           double checkTolerance = 1.0e-8;
    2002           if (!model_->factorization()->pivots())
    2003                checkTolerance = 1.0e-6;
    2004           if (model_->largestDualError() > checkTolerance)
    2005                tolerance *= model_->largestDualError() / checkTolerance;
    2006           // But cap
    2007           tolerance = CoinMin(1000.0, tolerance);
    2008      }
    2009      int number = model_->numberRows() + model_->numberColumns();
    2010      int iSequence;
    2011 
    2012      double * reducedCost = model_->djRegion();
    2013      int numberInfeasible = 0;
    2014           for (iSequence = 0; iSequence < number; iSequence++) {
    2015                double value = reducedCost[iSequence];
    2016                AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    2017 
    2018                switch(status) {
    2019 
    2020                case AbcSimplex::basic:
    2021                case AbcSimplex::isFixed:
    2022                     break;
    2023                case AbcSimplex::isFree:
    2024                case AbcSimplex::superBasic:
    2025                     if (fabs(value) > FREE_ACCEPT * tolerance)
    2026                          numberInfeasible++;
    2027                     break;
    2028                case AbcSimplex::atUpperBound:
    2029                     if (value > tolerance)
    2030                          numberInfeasible++;
    2031                     break;
    2032                case AbcSimplex::atLowerBound:
    2033                     if (value < -tolerance)
    2034                          numberInfeasible++;
    2035                }
    2036           }
    2037      return numberInfeasible == 0;
     1961  if (looksOptimal_)
     1962    return true; // user overrode
     1963  //**** THIS MUST MATCH the action coding above
     1964  double tolerance = model_->currentDualTolerance();
     1965  // we can't really trust infeasibilities if there is dual error
     1966  // this coding has to mimic coding in checkDualSolution
     1967  double error = CoinMin(1.0e-2, model_->largestDualError());
     1968  // allow tolerance at least slightly bigger than standard
     1969  tolerance = tolerance + error;
     1970  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
     1971    // we can't really trust infeasibilities if there is dual error
     1972    double checkTolerance = 1.0e-8;
     1973    if (!model_->factorization()->pivots())
     1974      checkTolerance = 1.0e-6;
     1975    if (model_->largestDualError() > checkTolerance)
     1976      tolerance *= model_->largestDualError() / checkTolerance;
     1977    // But cap
     1978    tolerance = CoinMin(1000.0, tolerance);
     1979  }
     1980  int number = model_->numberRows() + model_->numberColumns();
     1981  int iSequence;
     1982
     1983  double *reducedCost = model_->djRegion();
     1984  int numberInfeasible = 0;
     1985  for (iSequence = 0; iSequence < number; iSequence++) {
     1986    double value = reducedCost[iSequence];
     1987    AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     1988
     1989    switch (status) {
     1990
     1991    case AbcSimplex::basic:
     1992    case AbcSimplex::isFixed:
     1993      break;
     1994    case AbcSimplex::isFree:
     1995    case AbcSimplex::superBasic:
     1996      if (fabs(value) > FREE_ACCEPT * tolerance)
     1997        numberInfeasible++;
     1998      break;
     1999    case AbcSimplex::atUpperBound:
     2000      if (value > tolerance)
     2001        numberInfeasible++;
     2002      break;
     2003    case AbcSimplex::atLowerBound:
     2004      if (value < -tolerance)
     2005        numberInfeasible++;
     2006    }
     2007  }
     2008  return numberInfeasible == 0;
    20382009}
    20392010// Update djs doing partial pricing (dantzig)
    2040 int
    2041 AbcPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
    2042                                         int numberWanted,
    2043                                         int numberLook)
     2011int AbcPrimalColumnSteepest::partialPricing(CoinIndexedVector *updates,
     2012  int numberWanted,
     2013  int numberLook)
    20442014{
    2045      int number = 0;
    2046      int * index;
    2047      double * updateBy;
    2048      double * reducedCost;
    2049      double saveTolerance = model_->currentDualTolerance();
    2050      double tolerance = model_->currentDualTolerance();
    2051      // we can't really trust infeasibilities if there is dual error
    2052      // this coding has to mimic coding in checkDualSolution
    2053      double error = CoinMin(1.0e-2, model_->largestDualError());
    2054      // allow tolerance at least slightly bigger than standard
    2055      tolerance = tolerance + error;
    2056      if(model_->numberIterations() < model_->lastBadIteration() + 200) {
    2057           // we can't really trust infeasibilities if there is dual error
    2058           double checkTolerance = 1.0e-8;
    2059           if (!model_->factorization()->pivots())
    2060                checkTolerance = 1.0e-6;
    2061           if (model_->largestDualError() > checkTolerance)
    2062                tolerance *= model_->largestDualError() / checkTolerance;
    2063           // But cap
    2064           tolerance = CoinMin(1000.0, tolerance);
    2065      }
    2066      if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    2067        tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    2068      // So partial pricing can use
    2069      model_->setCurrentDualTolerance(tolerance);
    2070      model_->factorization()->updateColumnTranspose(*updates);
    2071      int numberColumns = model_->numberColumns();
    2072 
    2073      // Rows
    2074      reducedCost = model_->djRegion();
    2075 
    2076      number = updates->getNumElements();
    2077      index = updates->getIndices();
    2078      updateBy = updates->denseVector();
    2079      int j;
    2080      double * duals = model_->dualRowSolution();
    2081      for (j = 0; j < number; j++) {
    2082           int iSequence = index[j];
    2083           double value = duals[iSequence];
    2084           value -= updateBy[iSequence];
    2085           updateBy[iSequence] = 0.0;
    2086           duals[iSequence] = value;
    2087      }
    2088      //#define CLP_DEBUG
     2015  int number = 0;
     2016  int *index;
     2017  double *updateBy;
     2018  double *reducedCost;
     2019  double saveTolerance = model_->currentDualTolerance();
     2020  double tolerance = model_->currentDualTolerance();
     2021  // we can't really trust infeasibilities if there is dual error
     2022  // this coding has to mimic coding in checkDualSolution
     2023  double error = CoinMin(1.0e-2, model_->largestDualError());
     2024  // allow tolerance at least slightly bigger than standard
     2025  tolerance = tolerance + error;
     2026  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
     2027    // we can't really trust infeasibilities if there is dual error
     2028    double checkTolerance = 1.0e-8;
     2029    if (!model_->factorization()->pivots())
     2030      checkTolerance = 1.0e-6;
     2031    if (model_->largestDualError() > checkTolerance)
     2032      tolerance *= model_->largestDualError() / checkTolerance;
     2033    // But cap
     2034    tolerance = CoinMin(1000.0, tolerance);
     2035  }
     2036  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
     2037    tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
     2038  // So partial pricing can use
     2039  model_->setCurrentDualTolerance(tolerance);
     2040  model_->factorization()->updateColumnTranspose(*updates);
     2041  int numberColumns = model_->numberColumns();
     2042
     2043  // Rows
     2044  reducedCost = model_->djRegion();
     2045
     2046  number = updates->getNumElements();
     2047  index = updates->getIndices();
     2048  updateBy = updates->denseVector();
     2049  int j;
     2050  double *duals = model_->dualRowSolution();
     2051  for (j = 0; j < number; j++) {
     2052    int iSequence = index[j];
     2053    double value = duals[iSequence];
     2054    value -= updateBy[iSequence];
     2055    updateBy[iSequence] = 0.0;
     2056    duals[iSequence] = value;
     2057  }
     2058  //#define CLP_DEBUG
    20892059#ifdef CLP_DEBUG
    2090      // check duals
    2091      {
    2092           //work space
    2093           CoinIndexedVector arrayVector;
    2094           arrayVector.reserve(numberRows + 1000);
    2095           CoinIndexedVector workSpace;
    2096           workSpace.reserve(numberRows + 1000);
    2097 
    2098 
    2099           int iRow;
    2100           double * array = arrayVector.denseVector();
    2101           int * index = arrayVector.getIndices();
    2102           int number = 0;
    2103           int * pivotVariable = model_->pivotVariable();
    2104           double * cost = model_->costRegion();
    2105           for (iRow = 0; iRow < numberRows; iRow++) {
    2106                int iPivot = pivotVariable[iRow];
    2107                double value = cost[iPivot];
    2108                if (value) {
    2109                     array[iRow] = value;
    2110                     index[number++] = iRow;
    2111                }
     2060  // check duals
     2061  {
     2062    //work space
     2063    CoinIndexedVector arrayVector;
     2064    arrayVector.reserve(numberRows + 1000);
     2065    CoinIndexedVector workSpace;
     2066    workSpace.reserve(numberRows + 1000);
     2067
     2068    int iRow;
     2069    double *array = arrayVector.denseVector();
     2070    int *index = arrayVector.getIndices();
     2071    int number = 0;
     2072    int *pivotVariable = model_->pivotVariable();
     2073    double *cost = model_->costRegion();
     2074    for (iRow = 0; iRow < numberRows; iRow++) {
     2075      int iPivot = pivotVariable[iRow];
     2076      double value = cost[iPivot];
     2077      if (value) {
     2078        array[iRow] = value;
     2079        index[number++] = iRow;
     2080      }
     2081    }
     2082    arrayVector.setNumElements(number);
     2083
     2084    // Btran basic costs
     2085    model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
     2086
     2087    // now look at dual solution
     2088    for (iRow = 0; iRow < numberRows; iRow++) {
     2089      // slack
     2090      double value = array[iRow];
     2091      if (fabs(duals[iRow] - value) > 1.0e-3)
     2092        printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
     2093      //duals[iRow]=value;
     2094    }
     2095  }
     2096#endif
     2097#undef CLP_DEBUG
     2098  double bestDj = tolerance;
     2099  int bestSequence = -1;
     2100
     2101  const double *cost = model_->costRegion();
     2102
     2103  model_->abcMatrix()->setOriginalWanted(numberWanted);
     2104  model_->abcMatrix()->setCurrentWanted(numberWanted);
     2105  int iPassR = 0, iPassC = 0;
     2106  // Setup two passes
     2107  // This biases towards picking row variables
     2108  // This probably should be fixed
     2109  int startR[4];
     2110  const int *which = infeasible_->getIndices();
     2111  int nSlacks = infeasible_->getNumElements();
     2112  startR[1] = nSlacks;
     2113  startR[2] = 0;
     2114  double randomR = model_->randomNumberGenerator()->randomDouble();
     2115  double dstart = static_cast< double >(nSlacks) * randomR;
     2116  startR[0] = static_cast< int >(dstart);
     2117  startR[3] = startR[0];
     2118  double startC[4];
     2119  startC[1] = 1.0;
     2120  startC[2] = 0;
     2121  double randomC = model_->randomNumberGenerator()->randomDouble();
     2122  startC[0] = randomC;
     2123  startC[3] = randomC;
     2124  reducedCost = model_->djRegion();
     2125  int sequenceOut = model_->sequenceOut();
     2126  int chunk = CoinMin(1024, (numberColumns + nSlacks) / 32);
     2127#ifdef COIN_DETAIL
     2128  if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
     2129    printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
     2130    int i;
     2131    for (i = 0; i < 4; i++)
     2132      printf("start R %d C %g ", startR[i], startC[i]);
     2133    printf("\n");
     2134  }
     2135#endif
     2136  chunk = CoinMax(chunk, 256);
     2137  bool finishedR = false, finishedC = false;
     2138  bool doingR = randomR > randomC;
     2139  //doingR=false;
     2140  int saveNumberWanted = numberWanted;
     2141  while (!finishedR || !finishedC) {
     2142    if (finishedR)
     2143      doingR = false;
     2144    if (doingR) {
     2145      int saveSequence = bestSequence;
     2146      int start = startR[iPassR];
     2147      int end = CoinMin(startR[iPassR + 1], start + chunk / 2);
     2148      int jSequence;
     2149      for (jSequence = start; jSequence < end; jSequence++) {
     2150        int iSequence = which[jSequence];
     2151        if (iSequence != sequenceOut) {
     2152          double value;
     2153          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
     2154
     2155          switch (status) {
     2156
     2157          case AbcSimplex::basic:
     2158          case AbcSimplex::isFixed:
     2159            break;
     2160          case AbcSimplex::isFree:
     2161          case AbcSimplex::superBasic:
     2162            value = fabs(cost[iSequence] - duals[iSequence]);
     2163            if (value > FREE_ACCEPT * tolerance) {
     2164              numberWanted--;
     2165              // we are going to bias towards free (but only if reasonable)
     2166              value *= FREE_BIAS;
     2167              if (value > bestDj) {
     2168                // check flagged variable and correct dj
     2169                if (!model_->flagged(iSequence)) {
     2170                  bestDj = value;
     2171                  bestSequence = iSequence;
     2172                } else {
     2173                  // just to make sure we don't exit before got something
     2174                  numberWanted++;
     2175                }
     2176              }
     2177            }
     2178            break;
     2179          case AbcSimplex::atUpperBound:
     2180            value = cost[iSequence] - duals[iSequence];
     2181            if (value > tolerance) {
     2182              numberWanted--;
     2183              if (value > bestDj) {
     2184                // check flagged variable and correct dj
     2185                if (!model_->flagged(iSequence)) {
     2186                  bestDj = value;
     2187                  bestSequence = iSequence;
     2188                } else {
     2189                  // just to make sure we don't exit before got something
     2190                  numberWanted++;
     2191                }
     2192              }
     2193            }
     2194            break;
     2195          case AbcSimplex::atLowerBound:
     2196            value = -(cost[iSequence] - duals[iSequence]);
     2197            if (value > tolerance) {
     2198              numberWanted--;
     2199              if (value > bestDj) {
     2200                // check flagged variable and correct dj
     2201                if (!model_->flagged(iSequence)) {
     2202                  bestDj = value;
     2203                  bestSequence = iSequence;
     2204                } else {
     2205                  // just to make sure we don't exit before got something
     2206                  numberWanted++;
     2207                }
     2208              }
     2209            }
     2210            break;
    21122211          }
    2113           arrayVector.setNumElements(number);
    2114 
    2115           // Btran basic costs
    2116           model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
    2117 
    2118           // now look at dual solution
    2119           for (iRow = 0; iRow < numberRows; iRow++) {
    2120                // slack
    2121                double value = array[iRow];
    2122                if (fabs(duals[iRow] - value) > 1.0e-3)
    2123                     printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
    2124                //duals[iRow]=value;
    2125           }
    2126      }
    2127 #endif
    2128 #undef CLP_DEBUG
    2129      double bestDj = tolerance;
    2130      int bestSequence = -1;
    2131 
    2132      const double * cost = model_->costRegion();
    2133 
    2134      model_->abcMatrix()->setOriginalWanted(numberWanted);
    2135      model_->abcMatrix()->setCurrentWanted(numberWanted);
    2136      int iPassR = 0, iPassC = 0;
    2137      // Setup two passes
    2138      // This biases towards picking row variables
    2139      // This probably should be fixed
    2140      int startR[4];
    2141      const int * which = infeasible_->getIndices();
    2142      int nSlacks = infeasible_->getNumElements();
    2143      startR[1] = nSlacks;
    2144      startR[2] = 0;
    2145      double randomR = model_->randomNumberGenerator()->randomDouble();
    2146      double dstart = static_cast<double> (nSlacks) * randomR;
    2147      startR[0] = static_cast<int> (dstart);
    2148      startR[3] = startR[0];
    2149      double startC[4];
    2150      startC[1] = 1.0;
    2151      startC[2] = 0;
    2152      double randomC = model_->randomNumberGenerator()->randomDouble();
    2153      startC[0] = randomC;
    2154      startC[3] = randomC;
    2155      reducedCost = model_->djRegion();
    2156      int sequenceOut = model_->sequenceOut();
    2157      int chunk = CoinMin(1024, (numberColumns + nSlacks) / 32);
    2158 #ifdef COIN_DETAIL
    2159      if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
    2160           printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
    2161           int i;
    2162           for (i = 0; i < 4; i++)
    2163                printf("start R %d C %g ", startR[i], startC[i]);
    2164           printf("\n");
    2165      }
    2166 #endif
    2167      chunk = CoinMax(chunk, 256);
    2168      bool finishedR = false, finishedC = false;
    2169      bool doingR = randomR > randomC;
    2170      //doingR=false;
    2171      int saveNumberWanted = numberWanted;
    2172      while (!finishedR || !finishedC) {
    2173           if (finishedR)
    2174                doingR = false;
    2175           if (doingR) {
    2176                int saveSequence = bestSequence;
    2177                int start = startR[iPassR];
    2178                int end = CoinMin(startR[iPassR+1], start + chunk / 2);
    2179                int jSequence;
    2180                for (jSequence = start; jSequence < end; jSequence++) {
    2181                     int iSequence = which[jSequence];
    2182                     if (iSequence != sequenceOut) {
    2183                          double value;
    2184                          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
    2185 
    2186                          switch(status) {
    2187 
    2188                          case AbcSimplex::basic:
    2189                          case AbcSimplex::isFixed:
    2190                               break;
    2191                          case AbcSimplex::isFree:
    2192                          case AbcSimplex::superBasic:
    2193                               value = fabs(cost[iSequence] - duals[iSequence]);
    2194                               if (value > FREE_ACCEPT * tolerance) {
    2195                                    numberWanted--;
    2196                                    // we are going to bias towards free (but only if reasonable)
    2197                                    value *= FREE_BIAS;
    2198                                    if (value > bestDj) {
    2199                                         // check flagged variable and correct dj
    2200                                         if (!model_->flagged(iSequence)) {
    2201                                              bestDj = value;
    2202                                              bestSequence = iSequence;
    2203                                         } else {
    2204                                              // just to make sure we don't exit before got something
    2205                                              numberWanted++;
    2206                                         }
    2207                                    }
    2208                               }
    2209                               break;
    2210                          case AbcSimplex::atUpperBound:
    2211                               value = cost[iSequence] - duals[iSequence];
    2212                               if (value > tolerance) {
    2213                                    numberWanted--;
    2214                                    if (value > bestDj) {
    2215                                         // check flagged variable and correct dj
    2216                                         if (!model_->flagged(iSequence)) {
    2217                                              bestDj = value;
    2218                                              bestSequence = iSequence;
    2219                                         } else {
    2220                                              // just to make sure we don't exit before got something
    2221                                              numberWanted++;
    2222                                         }
    2223                                    }
    2224                               }
    2225                               break;
    2226                          case AbcSimplex::atLowerBound:
    2227                               value = -(cost[iSequence] - duals[iSequence]);
    2228                               if (value > tolerance) {
    2229                                    numberWanted--;
    2230                                    if (value > bestDj) {
    2231                                         // check flagged variable and correct dj
    2232                                         if (!model_->flagged(iSequence)) {
    2233                                              bestDj = value;
    2234                                              bestSequence = iSequence;
    2235                                         } else {
    2236                                              // just to make sure we don't exit before got something
    2237                                              numberWanted++;
    2238                                         }
    2239                                    }
    2240                               }
    2241                               break;
    2242                          }
    2243                     }
    2244                     if (!numberWanted)
    2245                          break;
    2246                }
    2247                numberLook -= (end - start);
    2248                if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
    2249                     numberWanted = 0; // give up
    2250                if (saveSequence != bestSequence) {
    2251                     // dj
    2252                     reducedCost[bestSequence] = cost[bestSequence] - duals[bestSequence];
    2253                     bestDj = fabs(reducedCost[bestSequence]);
    2254                     model_->abcMatrix()->setSavedBestSequence(bestSequence);
    2255                     model_->abcMatrix()->setSavedBestDj(reducedCost[bestSequence]);
    2256                }
    2257                model_->abcMatrix()->setCurrentWanted(numberWanted);
    2258                if (!numberWanted)
    2259                     break;
    2260                doingR = false;
    2261                // update start
    2262                startR[iPassR] = jSequence;
    2263                if (jSequence >= startR[iPassR+1]) {
    2264                     if (iPassR)
    2265                          finishedR = true;
    2266                     else
    2267                          iPassR = 2;
    2268                }
    2269           }
    2270           if (finishedC)
    2271                doingR = true;
    2272           if (!doingR) {
    2273             // temp
    2274             int saveSequence = bestSequence;
    2275             // Columns
    2276             double start = startC[iPassC];
    2277             // If we put this idea back then each function needs to update endFraction **
     2212        }
     2213        if (!numberWanted)
     2214          break;
     2215      }
     2216      numberLook -= (end - start);
     2217      if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
     2218        numberWanted = 0; // give up
     2219      if (saveSequence != bestSequence) {
     2220        // dj
     2221        reducedCost[bestSequence] = cost[bestSequence] - duals[bestSequence];
     2222        bestDj = fabs(reducedCost[bestSequence]);
     2223        model_->abcMatrix()->setSavedBestSequence(bestSequence);
     2224        model_->abcMatrix()->setSavedBestDj(reducedCost[bestSequence]);
     2225      }
     2226      model_->abcMatrix()->setCurrentWanted(numberWanted);
     2227      if (!numberWanted)
     2228        break;
     2229      doingR = false;
     2230      // update start
     2231      startR[iPassR] = jSequence;
     2232      if (jSequence >= startR[iPassR + 1]) {
     2233        if (iPassR)
     2234          finishedR = true;
     2235        else
     2236          iPassR = 2;
     2237      }
     2238    }
     2239    if (finishedC)
     2240      doingR = true;
     2241    if (!doingR) {
     2242      // temp
     2243      int saveSequence = bestSequence;
     2244      // Columns
     2245      double start = startC[iPassC];
     2246      // If we put this idea back then each function needs to update endFraction **
    22782247#if 0
    22792248            double dchunk = (static_cast<double> chunk) / (static_cast<double> numberColumns);
    22802249            double end = CoinMin(startC[iPassC+1], start + dchunk);;
    22812250#else
    2282             double end = startC[iPassC+1]; // force end
    2283 #endif
    2284             model_->abcMatrix()->partialPricing(start, end, bestSequence, numberWanted);
    2285             numberWanted = model_->abcMatrix()->currentWanted();
    2286             numberLook -= static_cast<int> ((end - start) * numberColumns);
    2287             if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
    2288               numberWanted = 0; // give up
    2289             if (bestSequence!=saveSequence) {
    2290               // dj
    2291               bestDj = fabs(reducedCost[bestSequence]);
    2292             }
    2293             if (!numberWanted)
    2294               break;
    2295             doingR = true;
    2296             // update start
    2297             startC[iPassC] = end;
    2298             if (end >= startC[iPassC+1] - 1.0e-8) {
    2299               if (iPassC)
    2300                 finishedC = true;
    2301               else
    2302                 iPassC = 2;
    2303             }
    2304           }
    2305      }
    2306      updates->setNumElements(0);
    2307 
    2308      // Restore tolerance
    2309      model_->setCurrentDualTolerance(saveTolerance);
     2251      double end = startC[iPassC + 1]; // force end
     2252#endif
     2253      model_->abcMatrix()->partialPricing(start, end, bestSequence, numberWanted);
     2254      numberWanted = model_->abcMatrix()->currentWanted();
     2255      numberLook -= static_cast< int >((end - start) * numberColumns);
     2256      if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
     2257        numberWanted = 0; // give up
     2258      if (bestSequence != saveSequence) {
     2259        // dj
     2260        bestDj = fabs(reducedCost[bestSequence]);
     2261      }
     2262      if (!numberWanted)
     2263        break;
     2264      doingR = true;
     2265      // update start
     2266      startC[iPassC] = end;
     2267      if (end >= startC[iPassC + 1] - 1.0e-8) {
     2268        if (iPassC)
     2269          finishedC = true;
     2270        else
     2271          iPassC = 2;
     2272      }
     2273    }
     2274  }
     2275  updates->setNumElements(0);
     2276
     2277  // Restore tolerance
     2278  model_->setCurrentDualTolerance(saveTolerance);
    23102279#ifndef NDEBUG
    2311      if (bestSequence >= 0) {
    2312           if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
    2313                assert(model_->reducedCost(bestSequence) < 0.0);
    2314           if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound)
    2315                assert(model_->reducedCost(bestSequence) > 0.0);
    2316      }
    2317 #endif
    2318      return bestSequence;
     2280  if (bestSequence >= 0) {
     2281    if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
     2282      assert(model_->reducedCost(bestSequence) < 0.0);
     2283    if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound)
     2284      assert(model_->reducedCost(bestSequence) > 0.0);
     2285  }
     2286#endif
     2287  return bestSequence;
    23192288}
     2289
     2290/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     2291*/
Note: See TracChangeset for help on using the changeset viewer.