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/ClpPrimalColumnSteepest.cpp

    r2373 r2385  
    2525#define ALT_UPDATE_WEIGHTS 2
    2626#endif
    27 #if ALT_UPDATE_WEIGHTS==1
    28 extern CoinIndexedVector * altVector[3];
    29 #endif
    30 static void debug1(int iSequence,double value,double weight)
     27#if ALT_UPDATE_WEIGHTS == 1
     28extern CoinIndexedVector *altVector[3];
     29#endif
     30static void debug1(int iSequence, double value, double weight)
    3131{
    3232  printf("xx %d inf %.20g wt %.20g\n",
    33          iSequence,value,weight);
     33    iSequence, value, weight);
    3434}
    3535//#define CLP_DEBUG
     
    4141// Default Constructor
    4242//-------------------------------------------------------------------
    43 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (int mode)
    44      : ClpPrimalColumnPivot(),
    45        devex_(0.0),
    46        weights_(NULL),
    47        infeasible_(NULL),
    48        alternateWeights_(NULL),
    49        savedWeights_(NULL),
    50        reference_(NULL),
    51        state_(-1),
    52        mode_(mode),
    53        infeasibilitiesState_(0),
    54        persistence_(normal),
    55        numberSwitched_(0),
    56        pivotSequence_(-1),
    57        savedPivotSequence_(-1),
    58        savedSequenceOut_(-1),
    59       sizeFactorization_(0)
     43ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(int mode)
     44  : ClpPrimalColumnPivot()
     45  , devex_(0.0)
     46  , weights_(NULL)
     47  , infeasible_(NULL)
     48  , alternateWeights_(NULL)
     49  , savedWeights_(NULL)
     50  , reference_(NULL)
     51  , state_(-1)
     52  , mode_(mode)
     53  , infeasibilitiesState_(0)
     54  , persistence_(normal)
     55  , numberSwitched_(0)
     56  , pivotSequence_(-1)
     57  , savedPivotSequence_(-1)
     58  , savedSequenceOut_(-1)
     59  , sizeFactorization_(0)
    6060{
    61      type_ = 2 + 64 * mode;
     61  type_ = 2 + 64 * mode;
    6262}
    6363//-------------------------------------------------------------------
    6464// Copy constructor
    6565//-------------------------------------------------------------------
    66 ClpPrimalColumnSteepest::ClpPrimalColumnSteepest (const ClpPrimalColumnSteepest & rhs)
    67      : ClpPrimalColumnPivot(rhs)
     66ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(const ClpPrimalColumnSteepest &rhs)
     67  : ClpPrimalColumnPivot(rhs)
    6868{
    69      state_ = rhs.state_;
    70      mode_ = rhs.mode_;
    71      infeasibilitiesState_ = rhs.infeasibilitiesState_;
    72      persistence_ = rhs.persistence_;
    73      numberSwitched_ = rhs.numberSwitched_;
    74      model_ = rhs.model_;
    75      pivotSequence_ = rhs.pivotSequence_;
    76      savedPivotSequence_ = rhs.savedPivotSequence_;
    77      savedSequenceOut_ = rhs.savedSequenceOut_;
    78      sizeFactorization_ = rhs.sizeFactorization_;
    79      devex_ = rhs.devex_;
    80      if ((model_ && model_->whatsChanged() & 1) != 0) {
    81           if (rhs.infeasible_) {
    82                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    83           } else {
    84                infeasible_ = NULL;
    85           }
    86           reference_ = NULL;
    87           if (rhs.weights_) {
    88                assert(model_);
    89                int number = model_->numberRows() + model_->numberColumns();
    90                assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
    91                weights_ = new double[number];
    92                CoinMemcpyN(rhs.weights_, number, weights_);
    93                savedWeights_ = new double[number];
    94                CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
    95                if (mode_ != 1) {
    96                     reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
    97                }
    98           } else {
    99                weights_ = NULL;
    100                savedWeights_ = NULL;
    101           }
    102           if (rhs.alternateWeights_) {
    103                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    104           } else {
    105                alternateWeights_ = NULL;
    106           }
    107      } else {
    108           infeasible_ = NULL;
    109           reference_ = NULL;
    110           weights_ = NULL;
    111           savedWeights_ = NULL;
    112           alternateWeights_ = NULL;
    113      }
     69  state_ = rhs.state_;
     70  mode_ = rhs.mode_;
     71  infeasibilitiesState_ = rhs.infeasibilitiesState_;
     72  persistence_ = rhs.persistence_;
     73  numberSwitched_ = rhs.numberSwitched_;
     74  model_ = rhs.model_;
     75  pivotSequence_ = rhs.pivotSequence_;
     76  savedPivotSequence_ = rhs.savedPivotSequence_;
     77  savedSequenceOut_ = rhs.savedSequenceOut_;
     78  sizeFactorization_ = rhs.sizeFactorization_;
     79  devex_ = rhs.devex_;
     80  if ((model_ && model_->whatsChanged() & 1) != 0) {
     81    if (rhs.infeasible_) {
     82      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     83    } else {
     84      infeasible_ = NULL;
     85    }
     86    reference_ = NULL;
     87    if (rhs.weights_) {
     88      assert(model_);
     89      int number = model_->numberRows() + model_->numberColumns();
     90      assert(number == rhs.model_->numberRows() + rhs.model_->numberColumns());
     91      weights_ = new double[number];
     92      CoinMemcpyN(rhs.weights_, number, weights_);
     93      savedWeights_ = new double[number];
     94      CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
     95      if (mode_ != 1) {
     96        reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
     97      }
     98    } else {
     99      weights_ = NULL;
     100      savedWeights_ = NULL;
     101    }
     102    if (rhs.alternateWeights_) {
     103      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     104    } else {
     105      alternateWeights_ = NULL;
     106    }
     107  } else {
     108    infeasible_ = NULL;
     109    reference_ = NULL;
     110    weights_ = NULL;
     111    savedWeights_ = NULL;
     112    alternateWeights_ = NULL;
     113  }
    114114}
    115115
     
    117117// Destructor
    118118//-------------------------------------------------------------------
    119 ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest ()
     119ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest()
    120120{
    121      delete [] weights_;
    122      delete infeasible_;
    123      delete alternateWeights_;
    124      delete [] savedWeights_;
    125      delete [] reference_;
     121  delete[] weights_;
     122  delete infeasible_;
     123  delete alternateWeights_;
     124  delete[] savedWeights_;
     125  delete[] reference_;
    126126}
    127127
     
    130130//-------------------------------------------------------------------
    131131ClpPrimalColumnSteepest &
    132 ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest& rhs)
     132ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest &rhs)
    133133{
    134      if (this != &rhs) {
    135           ClpPrimalColumnPivot::operator=(rhs);
    136           state_ = rhs.state_;
    137           mode_ = rhs.mode_;
    138           infeasibilitiesState_ = rhs.infeasibilitiesState_;
    139           persistence_ = rhs.persistence_;
    140           numberSwitched_ = rhs.numberSwitched_;
    141           model_ = rhs.model_;
    142           pivotSequence_ = rhs.pivotSequence_;
    143           savedPivotSequence_ = rhs.savedPivotSequence_;
    144           savedSequenceOut_ = rhs.savedSequenceOut_;
    145           sizeFactorization_ = rhs.sizeFactorization_;
    146           devex_ = rhs.devex_;
    147           delete [] weights_;
    148           delete [] reference_;
    149           reference_ = NULL;
    150           delete infeasible_;
    151           delete alternateWeights_;
    152           delete [] savedWeights_;
    153           savedWeights_ = NULL;
    154           if (rhs.infeasible_ != NULL) {
    155                infeasible_ = new CoinIndexedVector(rhs.infeasible_);
    156           } else {
    157                infeasible_ = NULL;
    158           }
    159           if (rhs.weights_ != NULL) {
    160                assert(model_);
    161                int number = model_->numberRows() + model_->numberColumns();
    162                assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
    163                weights_ = new double[number];
    164                CoinMemcpyN(rhs.weights_, number, weights_);
    165                savedWeights_ = new double[number];
    166                CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
    167                if (mode_ != 1) {
    168                     reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
    169                }
    170           } else {
    171                weights_ = NULL;
    172           }
    173           if (rhs.alternateWeights_ != NULL) {
    174                alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
    175           } else {
    176                alternateWeights_ = NULL;
    177           }
    178      }
    179      return *this;
     134  if (this != &rhs) {
     135    ClpPrimalColumnPivot::operator=(rhs);
     136    state_ = rhs.state_;
     137    mode_ = rhs.mode_;
     138    infeasibilitiesState_ = rhs.infeasibilitiesState_;
     139    persistence_ = rhs.persistence_;
     140    numberSwitched_ = rhs.numberSwitched_;
     141    model_ = rhs.model_;
     142    pivotSequence_ = rhs.pivotSequence_;
     143    savedPivotSequence_ = rhs.savedPivotSequence_;
     144    savedSequenceOut_ = rhs.savedSequenceOut_;
     145    sizeFactorization_ = rhs.sizeFactorization_;
     146    devex_ = rhs.devex_;
     147    delete[] weights_;
     148    delete[] reference_;
     149    reference_ = NULL;
     150    delete infeasible_;
     151    delete alternateWeights_;
     152    delete[] savedWeights_;
     153    savedWeights_ = NULL;
     154    if (rhs.infeasible_ != NULL) {
     155      infeasible_ = new CoinIndexedVector(rhs.infeasible_);
     156    } else {
     157      infeasible_ = NULL;
     158    }
     159    if (rhs.weights_ != NULL) {
     160      assert(model_);
     161      int number = model_->numberRows() + model_->numberColumns();
     162      assert(number == rhs.model_->numberRows() + rhs.model_->numberColumns());
     163      weights_ = new double[number];
     164      CoinMemcpyN(rhs.weights_, number, weights_);
     165      savedWeights_ = new double[number];
     166      CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
     167      if (mode_ != 1) {
     168        reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
     169      }
     170    } else {
     171      weights_ = NULL;
     172    }
     173    if (rhs.alternateWeights_ != NULL) {
     174      alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
     175    } else {
     176      alternateWeights_ = NULL;
     177    }
     178  }
     179  return *this;
    180180}
    181181// These have to match ClpPackedMatrix version
     
    183183#define ADD_ONE 1.0
    184184static void
    185 pivotColumnBit(clpTempInfo & info)
     185pivotColumnBit(clpTempInfo &info)
    186186{
    187   double * COIN_RESTRICT weights = const_cast<double *>(info.lower);
    188   const unsigned char * COIN_RESTRICT status = info.status;
    189   double tolerance=info.tolerance;
    190   double bestDj=info.primalRatio;
    191   int bestSequence=-1;
    192   double * COIN_RESTRICT infeas = const_cast<double *>(info.infeas);
    193   const int * COIN_RESTRICT start = info.which;
    194   const int * COIN_RESTRICT index = info.index;
     187  double *COIN_RESTRICT weights = const_cast< double * >(info.lower);
     188  const unsigned char *COIN_RESTRICT status = info.status;
     189  double tolerance = info.tolerance;
     190  double bestDj = info.primalRatio;
     191  int bestSequence = -1;
     192  double *COIN_RESTRICT infeas = const_cast< double * >(info.infeas);
     193  const int *COIN_RESTRICT start = info.which;
     194  const int *COIN_RESTRICT index = info.index;
    195195  for (int i = start[0]; i < start[1]; i++) {
    196196    int iSequence = index[i];
     
    199199    if (value > tolerance) {
    200200      if (value > bestDj * weight) {
    201         // check flagged variable and correct dj
    202         if ((status[iSequence]&64)==0) {
    203           bestDj = value / weight;
    204           bestSequence = iSequence;
    205         }
    206       }
    207     }
    208   }
    209   info.primalRatio=bestDj;
    210   info.numberAdded=bestSequence;
     201        // check flagged variable and correct dj
     202        if ((status[iSequence] & 64) == 0) {
     203          bestDj = value / weight;
     204          bestSequence = iSequence;
     205        }
     206      }
     207    }
     208  }
     209  info.primalRatio = bestDj;
     210  info.numberAdded = bestSequence;
    211211}
    212212// Returns pivot column, -1 if none
     
    214214        that is just +-weight where a feasibility changed.  It also has
    215215        reduced cost from last iteration in pivot row*/
    216 int
    217 ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector * updates,
    218                                      CoinIndexedVector * spareRow1,
    219                                      CoinIndexedVector * spareRow2,
    220                                      CoinIndexedVector * spareColumn1,
    221                                      CoinIndexedVector * spareColumn2)
     216int ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates,
     217  CoinIndexedVector *spareRow1,
     218  CoinIndexedVector *spareRow2,
     219  CoinIndexedVector *spareColumn1,
     220  CoinIndexedVector *spareColumn2)
    222221{
    223      assert(model_);
    224      if (model_->nonLinearCost()->lookBothWays() || model_->algorithm() == 2) {
    225           // Do old way
    226           updates->expand();
    227           return pivotColumnOldMethod(updates, spareRow1, spareRow2,
    228                                       spareColumn1, spareColumn2);
    229      }
    230      int number = 0;
    231      int * index;
    232      double tolerance = model_->currentDualTolerance();
    233      // we can't really trust infeasibilities if there is dual error
    234      // this coding has to mimic coding in checkDualSolution
    235      double error = CoinMin(1.0e-2, model_->largestDualError());
    236      // allow tolerance at least slightly bigger than standard
    237      tolerance = tolerance + error;
    238      int pivotRow = model_->pivotRow();
    239      int anyUpdates;
    240      double * infeas = infeasible_->denseVector();
    241 
    242      // Local copy of mode so can decide what to do
    243      int switchType;
    244      if (mode_ == 4)
    245           switchType = 5 - numberSwitched_;
    246      else if (mode_ >= 10)
    247           switchType = 3;
    248      else
    249           switchType = mode_;
    250      /* switchType -
     222  assert(model_);
     223  if (model_->nonLinearCost()->lookBothWays() || model_->algorithm() == 2) {
     224    // Do old way
     225    updates->expand();
     226    return pivotColumnOldMethod(updates, spareRow1, spareRow2,
     227      spareColumn1, spareColumn2);
     228  }
     229  int number = 0;
     230  int *index;
     231  double tolerance = model_->currentDualTolerance();
     232  // we can't really trust infeasibilities if there is dual error
     233  // this coding has to mimic coding in checkDualSolution
     234  double error = CoinMin(1.0e-2, model_->largestDualError());
     235  // allow tolerance at least slightly bigger than standard
     236  tolerance = tolerance + error;
     237  int pivotRow = model_->pivotRow();
     238  int anyUpdates;
     239  double *infeas = infeasible_->denseVector();
     240
     241  // Local copy of mode so can decide what to do
     242  int switchType;
     243  if (mode_ == 4)
     244    switchType = 5 - numberSwitched_;
     245  else if (mode_ >= 10)
     246    switchType = 3;
     247  else
     248    switchType = mode_;
     249    /* switchType -
    251250        0 - all exact devex
    252251        1 - all steepest
     
    257256        10 - can go to mini-sprint
    258257     */
    259      // Look at gub
     258    // Look at gub
    260259#if 1
    261      model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
     260  model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
    262261#else
    263      updates->clear();
    264      model_->computeDuals(NULL);
    265 #endif
    266      if (updates->getNumElements() > 1) {
    267           // would have to have two goes for devex, three for steepest
    268           anyUpdates = 2;
    269      } else if (updates->getNumElements()) {
    270           if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
    271                // reasonable size
    272                anyUpdates = 1;
    273                //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
    274                //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
    275           } else {
    276                // too small
    277                anyUpdates = 2;
     262  updates->clear();
     263  model_->computeDuals(NULL);
     264#endif
     265  if (updates->getNumElements() > 1) {
     266    // would have to have two goes for devex, three for steepest
     267    anyUpdates = 2;
     268  } else if (updates->getNumElements()) {
     269    if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
     270      // reasonable size
     271      anyUpdates = 1;
     272      //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
     273      //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
     274    } else {
     275      // too small
     276      anyUpdates = 2;
     277    }
     278  } else if (pivotSequence_ >= 0) {
     279    // just after re-factorization
     280    anyUpdates = -1;
     281  } else {
     282    // sub flip - nothing to do
     283    anyUpdates = 0;
     284  }
     285  int sequenceOut = model_->sequenceOut();
     286  if (switchType == 5) {
     287    // If known matrix then we will do partial pricing
     288    if (model_->clpMatrix()->canDoPartialPricing()) {
     289      pivotSequence_ = -1;
     290      pivotRow = -1;
     291      // See if to switch
     292      int numberRows = model_->numberRows();
     293      int numberWanted = 10;
     294      int numberColumns = model_->numberColumns();
     295      int numberHiddenRows = model_->clpMatrix()->hiddenRows();
     296      double ratio = static_cast< double >(sizeFactorization_ + numberHiddenRows) / static_cast< double >(numberRows + 2 * numberHiddenRows);
     297      // Number of dual infeasibilities at last invert
     298      int numberDual = model_->numberDualInfeasibilities();
     299      int numberLook = CoinMin(numberDual, numberColumns / 10);
     300      if (ratio < 1.0) {
     301        numberWanted = 100;
     302        numberLook /= 20;
     303        numberWanted = CoinMax(numberWanted, numberLook);
     304      } else if (ratio < 3.0) {
     305        numberWanted = 500;
     306        numberLook /= 15;
     307        numberWanted = CoinMax(numberWanted, numberLook);
     308      } else if (ratio < 4.0 || mode_ == 5) {
     309        numberWanted = 1000;
     310        numberLook /= 10;
     311        numberWanted = CoinMax(numberWanted, numberLook);
     312      } else if (mode_ != 5) {
     313        switchType = 4;
     314        // initialize
     315        numberSwitched_++;
     316        // Make sure will re-do
     317        delete[] weights_;
     318        weights_ = NULL;
     319        model_->computeDuals(NULL);
     320        saveWeights(model_, 4);
     321        anyUpdates = 0;
     322        COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
     323      }
     324      if (switchType == 5) {
     325        numberLook *= 5; // needs tuning for gub
     326        if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
     327          COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
     328            sizeFactorization_, ratio, numberWanted, numberLook));
     329        }
     330        // Update duals and row djs
     331        // Do partial pricing
     332        return partialPricing(updates, spareRow2,
     333          numberWanted, numberLook);
     334      }
     335    }
     336  }
     337  int bestSequence = -1;
     338  model_->spareIntArray_[3] = -3;
     339  if (switchType == 5) {
     340    if (anyUpdates > 0) {
     341      justDjs(updates, spareRow2,
     342        spareColumn1, spareColumn2);
     343    }
     344  } else if (anyUpdates == 1) {
     345    if (switchType < 4) {
     346      // exact etc when can use dj
     347      djsAndSteepest(updates, spareRow2,
     348        spareColumn1, spareColumn2);
     349      if (model_->spareIntArray_[3] > -2) {
     350        bestSequence = model_->spareIntArray_[3];
     351        infeasibilitiesState_ = 2;
     352      } else if (model_->spareIntArray_[3] == -2) {
     353        redoInfeasibilities();
     354      }
     355    } else {
     356      // devex etc when can use dj
     357      djsAndDevex(updates, spareRow2,
     358        spareColumn1, spareColumn2);
     359    }
     360  } else if (anyUpdates == -1) {
     361    if (switchType < 4) {
     362      // exact etc when djs okay
     363      justSteepest(updates, spareRow2,
     364        spareColumn1, spareColumn2);
     365    } else {
     366      // devex etc when djs okay
     367      justDevex(updates, spareRow2,
     368        spareColumn1, spareColumn2);
     369    }
     370  } else if (anyUpdates == 2) {
     371    if (switchType < 4) {
     372      // exact etc when have to use pivot
     373      djsAndSteepest2(updates, spareRow2,
     374        spareColumn1, spareColumn2);
     375    } else {
     376      // devex etc when have to use pivot
     377      djsAndDevex2(updates, spareRow2,
     378        spareColumn1, spareColumn2);
     379    }
     380  }
     381  // everything may have been done if vector copy
     382  if (infeasibilitiesState_ == 2) {
     383    // all done
     384    infeasibilitiesState_ = 1;
     385    model_->clpMatrix()->setSavedBestSequence(bestSequence);
     386    if (bestSequence >= 0)
     387      model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
     388    assert(sequenceOut != bestSequence);
     389    return bestSequence;
     390  } else if (infeasibilitiesState_ == 1) {
     391    // need to redo
     392    //infeasibilitiesState_ = 0;
     393    redoInfeasibilities();
     394  }
     395
     396#ifdef CLP_DEBUG
     397  alternateWeights_->checkClear();
     398#endif
     399  // make sure outgoing from last iteration okay
     400  if (sequenceOut >= 0) {
     401    ClpSimplex::Status status = model_->getStatus(sequenceOut);
     402    double value = model_->reducedCost(sequenceOut);
     403
     404    switch (status) {
     405
     406    case ClpSimplex::basic:
     407    case ClpSimplex::isFixed:
     408      break;
     409    case ClpSimplex::isFree:
     410    case ClpSimplex::superBasic:
     411      if (fabs(value) > FREE_ACCEPT * tolerance) {
     412        // we are going to bias towards free (but only if reasonable)
     413        value *= FREE_BIAS;
     414        // store square in list
     415        if (infeas[sequenceOut])
     416          infeas[sequenceOut] = value * value; // already there
     417        else
     418          infeasible_->quickAdd(sequenceOut, value * value);
     419      } else {
     420        infeasible_->zero(sequenceOut);
     421      }
     422      break;
     423    case ClpSimplex::atUpperBound:
     424      if (value > tolerance) {
     425        // store square in list
     426        if (infeas[sequenceOut])
     427          infeas[sequenceOut] = value * value; // already there
     428        else
     429          infeasible_->quickAdd(sequenceOut, value * value);
     430      } else {
     431        infeasible_->zero(sequenceOut);
     432      }
     433      break;
     434    case ClpSimplex::atLowerBound:
     435      if (value < -tolerance) {
     436        // store square in list
     437        if (infeas[sequenceOut])
     438          infeas[sequenceOut] = value * value; // already there
     439        else
     440          infeasible_->quickAdd(sequenceOut, value * value);
     441      } else {
     442        infeasible_->zero(sequenceOut);
     443      }
     444    }
     445  }
     446  // update of duals finished - now do pricing
     447  // See what sort of pricing
     448  int numberWanted = 10;
     449  number = infeasible_->getNumElements();
     450  int numberColumns = model_->numberColumns();
     451  if (switchType == 5) {
     452    pivotSequence_ = -1;
     453    pivotRow = -1;
     454    // See if to switch
     455    int numberRows = model_->numberRows();
     456    // ratio is done on number of columns here
     457    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
     458    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     459    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
     460    if (ratio < 1.0) {
     461      numberWanted = CoinMax(100, number / 200);
     462    } else if (ratio < 2.0 - 1.0) {
     463      numberWanted = CoinMax(500, number / 40);
     464    } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
     465      numberWanted = CoinMax(2000, number / 10);
     466      numberWanted = CoinMax(numberWanted, numberColumns / 30);
     467    } else if (mode_ != 5) {
     468      switchType = 4;
     469      // initialize
     470      numberSwitched_++;
     471      // Make sure will re-do
     472      delete[] weights_;
     473      weights_ = NULL;
     474      saveWeights(model_, 4);
     475      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
     476    }
     477    //if (model_->numberIterations()%1000==0)
     478    //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
     479  }
     480  int numberRows = model_->numberRows();
     481  // ratio is done on number of rows here
     482  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     483  if (switchType == 4) {
     484    // Still in devex mode
     485    // Go to steepest if lot of iterations?
     486    if (ratio < 5.0) {
     487      numberWanted = CoinMax(2000, number / 10);
     488      numberWanted = CoinMax(numberWanted, numberColumns / 20);
     489    } else if (ratio < 7.0) {
     490      numberWanted = CoinMax(2000, number / 5);
     491      numberWanted = CoinMax(numberWanted, numberColumns / 10);
     492    } else {
     493      // we can zero out
     494      updates->clear();
     495      spareColumn1->clear();
     496      switchType = 3;
     497      // initialize
     498      pivotSequence_ = -1;
     499      pivotRow = -1;
     500      numberSwitched_++;
     501      // Make sure will re-do
     502      delete[] weights_;
     503      weights_ = NULL;
     504      saveWeights(model_, 4);
     505      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
     506      updates->clear();
     507    }
     508    if (model_->numberIterations() % 1000 == 0)
     509      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
     510  }
     511  if (switchType < 4) {
     512    if (switchType < 2) {
     513      numberWanted = COIN_INT_MAX - 1;
     514    } else if (switchType == 2) {
     515      numberWanted = CoinMax(2000, number / 8);
     516    } else {
     517      if (ratio < 1.0) {
     518        numberWanted = CoinMax(2000, number / 20);
     519      } else if (ratio < 5.0) {
     520        numberWanted = CoinMax(2000, number / 10);
     521        numberWanted = CoinMax(numberWanted, numberColumns / 40);
     522      } else if (ratio < 10.0) {
     523        numberWanted = CoinMax(2000, number / 8);
     524        numberWanted = CoinMax(numberWanted, numberColumns / 20);
     525      } else {
     526        ratio = number * (ratio / 80.0);
     527        if (ratio > number) {
     528          numberWanted = number + 1;
     529        } else {
     530          numberWanted = CoinMax(2000, static_cast< int >(ratio));
     531          numberWanted = CoinMax(numberWanted, numberColumns / 10);
     532        }
     533      }
     534    }
     535    //if (model_->numberIterations()%1000==0)
     536    //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
     537    //switchType);
     538  }
     539
     540  double bestDj = 1.0e-30;
     541  bestSequence = -1;
     542#ifdef CLP_USER_DRIVEN
     543  // could go parallel?
     544  model_->eventHandler()->event(ClpEventHandler::beforeChooseIncoming);
     545#endif
     546
     547  int i, iSequence;
     548  index = infeasible_->getIndices();
     549  number = infeasible_->getNumElements();
     550  // Re-sort infeasible every 100 pivots
     551  // Not a good idea
     552  if (0 && model_->factorization()->pivots() > 0 && (model_->factorization()->pivots() % 100) == 0) {
     553    int nLook = model_->numberRows() + numberColumns;
     554    number = 0;
     555    for (i = 0; i < nLook; i++) {
     556      if (infeas[i]) {
     557        if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
     558          index[number++] = i;
     559        else
     560          infeas[i] = 0.0;
     561      }
     562    }
     563    infeasible_->setNumElements(number);
     564  }
     565  if (model_->numberIterations() < model_->lastBadIteration() + 200 && model_->factorization()->pivots() > 10) {
     566    // we can't really trust infeasibilities if there is dual error
     567    double checkTolerance = 1.0e-8;
     568    if (model_->largestDualError() > checkTolerance)
     569      tolerance *= model_->largestDualError() / checkTolerance;
     570    // But cap
     571    tolerance = CoinMin(1000.0, tolerance);
     572  }
     573#ifdef CLP_DEBUG
     574  if (model_->numberDualInfeasibilities() == 1)
     575    printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
     576      model_->largestDualError(), model_, model_->messageHandler(),
     577      number);
     578#endif
     579  // stop last one coming immediately
     580  double saveOutInfeasibility = 0.0;
     581  if (sequenceOut >= 0) {
     582    saveOutInfeasibility = infeas[sequenceOut];
     583    infeas[sequenceOut] = 0.0;
     584  }
     585  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
     586    tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
     587  tolerance *= tolerance; // as we are using squares
     588
     589  int iPass;
     590  // Setup two passes (unless all)
     591  if (mode_ > 1) {
     592    int start[4];
     593    start[1] = number;
     594    start[2] = 0;
     595    double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
     596    start[0] = static_cast< int >(dstart);
     597    start[3] = start[0];
     598    for (iPass = 0; iPass < 2; iPass++) {
     599      int end = start[2 * iPass + 1];
     600      if (switchType < 5) {
     601        for (i = start[2 * iPass]; i < end; i++) {
     602          iSequence = index[i];
     603          double value = infeas[iSequence];
     604          double weight = weights_[iSequence];
     605          //assert (weight);
     606          if (value > tolerance) {
     607            //weight=1.0;
     608            if (value > bestDj * weight) {
     609              // check flagged variable and correct dj
     610              if (!model_->flagged(iSequence)) {
     611                bestDj = value / weight;
     612                bestSequence = iSequence;
     613              } else {
     614                // just to make sure we don't exit before got something
     615                numberWanted++;
     616              }
     617            }
     618            numberWanted--;
    278619          }
    279      } else if (pivotSequence_ >= 0) {
    280           // just after re-factorization
    281           anyUpdates = -1;
    282      } else {
    283           // sub flip - nothing to do
    284           anyUpdates = 0;
    285      }
    286      int sequenceOut = model_->sequenceOut();
    287      if (switchType == 5) {
    288           // If known matrix then we will do partial pricing
    289           if (model_->clpMatrix()->canDoPartialPricing()) {
    290                pivotSequence_ = -1;
    291                pivotRow = -1;
    292                // See if to switch
    293                int numberRows = model_->numberRows();
    294                int numberWanted = 10;
    295                int numberColumns = model_->numberColumns();
    296                int numberHiddenRows = model_->clpMatrix()->hiddenRows();
    297                double ratio = static_cast<double> (sizeFactorization_ + numberHiddenRows) /
    298                               static_cast<double> (numberRows + 2 * numberHiddenRows);
    299                // Number of dual infeasibilities at last invert
    300                int numberDual = model_->numberDualInfeasibilities();
    301                int numberLook = CoinMin(numberDual, numberColumns / 10);
    302                if (ratio < 1.0) {
    303                     numberWanted = 100;
    304                     numberLook /= 20;
    305                     numberWanted = CoinMax(numberWanted, numberLook);
    306                } else if (ratio < 3.0) {
    307                     numberWanted = 500;
    308                     numberLook /= 15;
    309                     numberWanted = CoinMax(numberWanted, numberLook);
    310                } else if (ratio < 4.0 || mode_ == 5) {
    311                     numberWanted = 1000;
    312                     numberLook /= 10;
    313                     numberWanted = CoinMax(numberWanted, numberLook);
    314                } else if (mode_ != 5) {
    315                     switchType = 4;
    316                     // initialize
    317                     numberSwitched_++;
    318                     // Make sure will re-do
    319                     delete [] weights_;
    320                     weights_ = NULL;
    321                     model_->computeDuals(NULL);
    322                     saveWeights(model_, 4);
    323                     anyUpdates = 0;
    324                     COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
    325                }
    326                if (switchType == 5) {
    327                     numberLook *= 5; // needs tuning for gub
    328                     if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
    329                          COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
    330                                                   sizeFactorization_, ratio, numberWanted, numberLook));
    331                     }
    332                     // Update duals and row djs
    333                     // Do partial pricing
    334                     return partialPricing(updates, spareRow2,
    335                                           numberWanted, numberLook);
    336                }
     620          if (!numberWanted)
     621            break;
     622        }
     623      } else {
     624        // Dantzig
     625        for (i = start[2 * iPass]; i < end; i++) {
     626          iSequence = index[i];
     627          double value = infeas[iSequence];
     628          if (value > tolerance) {
     629            if (value > bestDj) {
     630              // check flagged variable and correct dj
     631              if (!model_->flagged(iSequence)) {
     632                bestDj = value;
     633                bestSequence = iSequence;
     634              } else {
     635                // just to make sure we don't exit before got something
     636                numberWanted++;
     637              }
     638            }
     639            numberWanted--;
    337640          }
    338      }
    339      int bestSequence = -1;
    340      model_->spareIntArray_[3] = -3;
    341      if (switchType == 5) {
    342           if (anyUpdates > 0) {
    343                justDjs(updates, spareRow2,
    344                        spareColumn1, spareColumn2);
    345           }
    346      } else if (anyUpdates == 1) {
    347           if (switchType < 4) {
    348                // exact etc when can use dj
    349                djsAndSteepest(updates, spareRow2,
    350                               spareColumn1, spareColumn2);
    351                if (model_->spareIntArray_[3]>-2) {
    352                  bestSequence = model_->spareIntArray_[3];
    353                  infeasibilitiesState_=2;
    354                } else if (model_->spareIntArray_[3]==-2) {
    355                  redoInfeasibilities();
    356                }
    357           } else {
    358                // devex etc when can use dj
    359                djsAndDevex(updates, spareRow2,
    360                            spareColumn1, spareColumn2);
    361           }
    362      } else if (anyUpdates == -1) {
    363           if (switchType < 4) {
    364                // exact etc when djs okay
    365                justSteepest(updates, spareRow2,
    366                             spareColumn1, spareColumn2);
    367           } else {
    368                // devex etc when djs okay
    369                justDevex(updates, spareRow2,
    370                          spareColumn1, spareColumn2);
    371           }
    372      } else if (anyUpdates == 2) {
    373           if (switchType < 4) {
    374                // exact etc when have to use pivot
    375                djsAndSteepest2(updates, spareRow2,
    376                                spareColumn1, spareColumn2);
    377           } else {
    378                // devex etc when have to use pivot
    379                djsAndDevex2(updates, spareRow2,
    380                             spareColumn1, spareColumn2);
    381           }
    382      }
    383      // everything may have been done if vector copy
    384      if (infeasibilitiesState_ == 2) {
    385        // all done
    386        infeasibilitiesState_ = 1;
    387        model_->clpMatrix()->setSavedBestSequence(bestSequence);
    388        if (bestSequence >= 0)
    389          model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
    390        assert (sequenceOut != bestSequence);
    391        return bestSequence;
    392      } else if (infeasibilitiesState_ == 1) {
    393        // need to redo
    394        //infeasibilitiesState_ = 0;
    395        redoInfeasibilities();
    396      }
    397        
    398 #ifdef CLP_DEBUG
    399      alternateWeights_->checkClear();
    400 #endif
    401      // make sure outgoing from last iteration okay
    402      if (sequenceOut >= 0) {
    403           ClpSimplex::Status status = model_->getStatus(sequenceOut);
    404           double value = model_->reducedCost(sequenceOut);
    405 
    406           switch(status) {
    407 
    408           case ClpSimplex::basic:
    409           case ClpSimplex::isFixed:
    410                break;
    411           case ClpSimplex::isFree:
    412           case ClpSimplex::superBasic:
    413                if (fabs(value) > FREE_ACCEPT * tolerance) {
    414                     // we are going to bias towards free (but only if reasonable)
    415                     value *= FREE_BIAS;
    416                     // store square in list
    417                     if (infeas[sequenceOut])
    418                          infeas[sequenceOut] = value * value; // already there
    419                     else
    420                          infeasible_->quickAdd(sequenceOut, value * value);
    421                } else {
    422                     infeasible_->zero(sequenceOut);
    423                }
    424                break;
    425           case ClpSimplex::atUpperBound:
    426                if (value > tolerance) {
    427                     // store square in list
    428                     if (infeas[sequenceOut])
    429                          infeas[sequenceOut] = value * value; // already there
    430                     else
    431                          infeasible_->quickAdd(sequenceOut, value * value);
    432                } else {
    433                     infeasible_->zero(sequenceOut);
    434                }
    435                break;
    436           case ClpSimplex::atLowerBound:
    437                if (value < -tolerance) {
    438                     // store square in list
    439                     if (infeas[sequenceOut])
    440                          infeas[sequenceOut] = value * value; // already there
    441                     else
    442                          infeasible_->quickAdd(sequenceOut, value * value);
    443                } else {
    444                     infeasible_->zero(sequenceOut);
    445                }
    446           }
    447      }
    448      // update of duals finished - now do pricing
    449      // See what sort of pricing
    450      int numberWanted = 10;
    451      number = infeasible_->getNumElements();
    452      int numberColumns = model_->numberColumns();
    453      if (switchType == 5) {
    454           pivotSequence_ = -1;
    455           pivotRow = -1;
    456           // See if to switch
    457           int numberRows = model_->numberRows();
    458           // ratio is done on number of columns here
    459           //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
    460           double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    461           //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
    462           if (ratio < 1.0) {
    463                numberWanted = CoinMax(100, number / 200);
    464           } else if (ratio < 2.0 - 1.0) {
    465                numberWanted = CoinMax(500, number / 40);
    466           } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
    467                numberWanted = CoinMax(2000, number / 10);
    468                numberWanted = CoinMax(numberWanted, numberColumns / 30);
    469           } else if (mode_ != 5) {
    470                switchType = 4;
    471                // initialize
    472                numberSwitched_++;
    473                // Make sure will re-do
    474                delete [] weights_;
    475                weights_ = NULL;
    476                saveWeights(model_, 4);
    477                COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
    478           }
    479           //if (model_->numberIterations()%1000==0)
    480           //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
    481      }
    482      int numberRows = model_->numberRows();
    483      // ratio is done on number of rows here
    484      double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    485      if(switchType == 4) {
    486           // Still in devex mode
    487           // Go to steepest if lot of iterations?
    488           if (ratio < 5.0) {
    489                numberWanted = CoinMax(2000, number / 10);
    490                numberWanted = CoinMax(numberWanted, numberColumns / 20);
    491           } else if (ratio < 7.0) {
    492                numberWanted = CoinMax(2000, number / 5);
    493                numberWanted = CoinMax(numberWanted, numberColumns / 10);
    494           } else {
    495                // we can zero out
    496                updates->clear();
    497                spareColumn1->clear();
    498                switchType = 3;
    499                // initialize
    500                pivotSequence_ = -1;
    501                pivotRow = -1;
    502                numberSwitched_++;
    503                // Make sure will re-do
    504                delete [] weights_;
    505                weights_ = NULL;
    506                saveWeights(model_, 4);
    507                COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
    508                updates->clear();
    509           }
    510           if (model_->numberIterations() % 1000 == 0)
    511             COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
    512      }
    513      if (switchType < 4) {
    514           if (switchType < 2 ) {
    515                numberWanted = COIN_INT_MAX-1;
    516           } else if (switchType == 2) {
    517                numberWanted = CoinMax(2000, number / 8);
    518           } else {
    519                if (ratio < 1.0) {
    520                     numberWanted = CoinMax(2000, number / 20);
    521                } else if (ratio < 5.0) {
    522                     numberWanted = CoinMax(2000, number / 10);
    523                     numberWanted = CoinMax(numberWanted, numberColumns / 40);
    524                } else if (ratio < 10.0) {
    525                     numberWanted = CoinMax(2000, number / 8);
    526                     numberWanted = CoinMax(numberWanted, numberColumns / 20);
    527                } else {
    528                     ratio = number * (ratio / 80.0);
    529                     if (ratio > number) {
    530                          numberWanted = number + 1;
    531                     } else {
    532                          numberWanted = CoinMax(2000, static_cast<int> (ratio));
    533                          numberWanted = CoinMax(numberWanted, numberColumns / 10);
    534                     }
    535                }
    536           }
    537           //if (model_->numberIterations()%1000==0)
    538           //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
    539           //switchType);
    540      }
    541 
    542 
    543      double bestDj = 1.0e-30;
    544      bestSequence = -1;
    545 #ifdef CLP_USER_DRIVEN
    546      // could go parallel?
    547      model_->eventHandler()->event(ClpEventHandler::beforeChooseIncoming);
    548 #endif
    549 
    550      int i, iSequence;
    551      index = infeasible_->getIndices();
    552      number = infeasible_->getNumElements();
    553      // Re-sort infeasible every 100 pivots
    554      // Not a good idea
    555      if (0 && model_->factorization()->pivots() > 0 &&
    556                (model_->factorization()->pivots() % 100) == 0) {
    557           int nLook = model_->numberRows() + numberColumns;
    558           number = 0;
    559           for (i = 0; i < nLook; i++) {
    560                if (infeas[i]) {
    561                     if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
    562                          index[number++] = i;
    563                     else
    564                          infeas[i] = 0.0;
    565                }
    566           }
    567           infeasible_->setNumElements(number);
    568      }
    569      if(model_->numberIterations() < model_->lastBadIteration() + 200 &&
    570                model_->factorization()->pivots() > 10) {
    571           // we can't really trust infeasibilities if there is dual error
    572           double checkTolerance = 1.0e-8;
    573           if (model_->largestDualError() > checkTolerance)
    574                tolerance *= model_->largestDualError() / checkTolerance;
    575           // But cap
    576           tolerance = CoinMin(1000.0, tolerance);
    577      }
    578 #ifdef CLP_DEBUG
    579      if (model_->numberDualInfeasibilities() == 1)
    580           printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
    581                  model_->largestDualError(), model_, model_->messageHandler(),
    582                  number);
    583 #endif
    584      // stop last one coming immediately
    585      double saveOutInfeasibility = 0.0;
    586      if (sequenceOut >= 0) {
    587           saveOutInfeasibility = infeas[sequenceOut];
    588           infeas[sequenceOut] = 0.0;
    589      }
    590      if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    591           tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    592      tolerance *= tolerance; // as we are using squares
    593 
    594      int iPass;
    595      // Setup two passes (unless all)
    596      if (mode_>1) {
    597        int start[4];
    598        start[1] = number;
    599        start[2] = 0;
    600        double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
    601        start[0] = static_cast<int> (dstart);
    602        start[3] = start[0];
    603        for (iPass = 0; iPass < 2; iPass++) {
    604          int end = start[2*iPass+1];
    605          if (switchType < 5) {
    606            for (i = start[2*iPass]; i < end; i++) {
    607              iSequence = index[i];
    608              double value = infeas[iSequence];
    609              double weight = weights_[iSequence];
    610              //assert (weight);
    611              if (value > tolerance) {
    612                //weight=1.0;
    613                if (value > bestDj * weight) {
    614                  // check flagged variable and correct dj
    615                  if (!model_->flagged(iSequence)) {
    616                    bestDj = value / weight;
    617                    bestSequence = iSequence;
    618                  } else {
    619                    // just to make sure we don't exit before got something
    620                    numberWanted++;
    621                  }
    622                }
    623                numberWanted--;
    624              }
    625              if (!numberWanted)
    626                break;
    627            }
    628          } else {
    629            // Dantzig
    630            for (i = start[2*iPass]; i < end; i++) {
    631              iSequence = index[i];
    632              double value = infeas[iSequence];
    633              if (value > tolerance) {
    634                if (value > bestDj) {
    635                  // check flagged variable and correct dj
    636                  if (!model_->flagged(iSequence)) {
    637                    bestDj = value;
    638                    bestSequence = iSequence;
    639                  } else {
    640                    // just to make sure we don't exit before got something
    641                    numberWanted++;
    642                  }
    643                }
    644                numberWanted--;
    645              }
    646              if (!numberWanted)
    647                break;
    648            }
    649          }
    650          if (!numberWanted)
    651            break;
    652        }
    653      } else {
     641          if (!numberWanted)
     642            break;
     643        }
     644      }
     645      if (!numberWanted)
     646        break;
     647    }
     648  } else {
    654649#if ABOCA_LITE
    655        int numberThreads=CoinMax(abcState(),1);
     650    int numberThreads = CoinMax(abcState(), 1);
    656651#define ABOCA_LITE_MAX ABOCA_LITE
    657652#else
    658        const int numberThreads=1;
     653    const int numberThreads = 1;
    659654#define ABOCA_LITE_MAX 1
    660655#endif
    661        if (0) {
    662         int iSequence;
    663         double value;
    664         double weight;
    665          iSequence=34841;
    666         value = infeas[iSequence];
    667         weight = weights_[iSequence];
    668          debug1(iSequence,value,weight);
    669          iSequence=34845;
    670         value = infeas[iSequence];
    671         weight = weights_[iSequence];
    672          debug1(iSequence,value,weight);
    673        }
    674        clpTempInfo info[ABOCA_LITE_MAX];
    675        int start_lite[2*ABOCA_LITE_MAX];
    676        int chunk0 = (number+numberThreads-1)/numberThreads;
    677        int n0=0;
    678        for (int i=0;i<numberThreads;i++) {
    679          int * startX=start_lite+2*i;
    680         info[i].primalRatio = bestDj;
    681         info[i].lower = weights_;
    682         info[i].infeas = infeas;
    683         info[i].index = index;
    684          info[i].status=const_cast<unsigned char *>(model_->statusArray());
    685          info[i].which=startX;
    686          info[i].tolerance=tolerance;
    687          startX[0]=n0;
    688          startX[1]=CoinMin(n0+chunk0,number);
    689         n0 += chunk0;
    690        }
    691        if (numberThreads==1) {
    692         pivotColumnBit(info[0]);
    693        } else {
    694          for (int i=0;i<numberThreads;i++) {
    695            cilk_spawn pivotColumnBit(info[i]);
    696         }
    697         cilk_sync;
    698        }
    699        for (int i=0;i<numberThreads;i++) {
    700         double bestDjX = info[i].primalRatio;
    701          if (bestDjX>bestDj) {
    702            bestDj=bestDjX;
    703            bestSequence=info[i].numberAdded;
    704         }
    705        }
    706      }
    707      //double largestWeight=0.0;
    708      //double smallestWeight=1.0e100;
    709      model_->clpMatrix()->setSavedBestSequence(bestSequence);
    710      if (bestSequence >= 0)
    711           model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
    712      if (sequenceOut >= 0) {
    713           infeas[sequenceOut] = saveOutInfeasibility;
    714      }
    715      /*if (model_->numberIterations()%100==0)
     656    if (0) {
     657      int iSequence;
     658      double value;
     659      double weight;
     660      iSequence = 34841;
     661      value = infeas[iSequence];
     662      weight = weights_[iSequence];
     663      debug1(iSequence, value, weight);
     664      iSequence = 34845;
     665      value = infeas[iSequence];
     666      weight = weights_[iSequence];
     667      debug1(iSequence, value, weight);
     668    }
     669    clpTempInfo info[ABOCA_LITE_MAX];
     670    int start_lite[2 * ABOCA_LITE_MAX];
     671    int chunk0 = (number + numberThreads - 1) / numberThreads;
     672    int n0 = 0;
     673    for (int i = 0; i < numberThreads; i++) {
     674      int *startX = start_lite + 2 * i;
     675      info[i].primalRatio = bestDj;
     676      info[i].lower = weights_;
     677      info[i].infeas = infeas;
     678      info[i].index = index;
     679      info[i].status = const_cast< unsigned char * >(model_->statusArray());
     680      info[i].which = startX;
     681      info[i].tolerance = tolerance;
     682      startX[0] = n0;
     683      startX[1] = CoinMin(n0 + chunk0, number);
     684      n0 += chunk0;
     685    }
     686    if (numberThreads == 1) {
     687      pivotColumnBit(info[0]);
     688    } else {
     689      for (int i = 0; i < numberThreads; i++) {
     690        cilk_spawn pivotColumnBit(info[i]);
     691      }
     692      cilk_sync;
     693    }
     694    for (int i = 0; i < numberThreads; i++) {
     695      double bestDjX = info[i].primalRatio;
     696      if (bestDjX > bestDj) {
     697        bestDj = bestDjX;
     698        bestSequence = info[i].numberAdded;
     699      }
     700    }
     701  }
     702  //double largestWeight=0.0;
     703  //double smallestWeight=1.0e100;
     704  model_->clpMatrix()->setSavedBestSequence(bestSequence);
     705  if (bestSequence >= 0)
     706    model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
     707  if (sequenceOut >= 0) {
     708    infeas[sequenceOut] = saveOutInfeasibility;
     709  }
     710  /*if (model_->numberIterations()%100==0)
    716711       printf("%d best %g\n",bestSequence,bestDj);*/
    717712
    718713#ifdef CLP_USER_DRIVEN
    719      // swap when working
    720      struct {int bestSequence;double bestDj;} tempInfo;
    721      tempInfo.bestDj=bestDj;tempInfo.bestSequence=bestSequence;
    722      model_->eventHandler()->eventWithInfo(ClpEventHandler::afterChooseIncoming,
    723                                    reinterpret_cast<void *>(&tempInfo));
     714  // swap when working
     715  struct {
     716    int bestSequence;
     717    double bestDj;
     718  } tempInfo;
     719  tempInfo.bestDj = bestDj;
     720  tempInfo.bestSequence = bestSequence;
     721  model_->eventHandler()->eventWithInfo(ClpEventHandler::afterChooseIncoming,
     722    reinterpret_cast< void * >(&tempInfo));
    724723#endif
    725724#ifndef NDEBUG
    726      if (bestSequence >= 0) {
    727           if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
    728                assert(model_->reducedCost(bestSequence) < 0.0);
    729           if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound) {
    730                assert(model_->reducedCost(bestSequence) > 0.0);
    731           }
    732      }
     725  if (bestSequence >= 0) {
     726    if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
     727      assert(model_->reducedCost(bestSequence) < 0.0);
     728    if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound) {
     729      assert(model_->reducedCost(bestSequence) > 0.0);
     730    }
     731  }
    733732#endif
    734733#ifdef ALT_UPDATE_WEIGHTSz
    735      printf("weights");
    736      for (int i=0;i<model_->numberColumns()+model_->numberRows();i++) {
    737        if (model_->getColumnStatus(i)==ClpSimplex::isFixed||
    738            model_->getColumnStatus(i)==ClpSimplex::basic)
    739          weights_[i]=1.0;
    740        printf(" %g",weights_[i]);
    741        if ((i%10)==9)
    742          printf("\n");
    743      }
    744      printf("\n");
     734  printf("weights");
     735  for (int i = 0; i < model_->numberColumns() + model_->numberRows(); i++) {
     736    if (model_->getColumnStatus(i) == ClpSimplex::isFixed || model_->getColumnStatus(i) == ClpSimplex::basic)
     737      weights_[i] = 1.0;
     738    printf(" %g", weights_[i]);
     739    if ((i % 10) == 9)
     740      printf("\n");
     741  }
     742  printf("\n");
    745743#endif
    746744#if 0
     
    750748       printf("column %d weight %g infeas %g\n",i,weights_[i],infeas[i]);
    751749#endif
    752      return bestSequence;
     750  return bestSequence;
    753751}
    754752// Just update djs
    755 void
    756 ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
    757                                  CoinIndexedVector * spareRow2,
    758                                  CoinIndexedVector * spareColumn1,
    759                                  CoinIndexedVector * spareColumn2)
     753void ClpPrimalColumnSteepest::justDjs(CoinIndexedVector *updates,
     754  CoinIndexedVector *spareRow2,
     755  CoinIndexedVector *spareColumn1,
     756  CoinIndexedVector *spareColumn2)
    760757{
    761      int iSection, j;
    762      int number = 0;
    763      int * index;
    764      double * updateBy;
    765      double * reducedCost;
    766      double tolerance = model_->currentDualTolerance();
    767      // we can't really trust infeasibilities if there is dual error
    768      // this coding has to mimic coding in checkDualSolution
    769      double error = CoinMin(1.0e-2, model_->largestDualError());
    770      // allow tolerance at least slightly bigger than standard
    771      tolerance = tolerance +  error;
    772      int pivotRow = model_->pivotRow();
    773      double * infeas = infeasible_->denseVector();
    774      //updates->scanAndPack();
    775      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    776 
    777      // put row of tableau in rowArray and columnArray (packed mode)
    778      model_->clpMatrix()->transposeTimes(model_, -1.0,
    779                                          updates, spareColumn2, spareColumn1);
    780      // normal
    781      for (iSection = 0; iSection < 2; iSection++) {
    782 
    783           reducedCost = model_->djRegion(iSection);
    784           int addSequence;
    785 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    786           double slack_multiplier;
    787 #endif
    788 
    789           if (!iSection) {
    790                number = updates->getNumElements();
    791                index = updates->getIndices();
    792                updateBy = updates->denseVector();
    793                addSequence = model_->numberColumns();
    794 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    795                slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
    796 #endif
    797           } else {
    798                number = spareColumn1->getNumElements();
    799                index = spareColumn1->getIndices();
    800                updateBy = spareColumn1->denseVector();
    801                addSequence = 0;
    802 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    803                slack_multiplier = 1.0;
    804 #endif
    805           }
    806 
    807           for (j = 0; j < number; j++) {
    808                int iSequence = index[j];
    809                double value = reducedCost[iSequence];
    810                value -= updateBy[j];
    811                updateBy[j] = 0.0;
    812                reducedCost[iSequence] = value;
    813                ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    814 
    815                switch(status) {
    816 
    817                case ClpSimplex::basic:
    818                     infeasible_->zero(iSequence + addSequence);
    819                case ClpSimplex::isFixed:
    820                     break;
    821                case ClpSimplex::isFree:
    822                case ClpSimplex::superBasic:
    823                     if (fabs(value) > FREE_ACCEPT * tolerance) {
    824                          // we are going to bias towards free (but only if reasonable)
    825                          value *= FREE_BIAS;
    826                          // store square in list
    827                          if (infeas[iSequence+addSequence])
    828                               infeas[iSequence+addSequence] = value * value; // already there
    829                          else
    830                               infeasible_->quickAdd(iSequence + addSequence, value * value);
    831                     } else {
    832                          infeasible_->zero(iSequence + addSequence);
    833                     }
    834                     break;
    835                case ClpSimplex::atUpperBound:
    836                     iSequence += addSequence;
    837                     if (value > tolerance) {
     758  int iSection, j;
     759  int number = 0;
     760  int *index;
     761  double *updateBy;
     762  double *reducedCost;
     763  double tolerance = model_->currentDualTolerance();
     764  // we can't really trust infeasibilities if there is dual error
     765  // this coding has to mimic coding in checkDualSolution
     766  double error = CoinMin(1.0e-2, model_->largestDualError());
     767  // allow tolerance at least slightly bigger than standard
     768  tolerance = tolerance + error;
     769  int pivotRow = model_->pivotRow();
     770  double *infeas = infeasible_->denseVector();
     771  //updates->scanAndPack();
     772  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     773
     774  // put row of tableau in rowArray and columnArray (packed mode)
     775  model_->clpMatrix()->transposeTimes(model_, -1.0,
     776    updates, spareColumn2, spareColumn1);
     777  // normal
     778  for (iSection = 0; iSection < 2; iSection++) {
     779
     780    reducedCost = model_->djRegion(iSection);
     781    int addSequence;
    838782#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    839                         value *= value*slack_multiplier;
     783    double slack_multiplier;
     784#endif
     785
     786    if (!iSection) {
     787      number = updates->getNumElements();
     788      index = updates->getIndices();
     789      updateBy = updates->denseVector();
     790      addSequence = model_->numberColumns();
     791#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     792      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
     793#endif
     794    } else {
     795      number = spareColumn1->getNumElements();
     796      index = spareColumn1->getIndices();
     797      updateBy = spareColumn1->denseVector();
     798      addSequence = 0;
     799#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     800      slack_multiplier = 1.0;
     801#endif
     802    }
     803
     804    for (j = 0; j < number; j++) {
     805      int iSequence = index[j];
     806      double value = reducedCost[iSequence];
     807      value -= updateBy[j];
     808      updateBy[j] = 0.0;
     809      reducedCost[iSequence] = value;
     810      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     811
     812      switch (status) {
     813
     814      case ClpSimplex::basic:
     815        infeasible_->zero(iSequence + addSequence);
     816      case ClpSimplex::isFixed:
     817        break;
     818      case ClpSimplex::isFree:
     819      case ClpSimplex::superBasic:
     820        if (fabs(value) > FREE_ACCEPT * tolerance) {
     821          // we are going to bias towards free (but only if reasonable)
     822          value *= FREE_BIAS;
     823          // store square in list
     824          if (infeas[iSequence + addSequence])
     825            infeas[iSequence + addSequence] = value * value; // already there
     826          else
     827            infeasible_->quickAdd(iSequence + addSequence, value * value);
     828        } else {
     829          infeasible_->zero(iSequence + addSequence);
     830        }
     831        break;
     832      case ClpSimplex::atUpperBound:
     833        iSequence += addSequence;
     834        if (value > tolerance) {
     835#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     836          value *= value * slack_multiplier;
    840837#else
    841                         value *= value;
    842 #endif
    843                          // store square in list
    844                          if (infeas[iSequence])
    845                               infeas[iSequence] = value; // already there
    846                          else
    847                               infeasible_->quickAdd(iSequence, value);
    848                     } else {
    849                          infeasible_->zero(iSequence);
    850                     }
    851                     break;
    852                case ClpSimplex::atLowerBound:
    853                     iSequence += addSequence;
    854                     if (value < -tolerance) {
     838          value *= value;
     839#endif
     840          // store square in list
     841          if (infeas[iSequence])
     842            infeas[iSequence] = value; // already there
     843          else
     844            infeasible_->quickAdd(iSequence, value);
     845        } else {
     846          infeasible_->zero(iSequence);
     847        }
     848        break;
     849      case ClpSimplex::atLowerBound:
     850        iSequence += addSequence;
     851        if (value < -tolerance) {
    855852#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    856                         value *= value*slack_multiplier;
     853          value *= value * slack_multiplier;
    857854#else
    858                         value *= value;
    859 #endif
    860                          // store square in list
    861                          if (infeas[iSequence])
    862                               infeas[iSequence] = value; // already there
    863                          else
    864                               infeasible_->quickAdd(iSequence, value);
    865                     } else {
    866                          infeasible_->zero(iSequence);
    867                     }
    868                }
    869           }
    870      }
    871      updates->setNumElements(0);
    872      spareColumn1->setNumElements(0);
    873      if (pivotRow >= 0) {
    874           // make sure infeasibility on incoming is 0.0
    875           int sequenceIn = model_->sequenceIn();
    876           infeasible_->zero(sequenceIn);
    877      }
     855          value *= value;
     856#endif
     857          // store square in list
     858          if (infeas[iSequence])
     859            infeas[iSequence] = value; // already there
     860          else
     861            infeasible_->quickAdd(iSequence, value);
     862        } else {
     863          infeasible_->zero(iSequence);
     864        }
     865      }
     866    }
     867  }
     868  updates->setNumElements(0);
     869  spareColumn1->setNumElements(0);
     870  if (pivotRow >= 0) {
     871    // make sure infeasibility on incoming is 0.0
     872    int sequenceIn = model_->sequenceIn();
     873    infeasible_->zero(sequenceIn);
     874  }
    878875}
    879876// Update djs, weights for Devex
    880 void
    881 ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
    882                                      CoinIndexedVector * spareRow2,
    883                                      CoinIndexedVector * spareColumn1,
    884                                      CoinIndexedVector * spareColumn2)
     877void ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector *updates,
     878  CoinIndexedVector *spareRow2,
     879  CoinIndexedVector *spareColumn1,
     880  CoinIndexedVector *spareColumn2)
    885881{
    886      int j;
    887      int number = 0;
    888      int * index;
    889      double * updateBy;
    890      double * reducedCost;
    891      double tolerance = model_->currentDualTolerance();
    892      // we can't really trust infeasibilities if there is dual error
    893      // this coding has to mimic coding in checkDualSolution
    894      double error = CoinMin(1.0e-2, model_->largestDualError());
    895      // allow tolerance at least slightly bigger than standard
    896      tolerance = tolerance +  error;
    897      // for weights update we use pivotSequence
    898      // unset in case sub flip
    899      assert (pivotSequence_ >= 0);
    900      assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
    901      pivotSequence_ = -1;
    902      double * infeas = infeasible_->denseVector();
    903      //updates->scanAndPack();
    904      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    905      // and we can see if reference
    906      //double referenceIn = 0.0;
    907      int sequenceIn = model_->sequenceIn();
    908      //if (mode_ != 1 && reference(sequenceIn))
    909      //   referenceIn = 1.0;
    910      // save outgoing weight round update
    911      double outgoingWeight = 0.0;
    912      int sequenceOut = model_->sequenceOut();
    913      if (sequenceOut >= 0)
    914           outgoingWeight = weights_[sequenceOut];
    915 
    916      double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
    917      // put row of tableau in rowArray and columnArray (packed mode)
    918      model_->clpMatrix()->transposeTimes(model_, -1.0,
    919                                          updates, spareColumn2, spareColumn1);
    920      // update weights
    921      double * weight;
    922      int numberColumns = model_->numberColumns();
    923      // rows
    924      reducedCost = model_->djRegion(0);
    925      int addSequence = model_->numberColumns();;
    926 
    927      number = updates->getNumElements();
    928      index = updates->getIndices();
    929      updateBy = updates->denseVector();
    930      weight = weights_ + numberColumns;
    931      // Devex
    932      for (j = 0; j < number; j++) {
    933           double thisWeight;
    934           double pivot;
    935           double value3;
    936           int iSequence = index[j];
    937           double value = reducedCost[iSequence];
    938           double value2 = updateBy[j];
    939           updateBy[j] = 0.0;
    940           value -= value2;
    941           reducedCost[iSequence] = value;
    942           ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    943 
    944           switch(status) {
    945 
    946           case ClpSimplex::basic:
    947                infeasible_->zero(iSequence + addSequence);
    948           case ClpSimplex::isFixed:
    949                break;
    950           case ClpSimplex::isFree:
    951           case ClpSimplex::superBasic:
    952                thisWeight = weight[iSequence];
    953                // row has -1
    954                pivot = value2 * scaleFactor;
    955                value3 = pivot * pivot * devex_;
    956                if (reference(iSequence + numberColumns))
    957                     value3 += 1.0;
    958                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    959                if (fabs(value) > FREE_ACCEPT * tolerance) {
    960                     // we are going to bias towards free (but only if reasonable)
    961                     value *= FREE_BIAS;
    962                     // store square in list
    963                     if (infeas[iSequence+addSequence])
    964                          infeas[iSequence+addSequence] = value * value; // already there
    965                     else
    966                          infeasible_->quickAdd(iSequence + addSequence, value * value);
    967                } else {
    968                     infeasible_->zero(iSequence + addSequence);
    969                }
    970                break;
    971           case ClpSimplex::atUpperBound:
    972                thisWeight = weight[iSequence];
    973                // row has -1
    974                pivot = value2 * scaleFactor;
    975                value3 = pivot * pivot * devex_;
    976                if (reference(iSequence + numberColumns))
    977                     value3 += 1.0;
    978                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    979                iSequence += addSequence;
    980                if (value > tolerance) {
    981                     // store square in list
     882  int j;
     883  int number = 0;
     884  int *index;
     885  double *updateBy;
     886  double *reducedCost;
     887  double tolerance = model_->currentDualTolerance();
     888  // we can't really trust infeasibilities if there is dual error
     889  // this coding has to mimic coding in checkDualSolution
     890  double error = CoinMin(1.0e-2, model_->largestDualError());
     891  // allow tolerance at least slightly bigger than standard
     892  tolerance = tolerance + error;
     893  // for weights update we use pivotSequence
     894  // unset in case sub flip
     895  assert(pivotSequence_ >= 0);
     896  assert(model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
     897  pivotSequence_ = -1;
     898  double *infeas = infeasible_->denseVector();
     899  //updates->scanAndPack();
     900  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     901  // and we can see if reference
     902  //double referenceIn = 0.0;
     903  int sequenceIn = model_->sequenceIn();
     904  //if (mode_ != 1 && reference(sequenceIn))
     905  //   referenceIn = 1.0;
     906  // save outgoing weight round update
     907  double outgoingWeight = 0.0;
     908  int sequenceOut = model_->sequenceOut();
     909  if (sequenceOut >= 0)
     910    outgoingWeight = weights_[sequenceOut];
     911
     912  double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
     913  // put row of tableau in rowArray and columnArray (packed mode)
     914  model_->clpMatrix()->transposeTimes(model_, -1.0,
     915    updates, spareColumn2, spareColumn1);
     916  // update weights
     917  double *weight;
     918  int numberColumns = model_->numberColumns();
     919  // rows
     920  reducedCost = model_->djRegion(0);
     921  int addSequence = model_->numberColumns();
     922  ;
     923
     924  number = updates->getNumElements();
     925  index = updates->getIndices();
     926  updateBy = updates->denseVector();
     927  weight = weights_ + numberColumns;
     928  // Devex
     929  for (j = 0; j < number; j++) {
     930    double thisWeight;
     931    double pivot;
     932    double value3;
     933    int iSequence = index[j];
     934    double value = reducedCost[iSequence];
     935    double value2 = updateBy[j];
     936    updateBy[j] = 0.0;
     937    value -= value2;
     938    reducedCost[iSequence] = value;
     939    ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     940
     941    switch (status) {
     942
     943    case ClpSimplex::basic:
     944      infeasible_->zero(iSequence + addSequence);
     945    case ClpSimplex::isFixed:
     946      break;
     947    case ClpSimplex::isFree:
     948    case ClpSimplex::superBasic:
     949      thisWeight = weight[iSequence];
     950      // row has -1
     951      pivot = value2 * scaleFactor;
     952      value3 = pivot * pivot * devex_;
     953      if (reference(iSequence + numberColumns))
     954        value3 += 1.0;
     955      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     956      if (fabs(value) > FREE_ACCEPT * tolerance) {
     957        // we are going to bias towards free (but only if reasonable)
     958        value *= FREE_BIAS;
     959        // store square in list
     960        if (infeas[iSequence + addSequence])
     961          infeas[iSequence + addSequence] = value * value; // already there
     962        else
     963          infeasible_->quickAdd(iSequence + addSequence, value * value);
     964      } else {
     965        infeasible_->zero(iSequence + addSequence);
     966      }
     967      break;
     968    case ClpSimplex::atUpperBound:
     969      thisWeight = weight[iSequence];
     970      // row has -1
     971      pivot = value2 * scaleFactor;
     972      value3 = pivot * pivot * devex_;
     973      if (reference(iSequence + numberColumns))
     974        value3 += 1.0;
     975      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     976      iSequence += addSequence;
     977      if (value > tolerance) {
     978        // store square in list
    982979#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    983                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
     980        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
    984981#else
    985                     value *= value;
    986 #endif
    987                     if (infeas[iSequence])
    988                          infeas[iSequence] = value; // already there
    989                     else
    990                          infeasible_->quickAdd(iSequence , value);
    991                } else {
    992                     infeasible_->zero(iSequence);
    993                }
    994                break;
    995           case ClpSimplex::atLowerBound:
    996                thisWeight = weight[iSequence];
    997                // row has -1
    998                pivot = value2 * scaleFactor;
    999                value3 = pivot * pivot * devex_;
    1000                if (reference(iSequence + numberColumns))
    1001                     value3 += 1.0;
    1002                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1003                iSequence += addSequence;
    1004                if (value < -tolerance) {
    1005                     // store square in list
     982        value *= value;
     983#endif
     984        if (infeas[iSequence])
     985          infeas[iSequence] = value; // already there
     986        else
     987          infeasible_->quickAdd(iSequence, value);
     988      } else {
     989        infeasible_->zero(iSequence);
     990      }
     991      break;
     992    case ClpSimplex::atLowerBound:
     993      thisWeight = weight[iSequence];
     994      // row has -1
     995      pivot = value2 * scaleFactor;
     996      value3 = pivot * pivot * devex_;
     997      if (reference(iSequence + numberColumns))
     998        value3 += 1.0;
     999      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1000      iSequence += addSequence;
     1001      if (value < -tolerance) {
     1002        // store square in list
    10061003#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1007                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
     1004        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
    10081005#else
    1009                     value *= value;
    1010 #endif
    1011                     if (infeas[iSequence])
    1012                          infeas[iSequence] = value; // already there
    1013                     else
    1014                          infeasible_->quickAdd(iSequence , value);
    1015                } else {
    1016                     infeasible_->zero(iSequence);
    1017                }
    1018           }
    1019      }
    1020 
    1021      // columns
    1022      weight = weights_;
    1023 
    1024      scaleFactor = -scaleFactor;
    1025      reducedCost = model_->djRegion(1);
    1026      number = spareColumn1->getNumElements();
    1027      index = spareColumn1->getIndices();
    1028      updateBy = spareColumn1->denseVector();
    1029 
    1030      // Devex
    1031 
    1032      for (j = 0; j < number; j++) {
    1033           double thisWeight;
    1034           double pivot;
    1035           double value3;
    1036           int iSequence = index[j];
    1037           double value = reducedCost[iSequence];
    1038           double value2 = updateBy[j];
    1039           value -= value2;
    1040           updateBy[j] = 0.0;
    1041           reducedCost[iSequence] = value;
    1042           ClpSimplex::Status status = model_->getStatus(iSequence);
    1043 
    1044           switch(status) {
    1045 
    1046           case ClpSimplex::basic:
    1047                infeasible_->zero(iSequence);
    1048           case ClpSimplex::isFixed:
    1049                break;
    1050           case ClpSimplex::isFree:
    1051           case ClpSimplex::superBasic:
    1052                thisWeight = weight[iSequence];
    1053                // row has -1
    1054                pivot = value2 * scaleFactor;
    1055                value3 = pivot * pivot * devex_;
    1056                if (reference(iSequence))
    1057                     value3 += 1.0;
    1058                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1059                if (fabs(value) > FREE_ACCEPT * tolerance) {
    1060                     // we are going to bias towards free (but only if reasonable)
    1061                     value *= FREE_BIAS;
    1062                     // store square in list
    1063                     if (infeas[iSequence])
    1064                          infeas[iSequence] = value * value; // already there
    1065                     else
    1066                          infeasible_->quickAdd(iSequence, value * value);
    1067                } else {
    1068                     infeasible_->zero(iSequence);
    1069                }
    1070                break;
    1071           case ClpSimplex::atUpperBound:
    1072                thisWeight = weight[iSequence];
    1073                // row has -1
    1074                pivot = value2 * scaleFactor;
    1075                value3 = pivot * pivot * devex_;
    1076                if (reference(iSequence))
    1077                     value3 += 1.0;
    1078                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1079                if (value > tolerance) {
    1080                     // store square in list
    1081                     if (infeas[iSequence])
    1082                          infeas[iSequence] = value * value; // already there
    1083                     else
    1084                          infeasible_->quickAdd(iSequence, value * value);
    1085                } else {
    1086                     infeasible_->zero(iSequence);
    1087                }
    1088                break;
    1089           case ClpSimplex::atLowerBound:
    1090                thisWeight = weight[iSequence];
    1091                // row has -1
    1092                pivot = value2 * scaleFactor;
    1093                value3 = pivot * pivot * devex_;
    1094                if (reference(iSequence))
    1095                     value3 += 1.0;
    1096                weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
    1097                if (value < -tolerance) {
    1098                     // store square in list
    1099                     if (infeas[iSequence])
    1100                          infeas[iSequence] = value * value; // already there
    1101                     else
    1102                          infeasible_->quickAdd(iSequence, value * value);
    1103                } else {
    1104                     infeasible_->zero(iSequence);
    1105                }
    1106           }
    1107      }
    1108      // restore outgoing weight
    1109      if (sequenceOut >= 0)
    1110           weights_[sequenceOut] = outgoingWeight;
    1111      // make sure infeasibility on incoming is 0.0
    1112      infeasible_->zero(sequenceIn);
    1113      spareRow2->setNumElements(0);
    1114      //#define SOME_DEBUG_1
     1006        value *= value;
     1007#endif
     1008        if (infeas[iSequence])
     1009          infeas[iSequence] = value; // already there
     1010        else
     1011          infeasible_->quickAdd(iSequence, value);
     1012      } else {
     1013        infeasible_->zero(iSequence);
     1014      }
     1015    }
     1016  }
     1017
     1018  // columns
     1019  weight = weights_;
     1020
     1021  scaleFactor = -scaleFactor;
     1022  reducedCost = model_->djRegion(1);
     1023  number = spareColumn1->getNumElements();
     1024  index = spareColumn1->getIndices();
     1025  updateBy = spareColumn1->denseVector();
     1026
     1027  // Devex
     1028
     1029  for (j = 0; j < number; j++) {
     1030    double thisWeight;
     1031    double pivot;
     1032    double value3;
     1033    int iSequence = index[j];
     1034    double value = reducedCost[iSequence];
     1035    double value2 = updateBy[j];
     1036    value -= value2;
     1037    updateBy[j] = 0.0;
     1038    reducedCost[iSequence] = value;
     1039    ClpSimplex::Status status = model_->getStatus(iSequence);
     1040
     1041    switch (status) {
     1042
     1043    case ClpSimplex::basic:
     1044      infeasible_->zero(iSequence);
     1045    case ClpSimplex::isFixed:
     1046      break;
     1047    case ClpSimplex::isFree:
     1048    case ClpSimplex::superBasic:
     1049      thisWeight = weight[iSequence];
     1050      // row has -1
     1051      pivot = value2 * scaleFactor;
     1052      value3 = pivot * pivot * devex_;
     1053      if (reference(iSequence))
     1054        value3 += 1.0;
     1055      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1056      if (fabs(value) > FREE_ACCEPT * tolerance) {
     1057        // we are going to bias towards free (but only if reasonable)
     1058        value *= FREE_BIAS;
     1059        // store square in list
     1060        if (infeas[iSequence])
     1061          infeas[iSequence] = value * value; // already there
     1062        else
     1063          infeasible_->quickAdd(iSequence, value * value);
     1064      } else {
     1065        infeasible_->zero(iSequence);
     1066      }
     1067      break;
     1068    case ClpSimplex::atUpperBound:
     1069      thisWeight = weight[iSequence];
     1070      // row has -1
     1071      pivot = value2 * scaleFactor;
     1072      value3 = pivot * pivot * devex_;
     1073      if (reference(iSequence))
     1074        value3 += 1.0;
     1075      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1076      if (value > tolerance) {
     1077        // store square in list
     1078        if (infeas[iSequence])
     1079          infeas[iSequence] = value * value; // already there
     1080        else
     1081          infeasible_->quickAdd(iSequence, value * value);
     1082      } else {
     1083        infeasible_->zero(iSequence);
     1084      }
     1085      break;
     1086    case ClpSimplex::atLowerBound:
     1087      thisWeight = weight[iSequence];
     1088      // row has -1
     1089      pivot = value2 * scaleFactor;
     1090      value3 = pivot * pivot * devex_;
     1091      if (reference(iSequence))
     1092        value3 += 1.0;
     1093      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
     1094      if (value < -tolerance) {
     1095        // store square in list
     1096        if (infeas[iSequence])
     1097          infeas[iSequence] = value * value; // already there
     1098        else
     1099          infeasible_->quickAdd(iSequence, value * value);
     1100      } else {
     1101        infeasible_->zero(iSequence);
     1102      }
     1103    }
     1104  }
     1105  // restore outgoing weight
     1106  if (sequenceOut >= 0)
     1107    weights_[sequenceOut] = outgoingWeight;
     1108  // make sure infeasibility on incoming is 0.0
     1109  infeasible_->zero(sequenceIn);
     1110  spareRow2->setNumElements(0);
     1111  //#define SOME_DEBUG_1
    11151112#ifdef SOME_DEBUG_1
    1116      // check for accuracy
    1117      int iCheck = 892;
    1118      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1119      int numberRows = model_->numberRows();
    1120      //int numberColumns = model_->numberColumns();
    1121      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1122           if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1123                     model_->getStatus(iCheck) != ClpSimplex::isFixed)
    1124                checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    1125      }
    1126 #endif
    1127      updates->setNumElements(0);
    1128      spareColumn1->setNumElements(0);
     1113  // check for accuracy
     1114  int iCheck = 892;
     1115  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1116  int numberRows = model_->numberRows();
     1117  //int numberColumns = model_->numberColumns();
     1118  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1119    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1120      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     1121  }
     1122#endif
     1123  updates->setNumElements(0);
     1124  spareColumn1->setNumElements(0);
    11291125}
    11301126// Update djs, weights for Steepest
    1131 void
    1132 ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates,
    1133                                         CoinIndexedVector * spareRow2,
    1134                                         CoinIndexedVector * spareColumn1,
    1135                                         CoinIndexedVector * spareColumn2)
     1127void ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector *updates,
     1128  CoinIndexedVector *spareRow2,
     1129  CoinIndexedVector *spareColumn1,
     1130  CoinIndexedVector *spareColumn2)
    11361131{
    1137      int j;
    1138      int number = 0;
    1139      int * index;
    1140      double * updateBy;
    1141      double * reducedCost;
    1142      double tolerance = model_->currentDualTolerance();
    1143      // we can't really trust infeasibilities if there is dual error
    1144      // this coding has to mimic coding in checkDualSolution
    1145      double error = CoinMin(1.0e-2, model_->largestDualError());
    1146      // allow tolerance at least slightly bigger than standard
    1147      tolerance = tolerance + error;
    1148      // for weights update we use pivotSequence
    1149      // unset in case sub flip
    1150      assert (pivotSequence_ >= 0);
    1151      assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
    1152      double * infeas = infeasible_->denseVector();
    1153      double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
    1154      assert (updates->getIndices()[0] == pivotSequence_);
    1155      pivotSequence_ = -1;
     1132  int j;
     1133  int number = 0;
     1134  int *index;
     1135  double *updateBy;
     1136  double *reducedCost;
     1137  double tolerance = model_->currentDualTolerance();
     1138  // we can't really trust infeasibilities if there is dual error
     1139  // this coding has to mimic coding in checkDualSolution
     1140  double error = CoinMin(1.0e-2, model_->largestDualError());
     1141  // allow tolerance at least slightly bigger than standard
     1142  tolerance = tolerance + error;
     1143  // for weights update we use pivotSequence
     1144  // unset in case sub flip
     1145  assert(pivotSequence_ >= 0);
     1146  assert(model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
     1147  double *infeas = infeasible_->denseVector();
     1148  double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
     1149  assert(updates->getIndices()[0] == pivotSequence_);
     1150  pivotSequence_ = -1;
    11561151#if 1
    1157      //updates->scanAndPack();
    1158      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    1159      //alternateWeights_->scanAndPack();
    1160 #if ALT_UPDATE_WEIGHTS !=2
    1161      model_->factorization()->updateColumnTranspose(spareRow2,
    1162                alternateWeights_);
    1163 #elif ALT_UPDATE_WEIGHTS==1
    1164      if (altVector[1]) {
    1165        int numberRows=model_->numberRows();
    1166        double * work1 = altVector[1]->denseVector();
    1167        double * worka = alternateWeights_->denseVector();
    1168        int iRow=-1;
    1169        double diff=1.0e-8;
    1170        for (int i=0;i<numberRows;i++) {
    1171          double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
    1172          double d=fabs(work1[i]-worka[i]);
    1173          if (dd>1.0e-6&&d>diff*dd) {
    1174            diff=d/dd;
    1175            iRow=i;
    1176         }
    1177        }
    1178        if (iRow>=0)
    1179         printf("largest difference of %g (%g,%g) on row %d\n",
    1180                 diff,work1[iRow],worka[iRow],iRow);
    1181      }
     1152  //updates->scanAndPack();
     1153  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     1154  //alternateWeights_->scanAndPack();
     1155#if ALT_UPDATE_WEIGHTS != 2
     1156  model_->factorization()->updateColumnTranspose(spareRow2,
     1157    alternateWeights_);
     1158#elif ALT_UPDATE_WEIGHTS == 1
     1159  if (altVector[1]) {
     1160    int numberRows = model_->numberRows();
     1161    double *work1 = altVector[1]->denseVector();
     1162    double *worka = alternateWeights_->denseVector();
     1163    int iRow = -1;
     1164    double diff = 1.0e-8;
     1165    for (int i = 0; i < numberRows; i++) {
     1166      double dd = CoinMax(fabs(work1[i]), fabs(worka[i]));
     1167      double d = fabs(work1[i] - worka[i]);
     1168      if (dd > 1.0e-6 && d > diff * dd) {
     1169        diff = d / dd;
     1170        iRow = i;
     1171      }
     1172    }
     1173    if (iRow >= 0)
     1174      printf("largest difference of %g (%g,%g) on row %d\n",
     1175        diff, work1[iRow], worka[iRow], iRow);
     1176  }
    11821177#endif
    11831178#else
    1184      model_->factorization()->updateTwoColumnsTranspose(spareRow2, updates,alternateWeights_);
    1185 #endif
    1186      // and we can see if reference
    1187      int sequenceIn = model_->sequenceIn();
    1188      double referenceIn;
    1189      if (mode_ != 1) {
    1190           if(reference(sequenceIn))
    1191                referenceIn = 1.0;
     1179  model_->factorization()->updateTwoColumnsTranspose(spareRow2, updates, alternateWeights_);
     1180#endif
     1181  // and we can see if reference
     1182  int sequenceIn = model_->sequenceIn();
     1183  double referenceIn;
     1184  if (mode_ != 1) {
     1185    if (reference(sequenceIn))
     1186      referenceIn = 1.0;
     1187    else
     1188      referenceIn = 0.0;
     1189  } else {
     1190    referenceIn = -1.0;
     1191  }
     1192  // save outgoing weight round update
     1193  double outgoingWeight = 0.0;
     1194  int sequenceOut = model_->sequenceOut();
     1195  if (sequenceOut >= 0)
     1196    outgoingWeight = weights_[sequenceOut];
     1197  // update row weights here so we can scale alternateWeights_
     1198  // update weights
     1199  double *weight;
     1200  double *other = alternateWeights_->denseVector();
     1201  int numberColumns = model_->numberColumns();
     1202  // if if (model_->clpMatrix()->canCombine(model_, pi1)) no need to index
     1203  // rows
     1204  reducedCost = model_->djRegion(0);
     1205  int addSequence = model_->numberColumns();
     1206  ;
     1207
     1208  number = updates->getNumElements();
     1209  index = updates->getIndices();
     1210  updateBy = updates->denseVector();
     1211  weight = weights_ + numberColumns;
     1212#ifdef CLP_USER_DRIVEN
     1213  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
     1214#endif
     1215
     1216  for (j = 0; j < number; j++) {
     1217    double thisWeight;
     1218    double pivot;
     1219    double modification;
     1220    double pivotSquared;
     1221    int iSequence = index[j];
     1222    double value2 = updateBy[j];
     1223    ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     1224    double value;
     1225
     1226    switch (status) {
     1227
     1228    case ClpSimplex::basic:
     1229      infeasible_->zero(iSequence + addSequence);
     1230      reducedCost[iSequence] = 0.0;
     1231    case ClpSimplex::isFixed:
     1232      break;
     1233    case ClpSimplex::isFree:
     1234    case ClpSimplex::superBasic:
     1235      value = reducedCost[iSequence] - value2;
     1236      modification = other[iSequence];
     1237      thisWeight = weight[iSequence];
     1238      // row has -1
     1239      pivot = value2 * scaleFactor;
     1240      pivotSquared = pivot * pivot;
     1241
     1242      thisWeight += pivotSquared * devex_ + pivot * modification;
     1243      reducedCost[iSequence] = value;
     1244      if (thisWeight < TRY_NORM) {
     1245        if (mode_ == 1) {
     1246          // steepest
     1247          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     1248        } else {
     1249          // exact
     1250          thisWeight = referenceIn * pivotSquared;
     1251          if (reference(iSequence + numberColumns))
     1252            thisWeight += 1.0;
     1253          thisWeight = CoinMax(thisWeight, TRY_NORM);
     1254        }
     1255      }
     1256      weight[iSequence] = thisWeight;
     1257      if (fabs(value) > FREE_ACCEPT * tolerance) {
     1258        // we are going to bias towards free (but only if reasonable)
     1259        value *= FREE_BIAS;
     1260        // store square in list
     1261        if (infeas[iSequence + addSequence])
     1262          infeas[iSequence + addSequence] = value * value; // already there
     1263        else
     1264          infeasible_->quickAdd(iSequence + addSequence, value * value);
     1265      } else {
     1266        infeasible_->zero(iSequence + addSequence);
     1267      }
     1268      break;
     1269    case ClpSimplex::atUpperBound:
     1270      value = reducedCost[iSequence] - value2;
     1271      modification = other[iSequence];
     1272      thisWeight = weight[iSequence];
     1273      // row has -1
     1274      pivot = value2 * scaleFactor;
     1275      pivotSquared = pivot * pivot;
     1276
     1277      thisWeight += pivotSquared * devex_ + pivot * modification;
     1278      reducedCost[iSequence] = value;
     1279      if (thisWeight < TRY_NORM) {
     1280        if (mode_ == 1) {
     1281          // steepest
     1282          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     1283        } else {
     1284          // exact
     1285          thisWeight = referenceIn * pivotSquared;
     1286          if (reference(iSequence + numberColumns))
     1287            thisWeight += 1.0;
     1288          thisWeight = CoinMax(thisWeight, TRY_NORM);
     1289        }
     1290      }
     1291      weight[iSequence] = thisWeight;
     1292      iSequence += addSequence;
     1293      if (value > tolerance) {
     1294        // store square in list
     1295#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1296        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
     1297#else
     1298        value *= value;
     1299#endif
     1300        if (infeas[iSequence])
     1301          infeas[iSequence] = value; // already there
     1302        else
     1303          infeasible_->quickAdd(iSequence, value);
     1304      } else {
     1305        infeasible_->zero(iSequence);
     1306      }
     1307      break;
     1308    case ClpSimplex::atLowerBound:
     1309      value = reducedCost[iSequence] - value2;
     1310      modification = other[iSequence];
     1311      thisWeight = weight[iSequence];
     1312      // row has -1
     1313      pivot = value2 * scaleFactor;
     1314      pivotSquared = pivot * pivot;
     1315
     1316      thisWeight += pivotSquared * devex_ + pivot * modification;
     1317      reducedCost[iSequence] = value;
     1318      if (thisWeight < TRY_NORM) {
     1319        if (mode_ == 1) {
     1320          // steepest
     1321          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     1322        } else {
     1323          // exact
     1324          thisWeight = referenceIn * pivotSquared;
     1325          if (reference(iSequence + numberColumns))
     1326            thisWeight += 1.0;
     1327          thisWeight = CoinMax(thisWeight, TRY_NORM);
     1328        }
     1329      }
     1330      weight[iSequence] = thisWeight;
     1331      iSequence += addSequence;
     1332      if (value < -tolerance) {
     1333        // store square in list
     1334#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1335        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
     1336#else
     1337        value *= value;
     1338#endif
     1339        if (infeas[iSequence])
     1340          infeas[iSequence] = value; // already there
     1341        else
     1342          infeasible_->quickAdd(iSequence, value);
     1343      } else {
     1344        infeasible_->zero(iSequence);
     1345      }
     1346    }
     1347  }
     1348#ifdef CLP_USER_DRIVEN
     1349  // could go parallel?
     1350  //model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
     1351#endif
     1352  // put row of tableau in rowArray and columnArray (packed)
     1353  // get subset which have nonzero tableau elements
     1354  int returnCode = transposeTimes2(updates, spareColumn1,
     1355    alternateWeights_, spareColumn2, spareRow2,
     1356    -scaleFactor);
     1357  // zero updateBy
     1358  CoinZeroN(updateBy, number);
     1359  alternateWeights_->clear();
     1360  // columns
     1361  assert(scaleFactor);
     1362  number = spareColumn1->getNumElements();
     1363  index = spareColumn1->getIndices();
     1364  updateBy = spareColumn1->denseVector();
     1365  if (returnCode != 2 && infeasibilitiesState_) {
     1366    //spareColumn1->clear();
     1367    redoInfeasibilities();
     1368  }
     1369  if (returnCode == 1) {
     1370    // most work already done
     1371    for (j = 0; j < number; j++) {
     1372      int iSequence = index[j];
     1373      double value = updateBy[j];
     1374      if (value) {
     1375        updateBy[j] = 0.0;
     1376        infeasible_->quickAdd(iSequence, value);
     1377      } else {
     1378        infeasible_->zero(iSequence);
     1379      }
     1380    }
     1381  } else if (returnCode == 0) {
     1382    reducedCost = model_->djRegion(1);
     1383
     1384    for (j = 0; j < number; j++) {
     1385      int iSequence = index[j];
     1386      double value = reducedCost[iSequence];
     1387      double value2 = updateBy[j];
     1388      updateBy[j] = 0.0;
     1389      value -= value2;
     1390      reducedCost[iSequence] = value;
     1391      ClpSimplex::Status status = model_->getStatus(iSequence);
     1392
     1393      switch (status) {
     1394
     1395      case ClpSimplex::basic:
     1396      case ClpSimplex::isFixed:
     1397        break;
     1398      case ClpSimplex::isFree:
     1399      case ClpSimplex::superBasic:
     1400        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1401          // we are going to bias towards free (but only if reasonable)
     1402          value *= FREE_BIAS;
     1403          // store square in list
     1404          if (infeas[iSequence])
     1405            infeas[iSequence] = value * value; // already there
    11921406          else
    1193                referenceIn = 0.0;
    1194      } else {
    1195           referenceIn = -1.0;
    1196      }
    1197      // save outgoing weight round update
    1198      double outgoingWeight = 0.0;
    1199      int sequenceOut = model_->sequenceOut();
    1200      if (sequenceOut >= 0)
    1201           outgoingWeight = weights_[sequenceOut];
    1202      // update row weights here so we can scale alternateWeights_
    1203      // update weights
    1204      double * weight;
    1205      double * other = alternateWeights_->denseVector();
    1206      int numberColumns = model_->numberColumns();
    1207      // if if (model_->clpMatrix()->canCombine(model_, pi1)) no need to index
    1208      // rows
    1209      reducedCost = model_->djRegion(0);
    1210      int addSequence = model_->numberColumns();;
    1211 
    1212      number = updates->getNumElements();
    1213      index = updates->getIndices();
    1214      updateBy = updates->denseVector();
    1215      weight = weights_ + numberColumns;
    1216 #ifdef CLP_USER_DRIVEN
    1217      model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
    1218 #endif
    1219 
    1220      for (j = 0; j < number; j++) {
    1221           double thisWeight;
    1222           double pivot;
    1223           double modification;
    1224           double pivotSquared;
    1225           int iSequence = index[j];
    1226           double value2 = updateBy[j];
    1227           ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    1228           double value;
    1229 
    1230           switch(status) {
    1231 
    1232           case ClpSimplex::basic:
    1233                infeasible_->zero(iSequence + addSequence);
    1234                reducedCost[iSequence] = 0.0;
    1235           case ClpSimplex::isFixed:
    1236                break;
    1237           case ClpSimplex::isFree:
    1238           case ClpSimplex::superBasic:
    1239                value = reducedCost[iSequence] - value2;
    1240                modification = other[iSequence];
    1241                thisWeight = weight[iSequence];
    1242                // row has -1
    1243                pivot = value2 * scaleFactor;
    1244                pivotSquared = pivot * pivot;
    1245 
    1246                thisWeight += pivotSquared * devex_ + pivot * modification;
    1247                reducedCost[iSequence] = value;
    1248                if (thisWeight < TRY_NORM) {
    1249                     if (mode_ == 1) {
    1250                          // steepest
    1251                          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    1252                     } else {
    1253                          // exact
    1254                          thisWeight = referenceIn * pivotSquared;
    1255                          if (reference(iSequence + numberColumns))
    1256                               thisWeight += 1.0;
    1257                          thisWeight = CoinMax(thisWeight, TRY_NORM);
    1258                     }
    1259                }
    1260                weight[iSequence] = thisWeight;
    1261                if (fabs(value) > FREE_ACCEPT * tolerance) {
    1262                     // we are going to bias towards free (but only if reasonable)
    1263                     value *= FREE_BIAS;
    1264                     // store square in list
    1265                     if (infeas[iSequence+addSequence])
    1266                          infeas[iSequence+addSequence] = value * value; // already there
    1267                     else
    1268                          infeasible_->quickAdd(iSequence + addSequence, value * value);
    1269                } else {
    1270                     infeasible_->zero(iSequence + addSequence);
    1271                }
    1272                break;
    1273           case ClpSimplex::atUpperBound:
    1274                value = reducedCost[iSequence] - value2;
    1275                modification = other[iSequence];
    1276                thisWeight = weight[iSequence];
    1277                // row has -1
    1278                pivot = value2 * scaleFactor;
    1279                pivotSquared = pivot * pivot;
    1280 
    1281                thisWeight += pivotSquared * devex_ + pivot * modification;
    1282                reducedCost[iSequence] = value;
    1283                if (thisWeight < TRY_NORM) {
    1284                     if (mode_ == 1) {
    1285                          // steepest
    1286                          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    1287                     } else {
    1288                          // exact
    1289                          thisWeight = referenceIn * pivotSquared;
    1290                          if (reference(iSequence + numberColumns))
    1291                               thisWeight += 1.0;
    1292                          thisWeight = CoinMax(thisWeight, TRY_NORM);
    1293                     }
    1294                }
    1295                weight[iSequence] = thisWeight;
    1296                iSequence += addSequence;
    1297                if (value > tolerance) {
    1298                     // store square in list
    1299 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1300                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
    1301 #else
    1302                     value *= value;
    1303 #endif
    1304                     if (infeas[iSequence])
    1305                          infeas[iSequence] = value; // already there
    1306                     else
    1307                          infeasible_->quickAdd(iSequence , value);
    1308                } else {
    1309                     infeasible_->zero(iSequence);
    1310                }
    1311                break;
    1312           case ClpSimplex::atLowerBound:
    1313                value = reducedCost[iSequence] - value2;
    1314                modification = other[iSequence];
    1315                thisWeight = weight[iSequence];
    1316                // row has -1
    1317                pivot = value2 * scaleFactor;
    1318                pivotSquared = pivot * pivot;
    1319 
    1320                thisWeight += pivotSquared * devex_ + pivot * modification;
    1321                reducedCost[iSequence] = value;
    1322                if (thisWeight < TRY_NORM) {
    1323                     if (mode_ == 1) {
    1324                          // steepest
    1325                          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    1326                     } else {
    1327                          // exact
    1328                          thisWeight = referenceIn * pivotSquared;
    1329                          if (reference(iSequence + numberColumns))
    1330                               thisWeight += 1.0;
    1331                          thisWeight = CoinMax(thisWeight, TRY_NORM);
    1332                     }
    1333                }
    1334                weight[iSequence] = thisWeight;
    1335                iSequence += addSequence;
    1336                if (value < -tolerance) {
    1337                     // store square in list
    1338 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1339                     value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
    1340 #else
    1341                     value *= value;
    1342 #endif
    1343                     if (infeas[iSequence])
    1344                          infeas[iSequence] = value; // already there
    1345                     else
    1346                          infeasible_->quickAdd(iSequence, value);
    1347                } else {
    1348                     infeasible_->zero(iSequence);
    1349                }
    1350           }
    1351      }
    1352 #ifdef CLP_USER_DRIVEN
    1353      // could go parallel?
    1354      //model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
    1355 #endif
    1356      // put row of tableau in rowArray and columnArray (packed)
    1357      // get subset which have nonzero tableau elements
    1358      int returnCode=transposeTimes2(updates, spareColumn1,
    1359                                     alternateWeights_, spareColumn2, spareRow2,
    1360                                     -scaleFactor);
    1361      // zero updateBy
    1362      CoinZeroN(updateBy, number);
    1363      alternateWeights_->clear();
    1364      // columns
    1365      assert (scaleFactor);
    1366      number = spareColumn1->getNumElements();
    1367      index = spareColumn1->getIndices();
    1368      updateBy = spareColumn1->denseVector();
    1369      if (returnCode != 2 && infeasibilitiesState_) {
    1370        //spareColumn1->clear();
    1371        redoInfeasibilities();
    1372      }
    1373      if (returnCode == 1) {
    1374        // most work already done
    1375        for (j = 0; j < number; j++) {
    1376          int iSequence = index[j];
    1377          double value = updateBy[j];
    1378          if (value) {
    1379            updateBy[j]=0.0;
    1380            infeasible_->quickAdd(iSequence, value);
    1381          } else {
    1382            infeasible_->zero(iSequence);
    1383          }
    1384        }
    1385      } else if (returnCode==0) {
    1386        reducedCost = model_->djRegion(1);
    1387 
    1388        for (j = 0; j < number; j++) {
    1389          int iSequence = index[j];
    1390          double value = reducedCost[iSequence];
    1391          double value2 = updateBy[j];
    1392          updateBy[j] = 0.0;
    1393          value -= value2;
    1394          reducedCost[iSequence] = value;
    1395          ClpSimplex::Status status = model_->getStatus(iSequence);
    1396          
    1397          switch(status) {
    1398            
    1399          case ClpSimplex::basic:
    1400          case ClpSimplex::isFixed:
    1401            break;
    1402          case ClpSimplex::isFree:
    1403          case ClpSimplex::superBasic:
    1404            if (fabs(value) > FREE_ACCEPT * tolerance) {
    1405              // we are going to bias towards free (but only if reasonable)
    1406              value *= FREE_BIAS;
    1407              // store square in list
    1408              if (infeas[iSequence])
    1409                infeas[iSequence] = value * value; // already there
    1410              else
    1411                infeasible_->quickAdd(iSequence, value * value);
    1412            } else {
    1413              infeasible_->zero(iSequence);
    1414            }
    1415            break;
    1416          case ClpSimplex::atUpperBound:
    1417            if (value > tolerance) {
    1418              // store square in list
    1419              if (infeas[iSequence])
    1420                infeas[iSequence] = value * value; // already there
    1421              else
    1422                infeasible_->quickAdd(iSequence, value * value);
    1423            } else {
    1424              infeasible_->zero(iSequence);
    1425            }
    1426            break;
    1427          case ClpSimplex::atLowerBound:
    1428            if (value < -tolerance) {
    1429              // store square in list
    1430              if (infeas[iSequence])
    1431                infeas[iSequence] = value * value; // already there
    1432              else
    1433                infeasible_->quickAdd(iSequence, value * value);
    1434            } else {
    1435              infeasible_->zero(iSequence);
    1436            }
    1437          }
    1438        }
    1439      } else {
    1440        //assert(!number);
    1441      }
    1442      // restore outgoing weight
    1443      if (sequenceOut >= 0) {
    1444        //#define GROW_REFERENCE
     1407            infeasible_->quickAdd(iSequence, value * value);
     1408        } else {
     1409          infeasible_->zero(iSequence);
     1410        }
     1411        break;
     1412      case ClpSimplex::atUpperBound:
     1413        if (value > tolerance) {
     1414          // store square in list
     1415          if (infeas[iSequence])
     1416            infeas[iSequence] = value * value; // already there
     1417          else
     1418            infeasible_->quickAdd(iSequence, value * value);
     1419        } else {
     1420          infeasible_->zero(iSequence);
     1421        }
     1422        break;
     1423      case ClpSimplex::atLowerBound:
     1424        if (value < -tolerance) {
     1425          // store square in list
     1426          if (infeas[iSequence])
     1427            infeas[iSequence] = value * value; // already there
     1428          else
     1429            infeasible_->quickAdd(iSequence, value * value);
     1430        } else {
     1431          infeasible_->zero(iSequence);
     1432        }
     1433      }
     1434    }
     1435  } else {
     1436    //assert(!number);
     1437  }
     1438  // restore outgoing weight
     1439  if (sequenceOut >= 0) {
     1440    //#define GROW_REFERENCE
    14451441#ifdef GROW_REFERENCE
    1446        if (!reference(sequenceOut)) {
    1447         outgoingWeight += 1.0;
    1448          setReference(sequenceOut,true);
    1449        }
    1450 #endif
    1451        weights_[sequenceOut] = outgoingWeight;
    1452        //if (model_->getStatus(sequenceOut) != ClpSimplex::basic &&
    1453        //  model_->getStatus(sequenceOut) != ClpSimplex::isFixed)
    1454        //checkAccuracy(sequenceOut, 1.0e-1, updates, spareRow2);
    1455      }
    1456      // make sure infeasibility on incoming is 0.0
    1457      infeasible_->zero(sequenceIn);
    1458      spareColumn2->setNumElements(0);
    1459      //#define SOME_DEBUG_1
     1442    if (!reference(sequenceOut)) {
     1443      outgoingWeight += 1.0;
     1444      setReference(sequenceOut, true);
     1445    }
     1446#endif
     1447    weights_[sequenceOut] = outgoingWeight;
     1448    //if (model_->getStatus(sequenceOut) != ClpSimplex::basic &&
     1449    //  model_->getStatus(sequenceOut) != ClpSimplex::isFixed)
     1450    //checkAccuracy(sequenceOut, 1.0e-1, updates, spareRow2);
     1451  }
     1452  // make sure infeasibility on incoming is 0.0
     1453  infeasible_->zero(sequenceIn);
     1454  spareColumn2->setNumElements(0);
     1455  //#define SOME_DEBUG_1
    14601456#ifdef SOME_DEBUG_1
    1461      // check for accuracy
    1462      int iCheck = 892;
    1463      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1464      int numberRows = model_->numberRows();
    1465      //int numberColumns = model_->numberColumns();
    1466      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1467           if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1468                     model_->getStatus(iCheck) != ClpSimplex::isFixed)
    1469                checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    1470      }
    1471 #endif
    1472      updates->setNumElements(0);
    1473      spareColumn1->setNumElements(0);
     1457  // check for accuracy
     1458  int iCheck = 892;
     1459  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1460  int numberRows = model_->numberRows();
     1461  //int numberColumns = model_->numberColumns();
     1462  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1463    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1464      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     1465  }
     1466#endif
     1467  updates->setNumElements(0);
     1468  spareColumn1->setNumElements(0);
    14741469}
    14751470// Update djs, weights for Devex
    1476 void
    1477 ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
    1478                                       CoinIndexedVector * spareRow2,
    1479                                       CoinIndexedVector * spareColumn1,
    1480                                       CoinIndexedVector * spareColumn2)
     1471void ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector *updates,
     1472  CoinIndexedVector *spareRow2,
     1473  CoinIndexedVector *spareColumn1,
     1474  CoinIndexedVector *spareColumn2)
    14811475{
    1482      int iSection, j;
    1483      int number = 0;
    1484      int * index;
    1485      double * updateBy;
    1486      double * reducedCost;
    1487      // dj could be very small (or even zero - take care)
    1488      double dj = model_->dualIn();
    1489      double tolerance = model_->currentDualTolerance();
    1490      // we can't really trust infeasibilities if there is dual error
    1491      // this coding has to mimic coding in checkDualSolution
    1492      double error = CoinMin(1.0e-2, model_->largestDualError());
    1493      // allow tolerance at least slightly bigger than standard
    1494      tolerance = tolerance +  error;
    1495      int pivotRow = model_->pivotRow();
    1496      double * infeas = infeasible_->denseVector();
    1497      //updates->scanAndPack();
    1498      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    1499 
    1500      // put row of tableau in rowArray and columnArray
    1501      model_->clpMatrix()->transposeTimes(model_, -1.0,
    1502                                          updates, spareColumn2, spareColumn1);
    1503      // normal
    1504      for (iSection = 0; iSection < 2; iSection++) {
    1505 
    1506           reducedCost = model_->djRegion(iSection);
    1507           int addSequence;
    1508 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1509           double slack_multiplier;
    1510 #endif
    1511 
    1512           if (!iSection) {
    1513                number = updates->getNumElements();
    1514                index = updates->getIndices();
    1515                updateBy = updates->denseVector();
    1516                addSequence = model_->numberColumns();
    1517 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1518                slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
    1519 #endif
    1520           } else {
    1521                number = spareColumn1->getNumElements();
    1522                index = spareColumn1->getIndices();
    1523                updateBy = spareColumn1->denseVector();
    1524                addSequence = 0;
    1525 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1526                slack_multiplier = 1.0;
    1527 #endif
    1528           }
    1529 
    1530           for (j = 0; j < number; j++) {
    1531                int iSequence = index[j];
    1532                double value = reducedCost[iSequence];
    1533                value -= updateBy[j];
    1534                updateBy[j] = 0.0;
    1535                reducedCost[iSequence] = value;
    1536                ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    1537 
    1538                switch(status) {
    1539 
    1540                case ClpSimplex::basic:
    1541                     infeasible_->zero(iSequence + addSequence);
    1542                case ClpSimplex::isFixed:
    1543                     break;
    1544                case ClpSimplex::isFree:
    1545                case ClpSimplex::superBasic:
    1546                     if (fabs(value) > FREE_ACCEPT * tolerance) {
    1547                          // we are going to bias towards free (but only if reasonable)
    1548                          value *= FREE_BIAS;
    1549                          // store square in list
    1550                          if (infeas[iSequence+addSequence])
    1551                               infeas[iSequence+addSequence] = value * value; // already there
    1552                          else
    1553                               infeasible_->quickAdd(iSequence + addSequence, value * value);
    1554                     } else {
    1555                          infeasible_->zero(iSequence + addSequence);
    1556                     }
    1557                     break;
    1558                case ClpSimplex::atUpperBound:
    1559                     iSequence += addSequence;
    1560                     if (value > tolerance) {
     1476  int iSection, j;
     1477  int number = 0;
     1478  int *index;
     1479  double *updateBy;
     1480  double *reducedCost;
     1481  // dj could be very small (or even zero - take care)
     1482  double dj = model_->dualIn();
     1483  double tolerance = model_->currentDualTolerance();
     1484  // we can't really trust infeasibilities if there is dual error
     1485  // this coding has to mimic coding in checkDualSolution
     1486  double error = CoinMin(1.0e-2, model_->largestDualError());
     1487  // allow tolerance at least slightly bigger than standard
     1488  tolerance = tolerance + error;
     1489  int pivotRow = model_->pivotRow();
     1490  double *infeas = infeasible_->denseVector();
     1491  //updates->scanAndPack();
     1492  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     1493
     1494  // put row of tableau in rowArray and columnArray
     1495  model_->clpMatrix()->transposeTimes(model_, -1.0,
     1496    updates, spareColumn2, spareColumn1);
     1497  // normal
     1498  for (iSection = 0; iSection < 2; iSection++) {
     1499
     1500    reducedCost = model_->djRegion(iSection);
     1501    int addSequence;
    15611502#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1562                         value *= value*slack_multiplier;
     1503    double slack_multiplier;
     1504#endif
     1505
     1506    if (!iSection) {
     1507      number = updates->getNumElements();
     1508      index = updates->getIndices();
     1509      updateBy = updates->denseVector();
     1510      addSequence = model_->numberColumns();
     1511#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1512      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
     1513#endif
     1514    } else {
     1515      number = spareColumn1->getNumElements();
     1516      index = spareColumn1->getIndices();
     1517      updateBy = spareColumn1->denseVector();
     1518      addSequence = 0;
     1519#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1520      slack_multiplier = 1.0;
     1521#endif
     1522    }
     1523
     1524    for (j = 0; j < number; j++) {
     1525      int iSequence = index[j];
     1526      double value = reducedCost[iSequence];
     1527      value -= updateBy[j];
     1528      updateBy[j] = 0.0;
     1529      reducedCost[iSequence] = value;
     1530      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     1531
     1532      switch (status) {
     1533
     1534      case ClpSimplex::basic:
     1535        infeasible_->zero(iSequence + addSequence);
     1536      case ClpSimplex::isFixed:
     1537        break;
     1538      case ClpSimplex::isFree:
     1539      case ClpSimplex::superBasic:
     1540        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1541          // we are going to bias towards free (but only if reasonable)
     1542          value *= FREE_BIAS;
     1543          // store square in list
     1544          if (infeas[iSequence + addSequence])
     1545            infeas[iSequence + addSequence] = value * value; // already there
     1546          else
     1547            infeasible_->quickAdd(iSequence + addSequence, value * value);
     1548        } else {
     1549          infeasible_->zero(iSequence + addSequence);
     1550        }
     1551        break;
     1552      case ClpSimplex::atUpperBound:
     1553        iSequence += addSequence;
     1554        if (value > tolerance) {
     1555#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1556          value *= value * slack_multiplier;
    15631557#else
    1564                         value *= value;
    1565 #endif
    1566                          // store square in list
    1567                          if (infeas[iSequence])
    1568                               infeas[iSequence] = value; // already there
    1569                          else
    1570                               infeasible_->quickAdd(iSequence, value);
    1571                     } else {
    1572                          infeasible_->zero(iSequence);
    1573                     }
    1574                     break;
    1575                case ClpSimplex::atLowerBound:
    1576                     iSequence += addSequence;
    1577                     if (value < -tolerance) {
     1558          value *= value;
     1559#endif
     1560          // store square in list
     1561          if (infeas[iSequence])
     1562            infeas[iSequence] = value; // already there
     1563          else
     1564            infeasible_->quickAdd(iSequence, value);
     1565        } else {
     1566          infeasible_->zero(iSequence);
     1567        }
     1568        break;
     1569      case ClpSimplex::atLowerBound:
     1570        iSequence += addSequence;
     1571        if (value < -tolerance) {
    15781572#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1579                         value *= value*slack_multiplier;
     1573          value *= value * slack_multiplier;
    15801574#else
    1581                         value *= value;
    1582 #endif
    1583                          // store square in list
    1584                          if (infeas[iSequence])
    1585                               infeas[iSequence] = value; // already there
    1586                          else
    1587                               infeasible_->quickAdd(iSequence, value);
    1588                     } else {
    1589                          infeasible_->zero(iSequence);
    1590                     }
    1591                }
    1592           }
    1593      }
    1594      // They are empty
    1595      updates->setNumElements(0);
    1596      spareColumn1->setNumElements(0);
    1597      // make sure infeasibility on incoming is 0.0
    1598      int sequenceIn = model_->sequenceIn();
    1599      infeasible_->zero(sequenceIn);
    1600      // for weights update we use pivotSequence
    1601      if (pivotSequence_ >= 0) {
    1602           pivotRow = pivotSequence_;
    1603           // unset in case sub flip
    1604           pivotSequence_ = -1;
    1605           // make sure infeasibility on incoming is 0.0
    1606           const int * pivotVariable = model_->pivotVariable();
    1607           sequenceIn = pivotVariable[pivotRow];
    1608           infeasible_->zero(sequenceIn);
    1609           // and we can see if reference
    1610           //double referenceIn = 0.0;
    1611           //if (mode_ != 1 && reference(sequenceIn))
    1612           //   referenceIn = 1.0;
    1613           // save outgoing weight round update
    1614           double outgoingWeight = 0.0;
    1615           int sequenceOut = model_->sequenceOut();
    1616           if (sequenceOut >= 0)
    1617                outgoingWeight = weights_[sequenceOut];
    1618           // update weights
    1619           updates->setNumElements(0);
    1620           spareColumn1->setNumElements(0);
    1621           // might as well set dj to 1
    1622           dj = 1.0;
    1623           updates->insert(pivotRow, -dj);
    1624           model_->factorization()->updateColumnTranspose(spareRow2, updates);
    1625           // put row of tableau in rowArray and columnArray
    1626           model_->clpMatrix()->transposeTimes(model_, -1.0,
    1627                                               updates, spareColumn2, spareColumn1);
    1628           double * weight;
    1629           int numberColumns = model_->numberColumns();
    1630           // rows
    1631           number = updates->getNumElements();
    1632           index = updates->getIndices();
    1633           updateBy = updates->denseVector();
    1634           weight = weights_ + numberColumns;
    1635 
    1636           assert (devex_ > 0.0);
    1637           for (j = 0; j < number; j++) {
    1638                int iSequence = index[j];
    1639                double thisWeight = weight[iSequence];
    1640                // row has -1
    1641                double pivot = - updateBy[iSequence];
    1642                updateBy[iSequence] = 0.0;
    1643                double value = pivot * pivot * devex_;
    1644                if (reference(iSequence + numberColumns))
    1645                     value += 1.0;
    1646                weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1647           }
    1648 
    1649           // columns
    1650           weight = weights_;
    1651 
    1652           number = spareColumn1->getNumElements();
    1653           index = spareColumn1->getIndices();
    1654           updateBy = spareColumn1->denseVector();
    1655           for (j = 0; j < number; j++) {
    1656                int iSequence = index[j];
    1657                double thisWeight = weight[iSequence];
    1658                // row has -1
    1659                double pivot = updateBy[iSequence];
    1660                updateBy[iSequence] = 0.0;
    1661                double value = pivot * pivot * devex_;
    1662                if (reference(iSequence))
    1663                     value += 1.0;
    1664                weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1665           }
    1666           // restore outgoing weight
    1667           if (sequenceOut >= 0)
    1668                weights_[sequenceOut] = outgoingWeight;
    1669           spareColumn2->setNumElements(0);
    1670           //#define SOME_DEBUG_1
     1575          value *= value;
     1576#endif
     1577          // store square in list
     1578          if (infeas[iSequence])
     1579            infeas[iSequence] = value; // already there
     1580          else
     1581            infeasible_->quickAdd(iSequence, value);
     1582        } else {
     1583          infeasible_->zero(iSequence);
     1584        }
     1585      }
     1586    }
     1587  }
     1588  // They are empty
     1589  updates->setNumElements(0);
     1590  spareColumn1->setNumElements(0);
     1591  // make sure infeasibility on incoming is 0.0
     1592  int sequenceIn = model_->sequenceIn();
     1593  infeasible_->zero(sequenceIn);
     1594  // for weights update we use pivotSequence
     1595  if (pivotSequence_ >= 0) {
     1596    pivotRow = pivotSequence_;
     1597    // unset in case sub flip
     1598    pivotSequence_ = -1;
     1599    // make sure infeasibility on incoming is 0.0
     1600    const int *pivotVariable = model_->pivotVariable();
     1601    sequenceIn = pivotVariable[pivotRow];
     1602    infeasible_->zero(sequenceIn);
     1603    // and we can see if reference
     1604    //double referenceIn = 0.0;
     1605    //if (mode_ != 1 && reference(sequenceIn))
     1606    //   referenceIn = 1.0;
     1607    // save outgoing weight round update
     1608    double outgoingWeight = 0.0;
     1609    int sequenceOut = model_->sequenceOut();
     1610    if (sequenceOut >= 0)
     1611      outgoingWeight = weights_[sequenceOut];
     1612    // update weights
     1613    updates->setNumElements(0);
     1614    spareColumn1->setNumElements(0);
     1615    // might as well set dj to 1
     1616    dj = 1.0;
     1617    updates->insert(pivotRow, -dj);
     1618    model_->factorization()->updateColumnTranspose(spareRow2, updates);
     1619    // put row of tableau in rowArray and columnArray
     1620    model_->clpMatrix()->transposeTimes(model_, -1.0,
     1621      updates, spareColumn2, spareColumn1);
     1622    double *weight;
     1623    int numberColumns = model_->numberColumns();
     1624    // rows
     1625    number = updates->getNumElements();
     1626    index = updates->getIndices();
     1627    updateBy = updates->denseVector();
     1628    weight = weights_ + numberColumns;
     1629
     1630    assert(devex_ > 0.0);
     1631    for (j = 0; j < number; j++) {
     1632      int iSequence = index[j];
     1633      double thisWeight = weight[iSequence];
     1634      // row has -1
     1635      double pivot = -updateBy[iSequence];
     1636      updateBy[iSequence] = 0.0;
     1637      double value = pivot * pivot * devex_;
     1638      if (reference(iSequence + numberColumns))
     1639        value += 1.0;
     1640      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1641    }
     1642
     1643    // columns
     1644    weight = weights_;
     1645
     1646    number = spareColumn1->getNumElements();
     1647    index = spareColumn1->getIndices();
     1648    updateBy = spareColumn1->denseVector();
     1649    for (j = 0; j < number; j++) {
     1650      int iSequence = index[j];
     1651      double thisWeight = weight[iSequence];
     1652      // row has -1
     1653      double pivot = updateBy[iSequence];
     1654      updateBy[iSequence] = 0.0;
     1655      double value = pivot * pivot * devex_;
     1656      if (reference(iSequence))
     1657        value += 1.0;
     1658      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1659    }
     1660    // restore outgoing weight
     1661    if (sequenceOut >= 0)
     1662      weights_[sequenceOut] = outgoingWeight;
     1663    spareColumn2->setNumElements(0);
     1664    //#define SOME_DEBUG_1
    16711665#ifdef SOME_DEBUG_1
    1672           // check for accuracy
    1673           int iCheck = 892;
    1674           //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1675           int numberRows = model_->numberRows();
    1676           //int numberColumns = model_->numberColumns();
    1677           for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1678                if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1679                          model_->getStatus(iCheck) != ClpSimplex::isFixed)
    1680                     checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    1681           }
    1682 #endif
    1683           updates->setNumElements(0);
    1684           spareColumn1->setNumElements(0);
    1685      }
     1666    // check for accuracy
     1667    int iCheck = 892;
     1668    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1669    int numberRows = model_->numberRows();
     1670    //int numberColumns = model_->numberColumns();
     1671    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1672      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1673        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     1674    }
     1675#endif
     1676    updates->setNumElements(0);
     1677    spareColumn1->setNumElements(0);
     1678  }
    16861679}
    16871680// Update djs, weights for Steepest
    1688 void
    1689 ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates,
    1690           CoinIndexedVector * spareRow2,
    1691           CoinIndexedVector * spareColumn1,
    1692           CoinIndexedVector * spareColumn2)
     1681void ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector *updates,
     1682  CoinIndexedVector *spareRow2,
     1683  CoinIndexedVector *spareColumn1,
     1684  CoinIndexedVector *spareColumn2)
    16931685{
    1694      int iSection, j;
    1695      int number = 0;
    1696      int * index;
    1697      double * updateBy;
    1698      double * reducedCost;
    1699      // dj could be very small (or even zero - take care)
    1700      double dj = model_->dualIn();
    1701      double tolerance = model_->currentDualTolerance();
    1702      // we can't really trust infeasibilities if there is dual error
    1703      // this coding has to mimic coding in checkDualSolution
    1704      double error = CoinMin(1.0e-2, model_->largestDualError());
    1705      // allow tolerance at least slightly bigger than standard
    1706      tolerance = tolerance + error;
    1707      int pivotRow = model_->pivotRow();
    1708      double * infeas = infeasible_->denseVector();
    1709      //updates->scanAndPack();
    1710      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    1711 
    1712      // put row of tableau in rowArray and columnArray
    1713      model_->clpMatrix()->transposeTimes(model_, -1.0,
    1714                                          updates, spareColumn2, spareColumn1);
     1686  int iSection, j;
     1687  int number = 0;
     1688  int *index;
     1689  double *updateBy;
     1690  double *reducedCost;
     1691  // dj could be very small (or even zero - take care)
     1692  double dj = model_->dualIn();
     1693  double tolerance = model_->currentDualTolerance();
     1694  // we can't really trust infeasibilities if there is dual error
     1695  // this coding has to mimic coding in checkDualSolution
     1696  double error = CoinMin(1.0e-2, model_->largestDualError());
     1697  // allow tolerance at least slightly bigger than standard
     1698  tolerance = tolerance + error;
     1699  int pivotRow = model_->pivotRow();
     1700  double *infeas = infeasible_->denseVector();
     1701  //updates->scanAndPack();
     1702  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     1703
     1704  // put row of tableau in rowArray and columnArray
     1705  model_->clpMatrix()->transposeTimes(model_, -1.0,
     1706    updates, spareColumn2, spareColumn1);
    17151707#ifdef CLP_USER_DRIVEN
    1716      model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
    1717 #endif
    1718      // normal
    1719      for (iSection = 0; iSection < 2; iSection++) {
    1720 
    1721           reducedCost = model_->djRegion(iSection);
    1722           int addSequence;
    1723 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1724           double slack_multiplier;
    1725 #endif
    1726 
    1727           if (!iSection) {
    1728                number = updates->getNumElements();
    1729                index = updates->getIndices();
    1730                updateBy = updates->denseVector();
    1731                addSequence = model_->numberColumns();
    1732 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1733                slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
    1734 #endif
     1708  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
     1709#endif
     1710  // normal
     1711  for (iSection = 0; iSection < 2; iSection++) {
     1712
     1713    reducedCost = model_->djRegion(iSection);
     1714    int addSequence;
     1715#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1716    double slack_multiplier;
     1717#endif
     1718
     1719    if (!iSection) {
     1720      number = updates->getNumElements();
     1721      index = updates->getIndices();
     1722      updateBy = updates->denseVector();
     1723      addSequence = model_->numberColumns();
     1724#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1725      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
     1726#endif
     1727    } else {
     1728      number = spareColumn1->getNumElements();
     1729      index = spareColumn1->getIndices();
     1730      updateBy = spareColumn1->denseVector();
     1731      addSequence = 0;
     1732#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1733      slack_multiplier = 1.0;
     1734#endif
     1735    }
     1736
     1737    for (j = 0; j < number; j++) {
     1738      int iSequence = index[j];
     1739      double value = reducedCost[iSequence];
     1740      value -= updateBy[j];
     1741      updateBy[j] = 0.0;
     1742      reducedCost[iSequence] = value;
     1743      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     1744
     1745      switch (status) {
     1746
     1747      case ClpSimplex::basic:
     1748        infeasible_->zero(iSequence + addSequence);
     1749      case ClpSimplex::isFixed:
     1750        break;
     1751      case ClpSimplex::isFree:
     1752      case ClpSimplex::superBasic:
     1753        if (fabs(value) > FREE_ACCEPT * tolerance) {
     1754          // we are going to bias towards free (but only if reasonable)
     1755          value *= FREE_BIAS;
     1756          // store square in list
     1757          if (infeas[iSequence + addSequence])
     1758            infeas[iSequence + addSequence] = value * value; // already there
     1759          else
     1760            infeasible_->quickAdd(iSequence + addSequence, value * value);
     1761        } else {
     1762          infeasible_->zero(iSequence + addSequence);
     1763        }
     1764        break;
     1765      case ClpSimplex::atUpperBound:
     1766        iSequence += addSequence;
     1767        if (value > tolerance) {
     1768#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1769          value *= value * slack_multiplier;
     1770#else
     1771          value *= value;
     1772#endif
     1773          // store square in list
     1774          if (infeas[iSequence])
     1775            infeas[iSequence] = value; // already there
     1776          else
     1777            infeasible_->quickAdd(iSequence, value);
     1778        } else {
     1779          infeasible_->zero(iSequence);
     1780        }
     1781        break;
     1782      case ClpSimplex::atLowerBound:
     1783        iSequence += addSequence;
     1784        if (value < -tolerance) {
     1785#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
     1786          value *= value * slack_multiplier;
     1787#else
     1788          value *= value;
     1789#endif
     1790          // store square in list
     1791          if (infeas[iSequence])
     1792            infeas[iSequence] = value; // already there
     1793          else
     1794            infeasible_->quickAdd(iSequence, value);
     1795        } else {
     1796          infeasible_->zero(iSequence);
     1797        }
     1798      }
     1799    }
     1800  }
     1801  // we can zero out as will have to get pivot row
     1802  // ***** move
     1803  updates->setNumElements(0);
     1804  spareColumn1->setNumElements(0);
     1805  if (pivotRow >= 0) {
     1806    // make sure infeasibility on incoming is 0.0
     1807    int sequenceIn = model_->sequenceIn();
     1808    infeasible_->zero(sequenceIn);
     1809  }
     1810  // for weights update we use pivotSequence
     1811  pivotRow = pivotSequence_;
     1812  // unset in case sub flip
     1813  pivotSequence_ = -1;
     1814  if (pivotRow >= 0) {
     1815    // make sure infeasibility on incoming is 0.0
     1816    const int *pivotVariable = model_->pivotVariable();
     1817    int sequenceIn = pivotVariable[pivotRow];
     1818    assert(sequenceIn == model_->sequenceIn());
     1819    infeasible_->zero(sequenceIn);
     1820    // and we can see if reference
     1821    double referenceIn;
     1822    if (mode_ != 1) {
     1823      if (reference(sequenceIn))
     1824        referenceIn = 1.0;
     1825      else
     1826        referenceIn = 0.0;
     1827    } else {
     1828      referenceIn = -1.0;
     1829    }
     1830    // save outgoing weight round update
     1831    double outgoingWeight = 0.0;
     1832    int sequenceOut = model_->sequenceOut();
     1833    if (sequenceOut >= 0)
     1834      outgoingWeight = weights_[sequenceOut];
     1835    // update weights
     1836    updates->setNumElements(0);
     1837    spareColumn1->setNumElements(0);
     1838    // might as well set dj to 1
     1839    dj = -1.0;
     1840    updates->createPacked(1, &pivotRow, &dj);
     1841    model_->factorization()->updateColumnTranspose(spareRow2, updates);
     1842    bool needSubset = (mode_ < 4 || numberSwitched_ > 1 || mode_ >= 10);
     1843
     1844    double *weight;
     1845    double *other = alternateWeights_->denseVector();
     1846    int numberColumns = model_->numberColumns();
     1847    // rows
     1848    number = updates->getNumElements();
     1849    index = updates->getIndices();
     1850    updateBy = updates->denseVector();
     1851    weight = weights_ + numberColumns;
     1852    if (needSubset) {
     1853#if ALT_UPDATE_WEIGHTS != 2
     1854      model_->factorization()->updateColumnTranspose(spareRow2,
     1855        alternateWeights_);
     1856#elif ALT_UPDATE_WEIGHTS == 1
     1857      if (altVector[1]) {
     1858        int numberRows = model_->numberRows();
     1859        double *work1 = altVector[1]->denseVector();
     1860        double *worka = alternateWeights_->denseVector();
     1861        int iRow = -1;
     1862        double diff = 1.0e-8;
     1863        for (int i = 0; i < numberRows; i++) {
     1864          double dd = CoinMax(fabs(work1[i]), fabs(worka[i]));
     1865          double d = fabs(work1[i] - worka[i]);
     1866          if (dd > 1.0e-6 && d > diff * dd) {
     1867            diff = d / dd;
     1868            iRow = i;
     1869          }
     1870        }
     1871        if (iRow >= 0)
     1872          printf("largest2 difference of %g (%g,%g) on row %d\n",
     1873            diff, work1[iRow], worka[iRow], iRow);
     1874      }
     1875#endif
     1876      // do alternateWeights_ here so can scale
     1877      for (j = 0; j < number; j++) {
     1878        int iSequence = index[j];
     1879        assert(iSequence >= 0 && iSequence < model_->numberRows());
     1880        double thisWeight = weight[iSequence];
     1881        // row has -1
     1882        double pivot = -updateBy[j];
     1883        double modification = other[iSequence];
     1884        double pivotSquared = pivot * pivot;
     1885
     1886        thisWeight += pivotSquared * devex_ + pivot * modification;
     1887        if (thisWeight < TRY_NORM) {
     1888          if (mode_ == 1) {
     1889            // steepest
     1890            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    17351891          } else {
    1736                number = spareColumn1->getNumElements();
    1737                index = spareColumn1->getIndices();
    1738                updateBy = spareColumn1->denseVector();
    1739                addSequence = 0;
    1740 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1741                slack_multiplier = 1.0;
    1742 #endif
     1892            // exact
     1893            thisWeight = referenceIn * pivotSquared;
     1894            if (reference(iSequence + numberColumns))
     1895              thisWeight += 1.0;
     1896            thisWeight = CoinMax(thisWeight, TRY_NORM);
    17431897          }
    1744 
    1745           for (j = 0; j < number; j++) {
    1746                int iSequence = index[j];
    1747                double value = reducedCost[iSequence];
    1748                value -= updateBy[j];
    1749                updateBy[j] = 0.0;
    1750                reducedCost[iSequence] = value;
    1751                ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    1752 
    1753                switch(status) {
    1754 
    1755                case ClpSimplex::basic:
    1756                     infeasible_->zero(iSequence + addSequence);
    1757                case ClpSimplex::isFixed:
    1758                     break;
    1759                case ClpSimplex::isFree:
    1760                case ClpSimplex::superBasic:
    1761                     if (fabs(value) > FREE_ACCEPT * tolerance) {
    1762                          // we are going to bias towards free (but only if reasonable)
    1763                          value *= FREE_BIAS;
    1764                          // store square in list
    1765                          if (infeas[iSequence+addSequence])
    1766                               infeas[iSequence+addSequence] = value * value; // already there
    1767                          else
    1768                               infeasible_->quickAdd(iSequence + addSequence, value * value);
    1769                     } else {
    1770                          infeasible_->zero(iSequence + addSequence);
    1771                     }
    1772                     break;
    1773                case ClpSimplex::atUpperBound:
    1774                     iSequence += addSequence;
    1775                     if (value > tolerance) {
    1776 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1777                         value *= value*slack_multiplier;
    1778 #else
    1779                         value *= value;
    1780 #endif
    1781                          // store square in list
    1782                          if (infeas[iSequence])
    1783                               infeas[iSequence] = value; // already there
    1784                          else
    1785                               infeasible_->quickAdd(iSequence, value);
    1786                     } else {
    1787                          infeasible_->zero(iSequence);
    1788                     }
    1789                     break;
    1790                case ClpSimplex::atLowerBound:
    1791                     iSequence += addSequence;
    1792                     if (value < -tolerance) {
    1793 #ifdef CLP_PRIMAL_SLACK_MULTIPLIER
    1794                         value *= value*slack_multiplier;
    1795 #else
    1796                         value *= value;
    1797 #endif
    1798                          // store square in list
    1799                          if (infeas[iSequence])
    1800                               infeas[iSequence] = value; // already there
    1801                          else
    1802                               infeasible_->quickAdd(iSequence, value);
    1803                     } else {
    1804                          infeasible_->zero(iSequence);
    1805                     }
    1806                }
    1807           }
    1808      }
    1809      // we can zero out as will have to get pivot row
    1810      // ***** move
    1811      updates->setNumElements(0);
    1812      spareColumn1->setNumElements(0);
    1813      if (pivotRow >= 0) {
    1814           // make sure infeasibility on incoming is 0.0
    1815           int sequenceIn = model_->sequenceIn();
    1816           infeasible_->zero(sequenceIn);
    1817      }
    1818      // for weights update we use pivotSequence
    1819      pivotRow = pivotSequence_;
    1820      // unset in case sub flip
    1821      pivotSequence_ = -1;
    1822      if (pivotRow >= 0) {
    1823           // make sure infeasibility on incoming is 0.0
    1824           const int * pivotVariable = model_->pivotVariable();
    1825           int sequenceIn = pivotVariable[pivotRow];
    1826           assert (sequenceIn == model_->sequenceIn());
    1827           infeasible_->zero(sequenceIn);
    1828           // and we can see if reference
    1829           double referenceIn;
    1830           if (mode_ != 1) {
    1831                if(reference(sequenceIn))
    1832                     referenceIn = 1.0;
    1833                else
    1834                     referenceIn = 0.0;
    1835           } else {
    1836                referenceIn = -1.0;
    1837           }
    1838           // save outgoing weight round update
    1839           double outgoingWeight = 0.0;
    1840           int sequenceOut = model_->sequenceOut();
    1841           if (sequenceOut >= 0)
    1842                outgoingWeight = weights_[sequenceOut];
    1843           // update weights
    1844           updates->setNumElements(0);
    1845           spareColumn1->setNumElements(0);
    1846           // might as well set dj to 1
    1847           dj = -1.0;
    1848           updates->createPacked(1, &pivotRow, &dj);
    1849           model_->factorization()->updateColumnTranspose(spareRow2, updates);
    1850           bool needSubset =  (mode_ < 4 || numberSwitched_ > 1 || mode_ >= 10);
    1851 
    1852           double * weight;
    1853           double * other = alternateWeights_->denseVector();
    1854           int numberColumns = model_->numberColumns();
    1855           // rows
    1856           number = updates->getNumElements();
    1857           index = updates->getIndices();
    1858           updateBy = updates->denseVector();
    1859           weight = weights_ + numberColumns;
    1860           if (needSubset) {
    1861 #if ALT_UPDATE_WEIGHTS !=2
    1862             model_->factorization()->updateColumnTranspose(spareRow2,
    1863                alternateWeights_);
    1864 #elif ALT_UPDATE_WEIGHTS==1
    1865      if (altVector[1]) {
    1866        int numberRows=model_->numberRows();
    1867        double * work1 = altVector[1]->denseVector();
    1868        double * worka = alternateWeights_->denseVector();
    1869        int iRow=-1;
    1870        double diff=1.0e-8;
    1871        for (int i=0;i<numberRows;i++) {
    1872          double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
    1873          double d=fabs(work1[i]-worka[i]);
    1874          if (dd>1.0e-6&&d>diff*dd) {
    1875            diff=d/dd;
    1876            iRow=i;
    1877          }
    1878        }
    1879        if (iRow>=0)
    1880          printf("largest2 difference of %g (%g,%g) on row %d\n",
    1881                 diff,work1[iRow],worka[iRow],iRow);
    1882      }
    1883 #endif
    1884                // do alternateWeights_ here so can scale
    1885                for (j = 0; j < number; j++) {
    1886                     int iSequence = index[j];
    1887                     assert (iSequence >= 0 && iSequence < model_->numberRows());
    1888                     double thisWeight = weight[iSequence];
    1889                     // row has -1
    1890                     double pivot = - updateBy[j];
    1891                     double modification = other[iSequence];
    1892                     double pivotSquared = pivot * pivot;
    1893 
    1894                     thisWeight += pivotSquared * devex_ + pivot * modification;
    1895                     if (thisWeight < TRY_NORM) {
    1896                          if (mode_ == 1) {
    1897                               // steepest
    1898                               thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    1899                          } else {
    1900                               // exact
    1901                               thisWeight = referenceIn * pivotSquared;
    1902                               if (reference(iSequence + numberColumns))
    1903                                    thisWeight += 1.0;
    1904                               thisWeight = CoinMax(thisWeight, TRY_NORM);
    1905                          }
    1906                     }
    1907                     weight[iSequence] = thisWeight;
    1908                }
    1909                transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2, 0.0);
    1910           } else {
    1911                // put row of tableau in rowArray and columnArray
    1912                model_->clpMatrix()->transposeTimes(model_, -1.0,
    1913                                                    updates, spareColumn2, spareColumn1);
    1914           }
    1915 
    1916           if (needSubset) {
    1917                CoinZeroN(updateBy, number);
    1918           } else if (mode_ == 4) {
    1919                // Devex
    1920                assert (devex_ > 0.0);
    1921                for (j = 0; j < number; j++) {
    1922                     int iSequence = index[j];
    1923                     double thisWeight = weight[iSequence];
    1924                     // row has -1
    1925                     double pivot = -updateBy[j];
    1926                     updateBy[j] = 0.0;
    1927                     double value = pivot * pivot * devex_;
    1928                     if (reference(iSequence + numberColumns))
    1929                          value += 1.0;
    1930                     weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1931                }
    1932           }
    1933 
    1934           // columns
    1935           weight = weights_;
    1936 
    1937           number = spareColumn1->getNumElements();
    1938           index = spareColumn1->getIndices();
    1939           updateBy = spareColumn1->denseVector();
    1940           if (needSubset) {
    1941                // Exact - already done
    1942           } else if (mode_ == 4) {
    1943                // Devex
    1944                for (j = 0; j < number; j++) {
    1945                     int iSequence = index[j];
    1946                     double thisWeight = weight[iSequence];
    1947                     // row has -1
    1948                     double pivot = updateBy[j];
    1949                     updateBy[j] = 0.0;
    1950                     double value = pivot * pivot * devex_;
    1951                     if (reference(iSequence))
    1952                          value += 1.0;
    1953                     weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    1954                }
    1955           }
    1956           // restore outgoing weight
    1957           if (sequenceOut >= 0)
    1958                weights_[sequenceOut] = outgoingWeight;
    1959           alternateWeights_->clear();
    1960           spareColumn2->setNumElements(0);
    1961           //#define SOME_DEBUG_1
     1898        }
     1899        weight[iSequence] = thisWeight;
     1900      }
     1901      transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2, 0.0);
     1902    } else {
     1903      // put row of tableau in rowArray and columnArray
     1904      model_->clpMatrix()->transposeTimes(model_, -1.0,
     1905        updates, spareColumn2, spareColumn1);
     1906    }
     1907
     1908    if (needSubset) {
     1909      CoinZeroN(updateBy, number);
     1910    } else if (mode_ == 4) {
     1911      // Devex
     1912      assert(devex_ > 0.0);
     1913      for (j = 0; j < number; j++) {
     1914        int iSequence = index[j];
     1915        double thisWeight = weight[iSequence];
     1916        // row has -1
     1917        double pivot = -updateBy[j];
     1918        updateBy[j] = 0.0;
     1919        double value = pivot * pivot * devex_;
     1920        if (reference(iSequence + numberColumns))
     1921          value += 1.0;
     1922        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1923      }
     1924    }
     1925
     1926    // columns
     1927    weight = weights_;
     1928
     1929    number = spareColumn1->getNumElements();
     1930    index = spareColumn1->getIndices();
     1931    updateBy = spareColumn1->denseVector();
     1932    if (needSubset) {
     1933      // Exact - already done
     1934    } else if (mode_ == 4) {
     1935      // Devex
     1936      for (j = 0; j < number; j++) {
     1937        int iSequence = index[j];
     1938        double thisWeight = weight[iSequence];
     1939        // row has -1
     1940        double pivot = updateBy[j];
     1941        updateBy[j] = 0.0;
     1942        double value = pivot * pivot * devex_;
     1943        if (reference(iSequence))
     1944          value += 1.0;
     1945        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     1946      }
     1947    }
     1948    // restore outgoing weight
     1949    if (sequenceOut >= 0)
     1950      weights_[sequenceOut] = outgoingWeight;
     1951    alternateWeights_->clear();
     1952    spareColumn2->setNumElements(0);
     1953    //#define SOME_DEBUG_1
    19621954#ifdef SOME_DEBUG_1
    1963           // check for accuracy
    1964           int iCheck = 892;
    1965           //printf("weight for iCheck is %g\n",weights_[iCheck]);
    1966           int numberRows = model_->numberRows();
    1967           //int numberColumns = model_->numberColumns();
    1968           for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    1969                if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    1970                          model_->getStatus(iCheck) != ClpSimplex::isFixed)
    1971                     checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    1972           }
    1973 #endif
    1974      }
    1975      updates->setNumElements(0);
    1976      spareColumn1->setNumElements(0);
     1955    // check for accuracy
     1956    int iCheck = 892;
     1957    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     1958    int numberRows = model_->numberRows();
     1959    //int numberColumns = model_->numberColumns();
     1960    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     1961      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     1962        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     1963    }
     1964#endif
     1965  }
     1966  updates->setNumElements(0);
     1967  spareColumn1->setNumElements(0);
    19771968}
    19781969// Updates two arrays for steepest
    1979 int
    1980 ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
    1981           const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
    1982           CoinIndexedVector * spare,
    1983           double scaleFactor)
     1970int ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
     1971  const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
     1972  CoinIndexedVector *spare,
     1973  double scaleFactor)
    19841974{
    1985      // see if reference
    1986      int sequenceIn = model_->sequenceIn();
    1987      double referenceIn;
    1988      if (mode_ != 1) {
    1989           if(reference(sequenceIn))
    1990                referenceIn = 1.0;
    1991           else
    1992                referenceIn = 0.0;
    1993      } else {
    1994           referenceIn = -1.0;
    1995      }
    1996      int returnCode=0;
    1997      if (model_->clpMatrix()->canCombine(model_, pi1)) {
    1998        double * infeas = scaleFactor ? infeasible_->denseVector() : NULL;
    1999        // put row of tableau in rowArray and columnArray
    2000        returnCode=model_->clpMatrix()->transposeTimes2(model_, pi1,
    2001                                                        dj1, pi2, spare,
    2002                                                        infeas,
    2003                                                        model_->djRegion(1),
    2004                                                referenceIn, devex_,
    2005                                                reference_,
    2006                                                weights_, scaleFactor);
    2007        if (model_->spareIntArray_[3]>-2)
    2008          returnCode=2;
    2009      } else {
    2010           // put row of tableau in rowArray and columnArray
    2011           model_->clpMatrix()->transposeTimes(model_, -1.0,
    2012                                               pi1, dj2, dj1);
    2013           // get subset which have nonzero tableau elements
    2014           model_->clpMatrix()->subsetTransposeTimes(model_, pi2, dj1, dj2);
    2015           bool killDjs = (scaleFactor == 0.0);
    2016           if (!scaleFactor)
    2017                scaleFactor = 1.0;
    2018           // columns
    2019           double * weight = weights_;
    2020 
    2021           int number = dj1->getNumElements();
    2022           const int * index = dj1->getIndices();
    2023           double * updateBy = dj1->denseVector();
    2024           double * updateBy2 = dj2->denseVector();
    2025 
    2026           for (int j = 0; j < number; j++) {
    2027                double thisWeight;
    2028                double pivot;
    2029                double pivotSquared;
    2030                int iSequence = index[j];
    2031                double value2 = updateBy[j];
    2032                if (killDjs)
    2033                     updateBy[j] = 0.0;
    2034                double modification = updateBy2[j];
    2035                updateBy2[j] = 0.0;
    2036                ClpSimplex::Status status = model_->getStatus(iSequence);
    2037 
    2038                if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
    2039                     thisWeight = weight[iSequence];
    2040                     pivot = value2 * scaleFactor;
    2041                     pivotSquared = pivot * pivot;
    2042 
    2043                     thisWeight += pivotSquared * devex_ + pivot * modification;
    2044                     if (thisWeight < TRY_NORM) {
    2045                          if (referenceIn < 0.0) {
    2046                               // steepest
    2047                               thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    2048                          } else {
    2049                               // exact
    2050                               thisWeight = referenceIn * pivotSquared;
    2051                               if (reference(iSequence))
    2052                                    thisWeight += 1.0;
    2053                               thisWeight = CoinMax(thisWeight, TRY_NORM);
    2054                          }
    2055                     }
    2056                     weight[iSequence] = thisWeight;
    2057                }
     1975  // see if reference
     1976  int sequenceIn = model_->sequenceIn();
     1977  double referenceIn;
     1978  if (mode_ != 1) {
     1979    if (reference(sequenceIn))
     1980      referenceIn = 1.0;
     1981    else
     1982      referenceIn = 0.0;
     1983  } else {
     1984    referenceIn = -1.0;
     1985  }
     1986  int returnCode = 0;
     1987  if (model_->clpMatrix()->canCombine(model_, pi1)) {
     1988    double *infeas = scaleFactor ? infeasible_->denseVector() : NULL;
     1989    // put row of tableau in rowArray and columnArray
     1990    returnCode = model_->clpMatrix()->transposeTimes2(model_, pi1,
     1991      dj1, pi2, spare,
     1992      infeas,
     1993      model_->djRegion(1),
     1994      referenceIn, devex_,
     1995      reference_,
     1996      weights_, scaleFactor);
     1997    if (model_->spareIntArray_[3] > -2)
     1998      returnCode = 2;
     1999  } else {
     2000    // put row of tableau in rowArray and columnArray
     2001    model_->clpMatrix()->transposeTimes(model_, -1.0,
     2002      pi1, dj2, dj1);
     2003    // get subset which have nonzero tableau elements
     2004    model_->clpMatrix()->subsetTransposeTimes(model_, pi2, dj1, dj2);
     2005    bool killDjs = (scaleFactor == 0.0);
     2006    if (!scaleFactor)
     2007      scaleFactor = 1.0;
     2008    // columns
     2009    double *weight = weights_;
     2010
     2011    int number = dj1->getNumElements();
     2012    const int *index = dj1->getIndices();
     2013    double *updateBy = dj1->denseVector();
     2014    double *updateBy2 = dj2->denseVector();
     2015
     2016    for (int j = 0; j < number; j++) {
     2017      double thisWeight;
     2018      double pivot;
     2019      double pivotSquared;
     2020      int iSequence = index[j];
     2021      double value2 = updateBy[j];
     2022      if (killDjs)
     2023        updateBy[j] = 0.0;
     2024      double modification = updateBy2[j];
     2025      updateBy2[j] = 0.0;
     2026      ClpSimplex::Status status = model_->getStatus(iSequence);
     2027
     2028      if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
     2029        thisWeight = weight[iSequence];
     2030        pivot = value2 * scaleFactor;
     2031        pivotSquared = pivot * pivot;
     2032
     2033        thisWeight += pivotSquared * devex_ + pivot * modification;
     2034        if (thisWeight < TRY_NORM) {
     2035          if (referenceIn < 0.0) {
     2036            // steepest
     2037            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     2038          } else {
     2039            // exact
     2040            thisWeight = referenceIn * pivotSquared;
     2041            if (reference(iSequence))
     2042              thisWeight += 1.0;
     2043            thisWeight = CoinMax(thisWeight, TRY_NORM);
    20582044          }
    2059      }
    2060      dj2->setNumElements(0);
    2061      return returnCode;
     2045        }
     2046        weight[iSequence] = thisWeight;
     2047      }
     2048    }
     2049  }
     2050  dj2->setNumElements(0);
     2051  return returnCode;
    20622052}
    20632053// Update weights for Devex
    2064 void
    2065 ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
    2066                                    CoinIndexedVector * spareRow2,
    2067                                    CoinIndexedVector * spareColumn1,
    2068                                    CoinIndexedVector * spareColumn2)
     2054void ClpPrimalColumnSteepest::justDevex(CoinIndexedVector *updates,
     2055  CoinIndexedVector *spareRow2,
     2056  CoinIndexedVector *spareColumn1,
     2057  CoinIndexedVector *spareColumn2)
    20692058{
    2070      int j;
    2071      int number = 0;
    2072      int * index;
    2073      double * updateBy;
    2074      // dj could be very small (or even zero - take care)
    2075      double dj = model_->dualIn();
    2076      double tolerance = model_->currentDualTolerance();
    2077      // we can't really trust infeasibilities if there is dual error
    2078      // this coding has to mimic coding in checkDualSolution
    2079      double error = CoinMin(1.0e-2, model_->largestDualError());
    2080      // allow tolerance at least slightly bigger than standard
    2081      tolerance = tolerance + error;
    2082      int pivotRow = model_->pivotRow();
    2083      // for weights update we use pivotSequence
    2084      pivotRow = pivotSequence_;
    2085      assert (pivotRow >= 0);
    2086      // make sure infeasibility on incoming is 0.0
    2087      const int * pivotVariable = model_->pivotVariable();
    2088      int sequenceIn = pivotVariable[pivotRow];
    2089      infeasible_->zero(sequenceIn);
    2090      // and we can see if reference
    2091      //double referenceIn = 0.0;
    2092      //if (mode_ != 1 && reference(sequenceIn))
    2093      //   referenceIn = 1.0;
    2094      // save outgoing weight round update
    2095      double outgoingWeight = 0.0;
    2096      int sequenceOut = model_->sequenceOut();
    2097      if (sequenceOut >= 0)
    2098           outgoingWeight = weights_[sequenceOut];
    2099      assert (!updates->getNumElements());
    2100      assert (!spareColumn1->getNumElements());
    2101      // unset in case sub flip
    2102      pivotSequence_ = -1;
    2103      // might as well set dj to 1
    2104      dj = -1.0;
    2105      updates->createPacked(1, &pivotRow, &dj);
    2106      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    2107      // put row of tableau in rowArray and columnArray
    2108      model_->clpMatrix()->transposeTimes(model_, -1.0,
    2109                                          updates, spareColumn2, spareColumn1);
    2110      double * weight;
    2111      int numberColumns = model_->numberColumns();
    2112      // rows
    2113      number = updates->getNumElements();
    2114      index = updates->getIndices();
    2115      updateBy = updates->denseVector();
    2116      weight = weights_ + numberColumns;
    2117 
    2118      // Devex
    2119      assert (devex_ > 0.0);
    2120      for (j = 0; j < number; j++) {
    2121           int iSequence = index[j];
    2122           double thisWeight = weight[iSequence];
    2123           // row has -1
    2124           double pivot = - updateBy[j];
    2125           updateBy[j] = 0.0;
    2126           double value = pivot * pivot * devex_;
    2127           if (reference(iSequence + numberColumns))
    2128                value += 1.0;
    2129           weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    2130      }
    2131 
    2132      // columns
    2133      weight = weights_;
    2134 
    2135      number = spareColumn1->getNumElements();
    2136      index = spareColumn1->getIndices();
    2137      updateBy = spareColumn1->denseVector();
    2138      // Devex
    2139      for (j = 0; j < number; j++) {
    2140           int iSequence = index[j];
    2141           double thisWeight = weight[iSequence];
    2142           // row has -1
    2143           double pivot = updateBy[j];
    2144           updateBy[j] = 0.0;
    2145           double value = pivot * pivot * devex_;
    2146           if (reference(iSequence))
    2147                value += 1.0;
    2148           weight[iSequence] = CoinMax(0.99 * thisWeight, value);
    2149      }
    2150      // restore outgoing weight
    2151      if (sequenceOut >= 0)
    2152           weights_[sequenceOut] = outgoingWeight;
    2153      spareColumn2->setNumElements(0);
    2154      //#define SOME_DEBUG_1
     2059  int j;
     2060  int number = 0;
     2061  int *index;
     2062  double *updateBy;
     2063  // dj could be very small (or even zero - take care)
     2064  double dj = model_->dualIn();
     2065  double tolerance = model_->currentDualTolerance();
     2066  // we can't really trust infeasibilities if there is dual error
     2067  // this coding has to mimic coding in checkDualSolution
     2068  double error = CoinMin(1.0e-2, model_->largestDualError());
     2069  // allow tolerance at least slightly bigger than standard
     2070  tolerance = tolerance + error;
     2071  int pivotRow = model_->pivotRow();
     2072  // for weights update we use pivotSequence
     2073  pivotRow = pivotSequence_;
     2074  assert(pivotRow >= 0);
     2075  // make sure infeasibility on incoming is 0.0
     2076  const int *pivotVariable = model_->pivotVariable();
     2077  int sequenceIn = pivotVariable[pivotRow];
     2078  infeasible_->zero(sequenceIn);
     2079  // and we can see if reference
     2080  //double referenceIn = 0.0;
     2081  //if (mode_ != 1 && reference(sequenceIn))
     2082  //   referenceIn = 1.0;
     2083  // save outgoing weight round update
     2084  double outgoingWeight = 0.0;
     2085  int sequenceOut = model_->sequenceOut();
     2086  if (sequenceOut >= 0)
     2087    outgoingWeight = weights_[sequenceOut];
     2088  assert(!updates->getNumElements());
     2089  assert(!spareColumn1->getNumElements());
     2090  // unset in case sub flip
     2091  pivotSequence_ = -1;
     2092  // might as well set dj to 1
     2093  dj = -1.0;
     2094  updates->createPacked(1, &pivotRow, &dj);
     2095  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     2096  // put row of tableau in rowArray and columnArray
     2097  model_->clpMatrix()->transposeTimes(model_, -1.0,
     2098    updates, spareColumn2, spareColumn1);
     2099  double *weight;
     2100  int numberColumns = model_->numberColumns();
     2101  // rows
     2102  number = updates->getNumElements();
     2103  index = updates->getIndices();
     2104  updateBy = updates->denseVector();
     2105  weight = weights_ + numberColumns;
     2106
     2107  // Devex
     2108  assert(devex_ > 0.0);
     2109  for (j = 0; j < number; j++) {
     2110    int iSequence = index[j];
     2111    double thisWeight = weight[iSequence];
     2112    // row has -1
     2113    double pivot = -updateBy[j];
     2114    updateBy[j] = 0.0;
     2115    double value = pivot * pivot * devex_;
     2116    if (reference(iSequence + numberColumns))
     2117      value += 1.0;
     2118    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     2119  }
     2120
     2121  // columns
     2122  weight = weights_;
     2123
     2124  number = spareColumn1->getNumElements();
     2125  index = spareColumn1->getIndices();
     2126  updateBy = spareColumn1->denseVector();
     2127  // Devex
     2128  for (j = 0; j < number; j++) {
     2129    int iSequence = index[j];
     2130    double thisWeight = weight[iSequence];
     2131    // row has -1
     2132    double pivot = updateBy[j];
     2133    updateBy[j] = 0.0;
     2134    double value = pivot * pivot * devex_;
     2135    if (reference(iSequence))
     2136      value += 1.0;
     2137    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     2138  }
     2139  // restore outgoing weight
     2140  if (sequenceOut >= 0)
     2141    weights_[sequenceOut] = outgoingWeight;
     2142  spareColumn2->setNumElements(0);
     2143  //#define SOME_DEBUG_1
    21552144#ifdef SOME_DEBUG_1
    2156      // check for accuracy
    2157      int iCheck = 892;
    2158      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    2159      int numberRows = model_->numberRows();
    2160      //int numberColumns = model_->numberColumns();
    2161      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    2162           if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    2163                     model_->getStatus(iCheck) != ClpSimplex::isFixed)
    2164                checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    2165      }
    2166 #endif
    2167      updates->setNumElements(0);
    2168      spareColumn1->setNumElements(0);
     2145  // check for accuracy
     2146  int iCheck = 892;
     2147  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     2148  int numberRows = model_->numberRows();
     2149  //int numberColumns = model_->numberColumns();
     2150  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     2151    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     2152      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     2153  }
     2154#endif
     2155  updates->setNumElements(0);
     2156  spareColumn1->setNumElements(0);
    21692157}
    21702158// Update weights for Steepest
    2171 void
    2172 ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates,
    2173                                       CoinIndexedVector * spareRow2,
    2174                                       CoinIndexedVector * spareColumn1,
    2175                                       CoinIndexedVector * spareColumn2)
     2159void ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector *updates,
     2160  CoinIndexedVector *spareRow2,
     2161  CoinIndexedVector *spareColumn1,
     2162  CoinIndexedVector *spareColumn2)
    21762163{
    2177      int j;
    2178      int number = 0;
    2179      int * index;
    2180      double * updateBy;
    2181      // dj could be very small (or even zero - take care)
    2182      double dj = model_->dualIn();
    2183      double tolerance = model_->currentDualTolerance();
    2184      // we can't really trust infeasibilities if there is dual error
    2185      // this coding has to mimic coding in checkDualSolution
    2186      double error = CoinMin(1.0e-2, model_->largestDualError());
    2187      // allow tolerance at least slightly bigger than standard
    2188      tolerance = tolerance + error;
    2189      int pivotRow = model_->pivotRow();
    2190      // for weights update we use pivotSequence
    2191      pivotRow = pivotSequence_;
    2192      // unset in case sub flip
    2193      pivotSequence_ = -1;
    2194      assert (pivotRow >= 0);
    2195      // make sure infeasibility on incoming is 0.0
    2196      const int * pivotVariable = model_->pivotVariable();
    2197      int sequenceIn = pivotVariable[pivotRow];
    2198      infeasible_->zero(sequenceIn);
    2199      // and we can see if reference
    2200      double referenceIn = 0.0;
    2201      if (mode_ != 1 && reference(sequenceIn))
    2202           referenceIn = 1.0;
    2203      // save outgoing weight round update
    2204      double outgoingWeight = 0.0;
    2205      int sequenceOut = model_->sequenceOut();
    2206      if (sequenceOut >= 0)
    2207           outgoingWeight = weights_[sequenceOut];
    2208      assert (!updates->getNumElements());
    2209      assert (!spareColumn1->getNumElements());
    2210      // update weights
    2211      //updates->setNumElements(0);
    2212      //spareColumn1->setNumElements(0);
    2213      // might as well set dj to 1
    2214      dj = -1.0;
    2215      updates->createPacked(1, &pivotRow, &dj);
    2216      model_->factorization()->updateColumnTranspose(spareRow2, updates);
    2217      // put row of tableau in rowArray and columnArray
    2218      model_->clpMatrix()->transposeTimes(model_, -1.0,
    2219                                          updates, spareColumn2, spareColumn1);
    2220      double * weight;
    2221      double * other = alternateWeights_->denseVector();
    2222      int numberColumns = model_->numberColumns();
    2223      // rows
    2224      number = updates->getNumElements();
    2225      index = updates->getIndices();
    2226      updateBy = updates->denseVector();
    2227      weight = weights_ + numberColumns;
    2228 
    2229      // Exact
    2230      // now update weight update array
    2231      //alternateWeights_->scanAndPack();
    2232 #if ALT_UPDATE_WEIGHTS !=2
    2233      model_->factorization()->updateColumnTranspose(spareRow2,
    2234                                                     alternateWeights_);
    2235 #elif ALT_UPDATE_WEIGHTS==1
    2236      if (altVector[1]) {
    2237        int numberRows=model_->numberRows();
    2238        double * work1 = altVector[1]->denseVector();
    2239        double * worka = alternateWeights_->denseVector();
    2240        int iRow=-1;
    2241        double diff=1.0e-8;
    2242        for (int i=0;i<numberRows;i++) {
    2243          double dd=CoinMax(fabs(work1[i]),fabs(worka[i]));
    2244          double d=fabs(work1[i]-worka[i]);
    2245          if (dd>1.0e-6&&d>diff*dd) {
    2246            diff=d/dd;
    2247            iRow=i;
    2248         }
    2249        }
    2250        if (iRow>=0)
    2251         printf("largest3 difference of %g (%g,%g) on row %d\n",
    2252                 diff,work1[iRow],worka[iRow],iRow);
    2253      }
    2254 #endif
    2255      // get subset which have nonzero tableau elements
    2256      model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
    2257                spareColumn1,
    2258                spareColumn2);
    2259      for (j = 0; j < number; j++) {
    2260           int iSequence = index[j];
    2261           double thisWeight = weight[iSequence];
    2262           // row has -1
    2263           double pivot = -updateBy[j];
    2264           updateBy[j] = 0.0;
    2265           double modification = other[iSequence];
    2266           double pivotSquared = pivot * pivot;
    2267 
    2268           thisWeight += pivotSquared * devex_ + pivot * modification;
    2269           if (thisWeight < TRY_NORM) {
    2270                if (mode_ == 1) {
    2271                     // steepest
    2272                     thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    2273                } else {
    2274                     // exact
    2275                     thisWeight = referenceIn * pivotSquared;
    2276                     if (reference(iSequence + numberColumns))
    2277                          thisWeight += 1.0;
    2278                     thisWeight = CoinMax(thisWeight, TRY_NORM);
    2279                }
    2280           }
    2281           weight[iSequence] = thisWeight;
    2282      }
    2283 
    2284      // columns
    2285      weight = weights_;
    2286 
    2287      number = spareColumn1->getNumElements();
    2288      index = spareColumn1->getIndices();
    2289      updateBy = spareColumn1->denseVector();
    2290      // Exact
    2291      double * updateBy2 = spareColumn2->denseVector();
    2292      for (j = 0; j < number; j++) {
    2293           int iSequence = index[j];
    2294           double thisWeight = weight[iSequence];
    2295           double pivot = updateBy[j];
    2296           updateBy[j] = 0.0;
    2297           double modification = updateBy2[j];
    2298           updateBy2[j] = 0.0;
    2299           double pivotSquared = pivot * pivot;
    2300 
    2301           thisWeight += pivotSquared * devex_ + pivot * modification;
    2302           if (thisWeight < TRY_NORM) {
    2303                if (mode_ == 1) {
    2304                     // steepest
    2305                     thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    2306                } else {
    2307                     // exact
    2308                     thisWeight = referenceIn * pivotSquared;
    2309                     if (reference(iSequence))
    2310                          thisWeight += 1.0;
    2311                     thisWeight = CoinMax(thisWeight, TRY_NORM);
    2312                }
    2313           }
    2314           weight[iSequence] = thisWeight;
    2315      }
    2316      // restore outgoing weight
    2317      if (sequenceOut >= 0)
    2318           weights_[sequenceOut] = outgoingWeight;
    2319      alternateWeights_->clear();
    2320      spareColumn2->setNumElements(0);
    2321      //#define SOME_DEBUG_1
     2164  int j;
     2165  int number = 0;
     2166  int *index;
     2167  double *updateBy;
     2168  // dj could be very small (or even zero - take care)
     2169  double dj = model_->dualIn();
     2170  double tolerance = model_->currentDualTolerance();
     2171  // we can't really trust infeasibilities if there is dual error
     2172  // this coding has to mimic coding in checkDualSolution
     2173  double error = CoinMin(1.0e-2, model_->largestDualError());
     2174  // allow tolerance at least slightly bigger than standard
     2175  tolerance = tolerance + error;
     2176  int pivotRow = model_->pivotRow();
     2177  // for weights update we use pivotSequence
     2178  pivotRow = pivotSequence_;
     2179  // unset in case sub flip
     2180  pivotSequence_ = -1;
     2181  assert(pivotRow >= 0);
     2182  // make sure infeasibility on incoming is 0.0
     2183  const int *pivotVariable = model_->pivotVariable();
     2184  int sequenceIn = pivotVariable[pivotRow];
     2185  infeasible_->zero(sequenceIn);
     2186  // and we can see if reference
     2187  double referenceIn = 0.0;
     2188  if (mode_ != 1 && reference(sequenceIn))
     2189    referenceIn = 1.0;
     2190  // save outgoing weight round update
     2191  double outgoingWeight = 0.0;
     2192  int sequenceOut = model_->sequenceOut();
     2193  if (sequenceOut >= 0)
     2194    outgoingWeight = weights_[sequenceOut];
     2195  assert(!updates->getNumElements());
     2196  assert(!spareColumn1->getNumElements());
     2197  // update weights
     2198  //updates->setNumElements(0);
     2199  //spareColumn1->setNumElements(0);
     2200  // might as well set dj to 1
     2201  dj = -1.0;
     2202  updates->createPacked(1, &pivotRow, &dj);
     2203  model_->factorization()->updateColumnTranspose(spareRow2, updates);
     2204  // put row of tableau in rowArray and columnArray
     2205  model_->clpMatrix()->transposeTimes(model_, -1.0,
     2206    updates, spareColumn2, spareColumn1);
     2207  double *weight;
     2208  double *other = alternateWeights_->denseVector();
     2209  int numberColumns = model_->numberColumns();
     2210  // rows
     2211  number = updates->getNumElements();
     2212  index = updates->getIndices();
     2213  updateBy = updates->denseVector();
     2214  weight = weights_ + numberColumns;
     2215
     2216  // Exact
     2217  // now update weight update array
     2218  //alternateWeights_->scanAndPack();
     2219#if ALT_UPDATE_WEIGHTS != 2
     2220  model_->factorization()->updateColumnTranspose(spareRow2,
     2221    alternateWeights_);
     2222#elif ALT_UPDATE_WEIGHTS == 1
     2223  if (altVector[1]) {
     2224    int numberRows = model_->numberRows();
     2225    double *work1 = altVector[1]->denseVector();
     2226    double *worka = alternateWeights_->denseVector();
     2227    int iRow = -1;
     2228    double diff = 1.0e-8;
     2229    for (int i = 0; i < numberRows; i++) {
     2230      double dd = CoinMax(fabs(work1[i]), fabs(worka[i]));
     2231      double d = fabs(work1[i] - worka[i]);
     2232      if (dd > 1.0e-6 && d > diff * dd) {
     2233        diff = d / dd;
     2234        iRow = i;
     2235      }
     2236    }
     2237    if (iRow >= 0)
     2238      printf("largest3 difference of %g (%g,%g) on row %d\n",
     2239        diff, work1[iRow], worka[iRow], iRow);
     2240  }
     2241#endif
     2242  // get subset which have nonzero tableau elements
     2243  model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
     2244    spareColumn1,
     2245    spareColumn2);
     2246  for (j = 0; j < number; j++) {
     2247    int iSequence = index[j];
     2248    double thisWeight = weight[iSequence];
     2249    // row has -1
     2250    double pivot = -updateBy[j];
     2251    updateBy[j] = 0.0;
     2252    double modification = other[iSequence];
     2253    double pivotSquared = pivot * pivot;
     2254
     2255    thisWeight += pivotSquared * devex_ + pivot * modification;
     2256    if (thisWeight < TRY_NORM) {
     2257      if (mode_ == 1) {
     2258        // steepest
     2259        thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     2260      } else {
     2261        // exact
     2262        thisWeight = referenceIn * pivotSquared;
     2263        if (reference(iSequence + numberColumns))
     2264          thisWeight += 1.0;
     2265        thisWeight = CoinMax(thisWeight, TRY_NORM);
     2266      }
     2267    }
     2268    weight[iSequence] = thisWeight;
     2269  }
     2270
     2271  // columns
     2272  weight = weights_;
     2273
     2274  number = spareColumn1->getNumElements();
     2275  index = spareColumn1->getIndices();
     2276  updateBy = spareColumn1->denseVector();
     2277  // Exact
     2278  double *updateBy2 = spareColumn2->denseVector();
     2279  for (j = 0; j < number; j++) {
     2280    int iSequence = index[j];
     2281    double thisWeight = weight[iSequence];
     2282    double pivot = updateBy[j];
     2283    updateBy[j] = 0.0;
     2284    double modification = updateBy2[j];
     2285    updateBy2[j] = 0.0;
     2286    double pivotSquared = pivot * pivot;
     2287
     2288    thisWeight += pivotSquared * devex_ + pivot * modification;
     2289    if (thisWeight < TRY_NORM) {
     2290      if (mode_ == 1) {
     2291        // steepest
     2292        thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
     2293      } else {
     2294        // exact
     2295        thisWeight = referenceIn * pivotSquared;
     2296        if (reference(iSequence))
     2297          thisWeight += 1.0;
     2298        thisWeight = CoinMax(thisWeight, TRY_NORM);
     2299      }
     2300    }
     2301    weight[iSequence] = thisWeight;
     2302  }
     2303  // restore outgoing weight
     2304  if (sequenceOut >= 0)
     2305    weights_[sequenceOut] = outgoingWeight;
     2306  alternateWeights_->clear();
     2307  spareColumn2->setNumElements(0);
     2308  //#define SOME_DEBUG_1
    23222309#ifdef SOME_DEBUG_1
    2323      // check for accuracy
    2324      int iCheck = 892;
    2325      //printf("weight for iCheck is %g\n",weights_[iCheck]);
    2326      int numberRows = model_->numberRows();
    2327      //int numberColumns = model_->numberColumns();
    2328      for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
    2329           if (model_->getStatus(iCheck) != ClpSimplex::basic &&
    2330                     model_->getStatus(iCheck) != ClpSimplex::isFixed)
    2331                checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
    2332      }
    2333 #endif
    2334      updates->setNumElements(0);
    2335      spareColumn1->setNumElements(0);
     2310  // check for accuracy
     2311  int iCheck = 892;
     2312  //printf("weight for iCheck is %g\n",weights_[iCheck]);
     2313  int numberRows = model_->numberRows();
     2314  //int numberColumns = model_->numberColumns();
     2315  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     2316    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     2317      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     2318  }
     2319#endif
     2320  updates->setNumElements(0);
     2321  spareColumn1->setNumElements(0);
    23362322}
    23372323// Returns pivot column, -1 if none
    2338 int
    2339 ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates,
    2340           CoinIndexedVector * ,
    2341           CoinIndexedVector * spareRow2,
    2342           CoinIndexedVector * spareColumn1,
    2343           CoinIndexedVector * spareColumn2)
     2324int ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector *updates,
     2325  CoinIndexedVector *,
     2326  CoinIndexedVector *spareRow2,
     2327  CoinIndexedVector *spareColumn1,
     2328  CoinIndexedVector *spareColumn2)
    23442329{
    2345      assert(model_);
    2346      int iSection, j;
    2347      int number = 0;
    2348      int * index;
    2349      double * updateBy;
    2350      double * reducedCost;
    2351      // dj could be very small (or even zero - take care)
    2352      double dj = model_->dualIn();
    2353      double tolerance = model_->currentDualTolerance();
    2354      // we can't really trust infeasibilities if there is dual error
    2355      // this coding has to mimic coding in checkDualSolution
    2356      double error = CoinMin(1.0e-2, model_->largestDualError());
    2357      // allow tolerance at least slightly bigger than standard
    2358      tolerance = tolerance + error;
    2359      int pivotRow = model_->pivotRow();
    2360      int anyUpdates;
    2361      double * infeas = infeasible_->denseVector();
    2362 
    2363      // Local copy of mode so can decide what to do
    2364      int switchType;
    2365      if (mode_ == 4)
    2366           switchType = 5 - numberSwitched_;
    2367      else if (mode_ >= 10)
    2368           switchType = 3;
    2369      else
    2370           switchType = mode_;
    2371      /* switchType -
     2330  assert(model_);
     2331  int iSection, j;
     2332  int number = 0;
     2333  int *index;
     2334  double *updateBy;
     2335  double *reducedCost;
     2336  // dj could be very small (or even zero - take care)
     2337  double dj = model_->dualIn();
     2338  double tolerance = model_->currentDualTolerance();
     2339  // we can't really trust infeasibilities if there is dual error
     2340  // this coding has to mimic coding in checkDualSolution
     2341  double error = CoinMin(1.0e-2, model_->largestDualError());
     2342  // allow tolerance at least slightly bigger than standard
     2343  tolerance = tolerance + error;
     2344  int pivotRow = model_->pivotRow();
     2345  int anyUpdates;
     2346  double *infeas = infeasible_->denseVector();
     2347
     2348  // Local copy of mode so can decide what to do
     2349  int switchType;
     2350  if (mode_ == 4)
     2351    switchType = 5 - numberSwitched_;
     2352  else if (mode_ >= 10)
     2353    switchType = 3;
     2354  else
     2355    switchType = mode_;
     2356  /* switchType -
    23722357        0 - all exact devex
    23732358        1 - all steepest
     
    23772362        5 - dantzig
    23782363     */
    2379      if (updates->getNumElements()) {
    2380           // would have to have two goes for devex, three for steepest
    2381           anyUpdates = 2;
    2382           // add in pivot contribution
    2383           if (pivotRow >= 0)
    2384                updates->add(pivotRow, -dj);
    2385      } else if (pivotRow >= 0) {
    2386           if (fabs(dj) > 1.0e-15) {
    2387                // some dj
    2388                updates->insert(pivotRow, -dj);
    2389                if (fabs(dj) > 1.0e-6) {
    2390                     // reasonable size
    2391                     anyUpdates = 1;
    2392                } else {
    2393                     // too small
    2394                     anyUpdates = 2;
    2395                }
    2396           } else {
    2397                // zero dj
    2398                anyUpdates = -1;
    2399           }
    2400      } else if (pivotSequence_ >= 0) {
    2401           // just after re-factorization
    2402           anyUpdates = -1;
    2403      } else {
    2404           // sub flip - nothing to do
    2405           anyUpdates = 0;
    2406      }
    2407 
    2408      if (anyUpdates > 0) {
    2409           model_->factorization()->updateColumnTranspose(spareRow2, updates);
    2410 
    2411           // put row of tableau in rowArray and columnArray
    2412           model_->clpMatrix()->transposeTimes(model_, -1.0,
    2413                                               updates, spareColumn2, spareColumn1);
    2414           // normal
    2415           for (iSection = 0; iSection < 2; iSection++) {
    2416 
    2417                reducedCost = model_->djRegion(iSection);
    2418                int addSequence;
    2419 
    2420                if (!iSection) {
    2421                     number = updates->getNumElements();
    2422                     index = updates->getIndices();
    2423                     updateBy = updates->denseVector();
    2424                     addSequence = model_->numberColumns();
    2425                } else {
    2426                     number = spareColumn1->getNumElements();
    2427                     index = spareColumn1->getIndices();
    2428                     updateBy = spareColumn1->denseVector();
    2429                     addSequence = 0;
    2430                }
    2431                if (!model_->nonLinearCost()->lookBothWays()) {
    2432 
    2433                     for (j = 0; j < number; j++) {
    2434                          int iSequence = index[j];
    2435                          double value = reducedCost[iSequence];
    2436                          value -= updateBy[iSequence];
    2437                          reducedCost[iSequence] = value;
    2438                          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    2439 
    2440                          switch(status) {
    2441 
    2442                          case ClpSimplex::basic:
    2443                               infeasible_->zero(iSequence + addSequence);
    2444                          case ClpSimplex::isFixed:
    2445                               break;
    2446                          case ClpSimplex::isFree:
    2447                          case ClpSimplex::superBasic:
    2448                               if (fabs(value) > FREE_ACCEPT * tolerance) {
    2449                                    // we are going to bias towards free (but only if reasonable)
    2450                                    value *= FREE_BIAS;
    2451                                    // store square in list
    2452                                    if (infeas[iSequence+addSequence])
    2453                                         infeas[iSequence+addSequence] = value * value; // already there
    2454                                    else
    2455                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2456                               } else {
    2457                                    infeasible_->zero(iSequence + addSequence);
    2458                               }
    2459                               break;
    2460                          case ClpSimplex::atUpperBound:
    2461                               if (value > tolerance) {
    2462                                    // store square in list
    2463                                    if (infeas[iSequence+addSequence])
    2464                                         infeas[iSequence+addSequence] = value * value; // already there
    2465                                    else
    2466                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2467                               } else {
    2468                                    infeasible_->zero(iSequence + addSequence);
    2469                               }
    2470                               break;
    2471                          case ClpSimplex::atLowerBound:
    2472                               if (value < -tolerance) {
    2473                                    // store square in list
    2474                                    if (infeas[iSequence+addSequence])
    2475                                         infeas[iSequence+addSequence] = value * value; // already there
    2476                                    else
    2477                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2478                               } else {
    2479                                    infeasible_->zero(iSequence + addSequence);
    2480                               }
    2481                          }
    2482                     }
    2483                } else {
    2484                     ClpNonLinearCost * nonLinear = model_->nonLinearCost();
    2485                     // We can go up OR down
    2486                     for (j = 0; j < number; j++) {
    2487                          int iSequence = index[j];
    2488                          double value = reducedCost[iSequence];
    2489                          value -= updateBy[iSequence];
    2490                          reducedCost[iSequence] = value;
    2491                          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    2492 
    2493                          switch(status) {
    2494 
    2495                          case ClpSimplex::basic:
    2496                               infeasible_->zero(iSequence + addSequence);
    2497                          case ClpSimplex::isFixed:
    2498                               break;
    2499                          case ClpSimplex::isFree:
    2500                          case ClpSimplex::superBasic:
    2501                               if (fabs(value) > FREE_ACCEPT * tolerance) {
    2502                                    // we are going to bias towards free (but only if reasonable)
    2503                                    value *= FREE_BIAS;
    2504                                    // store square in list
    2505                                    if (infeas[iSequence+addSequence])
    2506                                         infeas[iSequence+addSequence] = value * value; // already there
    2507                                    else
    2508                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2509                               } else {
    2510                                    infeasible_->zero(iSequence + addSequence);
    2511                               }
    2512                               break;
    2513                          case ClpSimplex::atUpperBound:
    2514                               if (value > tolerance) {
    2515                                    // store square in list
    2516                                    if (infeas[iSequence+addSequence])
    2517                                         infeas[iSequence+addSequence] = value * value; // already there
    2518                                    else
    2519                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2520                               } else {
    2521                                    // look other way - change up should be negative
    2522                                    value -= nonLinear->changeUpInCost(iSequence + addSequence);
    2523                                    if (value > -tolerance) {
    2524                                         infeasible_->zero(iSequence + addSequence);
    2525                                    } else {
    2526                                         // store square in list
    2527                                         if (infeas[iSequence+addSequence])
    2528                                              infeas[iSequence+addSequence] = value * value; // already there
    2529                                         else
    2530                                              infeasible_->quickAdd(iSequence + addSequence, value * value);
    2531                                    }
    2532                               }
    2533                               break;
    2534                          case ClpSimplex::atLowerBound:
    2535                               if (value < -tolerance) {
    2536                                    // store square in list
    2537                                    if (infeas[iSequence+addSequence])
    2538                                         infeas[iSequence+addSequence] = value * value; // already there
    2539                                    else
    2540                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2541                               } else {
    2542                                    // look other way - change down should be positive
    2543                                    value -= nonLinear->changeDownInCost(iSequence + addSequence);
    2544                                    if (value < tolerance) {
    2545                                         infeasible_->zero(iSequence + addSequence);
    2546                                    } else {
    2547                                         // store square in list
    2548                                         if (infeas[iSequence+addSequence])
    2549                                              infeas[iSequence+addSequence] = value * value; // already there
    2550                                         else
    2551                                              infeasible_->quickAdd(iSequence + addSequence, value * value);
    2552                                    }
    2553                               }
    2554                          }
    2555                     }
    2556                }
    2557           }
    2558           if (anyUpdates == 2) {
    2559                // we can zero out as will have to get pivot row
    2560                updates->clear();
    2561                spareColumn1->clear();
    2562           }
    2563           if (pivotRow >= 0) {
    2564                // make sure infeasibility on incoming is 0.0
    2565                int sequenceIn = model_->sequenceIn();
    2566                infeasible_->zero(sequenceIn);
    2567           }
    2568      }
    2569      // make sure outgoing from last iteration okay
    2570      int sequenceOut = model_->sequenceOut();
    2571      if (sequenceOut >= 0) {
    2572           ClpSimplex::Status status = model_->getStatus(sequenceOut);
    2573           double value = model_->reducedCost(sequenceOut);
    2574 
    2575           switch(status) {
     2364  if (updates->getNumElements()) {
     2365    // would have to have two goes for devex, three for steepest
     2366    anyUpdates = 2;
     2367    // add in pivot contribution
     2368    if (pivotRow >= 0)
     2369      updates->add(pivotRow, -dj);
     2370  } else if (pivotRow >= 0) {
     2371    if (fabs(dj) > 1.0e-15) {
     2372      // some dj
     2373      updates->insert(pivotRow, -dj);
     2374      if (fabs(dj) > 1.0e-6) {
     2375        // reasonable size
     2376        anyUpdates = 1;
     2377      } else {
     2378        // too small
     2379        anyUpdates = 2;
     2380      }
     2381    } else {
     2382      // zero dj
     2383      anyUpdates = -1;
     2384    }
     2385  } else if (pivotSequence_ >= 0) {
     2386    // just after re-factorization
     2387    anyUpdates = -1;
     2388  } else {
     2389    // sub flip - nothing to do
     2390    anyUpdates = 0;
     2391  }
     2392
     2393  if (anyUpdates > 0) {
     2394    model_->factorization()->updateColumnTranspose(spareRow2, updates);
     2395
     2396    // put row of tableau in rowArray and columnArray
     2397    model_->clpMatrix()->transposeTimes(model_, -1.0,
     2398      updates, spareColumn2, spareColumn1);
     2399    // normal
     2400    for (iSection = 0; iSection < 2; iSection++) {
     2401
     2402      reducedCost = model_->djRegion(iSection);
     2403      int addSequence;
     2404
     2405      if (!iSection) {
     2406        number = updates->getNumElements();
     2407        index = updates->getIndices();
     2408        updateBy = updates->denseVector();
     2409        addSequence = model_->numberColumns();
     2410      } else {
     2411        number = spareColumn1->getNumElements();
     2412        index = spareColumn1->getIndices();
     2413        updateBy = spareColumn1->denseVector();
     2414        addSequence = 0;
     2415      }
     2416      if (!model_->nonLinearCost()->lookBothWays()) {
     2417
     2418        for (j = 0; j < number; j++) {
     2419          int iSequence = index[j];
     2420          double value = reducedCost[iSequence];
     2421          value -= updateBy[iSequence];
     2422          reducedCost[iSequence] = value;
     2423          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     2424
     2425          switch (status) {
    25762426
    25772427          case ClpSimplex::basic:
     2428            infeasible_->zero(iSequence + addSequence);
    25782429          case ClpSimplex::isFixed:
    2579                break;
     2430            break;
    25802431          case ClpSimplex::isFree:
    25812432          case ClpSimplex::superBasic:
    2582                if (fabs(value) > FREE_ACCEPT * tolerance) {
    2583                     // we are going to bias towards free (but only if reasonable)
    2584                     value *= FREE_BIAS;
    2585                     // store square in list
    2586                     if (infeas[sequenceOut])
    2587                          infeas[sequenceOut] = value * value; // already there
    2588                     else
    2589                          infeasible_->quickAdd(sequenceOut, value * value);
    2590                } else {
    2591                     infeasible_->zero(sequenceOut);
    2592                }
    2593                break;
     2433            if (fabs(value) > FREE_ACCEPT * tolerance) {
     2434              // we are going to bias towards free (but only if reasonable)
     2435              value *= FREE_BIAS;
     2436              // store square in list
     2437              if (infeas[iSequence + addSequence])
     2438                infeas[iSequence + addSequence] = value * value; // already there
     2439              else
     2440                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2441            } else {
     2442              infeasible_->zero(iSequence + addSequence);
     2443            }
     2444            break;
    25942445          case ClpSimplex::atUpperBound:
    2595                if (value > tolerance) {
    2596                     // store square in list
    2597                     if (infeas[sequenceOut])
    2598                          infeas[sequenceOut] = value * value; // already there
    2599                     else
    2600                          infeasible_->quickAdd(sequenceOut, value * value);
    2601                } else {
    2602                     infeasible_->zero(sequenceOut);
    2603                }
    2604                break;
     2446            if (value > tolerance) {
     2447              // store square in list
     2448              if (infeas[iSequence + addSequence])
     2449                infeas[iSequence + addSequence] = value * value; // already there
     2450              else
     2451                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2452            } else {
     2453              infeasible_->zero(iSequence + addSequence);
     2454            }
     2455            break;
    26052456          case ClpSimplex::atLowerBound:
    2606                if (value < -tolerance) {
    2607                     // store square in list
    2608                     if (infeas[sequenceOut])
    2609                          infeas[sequenceOut] = value * value; // already there
    2610                     else
    2611                          infeasible_->quickAdd(sequenceOut, value * value);
    2612                } else {
    2613                     infeasible_->zero(sequenceOut);
    2614                }
     2457            if (value < -tolerance) {
     2458              // store square in list
     2459              if (infeas[iSequence + addSequence])
     2460                infeas[iSequence + addSequence] = value * value; // already there
     2461              else
     2462                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2463            } else {
     2464              infeasible_->zero(iSequence + addSequence);
     2465            }
    26152466          }
    2616      }
    2617 
    2618      // If in quadratic re-compute all
    2619      if (model_->algorithm() == 2) {
    2620           for (iSection = 0; iSection < 2; iSection++) {
    2621 
    2622                reducedCost = model_->djRegion(iSection);
    2623                int addSequence;
    2624                int iSequence;
    2625 
    2626                if (!iSection) {
    2627                     number = model_->numberRows();
    2628                     addSequence = model_->numberColumns();
    2629                } else {
    2630                     number = model_->numberColumns();
    2631                     addSequence = 0;
    2632                }
    2633 
    2634                if (!model_->nonLinearCost()->lookBothWays()) {
    2635                     for (iSequence = 0; iSequence < number; iSequence++) {
    2636                          double value = reducedCost[iSequence];
    2637                          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    2638 
    2639                          switch(status) {
    2640 
    2641                          case ClpSimplex::basic:
    2642                               infeasible_->zero(iSequence + addSequence);
    2643                          case ClpSimplex::isFixed:
    2644                               break;
    2645                          case ClpSimplex::isFree:
    2646                          case ClpSimplex::superBasic:
    2647                               if (fabs(value) > tolerance) {
    2648                                    // we are going to bias towards free (but only if reasonable)
    2649                                    value *= FREE_BIAS;
    2650                                    // store square in list
    2651                                    if (infeas[iSequence+addSequence])
    2652                                         infeas[iSequence+addSequence] = value * value; // already there
    2653                                    else
    2654                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2655                               } else {
    2656                                    infeasible_->zero(iSequence + addSequence);
    2657                               }
    2658                               break;
    2659                          case ClpSimplex::atUpperBound:
    2660                               if (value > tolerance) {
    2661                                    // store square in list
    2662                                    if (infeas[iSequence+addSequence])
    2663                                         infeas[iSequence+addSequence] = value * value; // already there
    2664                                    else
    2665                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2666                               } else {
    2667                                    infeasible_->zero(iSequence + addSequence);
    2668                               }
    2669                               break;
    2670                          case ClpSimplex::atLowerBound:
    2671                               if (value < -tolerance) {
    2672                                    // store square in list
    2673                                    if (infeas[iSequence+addSequence])
    2674                                         infeas[iSequence+addSequence] = value * value; // already there
    2675                                    else
    2676                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2677                               } else {
    2678                                    infeasible_->zero(iSequence + addSequence);
    2679                               }
    2680                          }
    2681                     }
    2682                } else {
    2683                     // we can go both ways
    2684                     ClpNonLinearCost * nonLinear = model_->nonLinearCost();
    2685                     for (iSequence = 0; iSequence < number; iSequence++) {
    2686                          double value = reducedCost[iSequence];
    2687                          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
    2688 
    2689                          switch(status) {
    2690 
    2691                          case ClpSimplex::basic:
    2692                               infeasible_->zero(iSequence + addSequence);
    2693                          case ClpSimplex::isFixed:
    2694                               break;
    2695                          case ClpSimplex::isFree:
    2696                          case ClpSimplex::superBasic:
    2697                               if (fabs(value) > tolerance) {
    2698                                    // we are going to bias towards free (but only if reasonable)
    2699                                    value *= FREE_BIAS;
    2700                                    // store square in list
    2701                                    if (infeas[iSequence+addSequence])
    2702                                         infeas[iSequence+addSequence] = value * value; // already there
    2703                                    else
    2704                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2705                               } else {
    2706                                    infeasible_->zero(iSequence + addSequence);
    2707                               }
    2708                               break;
    2709                          case ClpSimplex::atUpperBound:
    2710                               if (value > tolerance) {
    2711                                    // store square in list
    2712                                    if (infeas[iSequence+addSequence])
    2713                                         infeas[iSequence+addSequence] = value * value; // already there
    2714                                    else
    2715                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2716                               } else {
    2717                                    // look other way - change up should be negative
    2718                                    value -= nonLinear->changeUpInCost(iSequence + addSequence);
    2719                                    if (value > -tolerance) {
    2720                                         infeasible_->zero(iSequence + addSequence);
    2721                                    } else {
    2722                                         // store square in list
    2723                                         if (infeas[iSequence+addSequence])
    2724                                              infeas[iSequence+addSequence] = value * value; // already there
    2725                                         else
    2726                                              infeasible_->quickAdd(iSequence + addSequence, value * value);
    2727                                    }
    2728                               }
    2729                               break;
    2730                          case ClpSimplex::atLowerBound:
    2731                               if (value < -tolerance) {
    2732                                    // store square in list
    2733                                    if (infeas[iSequence+addSequence])
    2734                                         infeas[iSequence+addSequence] = value * value; // already there
    2735                                    else
    2736                                         infeasible_->quickAdd(iSequence + addSequence, value * value);
    2737                               } else {
    2738                                    // look other way - change down should be positive
    2739                                    value -= nonLinear->changeDownInCost(iSequence + addSequence);
    2740                                    if (value < tolerance) {
    2741                                         infeasible_->zero(iSequence + addSequence);
    2742                                    } else {
    2743                                         // store square in list
    2744                                         if (infeas[iSequence+addSequence])
    2745                                              infeas[iSequence+addSequence] = value * value; // already there
    2746                                         else
    2747                                              infeasible_->quickAdd(iSequence + addSequence, value * value);
    2748                                    }
    2749                               }
    2750                          }
    2751                     }
    2752                }
     2467        }
     2468      } else {
     2469        ClpNonLinearCost *nonLinear = model_->nonLinearCost();
     2470        // We can go up OR down
     2471        for (j = 0; j < number; j++) {
     2472          int iSequence = index[j];
     2473          double value = reducedCost[iSequence];
     2474          value -= updateBy[iSequence];
     2475          reducedCost[iSequence] = value;
     2476          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     2477
     2478          switch (status) {
     2479
     2480          case ClpSimplex::basic:
     2481            infeasible_->zero(iSequence + addSequence);
     2482          case ClpSimplex::isFixed:
     2483            break;
     2484          case ClpSimplex::isFree:
     2485          case ClpSimplex::superBasic:
     2486            if (fabs(value) > FREE_ACCEPT * tolerance) {
     2487              // we are going to bias towards free (but only if reasonable)
     2488              value *= FREE_BIAS;
     2489              // store square in list
     2490              if (infeas[iSequence + addSequence])
     2491                infeas[iSequence + addSequence] = value * value; // already there
     2492              else
     2493                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2494            } else {
     2495              infeasible_->zero(iSequence + addSequence);
     2496            }
     2497            break;
     2498          case ClpSimplex::atUpperBound:
     2499            if (value > tolerance) {
     2500              // store square in list
     2501              if (infeas[iSequence + addSequence])
     2502                infeas[iSequence + addSequence] = value * value; // already there
     2503              else
     2504                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2505            } else {
     2506              // look other way - change up should be negative
     2507              value -= nonLinear->changeUpInCost(iSequence + addSequence);
     2508              if (value > -tolerance) {
     2509                infeasible_->zero(iSequence + addSequence);
     2510              } else {
     2511                // store square in list
     2512                if (infeas[iSequence + addSequence])
     2513                  infeas[iSequence + addSequence] = value * value; // already there
     2514                else
     2515                  infeasible_->quickAdd(iSequence + addSequence, value * value);
     2516              }
     2517            }
     2518            break;
     2519          case ClpSimplex::atLowerBound:
     2520            if (value < -tolerance) {
     2521              // store square in list
     2522              if (infeas[iSequence + addSequence])
     2523                infeas[iSequence + addSequence] = value * value; // already there
     2524              else
     2525                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2526            } else {
     2527              // look other way - change down should be positive
     2528              value -= nonLinear->changeDownInCost(iSequence + addSequence);
     2529              if (value < tolerance) {
     2530                infeasible_->zero(iSequence + addSequence);
     2531              } else {
     2532                // store square in list
     2533                if (infeas[iSequence + addSequence])
     2534                  infeas[iSequence + addSequence] = value * value; // already there
     2535                else
     2536                  infeasible_->quickAdd(iSequence + addSequence, value * value);
     2537              }
     2538            }
    27532539          }
    2754      }
    2755      // See what sort of pricing
    2756      int numberWanted = 10;
    2757      number = infeasible_->getNumElements();
    2758      int numberColumns = model_->numberColumns();
    2759      if (switchType == 5) {
    2760           // we can zero out
    2761           updates->clear();
    2762           spareColumn1->clear();
    2763           pivotSequence_ = -1;
    2764           pivotRow = -1;
    2765           // See if to switch
    2766           int numberRows = model_->numberRows();
    2767           // ratio is done on number of columns here
    2768           //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
    2769           double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    2770           //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
    2771           if (ratio < 0.1) {
    2772                numberWanted = CoinMax(100, number / 200);
    2773           } else if (ratio < 0.3) {
    2774                numberWanted = CoinMax(500, number / 40);
    2775           } else if (ratio < 0.5 || mode_ == 5) {
    2776                numberWanted = CoinMax(2000, number / 10);
    2777                numberWanted = CoinMax(numberWanted, numberColumns / 30);
    2778           } else if (mode_ != 5) {
    2779                switchType = 4;
    2780                // initialize
    2781                numberSwitched_++;
    2782                // Make sure will re-do
    2783                delete [] weights_;
    2784                weights_ = NULL;
    2785                saveWeights(model_, 4);
    2786                COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
    2787                updates->clear();
     2540        }
     2541      }
     2542    }
     2543    if (anyUpdates == 2) {
     2544      // we can zero out as will have to get pivot row
     2545      updates->clear();
     2546      spareColumn1->clear();
     2547    }
     2548    if (pivotRow >= 0) {
     2549      // make sure infeasibility on incoming is 0.0
     2550      int sequenceIn = model_->sequenceIn();
     2551      infeasible_->zero(sequenceIn);
     2552    }
     2553  }
     2554  // make sure outgoing from last iteration okay
     2555  int sequenceOut = model_->sequenceOut();
     2556  if (sequenceOut >= 0) {
     2557    ClpSimplex::Status status = model_->getStatus(sequenceOut);
     2558    double value = model_->reducedCost(sequenceOut);
     2559
     2560    switch (status) {
     2561
     2562    case ClpSimplex::basic:
     2563    case ClpSimplex::isFixed:
     2564      break;
     2565    case ClpSimplex::isFree:
     2566    case ClpSimplex::superBasic:
     2567      if (fabs(value) > FREE_ACCEPT * tolerance) {
     2568        // we are going to bias towards free (but only if reasonable)
     2569        value *= FREE_BIAS;
     2570        // store square in list
     2571        if (infeas[sequenceOut])
     2572          infeas[sequenceOut] = value * value; // already there
     2573        else
     2574          infeasible_->quickAdd(sequenceOut, value * value);
     2575      } else {
     2576        infeasible_->zero(sequenceOut);
     2577      }
     2578      break;
     2579    case ClpSimplex::atUpperBound:
     2580      if (value > tolerance) {
     2581        // store square in list
     2582        if (infeas[sequenceOut])
     2583          infeas[sequenceOut] = value * value; // already there
     2584        else
     2585          infeasible_->quickAdd(sequenceOut, value * value);
     2586      } else {
     2587        infeasible_->zero(sequenceOut);
     2588      }
     2589      break;
     2590    case ClpSimplex::atLowerBound:
     2591      if (value < -tolerance) {
     2592        // store square in list
     2593        if (infeas[sequenceOut])
     2594          infeas[sequenceOut] = value * value; // already there
     2595        else
     2596          infeasible_->quickAdd(sequenceOut, value * value);
     2597      } else {
     2598        infeasible_->zero(sequenceOut);
     2599      }
     2600    }
     2601  }
     2602
     2603  // If in quadratic re-compute all
     2604  if (model_->algorithm() == 2) {
     2605    for (iSection = 0; iSection < 2; iSection++) {
     2606
     2607      reducedCost = model_->djRegion(iSection);
     2608      int addSequence;
     2609      int iSequence;
     2610
     2611      if (!iSection) {
     2612        number = model_->numberRows();
     2613        addSequence = model_->numberColumns();
     2614      } else {
     2615        number = model_->numberColumns();
     2616        addSequence = 0;
     2617      }
     2618
     2619      if (!model_->nonLinearCost()->lookBothWays()) {
     2620        for (iSequence = 0; iSequence < number; iSequence++) {
     2621          double value = reducedCost[iSequence];
     2622          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     2623
     2624          switch (status) {
     2625
     2626          case ClpSimplex::basic:
     2627            infeasible_->zero(iSequence + addSequence);
     2628          case ClpSimplex::isFixed:
     2629            break;
     2630          case ClpSimplex::isFree:
     2631          case ClpSimplex::superBasic:
     2632            if (fabs(value) > tolerance) {
     2633              // we are going to bias towards free (but only if reasonable)
     2634              value *= FREE_BIAS;
     2635              // store square in list
     2636              if (infeas[iSequence + addSequence])
     2637                infeas[iSequence + addSequence] = value * value; // already there
     2638              else
     2639                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2640            } else {
     2641              infeasible_->zero(iSequence + addSequence);
     2642            }
     2643            break;
     2644          case ClpSimplex::atUpperBound:
     2645            if (value > tolerance) {
     2646              // store square in list
     2647              if (infeas[iSequence + addSequence])
     2648                infeas[iSequence + addSequence] = value * value; // already there
     2649              else
     2650                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2651            } else {
     2652              infeasible_->zero(iSequence + addSequence);
     2653            }
     2654            break;
     2655          case ClpSimplex::atLowerBound:
     2656            if (value < -tolerance) {
     2657              // store square in list
     2658              if (infeas[iSequence + addSequence])
     2659                infeas[iSequence + addSequence] = value * value; // already there
     2660              else
     2661                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2662            } else {
     2663              infeasible_->zero(iSequence + addSequence);
     2664            }
    27882665          }
    2789           if (model_->numberIterations() % 1000 == 0)
    2790             COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
    2791      }
    2792      if(switchType == 4) {
    2793           // Still in devex mode
    2794           int numberRows = model_->numberRows();
    2795           // ratio is done on number of rows here
    2796           double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
    2797           // Go to steepest if lot of iterations?
    2798           if (ratio < 1.0) {
    2799                numberWanted = CoinMax(2000, number / 20);
    2800           } else if (ratio < 5.0) {
    2801                numberWanted = CoinMax(2000, number / 10);
    2802                numberWanted = CoinMax(numberWanted, numberColumns / 20);
     2666        }
     2667      } else {
     2668        // we can go both ways
     2669        ClpNonLinearCost *nonLinear = model_->nonLinearCost();
     2670        for (iSequence = 0; iSequence < number; iSequence++) {
     2671          double value = reducedCost[iSequence];
     2672          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
     2673
     2674          switch (status) {
     2675
     2676          case ClpSimplex::basic:
     2677            infeasible_->zero(iSequence + addSequence);
     2678          case ClpSimplex::isFixed:
     2679            break;
     2680          case ClpSimplex::isFree:
     2681          case ClpSimplex::superBasic:
     2682            if (fabs(value) > tolerance) {
     2683              // we are going to bias towards free (but only if reasonable)
     2684              value *= FREE_BIAS;
     2685              // store square in list
     2686              if (infeas[iSequence + addSequence])
     2687                infeas[iSequence + addSequence] = value * value; // already there
     2688              else
     2689                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2690            } else {
     2691              infeasible_->zero(iSequence + addSequence);
     2692            }
     2693            break;
     2694          case ClpSimplex::atUpperBound:
     2695            if (value > tolerance) {
     2696              // store square in list
     2697              if (infeas[iSequence + addSequence])
     2698                infeas[iSequence + addSequence] = value * value; // already there
     2699              else
     2700                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2701            } else {
     2702              // look other way - change up should be negative
     2703              value -= nonLinear->changeUpInCost(iSequence + addSequence);
     2704              if (value > -tolerance) {
     2705                infeasible_->zero(iSequence + addSequence);
     2706              } else {
     2707                // store square in list
     2708                if (infeas[iSequence + addSequence])
     2709                  infeas[iSequence + addSequence] = value * value; // already there
     2710                else
     2711                  infeasible_->quickAdd(iSequence + addSequence, value * value);
     2712              }
     2713            }
     2714            break;
     2715          case ClpSimplex::atLowerBound:
     2716            if (value < -tolerance) {
     2717              // store square in list
     2718              if (infeas[iSequence + addSequence])
     2719                infeas[iSequence + addSequence] = value * value; // already there
     2720              else
     2721                infeasible_->quickAdd(iSequence + addSequence, value * value);
     2722            } else {
     2723              // look other way - change down should be positive
     2724              value -= nonLinear->changeDownInCost(iSequence + addSequence);
     2725              if (value < tolerance) {
     2726                infeasible_->zero(iSequence + addSequence);
     2727              } else {
     2728                // store square in list
     2729                if (infeas[iSequence + addSequence])
     2730                  infeas[iSequence + addSequence] = value * value; // already there
     2731                else
     2732                  infeasible_->quickAdd(iSequence + addSequence, value * value);
     2733              }
     2734            }
     2735          }
     2736        }
     2737      }
     2738    }
     2739  }
     2740  // See what sort of pricing
     2741  int numberWanted = 10;
     2742  number = infeasible_->getNumElements();
     2743  int numberColumns = model_->numberColumns();
     2744  if (switchType == 5) {
     2745    // we can zero out
     2746    updates->clear();
     2747    spareColumn1->clear();
     2748    pivotSequence_ = -1;
     2749    pivotRow = -1;
     2750    // See if to switch
     2751    int numberRows = model_->numberRows();
     2752    // ratio is done on number of columns here
     2753    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
     2754    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     2755    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
     2756    if (ratio < 0.1) {
     2757      numberWanted = CoinMax(100, number / 200);
     2758    } else if (ratio < 0.3) {
     2759      numberWanted = CoinMax(500, number / 40);
     2760    } else if (ratio < 0.5 || mode_ == 5) {
     2761      numberWanted = CoinMax(2000, number / 10);
     2762      numberWanted = CoinMax(numberWanted, numberColumns / 30);
     2763    } else if (mode_ != 5) {
     2764      switchType = 4;
     2765      // initialize
     2766      numberSwitched_++;
     2767      // Make sure will re-do
     2768      delete[] weights_;
     2769      weights_ = NULL;
     2770      saveWeights(model_, 4);
     2771      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
     2772      updates->clear();
     2773    }
     2774    if (model_->numberIterations() % 1000 == 0)
     2775      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
     2776  }
     2777  if (switchType == 4) {
     2778    // Still in devex mode
     2779    int numberRows = model_->numberRows();
     2780    // ratio is done on number of rows here
     2781    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
     2782    // Go to steepest if lot of iterations?
     2783    if (ratio < 1.0) {
     2784      numberWanted = CoinMax(2000, number / 20);
     2785    } else if (ratio < 5.0) {
     2786      numberWanted = CoinMax(2000, number / 10);
     2787      numberWanted = CoinMax(numberWanted, numberColumns / 20);
     2788    } else {
     2789      // we can zero out
     2790      updates->clear();
     2791      spareColumn1->clear();
     2792      switchType = 3;
     2793      // initialize
     2794      pivotSequence_ = -1;
     2795      pivotRow = -1;
     2796      numberSwitched_++;
     2797      // Make sure will re-do
     2798      delete[] weights_;
     2799      weights_ = NULL;
     2800      saveWeights(model_, 4);
     2801      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
     2802      updates->clear();
     2803    }
     2804    if (model_->numberIterations() % 1000 == 0)
     2805      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
     2806  }
     2807  if (switchType < 4) {
     2808    if (switchType < 2) {
     2809      numberWanted = number + 1;
     2810    } else if (switchType == 2) {
     2811      numberWanted = CoinMax(2000, number / 8);
     2812    } else {
     2813      double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(model_->numberRows());
     2814      if (ratio < 1.0) {
     2815        numberWanted = CoinMax(2000, number / 20);
     2816      } else if (ratio < 5.0) {
     2817        numberWanted = CoinMax(2000, number / 10);
     2818        numberWanted = CoinMax(numberWanted, numberColumns / 20);
     2819      } else if (ratio < 10.0) {
     2820        numberWanted = CoinMax(2000, number / 8);
     2821        numberWanted = CoinMax(numberWanted, numberColumns / 20);
     2822      } else {
     2823        ratio = number * (ratio / 80.0);
     2824        if (ratio > number) {
     2825          numberWanted = number + 1;
     2826        } else {
     2827          numberWanted = CoinMax(2000, static_cast< int >(ratio));
     2828          numberWanted = CoinMax(numberWanted, numberColumns / 10);
     2829        }
     2830      }
     2831    }
     2832  }
     2833  // for weights update we use pivotSequence
     2834  pivotRow = pivotSequence_;
     2835  // unset in case sub flip
     2836  pivotSequence_ = -1;
     2837  if (pivotRow >= 0) {
     2838    // make sure infeasibility on incoming is 0.0
     2839    const int *pivotVariable = model_->pivotVariable();
     2840    int sequenceIn = pivotVariable[pivotRow];
     2841    infeasible_->zero(sequenceIn);
     2842    // and we can see if reference
     2843    double referenceIn = 0.0;
     2844    if (switchType != 1 && reference(sequenceIn))
     2845      referenceIn = 1.0;
     2846    // save outgoing weight round update
     2847    double outgoingWeight = 0.0;
     2848    if (sequenceOut >= 0)
     2849      outgoingWeight = weights_[sequenceOut];
     2850    // update weights
     2851    if (anyUpdates != 1) {
     2852      updates->setNumElements(0);
     2853      spareColumn1->setNumElements(0);
     2854      // might as well set dj to 1
     2855      dj = 1.0;
     2856      updates->insert(pivotRow, -dj);
     2857      model_->factorization()->updateColumnTranspose(spareRow2, updates);
     2858      // put row of tableau in rowArray and columnArray
     2859      model_->clpMatrix()->transposeTimes(model_, -1.0,
     2860        updates, spareColumn2, spareColumn1);
     2861    }
     2862    double *weight;
     2863    double *other = alternateWeights_->denseVector();
     2864    int numberColumns = model_->numberColumns();
     2865    double scaleFactor = -1.0 / dj; // as formula is with 1.0
     2866    // rows
     2867    number = updates->getNumElements();
     2868    index = updates->getIndices();
     2869    updateBy = updates->denseVector();
     2870    weight = weights_ + numberColumns;
     2871
     2872    if (switchType < 4) {
     2873      // Exact
     2874      // now update weight update array
     2875      model_->factorization()->updateColumnTranspose(spareRow2,
     2876        alternateWeights_);
     2877#ifdef ALT_UPDATE_WEIGHTS
     2878      abort();
     2879#endif
     2880      for (j = 0; j < number; j++) {
     2881        int iSequence = index[j];
     2882        double thisWeight = weight[iSequence];
     2883        // row has -1
     2884        double pivot = updateBy[iSequence] * scaleFactor;
     2885        updateBy[iSequence] = 0.0;
     2886        double modification = other[iSequence];
     2887        double pivotSquared = pivot * pivot;
     2888
     2889        thisWeight += pivotSquared * devex_ + pivot * modification;
     2890        if (thisWeight < TRY_NORM) {
     2891          if (switchType == 1) {
     2892            // steepest
     2893            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    28032894          } else {
    2804                // we can zero out
    2805                updates->clear();
    2806                spareColumn1->clear();
    2807                switchType = 3;
    2808                // initialize
    2809                pivotSequence_ = -1;
    2810                pivotRow = -1;
    2811                numberSwitched_++;
    2812                // Make sure will re-do
    2813                delete [] weights_;
    2814                weights_ = NULL;
    2815                saveWeights(model_, 4);
    2816                COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
    2817                updates->clear();
     2895            // exact
     2896            thisWeight = referenceIn * pivotSquared;
     2897            if (reference(iSequence + numberColumns))
     2898              thisWeight += 1.0;
     2899            thisWeight = CoinMax(thisWeight, TRY_NORM);
    28182900          }
    2819           if (model_->numberIterations() % 1000 == 0)
    2820             COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
    2821      }
    2822      if (switchType < 4) {
    2823           if (switchType < 2 ) {
    2824                numberWanted = number + 1;
    2825           } else if (switchType == 2) {
    2826                numberWanted = CoinMax(2000, number / 8);
     2901        }
     2902        weight[iSequence] = thisWeight;
     2903      }
     2904    } else if (switchType == 4) {
     2905      // Devex
     2906      assert(devex_ > 0.0);
     2907      for (j = 0; j < number; j++) {
     2908        int iSequence = index[j];
     2909        double thisWeight = weight[iSequence];
     2910        // row has -1
     2911        double pivot = updateBy[iSequence] * scaleFactor;
     2912        updateBy[iSequence] = 0.0;
     2913        double value = pivot * pivot * devex_;
     2914        if (reference(iSequence + numberColumns))
     2915          value += 1.0;
     2916        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     2917      }
     2918    }
     2919
     2920    // columns
     2921    weight = weights_;
     2922
     2923    scaleFactor = -scaleFactor;
     2924
     2925    number = spareColumn1->getNumElements();
     2926    index = spareColumn1->getIndices();
     2927    updateBy = spareColumn1->denseVector();
     2928    if (switchType < 4) {
     2929      // Exact
     2930      // get subset which have nonzero tableau elements
     2931      model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
     2932        spareColumn1,
     2933        spareColumn2);
     2934      double *updateBy2 = spareColumn2->denseVector();
     2935      for (j = 0; j < number; j++) {
     2936        int iSequence = index[j];
     2937        double thisWeight = weight[iSequence];
     2938        double pivot = updateBy[iSequence] * scaleFactor;
     2939        updateBy[iSequence] = 0.0;
     2940        double modification = updateBy2[j];
     2941        updateBy2[j] = 0.0;
     2942        double pivotSquared = pivot * pivot;
     2943
     2944        thisWeight += pivotSquared * devex_ + pivot * modification;
     2945        if (thisWeight < TRY_NORM) {
     2946          if (switchType == 1) {
     2947            // steepest
     2948            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    28272949          } else {
    2828                double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (model_->numberRows());
    2829                if (ratio < 1.0) {
    2830                     numberWanted = CoinMax(2000, number / 20);
    2831                } else if (ratio < 5.0) {
    2832                     numberWanted = CoinMax(2000, number / 10);
    2833                     numberWanted = CoinMax(numberWanted, numberColumns / 20);
    2834                } else if (ratio < 10.0) {
    2835                     numberWanted = CoinMax(2000, number / 8);
    2836                     numberWanted = CoinMax(numberWanted, numberColumns / 20);
    2837                } else {
    2838                     ratio = number * (ratio / 80.0);
    2839                     if (ratio > number) {
    2840                          numberWanted = number + 1;
    2841                     } else {
    2842                          numberWanted = CoinMax(2000, static_cast<int> (ratio));
    2843                          numberWanted = CoinMax(numberWanted, numberColumns / 10);
    2844                     }
    2845                }
     2950            // exact
     2951            thisWeight = referenceIn * pivotSquared;
     2952            if (reference(iSequence))
     2953              thisWeight += 1.0;
     2954            thisWeight = CoinMax(thisWeight, TRY_NORM);
    28462955          }
    2847      }
    2848      // for weights update we use pivotSequence
    2849      pivotRow = pivotSequence_;
    2850      // unset in case sub flip
    2851      pivotSequence_ = -1;
    2852      if (pivotRow >= 0) {
    2853           // make sure infeasibility on incoming is 0.0
    2854           const int * pivotVariable = model_->pivotVariable();
    2855           int sequenceIn = pivotVariable[pivotRow];
    2856           infeasible_->zero(sequenceIn);
    2857           // and we can see if reference
    2858           double referenceIn = 0.0;
    2859           if (switchType != 1 && reference(sequenceIn))
    2860                referenceIn = 1.0;
    2861           // save outgoing weight round update
    2862           double outgoingWeight = 0.0;
    2863           if (sequenceOut >= 0)
    2864                outgoingWeight = weights_[sequenceOut];
    2865           // update weights
    2866           if (anyUpdates != 1) {
    2867                updates->setNumElements(0);
    2868                spareColumn1->setNumElements(0);
    2869                // might as well set dj to 1
    2870                dj = 1.0;
    2871                updates->insert(pivotRow, -dj);
    2872                model_->factorization()->updateColumnTranspose(spareRow2, updates);
    2873                // put row of tableau in rowArray and columnArray
    2874                model_->clpMatrix()->transposeTimes(model_, -1.0,
    2875                                                    updates, spareColumn2, spareColumn1);
     2956        }
     2957        weight[iSequence] = thisWeight;
     2958      }
     2959    } else if (switchType == 4) {
     2960      // Devex
     2961      for (j = 0; j < number; j++) {
     2962        int iSequence = index[j];
     2963        double thisWeight = weight[iSequence];
     2964        // row has -1
     2965        double pivot = updateBy[iSequence] * scaleFactor;
     2966        updateBy[iSequence] = 0.0;
     2967        double value = pivot * pivot * devex_;
     2968        if (reference(iSequence))
     2969          value += 1.0;
     2970        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
     2971      }
     2972    }
     2973    // restore outgoing weight
     2974    if (sequenceOut >= 0)
     2975      weights_[sequenceOut] = outgoingWeight;
     2976    alternateWeights_->clear();
     2977    spareColumn2->setNumElements(0);
     2978    //#define SOME_DEBUG_1
     2979#ifdef SOME_DEBUG_1
     2980    // check for accuracy
     2981    int iCheck = 892;
     2982    //printf("weight for iCheck is %g\n",weights_[iCheck]);
     2983    int numberRows = model_->numberRows();
     2984    //int numberColumns = model_->numberColumns();
     2985    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
     2986      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
     2987        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
     2988    }
     2989#endif
     2990    updates->setNumElements(0);
     2991    spareColumn1->setNumElements(0);
     2992  }
     2993
     2994  // update of duals finished - now do pricing
     2995
     2996  double bestDj = 1.0e-30;
     2997  int bestSequence = -1;
     2998
     2999  int i, iSequence;
     3000  index = infeasible_->getIndices();
     3001  number = infeasible_->getNumElements();
     3002  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
     3003    // we can't really trust infeasibilities if there is dual error
     3004    double checkTolerance = 1.0e-8;
     3005    if (!model_->factorization()->pivots())
     3006      checkTolerance = 1.0e-6;
     3007    if (model_->largestDualError() > checkTolerance)
     3008      tolerance *= model_->largestDualError() / checkTolerance;
     3009    // But cap
     3010    tolerance = CoinMin(1000.0, tolerance);
     3011  }
     3012#ifdef CLP_DEBUG
     3013  if (model_->numberDualInfeasibilities() == 1)
     3014    printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
     3015      model_->largestDualError(), model_, model_->messageHandler(),
     3016      number);
     3017#endif
     3018  // stop last one coming immediately
     3019  double saveOutInfeasibility = 0.0;
     3020  if (sequenceOut >= 0) {
     3021    saveOutInfeasibility = infeas[sequenceOut];
     3022    infeas[sequenceOut] = 0.0;
     3023  }
     3024  tolerance *= tolerance; // as we are using squares
     3025
     3026  int iPass;
     3027  // Setup two passes
     3028  int start[4];
     3029  start[1] = number;
     3030  start[2] = 0;
     3031  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
     3032  start[0] = static_cast< int >(dstart);
     3033  start[3] = start[0];
     3034  //double largestWeight=0.0;
     3035  //double smallestWeight=1.0e100;
     3036  for (iPass = 0; iPass < 2; iPass++) {
     3037    int end = start[2 * iPass + 1];
     3038    if (switchType < 5) {
     3039      for (i = start[2 * iPass]; i < end; i++) {
     3040        iSequence = index[i];
     3041        double value = infeas[iSequence];
     3042        if (value > tolerance) {
     3043          double weight = weights_[iSequence];
     3044          //weight=1.0;
     3045          if (value > bestDj * weight) {
     3046            // check flagged variable and correct dj
     3047            if (!model_->flagged(iSequence)) {
     3048              bestDj = value / weight;
     3049              bestSequence = iSequence;
     3050            } else {
     3051              // just to make sure we don't exit before got something
     3052              numberWanted++;
     3053            }
    28763054          }
    2877           double * weight;
    2878           double * other = alternateWeights_->denseVector();
    2879           int numberColumns = model_->numberColumns();
    2880           double scaleFactor = -1.0 / dj; // as formula is with 1.0
    2881           // rows
    2882           number = updates->getNumElements();
    2883           index = updates->getIndices();
    2884           updateBy = updates->denseVector();
    2885           weight = weights_ + numberColumns;
    2886 
    2887           if (switchType < 4) {
    2888                // Exact
    2889                // now update weight update array
    2890                model_->factorization()->updateColumnTranspose(spareRow2,
    2891                          alternateWeights_);
    2892 #ifdef ALT_UPDATE_WEIGHTS
    2893                abort();
    2894 #endif
    2895                for (j = 0; j < number; j++) {
    2896                     int iSequence = index[j];
    2897                     double thisWeight = weight[iSequence];
    2898                     // row has -1
    2899                     double pivot = updateBy[iSequence] * scaleFactor;
    2900                     updateBy[iSequence] = 0.0;
    2901                     double modification = other[iSequence];
    2902                     double pivotSquared = pivot * pivot;
    2903 
    2904                     thisWeight += pivotSquared * devex_ + pivot * modification;
    2905                     if (thisWeight < TRY_NORM) {
    2906                          if (switchType == 1) {
    2907                               // steepest
    2908                               thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
    2909                          } else {
    2910                               // exact
    2911                               thisWeight = referenceIn * pivotSquared;
    2912                               if (reference(iSequence + numberColumns))
    2913                                    thisWeight += 1.0;
    2914                               thisWeight = CoinMax(thisWeight, TRY_NORM);
    2915                          }
    2916                     }
    2917                     weight[iSequence] = thisWeight;
    2918                }
    2919           } else if (switchType == 4) {
    2920                // Devex
    2921                assert (devex_ > 0.0);
    2922                for (j = 0; j < number; j++) {
    2923                     int iSequence = index[j];
    2924                     double thisWeight = weight[iSequence];
    2925<