source: trunk/Clp/src/ClpPrimalColumnSteepest.cpp @ 2470

Last change on this file since 2470 was 2444, checked in by forrest, 7 months ago

added assert in case #86 exists

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 145.8 KB
RevLine 
[1370]1/* $Id: ClpPrimalColumnSteepest.cpp 2444 2019-04-05 16:11:24Z stefan $ */
[2]2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
[1665]4// This code is licensed under the terms of the Eclipse Public License (EPL).
[2]5
[63]6#include "CoinPragma.hpp"
[2235]7#if ABOCA_LITE
8// 1 is not owner of abcState_
9#define ABCSTATE_LITE 1
10#endif
11//#define FAKE_CILK
[2]12
13#include "ClpSimplex.hpp"
14#include "ClpPrimalColumnSteepest.hpp"
[50]15#include "CoinIndexedVector.hpp"
[2]16#include "ClpFactorization.hpp"
[162]17#include "ClpNonLinearCost.hpp"
[2]18#include "ClpMessage.hpp"
[2235]19#include "ClpEventHandler.hpp"
20#include "ClpPackedMatrix.hpp"
[2]21#include "CoinHelperFunctions.hpp"
22#include <stdio.h>
[2235]23#if ABOCA_LITE
24#undef ALT_UPDATE_WEIGHTS
25#define ALT_UPDATE_WEIGHTS 2
26#endif
[2385]27#if ALT_UPDATE_WEIGHTS == 1
28extern CoinIndexedVector *altVector[3];
[2235]29#endif
[2385]30static void debug1(int iSequence, double value, double weight)
[2235]31{
32  printf("xx %d inf %.20g wt %.20g\n",
[2385]33    iSequence, value, weight);
[2235]34}
[341]35//#define CLP_DEBUG
[2]36//#############################################################################
37// Constructors / Destructor / Assignment
38//#############################################################################
39
40//-------------------------------------------------------------------
[1525]41// Default Constructor
[2]42//-------------------------------------------------------------------
[2385]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)
[2]60{
[2385]61  type_ = 2 + 64 * mode;
[2]62}
63//-------------------------------------------------------------------
[1525]64// Copy constructor
[2]65//-------------------------------------------------------------------
[2385]66ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(const ClpPrimalColumnSteepest &rhs)
67  : ClpPrimalColumnPivot(rhs)
[1525]68{
[2385]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  }
[2]114}
115
116//-------------------------------------------------------------------
[1525]117// Destructor
[2]118//-------------------------------------------------------------------
[2385]119ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest()
[2]120{
[2385]121  delete[] weights_;
122  delete infeasible_;
123  delete alternateWeights_;
124  delete[] savedWeights_;
125  delete[] reference_;
[2]126}
127
128//----------------------------------------------------------------
[1525]129// Assignment operator
[2]130//-------------------------------------------------------------------
131ClpPrimalColumnSteepest &
[2385]132ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest &rhs)
[2]133{
[2385]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;
[2]180}
[517]181// These have to match ClpPackedMatrix version
[2]182#define TRY_NORM 1.0e-4
183#define ADD_ONE 1.0
[2235]184static void
[2385]185pivotColumnBit(clpTempInfo &info)
[2235]186{
[2385]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;
[2235]195  for (int i = start[0]; i < start[1]; i++) {
196    int iSequence = index[i];
197    double value = infeas[iSequence];
198    double weight = weights[iSequence];
199    if (value > tolerance) {
200      if (value > bestDj * weight) {
[2385]201        // check flagged variable and correct dj
202        if ((status[iSequence] & 64) == 0) {
203          bestDj = value / weight;
204          bestSequence = iSequence;
205        }
[2235]206      }
207    }
208  }
[2385]209  info.primalRatio = bestDj;
210  info.numberAdded = bestSequence;
[2235]211}
[2]212// Returns pivot column, -1 if none
[259]213/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
[1525]214        that is just +-weight where a feasibility changed.  It also has
[259]215        reduced cost from last iteration in pivot row*/
[2385]216int ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates,
217  CoinIndexedVector *spareRow1,
218  CoinIndexedVector *spareRow2,
219  CoinIndexedVector *spareColumn1,
220  CoinIndexedVector *spareColumn2)
[2]221{
[2385]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();
[225]240
[2385]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 -
[1525]250        0 - all exact devex
251        1 - all steepest
252        2 - some exact devex
253        3 - auto some exact devex
254        4 - devex
255        5 - dantzig
256        10 - can go to mini-sprint
257     */
[2385]258    // Look at gub
[336]259#if 1
[2385]260  model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
[336]261#else
[2385]262  updates->clear();
263  model_->computeDuals(NULL);
[336]264#endif
[2385]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
[451]396#ifdef CLP_DEBUG
[2385]397  alternateWeights_->checkClear();
[451]398#endif
[2385]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);
[225]403
[2385]404    switch (status) {
[225]405
[2385]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  }
[225]539
[2385]540  double bestDj = 1.0e-30;
541  bestSequence = -1;
[2235]542#ifdef CLP_USER_DRIVEN
[2385]543  // could go parallel?
544  model_->eventHandler()->event(ClpEventHandler::beforeChooseIncoming);
[2235]545#endif
[1525]546
[2385]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  }
[225]573#ifdef CLP_DEBUG
[2385]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);
[225]578#endif
[2385]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
[225]588
[2385]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--;
619          }
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--;
640          }
641          if (!numberWanted)
642            break;
643        }
644      }
645      if (!numberWanted)
646        break;
647    }
648  } else {
[2235]649#if ABOCA_LITE
[2385]650    int numberThreads = CoinMax(abcState(), 1);
[2235]651#define ABOCA_LITE_MAX ABOCA_LITE
652#else
[2385]653    const int numberThreads = 1;
[2235]654#define ABOCA_LITE_MAX 1
655#endif
[2385]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)
[1525]711       printf("%d best %g\n",bestSequence,bestDj);*/
712
[2235]713#ifdef CLP_USER_DRIVEN
[2385]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));
[2235]723#endif
[225]724#ifndef NDEBUG
[2385]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  }
[225]732#endif
[2235]733#ifdef ALT_UPDATE_WEIGHTSz
[2385]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");
[2235]743#endif
[1878]744#if 0
745     for (int i=0;i<numberRows;i++)
746       printf("row %d weight %g infeas %g\n",i,weights_[i+numberColumns],infeas[i+numberColumns]);
747     for (int i=0;i<numberColumns;i++)
748       printf("column %d weight %g infeas %g\n",i,weights_[i],infeas[i]);
749#endif
[2385]750  return bestSequence;
[225]751}
752// Just update djs
[2385]753void ClpPrimalColumnSteepest::justDjs(CoinIndexedVector *updates,
754  CoinIndexedVector *spareRow2,
755  CoinIndexedVector *spareColumn1,
756  CoinIndexedVector *spareColumn2)
[225]757{
[2385]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);
[1197]773
[2385]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++) {
[1525]779
[2385]780    reducedCost = model_->djRegion(iSection);
781    int addSequence;
782#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
783    double slack_multiplier;
[1732]784#endif
[1525]785
[2385]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;
[1732]793#endif
[2385]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;
[1732]801#endif
[2385]802    }
[1525]803
[2385]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);
[1525]811
[2385]812      switch (status) {
[1525]813
[2385]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) {
[1732]835#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]836          value *= value * slack_multiplier;
[1732]837#else
[2385]838          value *= value;
[1732]839#endif
[2385]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) {
[1732]852#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]853          value *= value * slack_multiplier;
[1732]854#else
[2385]855          value *= value;
[1732]856#endif
[2385]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  }
[225]875}
876// Update djs, weights for Devex
[2385]877void ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector *updates,
878  CoinIndexedVector *spareRow2,
879  CoinIndexedVector *spareColumn1,
880  CoinIndexedVector *spareColumn2)
[225]881{
[2385]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];
[225]911
[2385]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  ;
[1525]923
[2385]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);
[1525]940
[2385]941    switch (status) {
[1525]942
[2385]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
[1732]979#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]980        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
[1732]981#else
[2385]982        value *= value;
[1732]983#endif
[2385]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
[1732]1003#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1004        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
[1732]1005#else
[2385]1006        value *= value;
[1732]1007#endif
[2385]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  }
[1525]1017
[2385]1018  // columns
1019  weight = weights_;
[1525]1020
[2385]1021  scaleFactor = -scaleFactor;
1022  reducedCost = model_->djRegion(1);
1023  number = spareColumn1->getNumElements();
1024  index = spareColumn1->getIndices();
1025  updateBy = spareColumn1->denseVector();
[1525]1026
[2385]1027  // Devex
[1525]1028
[2385]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);
[1525]1040
[2385]1041    switch (status) {
[1525]1042
[2385]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
[225]1112#ifdef SOME_DEBUG_1
[2385]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  }
[225]1122#endif
[2385]1123  updates->setNumElements(0);
1124  spareColumn1->setNumElements(0);
[225]1125}
1126// Update djs, weights for Steepest
[2385]1127void ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector *updates,
1128  CoinIndexedVector *spareRow2,
1129  CoinIndexedVector *spareColumn1,
1130  CoinIndexedVector *spareColumn2)
[225]1131{
[2385]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;
[2235]1151#if 1
[2385]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  }
[2235]1177#endif
1178#else
[2385]1179  model_->factorization()->updateTwoColumnsTranspose(spareRow2, updates, alternateWeights_);
[2235]1180#endif
[2385]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  ;
[225]1207
[2385]1208  number = updates->getNumElements();
1209  index = updates->getIndices();
1210  updateBy = updates->denseVector();
1211  weight = weights_ + numberColumns;
[2235]1212#ifdef CLP_USER_DRIVEN
[2385]1213  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
[2235]1214#endif
[1525]1215
[2385]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;
[1525]1225
[2385]1226    switch (status) {
[1525]1227
[2385]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;
[1525]1241
[2385]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;
[1525]1276
[2385]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
[1732]1295#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1296        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
[1732]1297#else
[2385]1298        value *= value;
[1732]1299#endif
[2385]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;
[1525]1315
[2385]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
[1732]1334#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1335        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
[1732]1336#else
[2385]1337        value *= value;
[1732]1338#endif
[2385]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  }
[2235]1348#ifdef CLP_USER_DRIVEN
[2385]1349  // could go parallel?
1350  //model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
[2235]1351#endif
[2385]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);
[1525]1383
[2385]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
1406          else
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
[2235]1441#ifdef GROW_REFERENCE
[2385]1442    if (!reference(sequenceOut)) {
1443      outgoingWeight += 1.0;
1444      setReference(sequenceOut, true);
1445    }
[2235]1446#endif
[2385]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
[225]1456#ifdef SOME_DEBUG_1
[2385]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  }
[225]1466#endif
[2385]1467  updates->setNumElements(0);
1468  spareColumn1->setNumElements(0);
[225]1469}
1470// Update djs, weights for Devex
[2385]1471void ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector *updates,
1472  CoinIndexedVector *spareRow2,
1473  CoinIndexedVector *spareColumn1,
1474  CoinIndexedVector *spareColumn2)
[225]1475{
[2385]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);
[1197]1493
[2385]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++) {
[1525]1499
[2385]1500    reducedCost = model_->djRegion(iSection);
1501    int addSequence;
1502#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1503    double slack_multiplier;
[1732]1504#endif
[1525]1505
[2385]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;
[1732]1513#endif
[2385]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;
[1732]1521#endif
[2385]1522    }
[1525]1523
[2385]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);
[1525]1531
[2385]1532      switch (status) {
[1525]1533
[2385]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) {
[1732]1555#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1556          value *= value * slack_multiplier;
[1732]1557#else
[2385]1558          value *= value;
[1732]1559#endif
[2385]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) {
[1732]1572#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1573          value *= value * slack_multiplier;
[1732]1574#else
[2385]1575          value *= value;
[1732]1576#endif
[2385]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;
[1525]1629
[2385]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    }
[1525]1642
[2385]1643    // columns
1644    weight = weights_;
[1525]1645
[2385]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
[225]1665#ifdef SOME_DEBUG_1
[2385]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    }
[225]1675#endif
[2385]1676    updates->setNumElements(0);
1677    spareColumn1->setNumElements(0);
1678  }
[225]1679}
1680// Update djs, weights for Steepest
[2385]1681void ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector *updates,
1682  CoinIndexedVector *spareRow2,
1683  CoinIndexedVector *spareColumn1,
1684  CoinIndexedVector *spareColumn2)
[225]1685{
[2385]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);
[1197]1703
[2385]1704  // put row of tableau in rowArray and columnArray
1705  model_->clpMatrix()->transposeTimes(model_, -1.0,
1706    updates, spareColumn2, spareColumn1);
[2235]1707#ifdef CLP_USER_DRIVEN
[2385]1708  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
[2235]1709#endif
[2385]1710  // normal
1711  for (iSection = 0; iSection < 2; iSection++) {
[517]1712
[2385]1713    reducedCost = model_->djRegion(iSection);
1714    int addSequence;
1715#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1716    double slack_multiplier;
[1732]1717#endif
[1197]1718
[2385]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;
[1732]1726#endif
[2385]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;
[1732]1734#endif
[2385]1735    }
[1525]1736
[2385]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);
[1525]1744
[2385]1745      switch (status) {
[1525]1746
[2385]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) {
[1732]1768#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1769          value *= value * slack_multiplier;
[1732]1770#else
[2385]1771          value *= value;
[1732]1772#endif
[2385]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) {
[1732]1785#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
[2385]1786          value *= value * slack_multiplier;
[1732]1787#else
[2385]1788          value *= value;
[1732]1789#endif
[2385]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;
[1525]1869          }
[2385]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      }
[2235]1875#endif
[2385]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;
[1525]1885
[2385]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);
[1525]1891          } else {
[2385]1892            // exact
1893            thisWeight = referenceIn * pivotSquared;
1894            if (reference(iSequence + numberColumns))
1895              thisWeight += 1.0;
1896            thisWeight = CoinMax(thisWeight, TRY_NORM);
[1525]1897          }
[2385]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    }
[1525]1907
[2385]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    }
[1525]1925
[2385]1926    // columns
1927    weight = weights_;
[1525]1928
[2385]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
[225]1954#ifdef SOME_DEBUG_1
[2385]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    }
[225]1964#endif
[2385]1965  }
1966  updates->setNumElements(0);
1967  spareColumn1->setNumElements(0);
[225]1968}
[517]1969// Updates two arrays for steepest
[2385]1970int ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
1971  const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
1972  CoinIndexedVector *spare,
1973  double scaleFactor)
[517]1974{
[2385]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_;
[1525]2010
[2385]2011    int number = dj1->getNumElements();
2012    const int *index = dj1->getIndices();
2013    double *updateBy = dj1->denseVector();
2014    double *updateBy2 = dj2->denseVector();
[1525]2015
[2385]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);
[1525]2027
[2385]2028      if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
2029        thisWeight = weight[iSequence];
2030        pivot = value2 * scaleFactor;
2031        pivotSquared = pivot * pivot;
[1525]2032
[2385]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);
[517]2044          }
[2385]2045        }
2046        weight[iSequence] = thisWeight;
2047      }
2048    }
2049  }
2050  dj2->setNumElements(0);
2051  return returnCode;
[517]2052}
[225]2053// Update weights for Devex
[2385]2054void ClpPrimalColumnSteepest::justDevex(CoinIndexedVector *updates,
2055  CoinIndexedVector *spareRow2,
2056  CoinIndexedVector *spareColumn1,
2057  CoinIndexedVector *spareColumn2)
[225]2058{
[2385]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;
[1525]2106
[2385]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  }
[1525]2120
[2385]2121  // columns
2122  weight = weights_;
[1525]2123
[2385]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
[225]2144#ifdef SOME_DEBUG_1
[2385]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  }
[225]2154#endif
[2385]2155  updates->setNumElements(0);
2156  spareColumn1->setNumElements(0);
[225]2157}
2158// Update weights for Steepest
[2385]2159void ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector *updates,
2160  CoinIndexedVector *spareRow2,
2161  CoinIndexedVector *spareColumn1,
2162  CoinIndexedVector *spareColumn2)
[225]2163{
[2385]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;
[1525]2215
[2385]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  }
[2235]2241#endif
[2385]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;
[1525]2254
[2385]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  }
[1525]2270
[2385]2271  // columns
2272  weight = weights_;
[1525]2273
[2385]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;
[1525]2287
[2385]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
[225]2309#ifdef SOME_DEBUG_1
[2385]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  }
[225]2319#endif
[2385]2320  updates->setNumElements(0);
2321  spareColumn1->setNumElements(0);
[225]2322}
2323// Returns pivot column, -1 if none
[2385]2324int ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector *updates,
2325  CoinIndexedVector *,
2326  CoinIndexedVector *spareRow2,
2327  CoinIndexedVector *spareColumn1,
2328  CoinIndexedVector *spareColumn2)
[225]2329{
[2385]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();
[2]2347
[2385]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 -
[1525]2357        0 - all exact devex
2358        1 - all steepest
2359        2 - some exact devex
2360        3 - auto some exact devex
2361        4 - devex
2362        5 - dantzig
2363     */
[2385]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  }
[2]2392
[2385]2393  if (anyUpdates > 0) {
2394    model_->factorization()->updateColumnTranspose(spareRow2, updates);
[225]2395
[2385]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++) {
[2]2401
[2385]2402      reducedCost = model_->djRegion(iSection);
2403      int addSequence;
[1197]2404
[2385]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()) {
[1197]2417
[2385]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);
[1525]2424
[2385]2425          switch (status) {
[1525]2426
[2385]2427          case ClpSimplex::basic:
2428            infeasible_->zero(iSequence + addSequence);
2429          case ClpSimplex::isFixed:
2430            break;
2431          case ClpSimplex::isFree:
2432          case ClpSimplex::superBasic:
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;
2445          case ClpSimplex::atUpperBound:
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;
2456          case ClpSimplex::atLowerBound:
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            }
[1525]2466          }
[2385]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);
[1525]2477
[2385]2478          switch (status) {
[1525]2479
2480          case ClpSimplex::basic:
[2385]2481            infeasible_->zero(iSequence + addSequence);
[1525]2482          case ClpSimplex::isFixed:
[2385]2483            break;
[1525]2484          case ClpSimplex::isFree:
2485          case ClpSimplex::superBasic:
[2385]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;
[1525]2498          case ClpSimplex::atUpperBound:
[2385]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;
[1525]2519          case ClpSimplex::atLowerBound:
[2385]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            }
[1525]2539          }
[2385]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);
[1525]2559
[2385]2560    switch (status) {
[1525]2561
[2385]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  }
[1525]2602
[2385]2603  // If in quadratic re-compute all
2604  if (model_->algorithm() == 2) {
2605    for (iSection = 0; iSection < 2; iSection++) {
[1525]2606
[2385]2607      reducedCost = model_->djRegion(iSection);
2608      int addSequence;
2609      int iSequence;
[1525]2610
[2385]2611      if (!iSection) {
2612        number = model_->numberRows();
2613        addSequence = model_->numberColumns();
2614      } else {
2615        number = model_->numberColumns();
2616        addSequence = 0;
2617      }
[1525]2618
[2385]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);
[1525]2623
[2385]2624          switch (status) {
[1525]2625
[2385]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            }
[1525]2665          }
[2385]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            }
[1525]2735          }
[2385]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;
[1525]2871
[2385]2872    if (switchType < 4) {
2873      // Exact
2874      // now update weight update array
2875      model_->factorization()->updateColumnTranspose(spareRow2,
2876        alternateWeights_);
[2235]2877#ifdef ALT_UPDATE_WEIGHTS
[2385]2878      abort();
[2235]2879#endif
[2385]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;
[1525]2888
[2385]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);
2894          } else {
2895            // exact
2896            thisWeight = referenceIn * pivotSquared;
2897            if (reference(iSequence + numberColumns))
2898              thisWeight += 1.0;
2899            thisWeight = CoinMax(thisWeight, TRY_NORM);
[1525]2900          }
[2385]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    }
[1525]2919
[2385]2920    // columns
2921    weight = weights_;
[1525]2922
[2385]2923    scaleFactor = -scaleFactor;
[1525]2924
[2385]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;
[1525]2943
[2385]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);
2949          } else {
2950            // exact
2951            thisWeight = referenceIn * pivotSquared;
2952            if (reference(iSequence))
2953              thisWeight += 1.0;
2954            thisWeight = CoinMax(thisWeight, TRY_NORM);
[1525]2955          }
[2385]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
[2]2979#ifdef SOME_DEBUG_1
[2385]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    }
[2]2989#endif
[2385]2990    updates->setNumElements(0);
2991    spareColumn1->setNumElements(0);
2992  }
[2]2993
[2385]2994  // update of duals finished - now do pricing
[2]2995
[2385]2996  double bestDj = 1.0e-30;
2997  int bestSequence = -1;
[2]2998
[2385]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  }
[225]3012#ifdef CLP_DEBUG
[2385]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);
[225]3017#endif
[2385]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
[225]3025
[2385]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            }
[1525]3054          }
[2385]3055        }
3056        numberWanted--;
3057        if (!numberWanted)
3058          break;
3059      }
3060    } else {
3061      // Dantzig
3062      for (i = start[2 * iPass]; i < end; i++) {
3063        iSequence = index[i];
3064        double value = infeas[iSequence];
3065        if (value > tolerance) {
3066          if (value > bestDj) {
3067            // check flagged variable and correct dj
3068            if (!model_->flagged(iSequence)) {
3069              bestDj = value;
3070              bestSequence = iSequence;
3071            } else {
3072              // just to make sure we don't exit before got something
3073              numberWanted++;
3074            }
3075          }
3076        }
3077        numberWanted--;
3078        if (!numberWanted)
3079          break;
3080      }
3081    }
3082    if (!numberWanted)
3083      break;
3084  }
3085  if (sequenceOut >= 0) {
3086    infeas[sequenceOut] = saveOutInfeasibility;
3087  }
3088  /*if (model_->numberIterations()%100==0)
[1525]3089       printf("%d best %g\n",bestSequence,bestDj);*/
[2385]3090  reducedCost = model_->djRegion();
3091  model_->clpMatrix()->setSavedBestSequence(bestSequence);
3092  if (bestSequence >= 0)
3093    model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
[1525]3094
[2]3095#ifdef CLP_DEBUG
[2385]3096  if (bestSequence >= 0) {
3097    if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
3098      assert(model_->reducedCost(bestSequence) < 0.0);
3099    if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound)
3100      assert(model_->reducedCost(bestSequence) > 0.0);
3101  }
[2]3102#endif
[2385]3103  return bestSequence;
[2]3104}
[618]3105// Called when maximum pivots changes
[2385]3106void ClpPrimalColumnSteepest::maximumPivotsChanged()
[618]3107{
[2385]3108  if (alternateWeights_ && alternateWeights_->capacity() != model_->numberRows() + model_->factorization()->maximumPivots()) {
3109    delete alternateWeights_;
3110    alternateWeights_ = new CoinIndexedVector();
3111    // enough space so can use it for factorization
3112    alternateWeights_->reserve(model_->numberRows() + model_->factorization()->maximumPivots());
3113  }
[618]3114}
[2385]3115void ClpPrimalColumnSteepest::redoInfeasibilities()
[2235]3116{
[2385]3117  double *COIN_RESTRICT infeas = infeasible_->denseVector();
3118  int *COIN_RESTRICT index = infeasible_->getIndices();
[2235]3119  double tolerance = model_->currentDualTolerance();
3120  // we can't really trust infeasibilities if there is dual error
3121  // this coding has to mimic coding in checkDualSolution
3122  double error = CoinMin(1.0e-2, model_->largestDualError());
3123  // allow tolerance at least slightly bigger than standard
[2385]3124  tolerance = tolerance + error;
[2235]3125  // reverse sign so test is cleaner
[2385]3126  tolerance = -tolerance;
[2235]3127  int number = model_->numberRows() + model_->numberColumns();
[2385]3128  int numberNonZero = 0;
3129  const double *COIN_RESTRICT reducedCost = model_->djRegion();
3130  const unsigned char *COIN_RESTRICT status = model_->statusArray();
3131  for (int iSequence = 0; iSequence < number; iSequence++) {
3132    unsigned char thisStatus = status[iSequence] & 7;
[2235]3133    double value = reducedCost[iSequence];
[2385]3134    infeas[iSequence] = 0.0;
3135    if (thisStatus == 3) {
3136    } else if ((thisStatus & 1) != 0) {
[2235]3137      // basic or fixed
[2385]3138      value = 0.0;
3139    } else if (thisStatus == 2) {
3140      value = -value;
[2235]3141    } else {
3142      // free or superbasic
3143      if (fabs(value) > FREE_ACCEPT * -tolerance) {
[2385]3144        // we are going to bias towards free (but only if reasonable)
3145        value = -fabs(value) * FREE_BIAS;
[2235]3146      } else {
[2385]3147        value = 0.0;
[2235]3148      }
3149    }
3150    if (value < tolerance) {
3151      // store square in list
3152      infeas[iSequence] = value * value;
3153      index[numberNonZero++] = iSequence;
3154    }
3155  }
3156  infeasible_->setNumElements(numberNonZero);
3157  infeasibilitiesState_ = 0;
3158}
[1525]3159/*
[2]3160   1) before factorization
3161   2) after factorization
3162   3) just redo infeasibilities
3163   4) restore weights
[1367]3164   5) at end of values pass (so need initialization)
[2]3165*/
[2385]3166void ClpPrimalColumnSteepest::saveWeights(ClpSimplex *model, int mode)
[2]3167{
[2385]3168  model_ = model;
3169  if (mode == 6) {
3170    // If incoming weight is 1.0 then return else as 5
3171    int sequenceIn = model_->sequenceIn();
3172    assert(sequenceIn >= 0 && sequenceIn < model_->numberRows() + model_->numberColumns());
[2444]3173    // possible weights_ was never set up // assert(weights_);
3174    if (weights_ && weights_[sequenceIn] == (mode_ != 1) ? 1.0 : 1.0 + ADD_ONE)
[2385]3175      return;
3176    else
3177      mode = 5;
3178  }
3179  if (mode_ == 4 || mode_ == 5) {
3180    if (mode == 1 && !weights_)
3181      numberSwitched_ = 0; // Reset
3182  }
3183  // alternateWeights_ is defined as indexed but is treated oddly
3184  // at times
3185  int numberRows = model_->numberRows();
3186  int numberColumns = model_->numberColumns();
3187  const int *pivotVariable = model_->pivotVariable();
3188  bool doInfeasibilities = true;
3189  if (mode == 1) {
3190    if (!model_->numberIterations())
3191      pivotSequence_ = -1;
[2235]3192#if ABOCA_LITE
[2385]3193    int numberThreads = abcState();
3194    if (numberThreads && mode_ > 1)
3195      mode_ = 0; // force exact
[2235]3196#endif
[2385]3197    if (weights_) {
3198      // Check if size has changed
3199      if (infeasible_->capacity() == numberRows + numberColumns && alternateWeights_->capacity() == numberRows + model_->factorization()->maximumPivots()) {
3200        //alternateWeights_->clear();
3201        if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
3202#if ALT_UPDATE_WEIGHTS != 2
3203          // save pivot order
3204          CoinMemcpyN(pivotVariable,
3205            numberRows, alternateWeights_->getIndices());
[2235]3206#endif
[2385]3207          // change from pivot row number to sequence number
3208          pivotSequence_ = pivotVariable[pivotSequence_];
3209        } else {
3210          pivotSequence_ = -1;
3211        }
3212        state_ = 1;
3213      } else {
3214        // size has changed - clear everything
3215        delete[] weights_;
3216        weights_ = NULL;
3217        delete infeasible_;
3218        infeasible_ = NULL;
3219        delete alternateWeights_;
3220        alternateWeights_ = NULL;
3221        delete[] savedWeights_;
3222        savedWeights_ = NULL;
3223        delete[] reference_;
3224        reference_ = NULL;
3225        state_ = -1;
3226        pivotSequence_ = -1;
3227      }
3228    }
3229  } else if (mode == 2 || mode == 4 || mode == 5) {
3230    // restore
3231    if (!weights_ || state_ == -1 || mode == 5) {
3232      // Partial is only allowed with certain types of matrix
3233      if ((mode_ != 4 && mode_ != 5) || numberSwitched_ || !model_->clpMatrix()->canDoPartialPricing()) {
3234        // initialize weights
3235        delete[] weights_;
3236        delete alternateWeights_;
3237        weights_ = new double[numberRows + numberColumns];
3238        alternateWeights_ = new CoinIndexedVector();
3239        // enough space so can use it for factorization
3240        alternateWeights_->reserve(numberRows + model_->factorization()->maximumPivots());
3241        initializeWeights();
3242        // create saved weights
3243        delete[] savedWeights_;
3244        savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
3245        // just do initialization
3246        mode = 3;
3247      } else {
3248        // Partial pricing
3249        // use region as somewhere to save non-fixed slacks
3250        // set up infeasibilities
3251        if (!infeasible_) {
3252          infeasible_ = new CoinIndexedVector();
3253          infeasible_->reserve(numberColumns + numberRows);
3254        }
3255        infeasible_->clear();
3256        int number = model_->numberRows() + model_->numberColumns();
3257        int iSequence;
3258        int numberLook = 0;
3259        int *which = infeasible_->getIndices();
3260        for (iSequence = model_->numberColumns(); iSequence < number; iSequence++) {
3261          ClpSimplex::Status status = model_->getStatus(iSequence);
3262          if (status != ClpSimplex::isFixed)
3263            which[numberLook++] = iSequence;
3264        }
3265        infeasible_->setNumElements(numberLook);
3266        doInfeasibilities = false;
3267      }
3268      savedPivotSequence_ = -2;
3269      savedSequenceOut_ = -2;
3270      if (pivotSequence_ < 0 || pivotSequence_ >= numberRows + numberColumns)
3271        pivotSequence_ = -1;
[1525]3272
[2385]3273    } else {
3274      if (mode != 4) {
3275        // save
3276        CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
3277        savedPivotSequence_ = pivotSequence_;
3278        savedSequenceOut_ = model_->sequenceOut();
3279      } else {
3280        // restore
3281        CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
3282        // was - but I think should not be
3283        //pivotSequence_= savedPivotSequence_;
3284        //model_->setSequenceOut(savedSequenceOut_);
3285        pivotSequence_ = -1;
3286        model_->setSequenceOut(-1);
3287        // indices are wrong so clear by hand
3288        //alternateWeights_->clear();
3289        CoinZeroN(alternateWeights_->denseVector(),
3290          alternateWeights_->capacity());
3291        alternateWeights_->setNumElements(0);
3292      }
3293    }
3294    state_ = 0;
3295    // set up infeasibilities
3296    if (!infeasible_) {
3297      infeasible_ = new CoinIndexedVector();
3298      infeasible_->reserve(numberColumns + numberRows);
3299    }
3300  }
3301  if (mode >= 2 && mode != 5) {
3302    if (mode != 3) {
3303      if (pivotSequence_ >= 0) {
3304#if ALT_UPDATE_WEIGHTS != 2
3305        // restore pivot row
3306        int iRow;
3307        // permute alternateWeights
3308        double *temp = model_->rowArray(3)->denseVector();
3309        ;
3310        double *work = alternateWeights_->denseVector();
3311        int *savePivotOrder = model_->rowArray(3)->getIndices();
3312        int *oldPivotOrder = alternateWeights_->getIndices();
3313        for (iRow = 0; iRow < numberRows; iRow++) {
3314          int iPivot = oldPivotOrder[iRow];
3315          temp[iPivot] = work[iRow];
3316          savePivotOrder[iRow] = iPivot;
3317        }
3318        int number = 0;
3319        int found = -1;
3320        int *which = oldPivotOrder;
3321        // find pivot row and re-create alternateWeights
3322        for (iRow = 0; iRow < numberRows; iRow++) {
3323          int iPivot = pivotVariable[iRow];
3324          if (iPivot == pivotSequence_)
3325            found = iRow;
3326          work[iRow] = temp[iPivot];
3327          if (work[iRow])
3328            which[number++] = iRow;
3329        }
3330        alternateWeights_->setNumElements(number);
[50]3331#ifdef CLP_DEBUG
[2385]3332        // Can happen but I should clean up
3333        assert(found >= 0);
[50]3334#endif
[2385]3335        pivotSequence_ = found;
3336        for (iRow = 0; iRow < numberRows; iRow++) {
3337          int iPivot = savePivotOrder[iRow];
3338          temp[iPivot] = 0.0;
3339        }
[2235]3340#else
[2385]3341        for (int iRow = 0; iRow < numberRows; iRow++) {
3342          int iPivot = pivotVariable[iRow];
3343          if (iPivot == pivotSequence_) {
3344            pivotSequence_ = iRow;
3345            break;
3346          }
3347        }
[2235]3348#endif
[2385]3349      } else {
3350        // Just clean up
3351        /* If this happens when alternateWeights_ is
[2281]3352                       in "save" mode then alternateWeights_->clear()
3353                       is disastrous.
3354                       As will be fairly dense anyway and this
3355                       rarely happens just zero out */
[2385]3356        if (alternateWeights_ && alternateWeights_->getNumElements()) {
3357          //alternateWeights_->clear();
3358          CoinZeroN(alternateWeights_->denseVector(),
3359            alternateWeights_->capacity());
3360          alternateWeights_->setNumElements(0);
3361        }
3362      }
3363    }
3364    // Save size of factorization
3365    if (!model_->factorization()->pivots())
3366      sizeFactorization_ = model_->factorization()->numberElements();
3367    if (!doInfeasibilities)
3368      return; // don't disturb infeasibilities
3369    double *COIN_RESTRICT infeas = infeasible_->denseVector();
3370    int *COIN_RESTRICT index = infeasible_->getIndices();
3371    int numberNonZero = 0;
3372    infeasibilitiesState_ = 0;
3373    double tolerance = model_->currentDualTolerance();
3374    int number = model_->numberRows() + model_->numberColumns();
3375    int iSequence;
[225]3376
[2385]3377    double *reducedCost = model_->djRegion();
3378    const double *lower = model_->lowerRegion();
3379    const double *upper = model_->upperRegion();
3380    const double *solution = model_->solutionRegion();
3381    double primalTolerance = model_->currentPrimalTolerance();
[225]3382
[2385]3383    if (!model_->nonLinearCost()->lookBothWays()) {
3384      const unsigned char *COIN_RESTRICT status = model_->statusArray();
3385#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
3386      for (iSequence = 0; iSequence < number; iSequence++) {
3387        double value = reducedCost[iSequence];
3388        infeas[iSequence] = 0.0;
3389        unsigned char thisStatus = status[iSequence] & 7;
3390        if (thisStatus == 3) {
3391        } else if ((thisStatus & 1) != 0) {
3392          // basic or fixed
3393          value = 0.0;
3394        } else if (thisStatus == 2) {
3395          value = -value;
3396        } else {
3397          // free or superbasic
3398          if (fabs(value) > FREE_ACCEPT * tolerance) {
3399            // we are going to bias towards free (but only if reasonable)
3400            value = -fabs(value) * FREE_BIAS;
3401          } else {
3402            value = 0.0;
3403          }
3404        }
3405        if (value < -tolerance) {
3406          infeas[iSequence] = value * value;
3407          index[numberNonZero++] = iSequence;
3408        }
3409      }
[1732]3410#else
[2385]3411      // Columns
3412      int numberColumns = model_->numberColumns();
3413      for (iSequence = 0; iSequence < numberColumns; iSequence++) {
3414        infeas[iSequence] = 0.0;
3415        double value = reducedCost[iSequence];
3416        unsigned char thisStatus = status[iSequence] & 7;
3417        if (thisStatus == 3) {
3418        } else if ((thisStatus & 1) != 0) {
3419          // basic or fixed
3420          value = 0.0;
3421        } else if (thisStatus == 2) {
3422          value = -value;
3423        } else {
3424          // free or superbasic
3425          if (fabs(value) > FREE_ACCEPT * tolerance) {
3426            // check hasn't slipped through
3427            if (solution[iSequence] < lower[iSequence] + primalTolerance) {
3428              model_->setStatus(iSequence, ClpSimplex::atLowerBound);
3429            } else if (solution[iSequence] > upper[iSequence] - primalTolerance) {
3430              model_->setStatus(iSequence, ClpSimplex::atUpperBound);
3431              value = -value;
3432            } else {
3433              // we are going to bias towards free (but only if reasonable)
3434              value = -fabs(value) * FREE_BIAS;
3435            }
3436          } else {
3437            value = 0.0;
3438          }
3439        }
3440        if (value < -tolerance) {
3441          infeas[iSequence] = value * value;
3442          index[numberNonZero++] = iSequence;
3443        }
3444      }
3445      // Rows
3446      for (; iSequence < number; iSequence++) {
3447        double value = reducedCost[iSequence];
3448        infeas[iSequence] = 0.0;
3449        unsigned char thisStatus = status[iSequence] & 7;
3450        if (thisStatus == 3) {
3451        } else if ((thisStatus & 1) != 0) {
3452          // basic or fixed
3453          value = 0.0;
3454        } else if (thisStatus == 2) {
3455          value = -value;
3456        } else {
3457          // free or superbasic
3458          if (fabs(value) > FREE_ACCEPT * tolerance) {
3459            // we are going to bias towards free (but only if reasonable)
3460            value = -fabs(value) * FREE_BIAS;
3461          } else {
3462            value = 0.0;
3463          }
3464        }
3465        if (value < -tolerance) {
3466          infeas[iSequence] = value * value;
3467          index[numberNonZero++] = iSequence;
3468        }
3469      }
[1732]3470#endif
[2385]3471      infeasible_->setNumElements(numberNonZero);
3472    } else {
3473      ClpNonLinearCost *nonLinear = model_->nonLinearCost();
3474      infeasible_->clear();
3475      // can go both ways
3476      for (iSequence = 0; iSequence < number; iSequence++) {
3477        double value = reducedCost[iSequence];
3478        ClpSimplex::Status status = model_->getStatus(iSequence);
[1525]3479
[2385]3480        switch (status) {
[1525]3481
[2385]3482        case ClpSimplex::basic:
3483        case ClpSimplex::isFixed:
3484          break;
3485        case ClpSimplex::isFree:
3486        case ClpSimplex::superBasic:
3487          if (fabs(value) > FREE_ACCEPT * tolerance) {
3488            // we are going to bias towards free (but only if reasonable)
3489            value *= FREE_BIAS;
3490            // store square in list
3491            infeasible_->quickAdd(iSequence, value * value);
[1525]3492          }
[2385]3493          break;
3494        case ClpSimplex::atUpperBound:
3495          if (value > tolerance) {
3496            infeasible_->quickAdd(iSequence, value * value);
3497          } else {
3498            // look other way - change up should be negative
3499            value -= nonLinear->changeUpInCost(iSequence);
3500            if (value < -tolerance) {
3501              // store square in list
3502              infeasible_->quickAdd(iSequence, value * value);
3503            }
3504          }
3505          break;
3506        case ClpSimplex::atLowerBound:
3507          if (value < -tolerance) {
3508            infeasible_->quickAdd(iSequence, value * value);
3509          } else {
3510            // look other way - change down should be positive
3511            value -= nonLinear->changeDownInCost(iSequence);
3512            if (value > tolerance) {
3513              // store square in list
3514              infeasible_->quickAdd(iSequence, value * value);
3515            }
3516          }
3517        }
3518      }