source: stable/1.17/Clp/src/ClpPrimalColumnSteepest.cpp @ 2431

Last change on this file since 2431 was 2431, checked in by stefan, 3 months ago

sync with trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 145.8 KB
Line 
1/* $Id: ClpPrimalColumnSteepest.cpp 2431 2019-03-15 15:56:51Z stefan $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#if ABOCA_LITE
8// 1 is not owner of abcState_
9#define ABCSTATE_LITE 1
10#endif
11//#define FAKE_CILK
12
13#include "ClpSimplex.hpp"
14#include "ClpPrimalColumnSteepest.hpp"
15#include "CoinIndexedVector.hpp"
16#include "ClpFactorization.hpp"
17#include "ClpNonLinearCost.hpp"
18#include "ClpMessage.hpp"
19#include "ClpEventHandler.hpp"
20#include "ClpPackedMatrix.hpp"
21#include "CoinHelperFunctions.hpp"
22#include <stdio.h>
23#if ABOCA_LITE
24#undef ALT_UPDATE_WEIGHTS
25#define ALT_UPDATE_WEIGHTS 2
26#endif
27#if ALT_UPDATE_WEIGHTS == 1
28extern CoinIndexedVector *altVector[3];
29#endif
30static void debug1(int iSequence, double value, double weight)
31{
32  printf("xx %d inf %.20g wt %.20g\n",
33    iSequence, value, weight);
34}
35//#define CLP_DEBUG
36//#############################################################################
37// Constructors / Destructor / Assignment
38//#############################################################################
39
40//-------------------------------------------------------------------
41// Default Constructor
42//-------------------------------------------------------------------
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)
60{
61  type_ = 2 + 64 * mode;
62}
63//-------------------------------------------------------------------
64// Copy constructor
65//-------------------------------------------------------------------
66ClpPrimalColumnSteepest::ClpPrimalColumnSteepest(const ClpPrimalColumnSteepest &rhs)
67  : ClpPrimalColumnPivot(rhs)
68{
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  }
114}
115
116//-------------------------------------------------------------------
117// Destructor
118//-------------------------------------------------------------------
119ClpPrimalColumnSteepest::~ClpPrimalColumnSteepest()
120{
121  delete[] weights_;
122  delete infeasible_;
123  delete alternateWeights_;
124  delete[] savedWeights_;
125  delete[] reference_;
126}
127
128//----------------------------------------------------------------
129// Assignment operator
130//-------------------------------------------------------------------
131ClpPrimalColumnSteepest &
132ClpPrimalColumnSteepest::operator=(const ClpPrimalColumnSteepest &rhs)
133{
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;
180}
181// These have to match ClpPackedMatrix version
182#define TRY_NORM 1.0e-4
183#define ADD_ONE 1.0
184static void
185pivotColumnBit(clpTempInfo &info)
186{
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;
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) {
201        // check flagged variable and correct dj
202        if ((status[iSequence] & 64) == 0) {
203          bestDj = value / weight;
204          bestSequence = iSequence;
205        }
206      }
207    }
208  }
209  info.primalRatio = bestDj;
210  info.numberAdded = bestSequence;
211}
212// Returns pivot column, -1 if none
213/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
214        that is just +-weight where a feasibility changed.  It also has
215        reduced cost from last iteration in pivot row*/
216int ClpPrimalColumnSteepest::pivotColumn(CoinIndexedVector *updates,
217  CoinIndexedVector *spareRow1,
218  CoinIndexedVector *spareRow2,
219  CoinIndexedVector *spareColumn1,
220  CoinIndexedVector *spareColumn2)
221{
222  assert(model_);
223  if (model_->nonLinearCost()->lookBothWays() || model_->algorithm() == 2) {
224    // Do old way
225    updates->expand();
226    return pivotColumnOldMethod(updates, spareRow1, spareRow2,
227      spareColumn1, spareColumn2);
228  }
229  int number = 0;
230  int *index;
231  double tolerance = model_->currentDualTolerance();
232  // we can't really trust infeasibilities if there is dual error
233  // this coding has to mimic coding in checkDualSolution
234  double error = CoinMin(1.0e-2, model_->largestDualError());
235  // allow tolerance at least slightly bigger than standard
236  tolerance = tolerance + error;
237  int pivotRow = model_->pivotRow();
238  int anyUpdates;
239  double *infeas = infeasible_->denseVector();
240
241  // Local copy of mode so can decide what to do
242  int switchType;
243  if (mode_ == 4)
244    switchType = 5 - numberSwitched_;
245  else if (mode_ >= 10)
246    switchType = 3;
247  else
248    switchType = mode_;
249    /* switchType -
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     */
258    // Look at gub
259#if 1
260  model_->clpMatrix()->dualExpanded(model_, updates, NULL, 4);
261#else
262  updates->clear();
263  model_->computeDuals(NULL);
264#endif
265  if (updates->getNumElements() > 1) {
266    // would have to have two goes for devex, three for steepest
267    anyUpdates = 2;
268  } else if (updates->getNumElements()) {
269    if (updates->getIndices()[0] == pivotRow && fabs(updates->denseVector()[0]) > 1.0e-6) {
270      // reasonable size
271      anyUpdates = 1;
272      //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
273      //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
274    } else {
275      // too small
276      anyUpdates = 2;
277    }
278  } else if (pivotSequence_ >= 0) {
279    // just after re-factorization
280    anyUpdates = -1;
281  } else {
282    // sub flip - nothing to do
283    anyUpdates = 0;
284  }
285  int sequenceOut = model_->sequenceOut();
286  if (switchType == 5) {
287    // If known matrix then we will do partial pricing
288    if (model_->clpMatrix()->canDoPartialPricing()) {
289      pivotSequence_ = -1;
290      pivotRow = -1;
291      // See if to switch
292      int numberRows = model_->numberRows();
293      int numberWanted = 10;
294      int numberColumns = model_->numberColumns();
295      int numberHiddenRows = model_->clpMatrix()->hiddenRows();
296      double ratio = static_cast< double >(sizeFactorization_ + numberHiddenRows) / static_cast< double >(numberRows + 2 * numberHiddenRows);
297      // Number of dual infeasibilities at last invert
298      int numberDual = model_->numberDualInfeasibilities();
299      int numberLook = CoinMin(numberDual, numberColumns / 10);
300      if (ratio < 1.0) {
301        numberWanted = 100;
302        numberLook /= 20;
303        numberWanted = CoinMax(numberWanted, numberLook);
304      } else if (ratio < 3.0) {
305        numberWanted = 500;
306        numberLook /= 15;
307        numberWanted = CoinMax(numberWanted, numberLook);
308      } else if (ratio < 4.0 || mode_ == 5) {
309        numberWanted = 1000;
310        numberLook /= 10;
311        numberWanted = CoinMax(numberWanted, numberLook);
312      } else if (mode_ != 5) {
313        switchType = 4;
314        // initialize
315        numberSwitched_++;
316        // Make sure will re-do
317        delete[] weights_;
318        weights_ = NULL;
319        model_->computeDuals(NULL);
320        saveWeights(model_, 4);
321        anyUpdates = 0;
322        COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
323      }
324      if (switchType == 5) {
325        numberLook *= 5; // needs tuning for gub
326        if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
327          COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
328            sizeFactorization_, ratio, numberWanted, numberLook));
329        }
330        // Update duals and row djs
331        // Do partial pricing
332        return partialPricing(updates, spareRow2,
333          numberWanted, numberLook);
334      }
335    }
336  }
337  int bestSequence = -1;
338  model_->spareIntArray_[3] = -3;
339  if (switchType == 5) {
340    if (anyUpdates > 0) {
341      justDjs(updates, spareRow2,
342        spareColumn1, spareColumn2);
343    }
344  } else if (anyUpdates == 1) {
345    if (switchType < 4) {
346      // exact etc when can use dj
347      djsAndSteepest(updates, spareRow2,
348        spareColumn1, spareColumn2);
349      if (model_->spareIntArray_[3] > -2) {
350        bestSequence = model_->spareIntArray_[3];
351        infeasibilitiesState_ = 2;
352      } else if (model_->spareIntArray_[3] == -2) {
353        redoInfeasibilities();
354      }
355    } else {
356      // devex etc when can use dj
357      djsAndDevex(updates, spareRow2,
358        spareColumn1, spareColumn2);
359    }
360  } else if (anyUpdates == -1) {
361    if (switchType < 4) {
362      // exact etc when djs okay
363      justSteepest(updates, spareRow2,
364        spareColumn1, spareColumn2);
365    } else {
366      // devex etc when djs okay
367      justDevex(updates, spareRow2,
368        spareColumn1, spareColumn2);
369    }
370  } else if (anyUpdates == 2) {
371    if (switchType < 4) {
372      // exact etc when have to use pivot
373      djsAndSteepest2(updates, spareRow2,
374        spareColumn1, spareColumn2);
375    } else {
376      // devex etc when have to use pivot
377      djsAndDevex2(updates, spareRow2,
378        spareColumn1, spareColumn2);
379    }
380  }
381  // everything may have been done if vector copy
382  if (infeasibilitiesState_ == 2) {
383    // all done
384    infeasibilitiesState_ = 1;
385    model_->clpMatrix()->setSavedBestSequence(bestSequence);
386    if (bestSequence >= 0)
387      model_->clpMatrix()->setSavedBestDj(model_->djRegion()[bestSequence]);
388    assert(sequenceOut != bestSequence);
389    return bestSequence;
390  } else if (infeasibilitiesState_ == 1) {
391    // need to redo
392    //infeasibilitiesState_ = 0;
393    redoInfeasibilities();
394  }
395
396#ifdef CLP_DEBUG
397  alternateWeights_->checkClear();
398#endif
399  // make sure outgoing from last iteration okay
400  if (sequenceOut >= 0) {
401    ClpSimplex::Status status = model_->getStatus(sequenceOut);
402    double value = model_->reducedCost(sequenceOut);
403
404    switch (status) {
405
406    case ClpSimplex::basic:
407    case ClpSimplex::isFixed:
408      break;
409    case ClpSimplex::isFree:
410    case ClpSimplex::superBasic:
411      if (fabs(value) > FREE_ACCEPT * tolerance) {
412        // we are going to bias towards free (but only if reasonable)
413        value *= FREE_BIAS;
414        // store square in list
415        if (infeas[sequenceOut])
416          infeas[sequenceOut] = value * value; // already there
417        else
418          infeasible_->quickAdd(sequenceOut, value * value);
419      } else {
420        infeasible_->zero(sequenceOut);
421      }
422      break;
423    case ClpSimplex::atUpperBound:
424      if (value > tolerance) {
425        // store square in list
426        if (infeas[sequenceOut])
427          infeas[sequenceOut] = value * value; // already there
428        else
429          infeasible_->quickAdd(sequenceOut, value * value);
430      } else {
431        infeasible_->zero(sequenceOut);
432      }
433      break;
434    case ClpSimplex::atLowerBound:
435      if (value < -tolerance) {
436        // store square in list
437        if (infeas[sequenceOut])
438          infeas[sequenceOut] = value * value; // already there
439        else
440          infeasible_->quickAdd(sequenceOut, value * value);
441      } else {
442        infeasible_->zero(sequenceOut);
443      }
444    }
445  }
446  // update of duals finished - now do pricing
447  // See what sort of pricing
448  int numberWanted = 10;
449  number = infeasible_->getNumElements();
450  int numberColumns = model_->numberColumns();
451  if (switchType == 5) {
452    pivotSequence_ = -1;
453    pivotRow = -1;
454    // See if to switch
455    int numberRows = model_->numberRows();
456    // ratio is done on number of columns here
457    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
458    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
459    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
460    if (ratio < 1.0) {
461      numberWanted = CoinMax(100, number / 200);
462    } else if (ratio < 2.0 - 1.0) {
463      numberWanted = CoinMax(500, number / 40);
464    } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
465      numberWanted = CoinMax(2000, number / 10);
466      numberWanted = CoinMax(numberWanted, numberColumns / 30);
467    } else if (mode_ != 5) {
468      switchType = 4;
469      // initialize
470      numberSwitched_++;
471      // Make sure will re-do
472      delete[] weights_;
473      weights_ = NULL;
474      saveWeights(model_, 4);
475      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
476    }
477    //if (model_->numberIterations()%1000==0)
478    //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
479  }
480  int numberRows = model_->numberRows();
481  // ratio is done on number of rows here
482  double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
483  if (switchType == 4) {
484    // Still in devex mode
485    // Go to steepest if lot of iterations?
486    if (ratio < 5.0) {
487      numberWanted = CoinMax(2000, number / 10);
488      numberWanted = CoinMax(numberWanted, numberColumns / 20);
489    } else if (ratio < 7.0) {
490      numberWanted = CoinMax(2000, number / 5);
491      numberWanted = CoinMax(numberWanted, numberColumns / 10);
492    } else {
493      // we can zero out
494      updates->clear();
495      spareColumn1->clear();
496      switchType = 3;
497      // initialize
498      pivotSequence_ = -1;
499      pivotRow = -1;
500      numberSwitched_++;
501      // Make sure will re-do
502      delete[] weights_;
503      weights_ = NULL;
504      saveWeights(model_, 4);
505      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
506      updates->clear();
507    }
508    if (model_->numberIterations() % 1000 == 0)
509      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
510  }
511  if (switchType < 4) {
512    if (switchType < 2) {
513      numberWanted = COIN_INT_MAX - 1;
514    } else if (switchType == 2) {
515      numberWanted = CoinMax(2000, number / 8);
516    } else {
517      if (ratio < 1.0) {
518        numberWanted = CoinMax(2000, number / 20);
519      } else if (ratio < 5.0) {
520        numberWanted = CoinMax(2000, number / 10);
521        numberWanted = CoinMax(numberWanted, numberColumns / 40);
522      } else if (ratio < 10.0) {
523        numberWanted = CoinMax(2000, number / 8);
524        numberWanted = CoinMax(numberWanted, numberColumns / 20);
525      } else {
526        ratio = number * (ratio / 80.0);
527        if (ratio > number) {
528          numberWanted = number + 1;
529        } else {
530          numberWanted = CoinMax(2000, static_cast< int >(ratio));
531          numberWanted = CoinMax(numberWanted, numberColumns / 10);
532        }
533      }
534    }
535    //if (model_->numberIterations()%1000==0)
536    //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
537    //switchType);
538  }
539
540  double bestDj = 1.0e-30;
541  bestSequence = -1;
542#ifdef CLP_USER_DRIVEN
543  // could go parallel?
544  model_->eventHandler()->event(ClpEventHandler::beforeChooseIncoming);
545#endif
546
547  int i, iSequence;
548  index = infeasible_->getIndices();
549  number = infeasible_->getNumElements();
550  // Re-sort infeasible every 100 pivots
551  // Not a good idea
552  if (0 && model_->factorization()->pivots() > 0 && (model_->factorization()->pivots() % 100) == 0) {
553    int nLook = model_->numberRows() + numberColumns;
554    number = 0;
555    for (i = 0; i < nLook; i++) {
556      if (infeas[i]) {
557        if (fabs(infeas[i]) > COIN_INDEXED_TINY_ELEMENT)
558          index[number++] = i;
559        else
560          infeas[i] = 0.0;
561      }
562    }
563    infeasible_->setNumElements(number);
564  }
565  if (model_->numberIterations() < model_->lastBadIteration() + 200 && model_->factorization()->pivots() > 10) {
566    // we can't really trust infeasibilities if there is dual error
567    double checkTolerance = 1.0e-8;
568    if (model_->largestDualError() > checkTolerance)
569      tolerance *= model_->largestDualError() / checkTolerance;
570    // But cap
571    tolerance = CoinMin(1000.0, tolerance);
572  }
573#ifdef CLP_DEBUG
574  if (model_->numberDualInfeasibilities() == 1)
575    printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
576      model_->largestDualError(), model_, model_->messageHandler(),
577      number);
578#endif
579  // stop last one coming immediately
580  double saveOutInfeasibility = 0.0;
581  if (sequenceOut >= 0) {
582    saveOutInfeasibility = infeas[sequenceOut];
583    infeas[sequenceOut] = 0.0;
584  }
585  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
586    tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
587  tolerance *= tolerance; // as we are using squares
588
589  int iPass;
590  // Setup two passes (unless all)
591  if (mode_ > 1) {
592    int start[4];
593    start[1] = number;
594    start[2] = 0;
595    double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
596    start[0] = static_cast< int >(dstart);
597    start[3] = start[0];
598    for (iPass = 0; iPass < 2; iPass++) {
599      int end = start[2 * iPass + 1];
600      if (switchType < 5) {
601        for (i = start[2 * iPass]; i < end; i++) {
602          iSequence = index[i];
603          double value = infeas[iSequence];
604          double weight = weights_[iSequence];
605          //assert (weight);
606          if (value > tolerance) {
607            //weight=1.0;
608            if (value > bestDj * weight) {
609              // check flagged variable and correct dj
610              if (!model_->flagged(iSequence)) {
611                bestDj = value / weight;
612                bestSequence = iSequence;
613              } else {
614                // just to make sure we don't exit before got something
615                numberWanted++;
616              }
617            }
618            numberWanted--;
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 {
649#if ABOCA_LITE
650    int numberThreads = CoinMax(abcState(), 1);
651#define ABOCA_LITE_MAX ABOCA_LITE
652#else
653    const int numberThreads = 1;
654#define ABOCA_LITE_MAX 1
655#endif
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)
711       printf("%d best %g\n",bestSequence,bestDj);*/
712
713#ifdef CLP_USER_DRIVEN
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));
723#endif
724#ifndef NDEBUG
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  }
732#endif
733#ifdef ALT_UPDATE_WEIGHTSz
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");
743#endif
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
750  return bestSequence;
751}
752// Just update djs
753void ClpPrimalColumnSteepest::justDjs(CoinIndexedVector *updates,
754  CoinIndexedVector *spareRow2,
755  CoinIndexedVector *spareColumn1,
756  CoinIndexedVector *spareColumn2)
757{
758  int iSection, j;
759  int number = 0;
760  int *index;
761  double *updateBy;
762  double *reducedCost;
763  double tolerance = model_->currentDualTolerance();
764  // we can't really trust infeasibilities if there is dual error
765  // this coding has to mimic coding in checkDualSolution
766  double error = CoinMin(1.0e-2, model_->largestDualError());
767  // allow tolerance at least slightly bigger than standard
768  tolerance = tolerance + error;
769  int pivotRow = model_->pivotRow();
770  double *infeas = infeasible_->denseVector();
771  //updates->scanAndPack();
772  model_->factorization()->updateColumnTranspose(spareRow2, updates);
773
774  // put row of tableau in rowArray and columnArray (packed mode)
775  model_->clpMatrix()->transposeTimes(model_, -1.0,
776    updates, spareColumn2, spareColumn1);
777  // normal
778  for (iSection = 0; iSection < 2; iSection++) {
779
780    reducedCost = model_->djRegion(iSection);
781    int addSequence;
782#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
783    double slack_multiplier;
784#endif
785
786    if (!iSection) {
787      number = updates->getNumElements();
788      index = updates->getIndices();
789      updateBy = updates->denseVector();
790      addSequence = model_->numberColumns();
791#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
792      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
793#endif
794    } else {
795      number = spareColumn1->getNumElements();
796      index = spareColumn1->getIndices();
797      updateBy = spareColumn1->denseVector();
798      addSequence = 0;
799#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
800      slack_multiplier = 1.0;
801#endif
802    }
803
804    for (j = 0; j < number; j++) {
805      int iSequence = index[j];
806      double value = reducedCost[iSequence];
807      value -= updateBy[j];
808      updateBy[j] = 0.0;
809      reducedCost[iSequence] = value;
810      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
811
812      switch (status) {
813
814      case ClpSimplex::basic:
815        infeasible_->zero(iSequence + addSequence);
816      case ClpSimplex::isFixed:
817        break;
818      case ClpSimplex::isFree:
819      case ClpSimplex::superBasic:
820        if (fabs(value) > FREE_ACCEPT * tolerance) {
821          // we are going to bias towards free (but only if reasonable)
822          value *= FREE_BIAS;
823          // store square in list
824          if (infeas[iSequence + addSequence])
825            infeas[iSequence + addSequence] = value * value; // already there
826          else
827            infeasible_->quickAdd(iSequence + addSequence, value * value);
828        } else {
829          infeasible_->zero(iSequence + addSequence);
830        }
831        break;
832      case ClpSimplex::atUpperBound:
833        iSequence += addSequence;
834        if (value > tolerance) {
835#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
836          value *= value * slack_multiplier;
837#else
838          value *= value;
839#endif
840          // store square in list
841          if (infeas[iSequence])
842            infeas[iSequence] = value; // already there
843          else
844            infeasible_->quickAdd(iSequence, value);
845        } else {
846          infeasible_->zero(iSequence);
847        }
848        break;
849      case ClpSimplex::atLowerBound:
850        iSequence += addSequence;
851        if (value < -tolerance) {
852#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
853          value *= value * slack_multiplier;
854#else
855          value *= value;
856#endif
857          // store square in list
858          if (infeas[iSequence])
859            infeas[iSequence] = value; // already there
860          else
861            infeasible_->quickAdd(iSequence, value);
862        } else {
863          infeasible_->zero(iSequence);
864        }
865      }
866    }
867  }
868  updates->setNumElements(0);
869  spareColumn1->setNumElements(0);
870  if (pivotRow >= 0) {
871    // make sure infeasibility on incoming is 0.0
872    int sequenceIn = model_->sequenceIn();
873    infeasible_->zero(sequenceIn);
874  }
875}
876// Update djs, weights for Devex
877void ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector *updates,
878  CoinIndexedVector *spareRow2,
879  CoinIndexedVector *spareColumn1,
880  CoinIndexedVector *spareColumn2)
881{
882  int j;
883  int number = 0;
884  int *index;
885  double *updateBy;
886  double *reducedCost;
887  double tolerance = model_->currentDualTolerance();
888  // we can't really trust infeasibilities if there is dual error
889  // this coding has to mimic coding in checkDualSolution
890  double error = CoinMin(1.0e-2, model_->largestDualError());
891  // allow tolerance at least slightly bigger than standard
892  tolerance = tolerance + error;
893  // for weights update we use pivotSequence
894  // unset in case sub flip
895  assert(pivotSequence_ >= 0);
896  assert(model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
897  pivotSequence_ = -1;
898  double *infeas = infeasible_->denseVector();
899  //updates->scanAndPack();
900  model_->factorization()->updateColumnTranspose(spareRow2, updates);
901  // and we can see if reference
902  //double referenceIn = 0.0;
903  int sequenceIn = model_->sequenceIn();
904  //if (mode_ != 1 && reference(sequenceIn))
905  //   referenceIn = 1.0;
906  // save outgoing weight round update
907  double outgoingWeight = 0.0;
908  int sequenceOut = model_->sequenceOut();
909  if (sequenceOut >= 0)
910    outgoingWeight = weights_[sequenceOut];
911
912  double scaleFactor = 1.0 / updates->denseVector()[0]; // as formula is with 1.0
913  // put row of tableau in rowArray and columnArray (packed mode)
914  model_->clpMatrix()->transposeTimes(model_, -1.0,
915    updates, spareColumn2, spareColumn1);
916  // update weights
917  double *weight;
918  int numberColumns = model_->numberColumns();
919  // rows
920  reducedCost = model_->djRegion(0);
921  int addSequence = model_->numberColumns();
922  ;
923
924  number = updates->getNumElements();
925  index = updates->getIndices();
926  updateBy = updates->denseVector();
927  weight = weights_ + numberColumns;
928  // Devex
929  for (j = 0; j < number; j++) {
930    double thisWeight;
931    double pivot;
932    double value3;
933    int iSequence = index[j];
934    double value = reducedCost[iSequence];
935    double value2 = updateBy[j];
936    updateBy[j] = 0.0;
937    value -= value2;
938    reducedCost[iSequence] = value;
939    ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
940
941    switch (status) {
942
943    case ClpSimplex::basic:
944      infeasible_->zero(iSequence + addSequence);
945    case ClpSimplex::isFixed:
946      break;
947    case ClpSimplex::isFree:
948    case ClpSimplex::superBasic:
949      thisWeight = weight[iSequence];
950      // row has -1
951      pivot = value2 * scaleFactor;
952      value3 = pivot * pivot * devex_;
953      if (reference(iSequence + numberColumns))
954        value3 += 1.0;
955      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
956      if (fabs(value) > FREE_ACCEPT * tolerance) {
957        // we are going to bias towards free (but only if reasonable)
958        value *= FREE_BIAS;
959        // store square in list
960        if (infeas[iSequence + addSequence])
961          infeas[iSequence + addSequence] = value * value; // already there
962        else
963          infeasible_->quickAdd(iSequence + addSequence, value * value);
964      } else {
965        infeasible_->zero(iSequence + addSequence);
966      }
967      break;
968    case ClpSimplex::atUpperBound:
969      thisWeight = weight[iSequence];
970      // row has -1
971      pivot = value2 * scaleFactor;
972      value3 = pivot * pivot * devex_;
973      if (reference(iSequence + numberColumns))
974        value3 += 1.0;
975      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
976      iSequence += addSequence;
977      if (value > tolerance) {
978        // store square in list
979#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
980        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
981#else
982        value *= value;
983#endif
984        if (infeas[iSequence])
985          infeas[iSequence] = value; // already there
986        else
987          infeasible_->quickAdd(iSequence, value);
988      } else {
989        infeasible_->zero(iSequence);
990      }
991      break;
992    case ClpSimplex::atLowerBound:
993      thisWeight = weight[iSequence];
994      // row has -1
995      pivot = value2 * scaleFactor;
996      value3 = pivot * pivot * devex_;
997      if (reference(iSequence + numberColumns))
998        value3 += 1.0;
999      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1000      iSequence += addSequence;
1001      if (value < -tolerance) {
1002        // store square in list
1003#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1004        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
1005#else
1006        value *= value;
1007#endif
1008        if (infeas[iSequence])
1009          infeas[iSequence] = value; // already there
1010        else
1011          infeasible_->quickAdd(iSequence, value);
1012      } else {
1013        infeasible_->zero(iSequence);
1014      }
1015    }
1016  }
1017
1018  // columns
1019  weight = weights_;
1020
1021  scaleFactor = -scaleFactor;
1022  reducedCost = model_->djRegion(1);
1023  number = spareColumn1->getNumElements();
1024  index = spareColumn1->getIndices();
1025  updateBy = spareColumn1->denseVector();
1026
1027  // Devex
1028
1029  for (j = 0; j < number; j++) {
1030    double thisWeight;
1031    double pivot;
1032    double value3;
1033    int iSequence = index[j];
1034    double value = reducedCost[iSequence];
1035    double value2 = updateBy[j];
1036    value -= value2;
1037    updateBy[j] = 0.0;
1038    reducedCost[iSequence] = value;
1039    ClpSimplex::Status status = model_->getStatus(iSequence);
1040
1041    switch (status) {
1042
1043    case ClpSimplex::basic:
1044      infeasible_->zero(iSequence);
1045    case ClpSimplex::isFixed:
1046      break;
1047    case ClpSimplex::isFree:
1048    case ClpSimplex::superBasic:
1049      thisWeight = weight[iSequence];
1050      // row has -1
1051      pivot = value2 * scaleFactor;
1052      value3 = pivot * pivot * devex_;
1053      if (reference(iSequence))
1054        value3 += 1.0;
1055      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1056      if (fabs(value) > FREE_ACCEPT * tolerance) {
1057        // we are going to bias towards free (but only if reasonable)
1058        value *= FREE_BIAS;
1059        // store square in list
1060        if (infeas[iSequence])
1061          infeas[iSequence] = value * value; // already there
1062        else
1063          infeasible_->quickAdd(iSequence, value * value);
1064      } else {
1065        infeasible_->zero(iSequence);
1066      }
1067      break;
1068    case ClpSimplex::atUpperBound:
1069      thisWeight = weight[iSequence];
1070      // row has -1
1071      pivot = value2 * scaleFactor;
1072      value3 = pivot * pivot * devex_;
1073      if (reference(iSequence))
1074        value3 += 1.0;
1075      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1076      if (value > tolerance) {
1077        // store square in list
1078        if (infeas[iSequence])
1079          infeas[iSequence] = value * value; // already there
1080        else
1081          infeasible_->quickAdd(iSequence, value * value);
1082      } else {
1083        infeasible_->zero(iSequence);
1084      }
1085      break;
1086    case ClpSimplex::atLowerBound:
1087      thisWeight = weight[iSequence];
1088      // row has -1
1089      pivot = value2 * scaleFactor;
1090      value3 = pivot * pivot * devex_;
1091      if (reference(iSequence))
1092        value3 += 1.0;
1093      weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1094      if (value < -tolerance) {
1095        // store square in list
1096        if (infeas[iSequence])
1097          infeas[iSequence] = value * value; // already there
1098        else
1099          infeasible_->quickAdd(iSequence, value * value);
1100      } else {
1101        infeasible_->zero(iSequence);
1102      }
1103    }
1104  }
1105  // restore outgoing weight
1106  if (sequenceOut >= 0)
1107    weights_[sequenceOut] = outgoingWeight;
1108  // make sure infeasibility on incoming is 0.0
1109  infeasible_->zero(sequenceIn);
1110  spareRow2->setNumElements(0);
1111  //#define SOME_DEBUG_1
1112#ifdef SOME_DEBUG_1
1113  // check for accuracy
1114  int iCheck = 892;
1115  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1116  int numberRows = model_->numberRows();
1117  //int numberColumns = model_->numberColumns();
1118  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1119    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
1120      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1121  }
1122#endif
1123  updates->setNumElements(0);
1124  spareColumn1->setNumElements(0);
1125}
1126// Update djs, weights for Steepest
1127void ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector *updates,
1128  CoinIndexedVector *spareRow2,
1129  CoinIndexedVector *spareColumn1,
1130  CoinIndexedVector *spareColumn2)
1131{
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;
1151#if 1
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  }
1177#endif
1178#else
1179  model_->factorization()->updateTwoColumnsTranspose(spareRow2, updates, alternateWeights_);
1180#endif
1181  // and we can see if reference
1182  int sequenceIn = model_->sequenceIn();
1183  double referenceIn;
1184  if (mode_ != 1) {
1185    if (reference(sequenceIn))
1186      referenceIn = 1.0;
1187    else
1188      referenceIn = 0.0;
1189  } else {
1190    referenceIn = -1.0;
1191  }
1192  // save outgoing weight round update
1193  double outgoingWeight = 0.0;
1194  int sequenceOut = model_->sequenceOut();
1195  if (sequenceOut >= 0)
1196    outgoingWeight = weights_[sequenceOut];
1197  // update row weights here so we can scale alternateWeights_
1198  // update weights
1199  double *weight;
1200  double *other = alternateWeights_->denseVector();
1201  int numberColumns = model_->numberColumns();
1202  // if if (model_->clpMatrix()->canCombine(model_, pi1)) no need to index
1203  // rows
1204  reducedCost = model_->djRegion(0);
1205  int addSequence = model_->numberColumns();
1206  ;
1207
1208  number = updates->getNumElements();
1209  index = updates->getIndices();
1210  updateBy = updates->denseVector();
1211  weight = weights_ + numberColumns;
1212#ifdef CLP_USER_DRIVEN
1213  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
1214#endif
1215
1216  for (j = 0; j < number; j++) {
1217    double thisWeight;
1218    double pivot;
1219    double modification;
1220    double pivotSquared;
1221    int iSequence = index[j];
1222    double value2 = updateBy[j];
1223    ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1224    double value;
1225
1226    switch (status) {
1227
1228    case ClpSimplex::basic:
1229      infeasible_->zero(iSequence + addSequence);
1230      reducedCost[iSequence] = 0.0;
1231    case ClpSimplex::isFixed:
1232      break;
1233    case ClpSimplex::isFree:
1234    case ClpSimplex::superBasic:
1235      value = reducedCost[iSequence] - value2;
1236      modification = other[iSequence];
1237      thisWeight = weight[iSequence];
1238      // row has -1
1239      pivot = value2 * scaleFactor;
1240      pivotSquared = pivot * pivot;
1241
1242      thisWeight += pivotSquared * devex_ + pivot * modification;
1243      reducedCost[iSequence] = value;
1244      if (thisWeight < TRY_NORM) {
1245        if (mode_ == 1) {
1246          // steepest
1247          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1248        } else {
1249          // exact
1250          thisWeight = referenceIn * pivotSquared;
1251          if (reference(iSequence + numberColumns))
1252            thisWeight += 1.0;
1253          thisWeight = CoinMax(thisWeight, TRY_NORM);
1254        }
1255      }
1256      weight[iSequence] = thisWeight;
1257      if (fabs(value) > FREE_ACCEPT * tolerance) {
1258        // we are going to bias towards free (but only if reasonable)
1259        value *= FREE_BIAS;
1260        // store square in list
1261        if (infeas[iSequence + addSequence])
1262          infeas[iSequence + addSequence] = value * value; // already there
1263        else
1264          infeasible_->quickAdd(iSequence + addSequence, value * value);
1265      } else {
1266        infeasible_->zero(iSequence + addSequence);
1267      }
1268      break;
1269    case ClpSimplex::atUpperBound:
1270      value = reducedCost[iSequence] - value2;
1271      modification = other[iSequence];
1272      thisWeight = weight[iSequence];
1273      // row has -1
1274      pivot = value2 * scaleFactor;
1275      pivotSquared = pivot * pivot;
1276
1277      thisWeight += pivotSquared * devex_ + pivot * modification;
1278      reducedCost[iSequence] = value;
1279      if (thisWeight < TRY_NORM) {
1280        if (mode_ == 1) {
1281          // steepest
1282          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1283        } else {
1284          // exact
1285          thisWeight = referenceIn * pivotSquared;
1286          if (reference(iSequence + numberColumns))
1287            thisWeight += 1.0;
1288          thisWeight = CoinMax(thisWeight, TRY_NORM);
1289        }
1290      }
1291      weight[iSequence] = thisWeight;
1292      iSequence += addSequence;
1293      if (value > tolerance) {
1294        // store square in list
1295#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1296        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
1297#else
1298        value *= value;
1299#endif
1300        if (infeas[iSequence])
1301          infeas[iSequence] = value; // already there
1302        else
1303          infeasible_->quickAdd(iSequence, value);
1304      } else {
1305        infeasible_->zero(iSequence);
1306      }
1307      break;
1308    case ClpSimplex::atLowerBound:
1309      value = reducedCost[iSequence] - value2;
1310      modification = other[iSequence];
1311      thisWeight = weight[iSequence];
1312      // row has -1
1313      pivot = value2 * scaleFactor;
1314      pivotSquared = pivot * pivot;
1315
1316      thisWeight += pivotSquared * devex_ + pivot * modification;
1317      reducedCost[iSequence] = value;
1318      if (thisWeight < TRY_NORM) {
1319        if (mode_ == 1) {
1320          // steepest
1321          thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1322        } else {
1323          // exact
1324          thisWeight = referenceIn * pivotSquared;
1325          if (reference(iSequence + numberColumns))
1326            thisWeight += 1.0;
1327          thisWeight = CoinMax(thisWeight, TRY_NORM);
1328        }
1329      }
1330      weight[iSequence] = thisWeight;
1331      iSequence += addSequence;
1332      if (value < -tolerance) {
1333        // store square in list
1334#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1335        value *= value * CLP_PRIMAL_SLACK_MULTIPLIER;
1336#else
1337        value *= value;
1338#endif
1339        if (infeas[iSequence])
1340          infeas[iSequence] = value; // already there
1341        else
1342          infeasible_->quickAdd(iSequence, value);
1343      } else {
1344        infeasible_->zero(iSequence);
1345      }
1346    }
1347  }
1348#ifdef CLP_USER_DRIVEN
1349  // could go parallel?
1350  //model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming,updates);
1351#endif
1352  // put row of tableau in rowArray and columnArray (packed)
1353  // get subset which have nonzero tableau elements
1354  int returnCode = transposeTimes2(updates, spareColumn1,
1355    alternateWeights_, spareColumn2, spareRow2,
1356    -scaleFactor);
1357  // zero updateBy
1358  CoinZeroN(updateBy, number);
1359  alternateWeights_->clear();
1360  // columns
1361  assert(scaleFactor);
1362  number = spareColumn1->getNumElements();
1363  index = spareColumn1->getIndices();
1364  updateBy = spareColumn1->denseVector();
1365  if (returnCode != 2 && infeasibilitiesState_) {
1366    //spareColumn1->clear();
1367    redoInfeasibilities();
1368  }
1369  if (returnCode == 1) {
1370    // most work already done
1371    for (j = 0; j < number; j++) {
1372      int iSequence = index[j];
1373      double value = updateBy[j];
1374      if (value) {
1375        updateBy[j] = 0.0;
1376        infeasible_->quickAdd(iSequence, value);
1377      } else {
1378        infeasible_->zero(iSequence);
1379      }
1380    }
1381  } else if (returnCode == 0) {
1382    reducedCost = model_->djRegion(1);
1383
1384    for (j = 0; j < number; j++) {
1385      int iSequence = index[j];
1386      double value = reducedCost[iSequence];
1387      double value2 = updateBy[j];
1388      updateBy[j] = 0.0;
1389      value -= value2;
1390      reducedCost[iSequence] = value;
1391      ClpSimplex::Status status = model_->getStatus(iSequence);
1392
1393      switch (status) {
1394
1395      case ClpSimplex::basic:
1396      case ClpSimplex::isFixed:
1397        break;
1398      case ClpSimplex::isFree:
1399      case ClpSimplex::superBasic:
1400        if (fabs(value) > FREE_ACCEPT * tolerance) {
1401          // we are going to bias towards free (but only if reasonable)
1402          value *= FREE_BIAS;
1403          // store square in list
1404          if (infeas[iSequence])
1405            infeas[iSequence] = value * value; // already there
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
1441#ifdef GROW_REFERENCE
1442    if (!reference(sequenceOut)) {
1443      outgoingWeight += 1.0;
1444      setReference(sequenceOut, true);
1445    }
1446#endif
1447    weights_[sequenceOut] = outgoingWeight;
1448    //if (model_->getStatus(sequenceOut) != ClpSimplex::basic &&
1449    //  model_->getStatus(sequenceOut) != ClpSimplex::isFixed)
1450    //checkAccuracy(sequenceOut, 1.0e-1, updates, spareRow2);
1451  }
1452  // make sure infeasibility on incoming is 0.0
1453  infeasible_->zero(sequenceIn);
1454  spareColumn2->setNumElements(0);
1455  //#define SOME_DEBUG_1
1456#ifdef SOME_DEBUG_1
1457  // check for accuracy
1458  int iCheck = 892;
1459  //printf("weight for iCheck is %g\n",weights_[iCheck]);
1460  int numberRows = model_->numberRows();
1461  //int numberColumns = model_->numberColumns();
1462  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1463    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
1464      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1465  }
1466#endif
1467  updates->setNumElements(0);
1468  spareColumn1->setNumElements(0);
1469}
1470// Update djs, weights for Devex
1471void ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector *updates,
1472  CoinIndexedVector *spareRow2,
1473  CoinIndexedVector *spareColumn1,
1474  CoinIndexedVector *spareColumn2)
1475{
1476  int iSection, j;
1477  int number = 0;
1478  int *index;
1479  double *updateBy;
1480  double *reducedCost;
1481  // dj could be very small (or even zero - take care)
1482  double dj = model_->dualIn();
1483  double tolerance = model_->currentDualTolerance();
1484  // we can't really trust infeasibilities if there is dual error
1485  // this coding has to mimic coding in checkDualSolution
1486  double error = CoinMin(1.0e-2, model_->largestDualError());
1487  // allow tolerance at least slightly bigger than standard
1488  tolerance = tolerance + error;
1489  int pivotRow = model_->pivotRow();
1490  double *infeas = infeasible_->denseVector();
1491  //updates->scanAndPack();
1492  model_->factorization()->updateColumnTranspose(spareRow2, updates);
1493
1494  // put row of tableau in rowArray and columnArray
1495  model_->clpMatrix()->transposeTimes(model_, -1.0,
1496    updates, spareColumn2, spareColumn1);
1497  // normal
1498  for (iSection = 0; iSection < 2; iSection++) {
1499
1500    reducedCost = model_->djRegion(iSection);
1501    int addSequence;
1502#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1503    double slack_multiplier;
1504#endif
1505
1506    if (!iSection) {
1507      number = updates->getNumElements();
1508      index = updates->getIndices();
1509      updateBy = updates->denseVector();
1510      addSequence = model_->numberColumns();
1511#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1512      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
1513#endif
1514    } else {
1515      number = spareColumn1->getNumElements();
1516      index = spareColumn1->getIndices();
1517      updateBy = spareColumn1->denseVector();
1518      addSequence = 0;
1519#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1520      slack_multiplier = 1.0;
1521#endif
1522    }
1523
1524    for (j = 0; j < number; j++) {
1525      int iSequence = index[j];
1526      double value = reducedCost[iSequence];
1527      value -= updateBy[j];
1528      updateBy[j] = 0.0;
1529      reducedCost[iSequence] = value;
1530      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1531
1532      switch (status) {
1533
1534      case ClpSimplex::basic:
1535        infeasible_->zero(iSequence + addSequence);
1536      case ClpSimplex::isFixed:
1537        break;
1538      case ClpSimplex::isFree:
1539      case ClpSimplex::superBasic:
1540        if (fabs(value) > FREE_ACCEPT * tolerance) {
1541          // we are going to bias towards free (but only if reasonable)
1542          value *= FREE_BIAS;
1543          // store square in list
1544          if (infeas[iSequence + addSequence])
1545            infeas[iSequence + addSequence] = value * value; // already there
1546          else
1547            infeasible_->quickAdd(iSequence + addSequence, value * value);
1548        } else {
1549          infeasible_->zero(iSequence + addSequence);
1550        }
1551        break;
1552      case ClpSimplex::atUpperBound:
1553        iSequence += addSequence;
1554        if (value > tolerance) {
1555#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1556          value *= value * slack_multiplier;
1557#else
1558          value *= value;
1559#endif
1560          // store square in list
1561          if (infeas[iSequence])
1562            infeas[iSequence] = value; // already there
1563          else
1564            infeasible_->quickAdd(iSequence, value);
1565        } else {
1566          infeasible_->zero(iSequence);
1567        }
1568        break;
1569      case ClpSimplex::atLowerBound:
1570        iSequence += addSequence;
1571        if (value < -tolerance) {
1572#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1573          value *= value * slack_multiplier;
1574#else
1575          value *= value;
1576#endif
1577          // store square in list
1578          if (infeas[iSequence])
1579            infeas[iSequence] = value; // already there
1580          else
1581            infeasible_->quickAdd(iSequence, value);
1582        } else {
1583          infeasible_->zero(iSequence);
1584        }
1585      }
1586    }
1587  }
1588  // They are empty
1589  updates->setNumElements(0);
1590  spareColumn1->setNumElements(0);
1591  // make sure infeasibility on incoming is 0.0
1592  int sequenceIn = model_->sequenceIn();
1593  infeasible_->zero(sequenceIn);
1594  // for weights update we use pivotSequence
1595  if (pivotSequence_ >= 0) {
1596    pivotRow = pivotSequence_;
1597    // unset in case sub flip
1598    pivotSequence_ = -1;
1599    // make sure infeasibility on incoming is 0.0
1600    const int *pivotVariable = model_->pivotVariable();
1601    sequenceIn = pivotVariable[pivotRow];
1602    infeasible_->zero(sequenceIn);
1603    // and we can see if reference
1604    //double referenceIn = 0.0;
1605    //if (mode_ != 1 && reference(sequenceIn))
1606    //   referenceIn = 1.0;
1607    // save outgoing weight round update
1608    double outgoingWeight = 0.0;
1609    int sequenceOut = model_->sequenceOut();
1610    if (sequenceOut >= 0)
1611      outgoingWeight = weights_[sequenceOut];
1612    // update weights
1613    updates->setNumElements(0);
1614    spareColumn1->setNumElements(0);
1615    // might as well set dj to 1
1616    dj = 1.0;
1617    updates->insert(pivotRow, -dj);
1618    model_->factorization()->updateColumnTranspose(spareRow2, updates);
1619    // put row of tableau in rowArray and columnArray
1620    model_->clpMatrix()->transposeTimes(model_, -1.0,
1621      updates, spareColumn2, spareColumn1);
1622    double *weight;
1623    int numberColumns = model_->numberColumns();
1624    // rows
1625    number = updates->getNumElements();
1626    index = updates->getIndices();
1627    updateBy = updates->denseVector();
1628    weight = weights_ + numberColumns;
1629
1630    assert(devex_ > 0.0);
1631    for (j = 0; j < number; j++) {
1632      int iSequence = index[j];
1633      double thisWeight = weight[iSequence];
1634      // row has -1
1635      double pivot = -updateBy[iSequence];
1636      updateBy[iSequence] = 0.0;
1637      double value = pivot * pivot * devex_;
1638      if (reference(iSequence + numberColumns))
1639        value += 1.0;
1640      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1641    }
1642
1643    // columns
1644    weight = weights_;
1645
1646    number = spareColumn1->getNumElements();
1647    index = spareColumn1->getIndices();
1648    updateBy = spareColumn1->denseVector();
1649    for (j = 0; j < number; j++) {
1650      int iSequence = index[j];
1651      double thisWeight = weight[iSequence];
1652      // row has -1
1653      double pivot = updateBy[iSequence];
1654      updateBy[iSequence] = 0.0;
1655      double value = pivot * pivot * devex_;
1656      if (reference(iSequence))
1657        value += 1.0;
1658      weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1659    }
1660    // restore outgoing weight
1661    if (sequenceOut >= 0)
1662      weights_[sequenceOut] = outgoingWeight;
1663    spareColumn2->setNumElements(0);
1664    //#define SOME_DEBUG_1
1665#ifdef SOME_DEBUG_1
1666    // check for accuracy
1667    int iCheck = 892;
1668    //printf("weight for iCheck is %g\n",weights_[iCheck]);
1669    int numberRows = model_->numberRows();
1670    //int numberColumns = model_->numberColumns();
1671    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1672      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
1673        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1674    }
1675#endif
1676    updates->setNumElements(0);
1677    spareColumn1->setNumElements(0);
1678  }
1679}
1680// Update djs, weights for Steepest
1681void ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector *updates,
1682  CoinIndexedVector *spareRow2,
1683  CoinIndexedVector *spareColumn1,
1684  CoinIndexedVector *spareColumn2)
1685{
1686  int iSection, j;
1687  int number = 0;
1688  int *index;
1689  double *updateBy;
1690  double *reducedCost;
1691  // dj could be very small (or even zero - take care)
1692  double dj = model_->dualIn();
1693  double tolerance = model_->currentDualTolerance();
1694  // we can't really trust infeasibilities if there is dual error
1695  // this coding has to mimic coding in checkDualSolution
1696  double error = CoinMin(1.0e-2, model_->largestDualError());
1697  // allow tolerance at least slightly bigger than standard
1698  tolerance = tolerance + error;
1699  int pivotRow = model_->pivotRow();
1700  double *infeas = infeasible_->denseVector();
1701  //updates->scanAndPack();
1702  model_->factorization()->updateColumnTranspose(spareRow2, updates);
1703
1704  // put row of tableau in rowArray and columnArray
1705  model_->clpMatrix()->transposeTimes(model_, -1.0,
1706    updates, spareColumn2, spareColumn1);
1707#ifdef CLP_USER_DRIVEN
1708  model_->eventHandler()->eventWithInfo(ClpEventHandler::beforeChooseIncoming, updates);
1709#endif
1710  // normal
1711  for (iSection = 0; iSection < 2; iSection++) {
1712
1713    reducedCost = model_->djRegion(iSection);
1714    int addSequence;
1715#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1716    double slack_multiplier;
1717#endif
1718
1719    if (!iSection) {
1720      number = updates->getNumElements();
1721      index = updates->getIndices();
1722      updateBy = updates->denseVector();
1723      addSequence = model_->numberColumns();
1724#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1725      slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
1726#endif
1727    } else {
1728      number = spareColumn1->getNumElements();
1729      index = spareColumn1->getIndices();
1730      updateBy = spareColumn1->denseVector();
1731      addSequence = 0;
1732#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1733      slack_multiplier = 1.0;
1734#endif
1735    }
1736
1737    for (j = 0; j < number; j++) {
1738      int iSequence = index[j];
1739      double value = reducedCost[iSequence];
1740      value -= updateBy[j];
1741      updateBy[j] = 0.0;
1742      reducedCost[iSequence] = value;
1743      ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
1744
1745      switch (status) {
1746
1747      case ClpSimplex::basic:
1748        infeasible_->zero(iSequence + addSequence);
1749      case ClpSimplex::isFixed:
1750        break;
1751      case ClpSimplex::isFree:
1752      case ClpSimplex::superBasic:
1753        if (fabs(value) > FREE_ACCEPT * tolerance) {
1754          // we are going to bias towards free (but only if reasonable)
1755          value *= FREE_BIAS;
1756          // store square in list
1757          if (infeas[iSequence + addSequence])
1758            infeas[iSequence + addSequence] = value * value; // already there
1759          else
1760            infeasible_->quickAdd(iSequence + addSequence, value * value);
1761        } else {
1762          infeasible_->zero(iSequence + addSequence);
1763        }
1764        break;
1765      case ClpSimplex::atUpperBound:
1766        iSequence += addSequence;
1767        if (value > tolerance) {
1768#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1769          value *= value * slack_multiplier;
1770#else
1771          value *= value;
1772#endif
1773          // store square in list
1774          if (infeas[iSequence])
1775            infeas[iSequence] = value; // already there
1776          else
1777            infeasible_->quickAdd(iSequence, value);
1778        } else {
1779          infeasible_->zero(iSequence);
1780        }
1781        break;
1782      case ClpSimplex::atLowerBound:
1783        iSequence += addSequence;
1784        if (value < -tolerance) {
1785#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1786          value *= value * slack_multiplier;
1787#else
1788          value *= value;
1789#endif
1790          // store square in list
1791          if (infeas[iSequence])
1792            infeas[iSequence] = value; // already there
1793          else
1794            infeasible_->quickAdd(iSequence, value);
1795        } else {
1796          infeasible_->zero(iSequence);
1797        }
1798      }
1799    }
1800  }
1801  // we can zero out as will have to get pivot row
1802  // ***** move
1803  updates->setNumElements(0);
1804  spareColumn1->setNumElements(0);
1805  if (pivotRow >= 0) {
1806    // make sure infeasibility on incoming is 0.0
1807    int sequenceIn = model_->sequenceIn();
1808    infeasible_->zero(sequenceIn);
1809  }
1810  // for weights update we use pivotSequence
1811  pivotRow = pivotSequence_;
1812  // unset in case sub flip
1813  pivotSequence_ = -1;
1814  if (pivotRow >= 0) {
1815    // make sure infeasibility on incoming is 0.0
1816    const int *pivotVariable = model_->pivotVariable();
1817    int sequenceIn = pivotVariable[pivotRow];
1818    assert(sequenceIn == model_->sequenceIn());
1819    infeasible_->zero(sequenceIn);
1820    // and we can see if reference
1821    double referenceIn;
1822    if (mode_ != 1) {
1823      if (reference(sequenceIn))
1824        referenceIn = 1.0;
1825      else
1826        referenceIn = 0.0;
1827    } else {
1828      referenceIn = -1.0;
1829    }
1830    // save outgoing weight round update
1831    double outgoingWeight = 0.0;
1832    int sequenceOut = model_->sequenceOut();
1833    if (sequenceOut >= 0)
1834      outgoingWeight = weights_[sequenceOut];
1835    // update weights
1836    updates->setNumElements(0);
1837    spareColumn1->setNumElements(0);
1838    // might as well set dj to 1
1839    dj = -1.0;
1840    updates->createPacked(1, &pivotRow, &dj);
1841    model_->factorization()->updateColumnTranspose(spareRow2, updates);
1842    bool needSubset = (mode_ < 4 || numberSwitched_ > 1 || mode_ >= 10);
1843
1844    double *weight;
1845    double *other = alternateWeights_->denseVector();
1846    int numberColumns = model_->numberColumns();
1847    // rows
1848    number = updates->getNumElements();
1849    index = updates->getIndices();
1850    updateBy = updates->denseVector();
1851    weight = weights_ + numberColumns;
1852    if (needSubset) {
1853#if ALT_UPDATE_WEIGHTS != 2
1854      model_->factorization()->updateColumnTranspose(spareRow2,
1855        alternateWeights_);
1856#elif ALT_UPDATE_WEIGHTS == 1
1857      if (altVector[1]) {
1858        int numberRows = model_->numberRows();
1859        double *work1 = altVector[1]->denseVector();
1860        double *worka = alternateWeights_->denseVector();
1861        int iRow = -1;
1862        double diff = 1.0e-8;
1863        for (int i = 0; i < numberRows; i++) {
1864          double dd = CoinMax(fabs(work1[i]), fabs(worka[i]));
1865          double d = fabs(work1[i] - worka[i]);
1866          if (dd > 1.0e-6 && d > diff * dd) {
1867            diff = d / dd;
1868            iRow = i;
1869          }
1870        }
1871        if (iRow >= 0)
1872          printf("largest2 difference of %g (%g,%g) on row %d\n",
1873            diff, work1[iRow], worka[iRow], iRow);
1874      }
1875#endif
1876      // do alternateWeights_ here so can scale
1877      for (j = 0; j < number; j++) {
1878        int iSequence = index[j];
1879        assert(iSequence >= 0 && iSequence < model_->numberRows());
1880        double thisWeight = weight[iSequence];
1881        // row has -1
1882        double pivot = -updateBy[j];
1883        double modification = other[iSequence];
1884        double pivotSquared = pivot * pivot;
1885
1886        thisWeight += pivotSquared * devex_ + pivot * modification;
1887        if (thisWeight < TRY_NORM) {
1888          if (mode_ == 1) {
1889            // steepest
1890            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
1891          } else {
1892            // exact
1893            thisWeight = referenceIn * pivotSquared;
1894            if (reference(iSequence + numberColumns))
1895              thisWeight += 1.0;
1896            thisWeight = CoinMax(thisWeight, TRY_NORM);
1897          }
1898        }
1899        weight[iSequence] = thisWeight;
1900      }
1901      transposeTimes2(updates, spareColumn1, alternateWeights_, spareColumn2, spareRow2, 0.0);
1902    } else {
1903      // put row of tableau in rowArray and columnArray
1904      model_->clpMatrix()->transposeTimes(model_, -1.0,
1905        updates, spareColumn2, spareColumn1);
1906    }
1907
1908    if (needSubset) {
1909      CoinZeroN(updateBy, number);
1910    } else if (mode_ == 4) {
1911      // Devex
1912      assert(devex_ > 0.0);
1913      for (j = 0; j < number; j++) {
1914        int iSequence = index[j];
1915        double thisWeight = weight[iSequence];
1916        // row has -1
1917        double pivot = -updateBy[j];
1918        updateBy[j] = 0.0;
1919        double value = pivot * pivot * devex_;
1920        if (reference(iSequence + numberColumns))
1921          value += 1.0;
1922        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1923      }
1924    }
1925
1926    // columns
1927    weight = weights_;
1928
1929    number = spareColumn1->getNumElements();
1930    index = spareColumn1->getIndices();
1931    updateBy = spareColumn1->denseVector();
1932    if (needSubset) {
1933      // Exact - already done
1934    } else if (mode_ == 4) {
1935      // Devex
1936      for (j = 0; j < number; j++) {
1937        int iSequence = index[j];
1938        double thisWeight = weight[iSequence];
1939        // row has -1
1940        double pivot = updateBy[j];
1941        updateBy[j] = 0.0;
1942        double value = pivot * pivot * devex_;
1943        if (reference(iSequence))
1944          value += 1.0;
1945        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1946      }
1947    }
1948    // restore outgoing weight
1949    if (sequenceOut >= 0)
1950      weights_[sequenceOut] = outgoingWeight;
1951    alternateWeights_->clear();
1952    spareColumn2->setNumElements(0);
1953    //#define SOME_DEBUG_1
1954#ifdef SOME_DEBUG_1
1955    // check for accuracy
1956    int iCheck = 892;
1957    //printf("weight for iCheck is %g\n",weights_[iCheck]);
1958    int numberRows = model_->numberRows();
1959    //int numberColumns = model_->numberColumns();
1960    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1961      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
1962        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
1963    }
1964#endif
1965  }
1966  updates->setNumElements(0);
1967  spareColumn1->setNumElements(0);
1968}
1969// Updates two arrays for steepest
1970int ClpPrimalColumnSteepest::transposeTimes2(const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
1971  const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
1972  CoinIndexedVector *spare,
1973  double scaleFactor)
1974{
1975  // see if reference
1976  int sequenceIn = model_->sequenceIn();
1977  double referenceIn;
1978  if (mode_ != 1) {
1979    if (reference(sequenceIn))
1980      referenceIn = 1.0;
1981    else
1982      referenceIn = 0.0;
1983  } else {
1984    referenceIn = -1.0;
1985  }
1986  int returnCode = 0;
1987  if (model_->clpMatrix()->canCombine(model_, pi1)) {
1988    double *infeas = scaleFactor ? infeasible_->denseVector() : NULL;
1989    // put row of tableau in rowArray and columnArray
1990    returnCode = model_->clpMatrix()->transposeTimes2(model_, pi1,
1991      dj1, pi2, spare,
1992      infeas,
1993      model_->djRegion(1),
1994      referenceIn, devex_,
1995      reference_,
1996      weights_, scaleFactor);
1997    if (model_->spareIntArray_[3] > -2)
1998      returnCode = 2;
1999  } else {
2000    // put row of tableau in rowArray and columnArray
2001    model_->clpMatrix()->transposeTimes(model_, -1.0,
2002      pi1, dj2, dj1);
2003    // get subset which have nonzero tableau elements
2004    model_->clpMatrix()->subsetTransposeTimes(model_, pi2, dj1, dj2);
2005    bool killDjs = (scaleFactor == 0.0);
2006    if (!scaleFactor)
2007      scaleFactor = 1.0;
2008    // columns
2009    double *weight = weights_;
2010
2011    int number = dj1->getNumElements();
2012    const int *index = dj1->getIndices();
2013    double *updateBy = dj1->denseVector();
2014    double *updateBy2 = dj2->denseVector();
2015
2016    for (int j = 0; j < number; j++) {
2017      double thisWeight;
2018      double pivot;
2019      double pivotSquared;
2020      int iSequence = index[j];
2021      double value2 = updateBy[j];
2022      if (killDjs)
2023        updateBy[j] = 0.0;
2024      double modification = updateBy2[j];
2025      updateBy2[j] = 0.0;
2026      ClpSimplex::Status status = model_->getStatus(iSequence);
2027
2028      if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) {
2029        thisWeight = weight[iSequence];
2030        pivot = value2 * scaleFactor;
2031        pivotSquared = pivot * pivot;
2032
2033        thisWeight += pivotSquared * devex_ + pivot * modification;
2034        if (thisWeight < TRY_NORM) {
2035          if (referenceIn < 0.0) {
2036            // steepest
2037            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2038          } else {
2039            // exact
2040            thisWeight = referenceIn * pivotSquared;
2041            if (reference(iSequence))
2042              thisWeight += 1.0;
2043            thisWeight = CoinMax(thisWeight, TRY_NORM);
2044          }
2045        }
2046        weight[iSequence] = thisWeight;
2047      }
2048    }
2049  }
2050  dj2->setNumElements(0);
2051  return returnCode;
2052}
2053// Update weights for Devex
2054void ClpPrimalColumnSteepest::justDevex(CoinIndexedVector *updates,
2055  CoinIndexedVector *spareRow2,
2056  CoinIndexedVector *spareColumn1,
2057  CoinIndexedVector *spareColumn2)
2058{
2059  int j;
2060  int number = 0;
2061  int *index;
2062  double *updateBy;
2063  // dj could be very small (or even zero - take care)
2064  double dj = model_->dualIn();
2065  double tolerance = model_->currentDualTolerance();
2066  // we can't really trust infeasibilities if there is dual error
2067  // this coding has to mimic coding in checkDualSolution
2068  double error = CoinMin(1.0e-2, model_->largestDualError());
2069  // allow tolerance at least slightly bigger than standard
2070  tolerance = tolerance + error;
2071  int pivotRow = model_->pivotRow();
2072  // for weights update we use pivotSequence
2073  pivotRow = pivotSequence_;
2074  assert(pivotRow >= 0);
2075  // make sure infeasibility on incoming is 0.0
2076  const int *pivotVariable = model_->pivotVariable();
2077  int sequenceIn = pivotVariable[pivotRow];
2078  infeasible_->zero(sequenceIn);
2079  // and we can see if reference
2080  //double referenceIn = 0.0;
2081  //if (mode_ != 1 && reference(sequenceIn))
2082  //   referenceIn = 1.0;
2083  // save outgoing weight round update
2084  double outgoingWeight = 0.0;
2085  int sequenceOut = model_->sequenceOut();
2086  if (sequenceOut >= 0)
2087    outgoingWeight = weights_[sequenceOut];
2088  assert(!updates->getNumElements());
2089  assert(!spareColumn1->getNumElements());
2090  // unset in case sub flip
2091  pivotSequence_ = -1;
2092  // might as well set dj to 1
2093  dj = -1.0;
2094  updates->createPacked(1, &pivotRow, &dj);
2095  model_->factorization()->updateColumnTranspose(spareRow2, updates);
2096  // put row of tableau in rowArray and columnArray
2097  model_->clpMatrix()->transposeTimes(model_, -1.0,
2098    updates, spareColumn2, spareColumn1);
2099  double *weight;
2100  int numberColumns = model_->numberColumns();
2101  // rows
2102  number = updates->getNumElements();
2103  index = updates->getIndices();
2104  updateBy = updates->denseVector();
2105  weight = weights_ + numberColumns;
2106
2107  // Devex
2108  assert(devex_ > 0.0);
2109  for (j = 0; j < number; j++) {
2110    int iSequence = index[j];
2111    double thisWeight = weight[iSequence];
2112    // row has -1
2113    double pivot = -updateBy[j];
2114    updateBy[j] = 0.0;
2115    double value = pivot * pivot * devex_;
2116    if (reference(iSequence + numberColumns))
2117      value += 1.0;
2118    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2119  }
2120
2121  // columns
2122  weight = weights_;
2123
2124  number = spareColumn1->getNumElements();
2125  index = spareColumn1->getIndices();
2126  updateBy = spareColumn1->denseVector();
2127  // Devex
2128  for (j = 0; j < number; j++) {
2129    int iSequence = index[j];
2130    double thisWeight = weight[iSequence];
2131    // row has -1
2132    double pivot = updateBy[j];
2133    updateBy[j] = 0.0;
2134    double value = pivot * pivot * devex_;
2135    if (reference(iSequence))
2136      value += 1.0;
2137    weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2138  }
2139  // restore outgoing weight
2140  if (sequenceOut >= 0)
2141    weights_[sequenceOut] = outgoingWeight;
2142  spareColumn2->setNumElements(0);
2143  //#define SOME_DEBUG_1
2144#ifdef SOME_DEBUG_1
2145  // check for accuracy
2146  int iCheck = 892;
2147  //printf("weight for iCheck is %g\n",weights_[iCheck]);
2148  int numberRows = model_->numberRows();
2149  //int numberColumns = model_->numberColumns();
2150  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
2151    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
2152      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
2153  }
2154#endif
2155  updates->setNumElements(0);
2156  spareColumn1->setNumElements(0);
2157}
2158// Update weights for Steepest
2159void ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector *updates,
2160  CoinIndexedVector *spareRow2,
2161  CoinIndexedVector *spareColumn1,
2162  CoinIndexedVector *spareColumn2)
2163{
2164  int j;
2165  int number = 0;
2166  int *index;
2167  double *updateBy;
2168  // dj could be very small (or even zero - take care)
2169  double dj = model_->dualIn();
2170  double tolerance = model_->currentDualTolerance();
2171  // we can't really trust infeasibilities if there is dual error
2172  // this coding has to mimic coding in checkDualSolution
2173  double error = CoinMin(1.0e-2, model_->largestDualError());
2174  // allow tolerance at least slightly bigger than standard
2175  tolerance = tolerance + error;
2176  int pivotRow = model_->pivotRow();
2177  // for weights update we use pivotSequence
2178  pivotRow = pivotSequence_;
2179  // unset in case sub flip
2180  pivotSequence_ = -1;
2181  assert(pivotRow >= 0);
2182  // make sure infeasibility on incoming is 0.0
2183  const int *pivotVariable = model_->pivotVariable();
2184  int sequenceIn = pivotVariable[pivotRow];
2185  infeasible_->zero(sequenceIn);
2186  // and we can see if reference
2187  double referenceIn = 0.0;
2188  if (mode_ != 1 && reference(sequenceIn))
2189    referenceIn = 1.0;
2190  // save outgoing weight round update
2191  double outgoingWeight = 0.0;
2192  int sequenceOut = model_->sequenceOut();
2193  if (sequenceOut >= 0)
2194    outgoingWeight = weights_[sequenceOut];
2195  assert(!updates->getNumElements());
2196  assert(!spareColumn1->getNumElements());
2197  // update weights
2198  //updates->setNumElements(0);
2199  //spareColumn1->setNumElements(0);
2200  // might as well set dj to 1
2201  dj = -1.0;
2202  updates->createPacked(1, &pivotRow, &dj);
2203  model_->factorization()->updateColumnTranspose(spareRow2, updates);
2204  // put row of tableau in rowArray and columnArray
2205  model_->clpMatrix()->transposeTimes(model_, -1.0,
2206    updates, spareColumn2, spareColumn1);
2207  double *weight;
2208  double *other = alternateWeights_->denseVector();
2209  int numberColumns = model_->numberColumns();
2210  // rows
2211  number = updates->getNumElements();
2212  index = updates->getIndices();
2213  updateBy = updates->denseVector();
2214  weight = weights_ + numberColumns;
2215
2216  // Exact
2217  // now update weight update array
2218  //alternateWeights_->scanAndPack();
2219#if ALT_UPDATE_WEIGHTS != 2
2220  model_->factorization()->updateColumnTranspose(spareRow2,
2221    alternateWeights_);
2222#elif ALT_UPDATE_WEIGHTS == 1
2223  if (altVector[1]) {
2224    int numberRows = model_->numberRows();
2225    double *work1 = altVector[1]->denseVector();
2226    double *worka = alternateWeights_->denseVector();
2227    int iRow = -1;
2228    double diff = 1.0e-8;
2229    for (int i = 0; i < numberRows; i++) {
2230      double dd = CoinMax(fabs(work1[i]), fabs(worka[i]));
2231      double d = fabs(work1[i] - worka[i]);
2232      if (dd > 1.0e-6 && d > diff * dd) {
2233        diff = d / dd;
2234        iRow = i;
2235      }
2236    }
2237    if (iRow >= 0)
2238      printf("largest3 difference of %g (%g,%g) on row %d\n",
2239        diff, work1[iRow], worka[iRow], iRow);
2240  }
2241#endif
2242  // get subset which have nonzero tableau elements
2243  model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
2244    spareColumn1,
2245    spareColumn2);
2246  for (j = 0; j < number; j++) {
2247    int iSequence = index[j];
2248    double thisWeight = weight[iSequence];
2249    // row has -1
2250    double pivot = -updateBy[j];
2251    updateBy[j] = 0.0;
2252    double modification = other[iSequence];
2253    double pivotSquared = pivot * pivot;
2254
2255    thisWeight += pivotSquared * devex_ + pivot * modification;
2256    if (thisWeight < TRY_NORM) {
2257      if (mode_ == 1) {
2258        // steepest
2259        thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2260      } else {
2261        // exact
2262        thisWeight = referenceIn * pivotSquared;
2263        if (reference(iSequence + numberColumns))
2264          thisWeight += 1.0;
2265        thisWeight = CoinMax(thisWeight, TRY_NORM);
2266      }
2267    }
2268    weight[iSequence] = thisWeight;
2269  }
2270
2271  // columns
2272  weight = weights_;
2273
2274  number = spareColumn1->getNumElements();
2275  index = spareColumn1->getIndices();
2276  updateBy = spareColumn1->denseVector();
2277  // Exact
2278  double *updateBy2 = spareColumn2->denseVector();
2279  for (j = 0; j < number; j++) {
2280    int iSequence = index[j];
2281    double thisWeight = weight[iSequence];
2282    double pivot = updateBy[j];
2283    updateBy[j] = 0.0;
2284    double modification = updateBy2[j];
2285    updateBy2[j] = 0.0;
2286    double pivotSquared = pivot * pivot;
2287
2288    thisWeight += pivotSquared * devex_ + pivot * modification;
2289    if (thisWeight < TRY_NORM) {
2290      if (mode_ == 1) {
2291        // steepest
2292        thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2293      } else {
2294        // exact
2295        thisWeight = referenceIn * pivotSquared;
2296        if (reference(iSequence))
2297          thisWeight += 1.0;
2298        thisWeight = CoinMax(thisWeight, TRY_NORM);
2299      }
2300    }
2301    weight[iSequence] = thisWeight;
2302  }
2303  // restore outgoing weight
2304  if (sequenceOut >= 0)
2305    weights_[sequenceOut] = outgoingWeight;
2306  alternateWeights_->clear();
2307  spareColumn2->setNumElements(0);
2308  //#define SOME_DEBUG_1
2309#ifdef SOME_DEBUG_1
2310  // check for accuracy
2311  int iCheck = 892;
2312  //printf("weight for iCheck is %g\n",weights_[iCheck]);
2313  int numberRows = model_->numberRows();
2314  //int numberColumns = model_->numberColumns();
2315  for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
2316    if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
2317      checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
2318  }
2319#endif
2320  updates->setNumElements(0);
2321  spareColumn1->setNumElements(0);
2322}
2323// Returns pivot column, -1 if none
2324int ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector *updates,
2325  CoinIndexedVector *,
2326  CoinIndexedVector *spareRow2,
2327  CoinIndexedVector *spareColumn1,
2328  CoinIndexedVector *spareColumn2)
2329{
2330  assert(model_);
2331  int iSection, j;
2332  int number = 0;
2333  int *index;
2334  double *updateBy;
2335  double *reducedCost;
2336  // dj could be very small (or even zero - take care)
2337  double dj = model_->dualIn();
2338  double tolerance = model_->currentDualTolerance();
2339  // we can't really trust infeasibilities if there is dual error
2340  // this coding has to mimic coding in checkDualSolution
2341  double error = CoinMin(1.0e-2, model_->largestDualError());
2342  // allow tolerance at least slightly bigger than standard
2343  tolerance = tolerance + error;
2344  int pivotRow = model_->pivotRow();
2345  int anyUpdates;
2346  double *infeas = infeasible_->denseVector();
2347
2348  // Local copy of mode so can decide what to do
2349  int switchType;
2350  if (mode_ == 4)
2351    switchType = 5 - numberSwitched_;
2352  else if (mode_ >= 10)
2353    switchType = 3;
2354  else
2355    switchType = mode_;
2356  /* switchType -
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     */
2364  if (updates->getNumElements()) {
2365    // would have to have two goes for devex, three for steepest
2366    anyUpdates = 2;
2367    // add in pivot contribution
2368    if (pivotRow >= 0)
2369      updates->add(pivotRow, -dj);
2370  } else if (pivotRow >= 0) {
2371    if (fabs(dj) > 1.0e-15) {
2372      // some dj
2373      updates->insert(pivotRow, -dj);
2374      if (fabs(dj) > 1.0e-6) {
2375        // reasonable size
2376        anyUpdates = 1;
2377      } else {
2378        // too small
2379        anyUpdates = 2;
2380      }
2381    } else {
2382      // zero dj
2383      anyUpdates = -1;
2384    }
2385  } else if (pivotSequence_ >= 0) {
2386    // just after re-factorization
2387    anyUpdates = -1;
2388  } else {
2389    // sub flip - nothing to do
2390    anyUpdates = 0;
2391  }
2392
2393  if (anyUpdates > 0) {
2394    model_->factorization()->updateColumnTranspose(spareRow2, updates);
2395
2396    // put row of tableau in rowArray and columnArray
2397    model_->clpMatrix()->transposeTimes(model_, -1.0,
2398      updates, spareColumn2, spareColumn1);
2399    // normal
2400    for (iSection = 0; iSection < 2; iSection++) {
2401
2402      reducedCost = model_->djRegion(iSection);
2403      int addSequence;
2404
2405      if (!iSection) {
2406        number = updates->getNumElements();
2407        index = updates->getIndices();
2408        updateBy = updates->denseVector();
2409        addSequence = model_->numberColumns();
2410      } else {
2411        number = spareColumn1->getNumElements();
2412        index = spareColumn1->getIndices();
2413        updateBy = spareColumn1->denseVector();
2414        addSequence = 0;
2415      }
2416      if (!model_->nonLinearCost()->lookBothWays()) {
2417
2418        for (j = 0; j < number; j++) {
2419          int iSequence = index[j];
2420          double value = reducedCost[iSequence];
2421          value -= updateBy[iSequence];
2422          reducedCost[iSequence] = value;
2423          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2424
2425          switch (status) {
2426
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            }
2466          }
2467        }
2468      } else {
2469        ClpNonLinearCost *nonLinear = model_->nonLinearCost();
2470        // We can go up OR down
2471        for (j = 0; j < number; j++) {
2472          int iSequence = index[j];
2473          double value = reducedCost[iSequence];
2474          value -= updateBy[iSequence];
2475          reducedCost[iSequence] = value;
2476          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2477
2478          switch (status) {
2479
2480          case ClpSimplex::basic:
2481            infeasible_->zero(iSequence + addSequence);
2482          case ClpSimplex::isFixed:
2483            break;
2484          case ClpSimplex::isFree:
2485          case ClpSimplex::superBasic:
2486            if (fabs(value) > FREE_ACCEPT * tolerance) {
2487              // we are going to bias towards free (but only if reasonable)
2488              value *= FREE_BIAS;
2489              // store square in list
2490              if (infeas[iSequence + addSequence])
2491                infeas[iSequence + addSequence] = value * value; // already there
2492              else
2493                infeasible_->quickAdd(iSequence + addSequence, value * value);
2494            } else {
2495              infeasible_->zero(iSequence + addSequence);
2496            }
2497            break;
2498          case ClpSimplex::atUpperBound:
2499            if (value > tolerance) {
2500              // store square in list
2501              if (infeas[iSequence + addSequence])
2502                infeas[iSequence + addSequence] = value * value; // already there
2503              else
2504                infeasible_->quickAdd(iSequence + addSequence, value * value);
2505            } else {
2506              // look other way - change up should be negative
2507              value -= nonLinear->changeUpInCost(iSequence + addSequence);
2508              if (value > -tolerance) {
2509                infeasible_->zero(iSequence + addSequence);
2510              } else {
2511                // store square in list
2512                if (infeas[iSequence + addSequence])
2513                  infeas[iSequence + addSequence] = value * value; // already there
2514                else
2515                  infeasible_->quickAdd(iSequence + addSequence, value * value);
2516              }
2517            }
2518            break;
2519          case ClpSimplex::atLowerBound:
2520            if (value < -tolerance) {
2521              // store square in list
2522              if (infeas[iSequence + addSequence])
2523                infeas[iSequence + addSequence] = value * value; // already there
2524              else
2525                infeasible_->quickAdd(iSequence + addSequence, value * value);
2526            } else {
2527              // look other way - change down should be positive
2528              value -= nonLinear->changeDownInCost(iSequence + addSequence);
2529              if (value < tolerance) {
2530                infeasible_->zero(iSequence + addSequence);
2531              } else {
2532                // store square in list
2533                if (infeas[iSequence + addSequence])
2534                  infeas[iSequence + addSequence] = value * value; // already there
2535                else
2536                  infeasible_->quickAdd(iSequence + addSequence, value * value);
2537              }
2538            }
2539          }
2540        }
2541      }
2542    }
2543    if (anyUpdates == 2) {
2544      // we can zero out as will have to get pivot row
2545      updates->clear();
2546      spareColumn1->clear();
2547    }
2548    if (pivotRow >= 0) {
2549      // make sure infeasibility on incoming is 0.0
2550      int sequenceIn = model_->sequenceIn();
2551      infeasible_->zero(sequenceIn);
2552    }
2553  }
2554  // make sure outgoing from last iteration okay
2555  int sequenceOut = model_->sequenceOut();
2556  if (sequenceOut >= 0) {
2557    ClpSimplex::Status status = model_->getStatus(sequenceOut);
2558    double value = model_->reducedCost(sequenceOut);
2559
2560    switch (status) {
2561
2562    case ClpSimplex::basic:
2563    case ClpSimplex::isFixed:
2564      break;
2565    case ClpSimplex::isFree:
2566    case ClpSimplex::superBasic:
2567      if (fabs(value) > FREE_ACCEPT * tolerance) {
2568        // we are going to bias towards free (but only if reasonable)
2569        value *= FREE_BIAS;
2570        // store square in list
2571        if (infeas[sequenceOut])
2572          infeas[sequenceOut] = value * value; // already there
2573        else
2574          infeasible_->quickAdd(sequenceOut, value * value);
2575      } else {
2576        infeasible_->zero(sequenceOut);
2577      }
2578      break;
2579    case ClpSimplex::atUpperBound:
2580      if (value > tolerance) {
2581        // store square in list
2582        if (infeas[sequenceOut])
2583          infeas[sequenceOut] = value * value; // already there
2584        else
2585          infeasible_->quickAdd(sequenceOut, value * value);
2586      } else {
2587        infeasible_->zero(sequenceOut);
2588      }
2589      break;
2590    case ClpSimplex::atLowerBound:
2591      if (value < -tolerance) {
2592        // store square in list
2593        if (infeas[sequenceOut])
2594          infeas[sequenceOut] = value * value; // already there
2595        else
2596          infeasible_->quickAdd(sequenceOut, value * value);
2597      } else {
2598        infeasible_->zero(sequenceOut);
2599      }
2600    }
2601  }
2602
2603  // If in quadratic re-compute all
2604  if (model_->algorithm() == 2) {
2605    for (iSection = 0; iSection < 2; iSection++) {
2606
2607      reducedCost = model_->djRegion(iSection);
2608      int addSequence;
2609      int iSequence;
2610
2611      if (!iSection) {
2612        number = model_->numberRows();
2613        addSequence = model_->numberColumns();
2614      } else {
2615        number = model_->numberColumns();
2616        addSequence = 0;
2617      }
2618
2619      if (!model_->nonLinearCost()->lookBothWays()) {
2620        for (iSequence = 0; iSequence < number; iSequence++) {
2621          double value = reducedCost[iSequence];
2622          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2623
2624          switch (status) {
2625
2626          case ClpSimplex::basic:
2627            infeasible_->zero(iSequence + addSequence);
2628          case ClpSimplex::isFixed:
2629            break;
2630          case ClpSimplex::isFree:
2631          case ClpSimplex::superBasic:
2632            if (fabs(value) > tolerance) {
2633              // we are going to bias towards free (but only if reasonable)
2634              value *= FREE_BIAS;
2635              // store square in list
2636              if (infeas[iSequence + addSequence])
2637                infeas[iSequence + addSequence] = value * value; // already there
2638              else
2639                infeasible_->quickAdd(iSequence + addSequence, value * value);
2640            } else {
2641              infeasible_->zero(iSequence + addSequence);
2642            }
2643            break;
2644          case ClpSimplex::atUpperBound:
2645            if (value > tolerance) {
2646              // store square in list
2647              if (infeas[iSequence + addSequence])
2648                infeas[iSequence + addSequence] = value * value; // already there
2649              else
2650                infeasible_->quickAdd(iSequence + addSequence, value * value);
2651            } else {
2652              infeasible_->zero(iSequence + addSequence);
2653            }
2654            break;
2655          case ClpSimplex::atLowerBound:
2656            if (value < -tolerance) {
2657              // store square in list
2658              if (infeas[iSequence + addSequence])
2659                infeas[iSequence + addSequence] = value * value; // already there
2660              else
2661                infeasible_->quickAdd(iSequence + addSequence, value * value);
2662            } else {
2663              infeasible_->zero(iSequence + addSequence);
2664            }
2665          }
2666        }
2667      } else {
2668        // we can go both ways
2669        ClpNonLinearCost *nonLinear = model_->nonLinearCost();
2670        for (iSequence = 0; iSequence < number; iSequence++) {
2671          double value = reducedCost[iSequence];
2672          ClpSimplex::Status status = model_->getStatus(iSequence + addSequence);
2673
2674          switch (status) {
2675
2676          case ClpSimplex::basic:
2677            infeasible_->zero(iSequence + addSequence);
2678          case ClpSimplex::isFixed:
2679            break;
2680          case ClpSimplex::isFree:
2681          case ClpSimplex::superBasic:
2682            if (fabs(value) > tolerance) {
2683              // we are going to bias towards free (but only if reasonable)
2684              value *= FREE_BIAS;
2685              // store square in list
2686              if (infeas[iSequence + addSequence])
2687                infeas[iSequence + addSequence] = value * value; // already there
2688              else
2689                infeasible_->quickAdd(iSequence + addSequence, value * value);
2690            } else {
2691              infeasible_->zero(iSequence + addSequence);
2692            }
2693            break;
2694          case ClpSimplex::atUpperBound:
2695            if (value > tolerance) {
2696              // store square in list
2697              if (infeas[iSequence + addSequence])
2698                infeas[iSequence + addSequence] = value * value; // already there
2699              else
2700                infeasible_->quickAdd(iSequence + addSequence, value * value);
2701            } else {
2702              // look other way - change up should be negative
2703              value -= nonLinear->changeUpInCost(iSequence + addSequence);
2704              if (value > -tolerance) {
2705                infeasible_->zero(iSequence + addSequence);
2706              } else {
2707                // store square in list
2708                if (infeas[iSequence + addSequence])
2709                  infeas[iSequence + addSequence] = value * value; // already there
2710                else
2711                  infeasible_->quickAdd(iSequence + addSequence, value * value);
2712              }
2713            }
2714            break;
2715          case ClpSimplex::atLowerBound:
2716            if (value < -tolerance) {
2717              // store square in list
2718              if (infeas[iSequence + addSequence])
2719                infeas[iSequence + addSequence] = value * value; // already there
2720              else
2721                infeasible_->quickAdd(iSequence + addSequence, value * value);
2722            } else {
2723              // look other way - change down should be positive
2724              value -= nonLinear->changeDownInCost(iSequence + addSequence);
2725              if (value < tolerance) {
2726                infeasible_->zero(iSequence + addSequence);
2727              } else {
2728                // store square in list
2729                if (infeas[iSequence + addSequence])
2730                  infeas[iSequence + addSequence] = value * value; // already there
2731                else
2732                  infeasible_->quickAdd(iSequence + addSequence, value * value);
2733              }
2734            }
2735          }
2736        }
2737      }
2738    }
2739  }
2740  // See what sort of pricing
2741  int numberWanted = 10;
2742  number = infeasible_->getNumElements();
2743  int numberColumns = model_->numberColumns();
2744  if (switchType == 5) {
2745    // we can zero out
2746    updates->clear();
2747    spareColumn1->clear();
2748    pivotSequence_ = -1;
2749    pivotRow = -1;
2750    // See if to switch
2751    int numberRows = model_->numberRows();
2752    // ratio is done on number of columns here
2753    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
2754    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
2755    //double ratio = static_cast<double> sizeFactorization_/static_cast<double> model_->clpMatrix()->getNumElements();
2756    if (ratio < 0.1) {
2757      numberWanted = CoinMax(100, number / 200);
2758    } else if (ratio < 0.3) {
2759      numberWanted = CoinMax(500, number / 40);
2760    } else if (ratio < 0.5 || mode_ == 5) {
2761      numberWanted = CoinMax(2000, number / 10);
2762      numberWanted = CoinMax(numberWanted, numberColumns / 30);
2763    } else if (mode_ != 5) {
2764      switchType = 4;
2765      // initialize
2766      numberSwitched_++;
2767      // Make sure will re-do
2768      delete[] weights_;
2769      weights_ = NULL;
2770      saveWeights(model_, 4);
2771      COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
2772      updates->clear();
2773    }
2774    if (model_->numberIterations() % 1000 == 0)
2775      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
2776  }
2777  if (switchType == 4) {
2778    // Still in devex mode
2779    int numberRows = model_->numberRows();
2780    // ratio is done on number of rows here
2781    double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(numberRows);
2782    // Go to steepest if lot of iterations?
2783    if (ratio < 1.0) {
2784      numberWanted = CoinMax(2000, number / 20);
2785    } else if (ratio < 5.0) {
2786      numberWanted = CoinMax(2000, number / 10);
2787      numberWanted = CoinMax(numberWanted, numberColumns / 20);
2788    } else {
2789      // we can zero out
2790      updates->clear();
2791      spareColumn1->clear();
2792      switchType = 3;
2793      // initialize
2794      pivotSequence_ = -1;
2795      pivotRow = -1;
2796      numberSwitched_++;
2797      // Make sure will re-do
2798      delete[] weights_;
2799      weights_ = NULL;
2800      saveWeights(model_, 4);
2801      COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
2802      updates->clear();
2803    }
2804    if (model_->numberIterations() % 1000 == 0)
2805      COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d\n", sizeFactorization_, ratio, numberWanted));
2806  }
2807  if (switchType < 4) {
2808    if (switchType < 2) {
2809      numberWanted = number + 1;
2810    } else if (switchType == 2) {
2811      numberWanted = CoinMax(2000, number / 8);
2812    } else {
2813      double ratio = static_cast< double >(sizeFactorization_) / static_cast< double >(model_->numberRows());
2814      if (ratio < 1.0) {
2815        numberWanted = CoinMax(2000, number / 20);
2816      } else if (ratio < 5.0) {
2817        numberWanted = CoinMax(2000, number / 10);
2818        numberWanted = CoinMax(numberWanted, numberColumns / 20);
2819      } else if (ratio < 10.0) {
2820        numberWanted = CoinMax(2000, number / 8);
2821        numberWanted = CoinMax(numberWanted, numberColumns / 20);
2822      } else {
2823        ratio = number * (ratio / 80.0);
2824        if (ratio > number) {
2825          numberWanted = number + 1;
2826        } else {
2827          numberWanted = CoinMax(2000, static_cast< int >(ratio));
2828          numberWanted = CoinMax(numberWanted, numberColumns / 10);
2829        }
2830      }
2831    }
2832  }
2833  // for weights update we use pivotSequence
2834  pivotRow = pivotSequence_;
2835  // unset in case sub flip
2836  pivotSequence_ = -1;
2837  if (pivotRow >= 0) {
2838    // make sure infeasibility on incoming is 0.0
2839    const int *pivotVariable = model_->pivotVariable();
2840    int sequenceIn = pivotVariable[pivotRow];
2841    infeasible_->zero(sequenceIn);
2842    // and we can see if reference
2843    double referenceIn = 0.0;
2844    if (switchType != 1 && reference(sequenceIn))
2845      referenceIn = 1.0;
2846    // save outgoing weight round update
2847    double outgoingWeight = 0.0;
2848    if (sequenceOut >= 0)
2849      outgoingWeight = weights_[sequenceOut];
2850    // update weights
2851    if (anyUpdates != 1) {
2852      updates->setNumElements(0);
2853      spareColumn1->setNumElements(0);
2854      // might as well set dj to 1
2855      dj = 1.0;
2856      updates->insert(pivotRow, -dj);
2857      model_->factorization()->updateColumnTranspose(spareRow2, updates);
2858      // put row of tableau in rowArray and columnArray
2859      model_->clpMatrix()->transposeTimes(model_, -1.0,
2860        updates, spareColumn2, spareColumn1);
2861    }
2862    double *weight;
2863    double *other = alternateWeights_->denseVector();
2864    int numberColumns = model_->numberColumns();
2865    double scaleFactor = -1.0 / dj; // as formula is with 1.0
2866    // rows
2867    number = updates->getNumElements();
2868    index = updates->getIndices();
2869    updateBy = updates->denseVector();
2870    weight = weights_ + numberColumns;
2871
2872    if (switchType < 4) {
2873      // Exact
2874      // now update weight update array
2875      model_->factorization()->updateColumnTranspose(spareRow2,
2876        alternateWeights_);
2877#ifdef ALT_UPDATE_WEIGHTS
2878      abort();
2879#endif
2880      for (j = 0; j < number; j++) {
2881        int iSequence = index[j];
2882        double thisWeight = weight[iSequence];
2883        // row has -1
2884        double pivot = updateBy[iSequence] * scaleFactor;
2885        updateBy[iSequence] = 0.0;
2886        double modification = other[iSequence];
2887        double pivotSquared = pivot * pivot;
2888
2889        thisWeight += pivotSquared * devex_ + pivot * modification;
2890        if (thisWeight < TRY_NORM) {
2891          if (switchType == 1) {
2892            // steepest
2893            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2894          } else {
2895            // exact
2896            thisWeight = referenceIn * pivotSquared;
2897            if (reference(iSequence + numberColumns))
2898              thisWeight += 1.0;
2899            thisWeight = CoinMax(thisWeight, TRY_NORM);
2900          }
2901        }
2902        weight[iSequence] = thisWeight;
2903      }
2904    } else if (switchType == 4) {
2905      // Devex
2906      assert(devex_ > 0.0);
2907      for (j = 0; j < number; j++) {
2908        int iSequence = index[j];
2909        double thisWeight = weight[iSequence];
2910        // row has -1
2911        double pivot = updateBy[iSequence] * scaleFactor;
2912        updateBy[iSequence] = 0.0;
2913        double value = pivot * pivot * devex_;
2914        if (reference(iSequence + numberColumns))
2915          value += 1.0;
2916        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2917      }
2918    }
2919
2920    // columns
2921    weight = weights_;
2922
2923    scaleFactor = -scaleFactor;
2924
2925    number = spareColumn1->getNumElements();
2926    index = spareColumn1->getIndices();
2927    updateBy = spareColumn1->denseVector();
2928    if (switchType < 4) {
2929      // Exact
2930      // get subset which have nonzero tableau elements
2931      model_->clpMatrix()->subsetTransposeTimes(model_, alternateWeights_,
2932        spareColumn1,
2933        spareColumn2);
2934      double *updateBy2 = spareColumn2->denseVector();
2935      for (j = 0; j < number; j++) {
2936        int iSequence = index[j];
2937        double thisWeight = weight[iSequence];
2938        double pivot = updateBy[iSequence] * scaleFactor;
2939        updateBy[iSequence] = 0.0;
2940        double modification = updateBy2[j];
2941        updateBy2[j] = 0.0;
2942        double pivotSquared = pivot * pivot;
2943
2944        thisWeight += pivotSquared * devex_ + pivot * modification;
2945        if (thisWeight < TRY_NORM) {
2946          if (switchType == 1) {
2947            // steepest
2948            thisWeight = CoinMax(TRY_NORM, ADD_ONE + pivotSquared);
2949          } else {
2950            // exact
2951            thisWeight = referenceIn * pivotSquared;
2952            if (reference(iSequence))
2953              thisWeight += 1.0;
2954            thisWeight = CoinMax(thisWeight, TRY_NORM);
2955          }
2956        }
2957        weight[iSequence] = thisWeight;
2958      }
2959    } else if (switchType == 4) {
2960      // Devex
2961      for (j = 0; j < number; j++) {
2962        int iSequence = index[j];
2963        double thisWeight = weight[iSequence];
2964        // row has -1
2965        double pivot = updateBy[iSequence] * scaleFactor;
2966        updateBy[iSequence] = 0.0;
2967        double value = pivot * pivot * devex_;
2968        if (reference(iSequence))
2969          value += 1.0;
2970        weight[iSequence] = CoinMax(0.99 * thisWeight, value);
2971      }
2972    }
2973    // restore outgoing weight
2974    if (sequenceOut >= 0)
2975      weights_[sequenceOut] = outgoingWeight;
2976    alternateWeights_->clear();
2977    spareColumn2->setNumElements(0);
2978    //#define SOME_DEBUG_1
2979#ifdef SOME_DEBUG_1
2980    // check for accuracy
2981    int iCheck = 892;
2982    //printf("weight for iCheck is %g\n",weights_[iCheck]);
2983    int numberRows = model_->numberRows();
2984    //int numberColumns = model_->numberColumns();
2985    for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
2986      if (model_->getStatus(iCheck) != ClpSimplex::basic && model_->getStatus(iCheck) != ClpSimplex::isFixed)
2987        checkAccuracy(iCheck, 1.0e-1, updates, spareRow2);
2988    }
2989#endif
2990    updates->setNumElements(0);
2991    spareColumn1->setNumElements(0);
2992  }
2993
2994  // update of duals finished - now do pricing
2995
2996  double bestDj = 1.0e-30;
2997  int bestSequence = -1;
2998
2999  int i, iSequence;
3000  index = infeasible_->getIndices();
3001  number = infeasible_->getNumElements();
3002  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
3003    // we can't really trust infeasibilities if there is dual error
3004    double checkTolerance = 1.0e-8;
3005    if (!model_->factorization()->pivots())
3006      checkTolerance = 1.0e-6;
3007    if (model_->largestDualError() > checkTolerance)
3008      tolerance *= model_->largestDualError() / checkTolerance;
3009    // But cap
3010    tolerance = CoinMin(1000.0, tolerance);
3011  }
3012#ifdef CLP_DEBUG
3013  if (model_->numberDualInfeasibilities() == 1)
3014    printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
3015      model_->largestDualError(), model_, model_->messageHandler(),
3016      number);
3017#endif
3018  // stop last one coming immediately
3019  double saveOutInfeasibility = 0.0;
3020  if (sequenceOut >= 0) {
3021    saveOutInfeasibility = infeas[sequenceOut];
3022    infeas[sequenceOut] = 0.0;
3023  }
3024  tolerance *= tolerance; // as we are using squares
3025
3026  int iPass;
3027  // Setup two passes
3028  int start[4];
3029  start[1] = number;
3030  start[2] = 0;
3031  double dstart = static_cast< double >(number) * model_->randomNumberGenerator()->randomDouble();
3032  start[0] = static_cast< int >(dstart);
3033  start[3] = start[0];
3034  //double largestWeight=0.0;
3035  //double smallestWeight=1.0e100;
3036  for (iPass = 0; iPass < 2; iPass++) {
3037    int end = start[2 * iPass + 1];
3038    if (switchType < 5) {
3039      for (i = start[2 * iPass]; i < end; i++) {
3040        iSequence = index[i];
3041        double value = infeas[iSequence];
3042        if (value > tolerance) {
3043          double weight = weights_[iSequence];
3044          //weight=1.0;
3045          if (value > bestDj * weight) {
3046            // check flagged variable and correct dj
3047            if (!model_->flagged(iSequence)) {
3048              bestDj = value / weight;
3049              bestSequence = iSequence;
3050            } else {
3051              // just to make sure we don't exit before got something
3052              numberWanted++;
3053            }
3054          }
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)
3089       printf("%d best %g\n",bestSequence,bestDj);*/
3090  reducedCost = model_->djRegion();
3091  model_->clpMatrix()->setSavedBestSequence(bestSequence);
3092  if (bestSequence >= 0)
3093    model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
3094
3095#ifdef CLP_DEBUG
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  }
3102#endif
3103  return bestSequence;
3104}
3105// Called when maximum pivots changes
3106void ClpPrimalColumnSteepest::maximumPivotsChanged()
3107{
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  }
3114}
3115void ClpPrimalColumnSteepest::redoInfeasibilities()
3116{
3117  double *COIN_RESTRICT infeas = infeasible_->denseVector();
3118  int *COIN_RESTRICT index = infeasible_->getIndices();
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
3124  tolerance = tolerance + error;
3125  // reverse sign so test is cleaner
3126  tolerance = -tolerance;
3127  int number = model_->numberRows() + model_->numberColumns();
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;
3133    double value = reducedCost[iSequence];
3134    infeas[iSequence] = 0.0;
3135    if (thisStatus == 3) {
3136    } else if ((thisStatus & 1) != 0) {
3137      // basic or fixed
3138      value = 0.0;
3139    } else if (thisStatus == 2) {
3140      value = -value;
3141    } else {
3142      // free or superbasic
3143      if (fabs(value) > FREE_ACCEPT * -tolerance) {
3144        // we are going to bias towards free (but only if reasonable)
3145        value = -fabs(value) * FREE_BIAS;
3146      } else {
3147        value = 0.0;
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}
3159/*
3160   1) before factorization
3161   2) after factorization
3162   3) just redo infeasibilities
3163   4) restore weights
3164   5) at end of values pass (so need initialization)
3165*/
3166void ClpPrimalColumnSteepest::saveWeights(ClpSimplex *model, int mode)
3167{
3168  model_ = model;
3169  if (mode == 6) {
3170    // If incoming weight is 1.0 then return else as 5
3171    assert(weights_);
3172    int sequenceIn = model_->sequenceIn();
3173    assert(sequenceIn >= 0 && sequenceIn < model_->numberRows() + model_->numberColumns());
3174    if (weights_[sequenceIn] == (mode_ != 1) ? 1.0 : 1.0 + ADD_ONE)
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;
3192#if ABOCA_LITE
3193    int numberThreads = abcState();
3194    if (numberThreads && mode_ > 1)
3195      mode_ = 0; // force exact
3196#endif
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());
3206#endif
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;
3272
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);
3331#ifdef CLP_DEBUG
3332        // Can happen but I should clean up
3333        assert(found >= 0);
3334#endif
3335        pivotSequence_ = found;
3336        for (iRow = 0; iRow < numberRows; iRow++) {
3337          int iPivot = savePivotOrder[iRow];
3338          temp[iPivot] = 0.0;
3339        }
3340#else
3341        for (int iRow = 0; iRow < numberRows; iRow++) {
3342          int iPivot = pivotVariable[iRow];
3343          if (iPivot == pivotSequence_) {
3344            pivotSequence_ = iRow;
3345            break;
3346          }
3347        }
3348#endif
3349      } else {
3350        // Just clean up
3351        /* If this happens when alternateWeights_ is
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 */
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;
3376
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();
3382
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      }
3410#else
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      }
3470#endif
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);
3479
3480        switch (status) {
3481
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);
3492          }
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      }
3519    }
3520  }
3521}
3522// Gets rid of last update
3523void ClpPrimalColumnSteepest::unrollWeights()
3524{
3525  if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3526    return;
3527  double *saved = alternateWeights_->denseVector();
3528  int number = alternateWeights_->getNumElements();
3529  int *which = alternateWeights_->getIndices();
3530  int i;
3531  for (i = 0; i < number; i++) {
3532    int iRow = which[i];
3533    weights_[iRow] = saved[iRow];
3534    saved[iRow] = 0.0;
3535  }
3536  alternateWeights_->setNumElements(0);
3537}
3538
3539//-------------------------------------------------------------------
3540// Clone
3541//-------------------------------------------------------------------
3542ClpPrimalColumnPivot *ClpPrimalColumnSteepest::clone(bool CopyData) const
3543{
3544  if (CopyData) {
3545    return new ClpPrimalColumnSteepest(*this);
3546  } else {
3547    return new ClpPrimalColumnSteepest();
3548  }
3549}
3550void ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector *input)
3551{
3552  // Local copy of mode so can decide what to do
3553  int switchType = mode_;
3554  if (mode_ == 4 && numberSwitched_)
3555    switchType = 3;
3556  else if (mode_ == 4 || mode_ == 5)
3557    return;
3558  int number = input->getNumElements();
3559  int *which = input->getIndices();
3560  double *work = input->denseVector();
3561#if ALT_UPDATE_WEIGHTS != 2
3562  int newNumber = 0;
3563  int *newWhich = alternateWeights_->getIndices();
3564  double *newWork = alternateWeights_->denseVector();
3565#endif
3566#ifdef ALT_UPDATE_WEIGHTSz
3567  {
3568    //int newNumber2 = 0;
3569    if (!altVector[0]) {
3570      altVector[0] = new CoinIndexedVector(2000);
3571      altVector[1] = new CoinIndexedVector(2000);
3572      altVector[2] = new CoinIndexedVector(2000);
3573    }
3574    altVector[0]->clear();
3575    // get updated pivot row
3576    int pivotRow = model_->pivotRow();
3577    // should it be - or what
3578    altVector[0]->quickAdd(pivotRow, model_->dualIn());
3579    model_->factorization()->updateColumnTranspose(altVector[2],
3580      altVector[0]);
3581    double *work2 = altVector[0]->denseVector();
3582    //altVector[1]->clear();
3583    int *newWhich2 = altVector[1]->getIndices();
3584    double *newWork2 = altVector[1]->denseVector();
3585    int number2 = altVector[1]->getNumElements();
3586    int nRow = model_->numberRows();
3587    int nCol = model_->numberColumns();
3588    int nTotal = nRow + nCol;
3589    double *temp = new double[2 * nTotal + nRow];
3590    memset(temp, 0, (2 * nTotal + nRow) * sizeof(double));
3591    double *pivRow = temp + nTotal;
3592    double *temp2 = temp + nCol;
3593    double *temp2P = pivRow + nCol;
3594    double *piU = pivRow + nTotal;
3595    double devex = 0.0;
3596    double scaleFactor = 1.0 / model_->dualIn();
3597    const int *pivotVariable = model_->pivotVariable();
3598    for (int i = 0; i < number; i++) {
3599      int iRow = which[i];
3600      int iPivot = pivotVariable[iRow];
3601      if (reference(iPivot)) {
3602        devex += work[iRow] * work[iRow];
3603      }
3604    }
3605    int sequenceIn = model_->sequenceIn();
3606    int sequenceOut = model_->sequenceOut();
3607    for (int i = 0; i < number2; i++) {
3608      int iRow = newWhich2[i];
3609      temp2[iRow] = newWork2[iRow];
3610    }
3611    //if (!input->packedMode()) {
3612    for (int i = 0; i < number; i++) {
3613      int iRow = which[i];
3614      piU[iRow] = work2[iRow];
3615      temp2P[iRow] = work2[iRow];
3616    }
3617    double alpha = model_->alpha();
3618    const int *row = model_->matrix()->getIndices();
3619    const CoinBigIndex *columnStart = model_->matrix()->getVectorStarts();
3620    const int *columnLength = model_->matrix()->getVectorLengths();
3621    const double *element = model_->matrix()->getElements();
3622    for (int i = 0; i < nCol; i++) {
3623      CoinBigIndex start = columnStart[i];
3624      CoinBigIndex end = start + columnLength[i];
3625      double value = 0.0;
3626      double value2 = 0.0;
3627      for (CoinBigIndex j = start; j < end; j++) {
3628        int iRow = row[j];
3629        value -= piU[iRow] * element[j];
3630        value2 -= newWork2[iRow] * element[j];
3631      }
3632      pivRow[i] = value;
3633      temp[i] = value2;
3634    }
3635    const unsigned char *COIN_RESTRICT statusArray = model_->statusArray();
3636    for (int i = 0; i < nTotal; i++) {
3637      unsigned char thisStatus = statusArray[i] & 7;
3638#if 0
3639        if (thisStatus==3) {
3640          // lower
3641        } else if ((thisStatus&1)!=0) {
3642          // basic or fixed
3643          //value=0.0;
3644        } else if (thisStatus==2) {
3645          // uppervalue=-value;
3646        } else {
3647          // free or superbasic
3648          //if (fabs(value) > FREE_ACCEPT * -dualTolerance) {
3649            // we are going to bias towards free (but only if reasonable)
3650            //value = -fabs(value)*FREE_BIAS;
3651          //} else {
3652          //value=0.0;
3653          //}
3654        }
3655#else
3656      if (thisStatus != 1 && thisStatus != 5) {
3657        double pivot = pivRow[i] * scaleFactor;
3658        double modification = temp[i];
3659        double thisWeight = weights_[i];
3660        double pivotSquared = pivot * pivot;
3661        double newWeight = thisWeight + pivotSquared * devexA - 2.0 * pivot * modification;
3662        temp[i] = newWeight;
3663      } else {
3664        temp[i] = 1.0;
3665      }
3666#endif
3667    }
3668    temp[sequenceOut] = devexA / (alpha * alpha);
3669    // to keep clean for debug
3670#ifndef NDEBUG
3671    {
3672      if (sequenceOut < nCol) {
3673        if (model_->columnLower()[sequenceOut] == model_->columnUpper()[sequenceOut])
3674          temp[sequenceOut] = 1.0;
3675      } else {
3676        if (model_->rowLower()[sequenceOut - nCol] == model_->rowUpper()[sequenceOut - nCol])
3677          temp[sequenceOut] = 1.0;
3678      }
3679    }
3680#endif
3681    temp[sequenceIn] = 1.0;
3682    for (int i = 0; i < nTotal; i++) {
3683      printf("%g ", temp[i]);
3684      if ((i % 10) == 9)
3685        printf("\n");
3686    }
3687    if (((nTotal - 1) % 10) != 9)
3688      printf("\n");
3689    delete[] temp;
3690  }
3691#endif
3692  int i;
3693  int sequenceIn = model_->sequenceIn();
3694  int sequenceOut = model_->sequenceOut();
3695  const int *pivotVariable = model_->pivotVariable();
3696
3697  int pivotRow = model_->pivotRow();
3698  pivotSequence_ = pivotRow;
3699
3700  devex_ = 0.0;
3701  // Can't create alternateWeights_ as packed as needed unpacked
3702  if (!input->packedMode()) {
3703    if (pivotRow >= 0) {
3704      if (switchType == 1) {
3705        for (i = 0; i < number; i++) {
3706          int iRow = which[i];
3707          devex_ += work[iRow] * work[iRow];
3708#if ALT_UPDATE_WEIGHTS != 2
3709          newWork[iRow] = -2.0 * work[iRow];
3710#endif
3711        }
3712#if ALT_UPDATE_WEIGHTS != 2
3713        newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3714#endif
3715        devex_ += ADD_ONE;
3716        weights_[sequenceOut] = 1.0 + ADD_ONE;
3717#if ALT_UPDATE_WEIGHTS != 2
3718        CoinMemcpyN(which, number, newWhich);
3719        alternateWeights_->setNumElements(number);
3720#endif
3721      } else {
3722        if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3723          for (i = 0; i < number; i++) {
3724            int iRow = which[i];
3725            int iPivot = pivotVariable[iRow];
3726            if (reference(iPivot)) {
3727              devex_ += work[iRow] * work[iRow];
3728#if ALT_UPDATE_WEIGHTS != 2
3729              newWork[iRow] = -2.0 * work[iRow];
3730              newWhich[newNumber++] = iRow;
3731#endif
3732            }
3733          }
3734#if ALT_UPDATE_WEIGHTS != 2
3735          if (!newWork[pivotRow] && devex_ > 0.0)
3736            newWhich[newNumber++] = pivotRow; // add if not already in
3737          newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3738#endif
3739        } else {
3740          for (i = 0; i < number; i++) {
3741            int iRow = which[i];
3742            int iPivot = pivotVariable[iRow];
3743            if (reference(iPivot))
3744              devex_ += work[iRow] * work[iRow];
3745          }
3746        }
3747        if (reference(sequenceIn)) {
3748          devex_ += 1.0;
3749        } else {
3750        }
3751        if (reference(sequenceOut)) {
3752          weights_[sequenceOut] = 1.0 + 1.0;
3753        } else {
3754          weights_[sequenceOut] = 1.0;
3755        }
3756#if ALT_UPDATE_WEIGHTS != 2
3757        alternateWeights_->setNumElements(newNumber);
3758#endif
3759      }
3760    } else {
3761      if (switchType == 1) {
3762        for (i = 0; i < number; i++) {
3763          int iRow = which[i];
3764          devex_ += work[iRow] * work[iRow];
3765        }
3766        devex_ += ADD_ONE;
3767      } else {
3768        for (i = 0; i < number; i++) {
3769          int iRow = which[i];
3770          int iPivot = pivotVariable[iRow];
3771          if (reference(iPivot)) {
3772            devex_ += work[iRow] * work[iRow];
3773          }
3774        }
3775        if (reference(sequenceIn))
3776          devex_ += 1.0;
3777      }
3778    }
3779  } else {
3780    // packed input
3781    if (pivotRow >= 0) {
3782      if (switchType == 1) {
3783        for (i = 0; i < number; i++) {
3784#if ALT_UPDATE_WEIGHTS != 2
3785          int iRow = which[i];
3786#endif
3787          devex_ += work[i] * work[i];
3788#if ALT_UPDATE_WEIGHTS != 2
3789          newWork[iRow] = -2.0 * work[i];
3790#endif
3791        }
3792#if ALT_UPDATE_WEIGHTS != 2
3793        newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3794#endif
3795        devex_ += ADD_ONE;
3796        weights_[sequenceOut] = 1.0 + ADD_ONE;
3797#if ALT_UPDATE_WEIGHTS != 2
3798        CoinMemcpyN(which, number, newWhich);
3799        alternateWeights_->setNumElements(number);
3800#endif
3801      } else {
3802        if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3803          for (i = 0; i < number; i++) {
3804            int iRow = which[i];
3805            int iPivot = pivotVariable[iRow];
3806            if (reference(iPivot)) {
3807              devex_ += work[i] * work[i];
3808#if ALT_UPDATE_WEIGHTS != 2
3809              newWork[iRow] = -2.0 * work[i];
3810              newWhich[newNumber++] = iRow;
3811#endif
3812            }
3813          }
3814#if ALT_UPDATE_WEIGHTS != 2
3815          if (!newWork[pivotRow] && devex_ > 0.0)
3816            newWhich[newNumber++] = pivotRow; // add if not already in
3817          newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3818#endif
3819        } else {
3820          for (i = 0; i < number; i++) {
3821            int iRow = which[i];
3822            int iPivot = pivotVariable[iRow];
3823            if (reference(iPivot))
3824              devex_ += work[i] * work[i];
3825          }
3826        }
3827        if (reference(sequenceIn)) {
3828          devex_ += 1.0;
3829        } else {
3830        }
3831        if (reference(sequenceOut)) {
3832          weights_[sequenceOut] = 1.0 + 1.0;
3833        } else {
3834          weights_[sequenceOut] = 1.0;
3835        }
3836#if ALT_UPDATE_WEIGHTS != 2
3837        alternateWeights_->setNumElements(newNumber);
3838#endif
3839      }
3840    } else {
3841      if (switchType == 1) {
3842        for (i = 0; i < number; i++) {
3843          devex_ += work[i] * work[i];
3844        }
3845        devex_ += ADD_ONE;
3846      } else {
3847        for (i = 0; i < number; i++) {
3848          int iRow = which[i];
3849          int iPivot = pivotVariable[iRow];
3850          if (reference(iPivot)) {
3851            devex_ += work[i] * work[i];
3852          }
3853        }
3854        if (reference(sequenceIn))
3855          devex_ += 1.0;
3856      }
3857    }
3858  }
3859  if (devex_ < 1.001e-30) {
3860    COIN_DETAIL_PRINT(printf("devex of incoming tiny %d %g\n", sequenceIn, devex_));
3861    devex_ = 1.0e-30;
3862  }
3863  double oldDevex = weights_[sequenceIn];
3864#ifdef CLP_DEBUG
3865  if ((model_->messageHandler()->logLevel() & 32))
3866    printf("old weight %g, new %g\n", oldDevex, devex_);
3867#endif
3868  double check = CoinMax(devex_, oldDevex) + 0.1;
3869  weights_[sequenceIn] = devex_;
3870  double testValue = 0.1;
3871  if (mode_ == 4 && numberSwitched_ == 1)
3872    testValue = 0.5;
3873  if (fabs(devex_ - oldDevex) > testValue * check) {
3874#ifdef CLP_DEBUG
3875    if ((model_->messageHandler()->logLevel() & 48) == 16)
3876      printf("old weight %g, new %g\n", oldDevex, devex_);
3877#endif
3878    //printf("old weight %g, new %g\n",oldDevex,devex_);
3879    testValue = 0.99;
3880    if (mode_ == 1)
3881      testValue = 1.01e1; // make unlikely to do if steepest
3882    else if (mode_ == 4 && numberSwitched_ == 1)
3883      testValue = 0.9;
3884    double difference = fabs(devex_ - oldDevex);
3885    if (difference > testValue * check) {
3886      // need to redo
3887      model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
3888        *model_->messagesPointer())
3889        << oldDevex << devex_
3890        << CoinMessageEol;
3891      initializeWeights();
3892      // redo devex_
3893      if (pivotRow >= 0)
3894        devex_ = 1.0;
3895    }
3896  }
3897  if (pivotRow >= 0) {
3898    // set outgoing weight here
3899    double alpha = model_->alpha();
3900    if (fabs(alpha) > 1.0e15) {
3901      COIN_DETAIL_PRINT(printf("alpha %g for %d !!\n", alpha, model_->sequenceOut()));
3902      alpha = 1.0e15;
3903    }
3904    weights_[model_->sequenceOut()] = devex_ / (alpha * alpha);
3905  }
3906}
3907// Checks accuracy - just for debug
3908void ClpPrimalColumnSteepest::checkAccuracy(int sequence,
3909  double relativeTolerance,
3910  CoinIndexedVector *rowArray1,
3911  CoinIndexedVector *rowArray2)
3912{
3913  if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3914    return;
3915  model_->unpack(rowArray1, sequence);
3916  model_->factorization()->updateColumn(rowArray2, rowArray1);
3917  int number = rowArray1->getNumElements();
3918  int *which = rowArray1->getIndices();
3919  double *work = rowArray1->denseVector();
3920  const int *pivotVariable = model_->pivotVariable();
3921
3922  double devex = 0.0;
3923  int i;
3924
3925  if (mode_ == 1) {
3926    for (i = 0; i < number; i++) {
3927      int iRow = which[i];
3928      devex += work[iRow] * work[iRow];
3929      work[iRow] = 0.0;
3930    }
3931    devex += ADD_ONE;
3932  } else {
3933    for (i = 0; i < number; i++) {
3934      int iRow = which[i];
3935      int iPivot = pivotVariable[iRow];
3936      if (reference(iPivot)) {
3937        devex += work[iRow] * work[iRow];
3938      }
3939      work[iRow] = 0.0;
3940    }
3941    if (reference(sequence))
3942      devex += 1.0;
3943  }
3944
3945  double oldDevex = CoinMax(weights_[sequence], 1.0e-4);
3946  devex = CoinMax(devex, 1.0e-4);
3947  double check = CoinMax(devex, oldDevex);
3948  ;
3949  rowArray1->setNumElements(0);
3950  if (fabs(devex - oldDevex) > relativeTolerance * check) {
3951    //COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
3952    printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex);
3953    if (mode_ == 0) {
3954      rowArray1->setNumElements(0);
3955      model_->unpack(rowArray1, sequence);
3956      number = rowArray1->getNumElements();
3957      for (i = 0; i < number; i++)
3958        printf("(%d,%g) ", which[i], work[which[i]]);
3959      printf("\n");
3960      model_->factorization()->updateColumn(rowArray2, rowArray1);
3961      number = rowArray1->getNumElements();
3962      for (i = 0; i < number; i++)
3963        printf("(%d,%g) ", which[i], work[which[i]]);
3964      printf("\n");
3965      devex = 0.0;
3966      for (i = 0; i < number; i++) {
3967        int iRow = which[i];
3968        int iPivot = pivotVariable[iRow];
3969        if (reference(iPivot)) {
3970          devex += work[iRow] * work[iRow];
3971        }
3972        work[iRow] = 0.0;
3973      }
3974      if (reference(sequence))
3975        devex += 1.0;
3976    }
3977    // update so won't print again
3978    weights_[sequence] = devex;
3979  }
3980}
3981
3982// Initialize weights
3983void ClpPrimalColumnSteepest::initializeWeights()
3984{
3985  int numberRows = model_->numberRows();
3986  int numberColumns = model_->numberColumns();
3987  int number = numberRows + numberColumns;
3988  int iSequence;
3989  if (mode_ != 1) {
3990    // initialize to 1.0
3991    // and set reference framework
3992    if (!reference_) {
3993      int nWords = (number + 31) >> 5;
3994      reference_ = new unsigned int[nWords];
3995      CoinZeroN(reference_, nWords);
3996    }
3997
3998    for (iSequence = 0; iSequence < number; iSequence++) {
3999      weights_[iSequence] = 1.0;
4000      if (model_->getStatus(iSequence) == ClpSimplex::basic) {
4001        setReference(iSequence, false);
4002      } else {
4003        setReference(iSequence, true);
4004      }
4005    }
4006  } else {
4007    CoinIndexedVector *temp = new CoinIndexedVector();
4008    temp->reserve(numberRows + model_->factorization()->maximumPivots());
4009    double *array = alternateWeights_->denseVector();
4010    int *which = alternateWeights_->getIndices();
4011
4012    for (iSequence = 0; iSequence < number; iSequence++) {
4013      weights_[iSequence] = 1.0 + ADD_ONE;
4014      if (model_->getStatus(iSequence) != ClpSimplex::basic && model_->getStatus(iSequence) != ClpSimplex::isFixed) {
4015        model_->unpack(alternateWeights_, iSequence);
4016        double value = ADD_ONE;
4017        model_->factorization()->updateColumn(temp, alternateWeights_);
4018        int number = alternateWeights_->getNumElements();
4019        int j;
4020        for (j = 0; j < number; j++) {
4021          int iRow = which[j];
4022          value += array[iRow] * array[iRow];
4023          array[iRow] = 0.0;
4024        }
4025        alternateWeights_->setNumElements(0);
4026        weights_[iSequence] = value;
4027      }
4028    }
4029    delete temp;
4030  }
4031}
4032// Gets rid of all arrays
4033void ClpPrimalColumnSteepest::clearArrays()
4034{
4035  if (persistence_ == normal) {
4036    delete[] weights_;
4037    weights_ = NULL;
4038    delete infeasible_;
4039    infeasible_ = NULL;
4040    delete alternateWeights_;
4041    alternateWeights_ = NULL;
4042    delete[] savedWeights_;
4043    savedWeights_ = NULL;
4044    delete[] reference_;
4045    reference_ = NULL;
4046  }
4047  pivotSequence_ = -1;
4048  state_ = -1;
4049  savedPivotSequence_ = -1;
4050  savedSequenceOut_ = -1;
4051  devex_ = 0.0;
4052}
4053// Returns true if would not find any column
4054bool ClpPrimalColumnSteepest::looksOptimal() const
4055{
4056  if (looksOptimal_)
4057    return true; // user overrode
4058  //**** THIS MUST MATCH the action coding above
4059  double tolerance = model_->currentDualTolerance();
4060  // we can't really trust infeasibilities if there is dual error
4061  // this coding has to mimic coding in checkDualSolution
4062  double error = CoinMin(1.0e-2, model_->largestDualError());
4063  // allow tolerance at least slightly bigger than standard
4064  tolerance = tolerance + error;
4065  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
4066    // we can't really trust infeasibilities if there is dual error
4067    double checkTolerance = 1.0e-8;
4068    if (!model_->factorization()->pivots())
4069      checkTolerance = 1.0e-6;
4070    if (model_->largestDualError() > checkTolerance)
4071      tolerance *= model_->largestDualError() / checkTolerance;
4072    // But cap
4073    tolerance = CoinMin(1000.0, tolerance);
4074  }
4075  int number = model_->numberRows() + model_->numberColumns();
4076  int iSequence;
4077
4078  double *reducedCost = model_->djRegion();
4079  int numberInfeasible = 0;
4080  if (!model_->nonLinearCost()->lookBothWays()) {
4081    for (iSequence = 0; iSequence < number; iSequence++) {
4082      double value = reducedCost[iSequence];
4083      ClpSimplex::Status status = model_->getStatus(iSequence);
4084
4085      switch (status) {
4086
4087      case ClpSimplex::basic:
4088      case ClpSimplex::isFixed:
4089        break;
4090      case ClpSimplex::isFree:
4091      case ClpSimplex::superBasic:
4092        if (fabs(value) > FREE_ACCEPT * tolerance)
4093          numberInfeasible++;
4094        break;
4095      case ClpSimplex::atUpperBound:
4096        if (value > tolerance)
4097          numberInfeasible++;
4098        break;
4099      case ClpSimplex::atLowerBound:
4100        if (value < -tolerance)
4101          numberInfeasible++;
4102      }
4103    }
4104  } else {
4105    ClpNonLinearCost *nonLinear = model_->nonLinearCost();
4106    // can go both ways
4107    for (iSequence = 0; iSequence < number; iSequence++) {
4108      double value = reducedCost[iSequence];
4109      ClpSimplex::Status status = model_->getStatus(iSequence);
4110
4111      switch (status) {
4112
4113      case ClpSimplex::basic:
4114      case ClpSimplex::isFixed:
4115        break;
4116      case ClpSimplex::isFree:
4117      case ClpSimplex::superBasic:
4118        if (fabs(value) > FREE_ACCEPT * tolerance)
4119          numberInfeasible++;
4120        break;
4121      case ClpSimplex::atUpperBound:
4122        if (value > tolerance) {
4123          numberInfeasible++;
4124        } else {
4125          // look other way - change up should be negative
4126          value -= nonLinear->changeUpInCost(iSequence);
4127          if (value < -tolerance)
4128            numberInfeasible++;
4129        }
4130        break;
4131      case ClpSimplex::atLowerBound:
4132        if (value < -tolerance) {
4133          numberInfeasible++;
4134        } else {
4135          // look other way - change down should be positive
4136          value -= nonLinear->changeDownInCost(iSequence);
4137          if (value > tolerance)
4138            numberInfeasible++;
4139        }
4140      }
4141    }
4142  }
4143  return numberInfeasible == 0;
4144}
4145/* Returns number of extra columns for sprint algorithm - 0 means off.
4146   Also number of iterations before recompute
4147*/
4148int ClpPrimalColumnSteepest::numberSprintColumns(int &numberIterations) const
4149{
4150  numberIterations = 0;
4151  int numberAdd = 0;
4152  if (!numberSwitched_ && mode_ >= 10) {
4153    numberIterations = CoinMin(2000, model_->numberRows() / 5);
4154    numberIterations = CoinMax(numberIterations, model_->factorizationFrequency());
4155    numberIterations = CoinMax(numberIterations, 500);
4156    if (mode_ == 10) {
4157      numberAdd = CoinMax(300, model_->numberColumns() / 10);
4158      numberAdd = CoinMax(numberAdd, model_->numberRows() / 5);
4159      // fake all
4160      //numberAdd=1000000;
4161      numberAdd = CoinMin(numberAdd, model_->numberColumns());
4162    } else {
4163      abort();
4164    }
4165  }
4166  return numberAdd;
4167}
4168// Switch off sprint idea
4169void ClpPrimalColumnSteepest::switchOffSprint()
4170{
4171  numberSwitched_ = 10;
4172}
4173// Update djs doing partial pricing (dantzig)
4174int ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector *updates,
4175  CoinIndexedVector *spareRow2,
4176  int numberWanted,
4177  int numberLook)
4178{
4179  int number = 0;
4180  int *index;
4181  double *updateBy;
4182  double *reducedCost;
4183  double saveTolerance = model_->currentDualTolerance();
4184  double tolerance = model_->currentDualTolerance();
4185  // we can't really trust infeasibilities if there is dual error
4186  // this coding has to mimic coding in checkDualSolution
4187  double error = CoinMin(1.0e-2, model_->largestDualError());
4188  // allow tolerance at least slightly bigger than standard
4189  tolerance = tolerance + error;
4190  if (model_->numberIterations() < model_->lastBadIteration() + 200) {
4191    // we can't really trust infeasibilities if there is dual error
4192    double checkTolerance = 1.0e-8;
4193    if (!model_->factorization()->pivots())
4194      checkTolerance = 1.0e-6;
4195    if (model_->largestDualError() > checkTolerance)
4196      tolerance *= model_->largestDualError() / checkTolerance;
4197    // But cap
4198    tolerance = CoinMin(1000.0, tolerance);
4199  }
4200  if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
4201    tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
4202  // So partial pricing can use
4203  model_->setCurrentDualTolerance(tolerance);
4204  model_->factorization()->updateColumnTranspose(spareRow2, updates);
4205  int numberColumns = model_->numberColumns();
4206
4207  // Rows
4208  reducedCost = model_->djRegion(0);
4209
4210  number = updates->getNumElements();
4211  index = updates->getIndices();
4212  updateBy = updates->denseVector();
4213  int j;
4214  double *duals = model_->dualRowSolution();
4215  for (j = 0; j < number; j++) {
4216    int iSequence = index[j];
4217    double value = duals[iSequence];
4218    value -= updateBy[j];
4219    updateBy[j] = 0.0;
4220    duals[iSequence] = value;
4221  }
4222  //#define CLP_DEBUG
4223#ifdef CLP_DEBUG
4224  // check duals
4225  {
4226    int numberRows = model_->numberRows();
4227    //work space
4228    CoinIndexedVector arrayVector;
4229    arrayVector.reserve(numberRows + 1000);
4230    CoinIndexedVector workSpace;
4231    workSpace.reserve(numberRows + 1000);
4232
4233    int iRow;
4234    double *array = arrayVector.denseVector();
4235    int *index = arrayVector.getIndices();
4236    int number = 0;
4237    int *pivotVariable = model_->pivotVariable();
4238    double *cost = model_->costRegion();
4239    for (iRow = 0; iRow < numberRows; iRow++) {
4240      int iPivot = pivotVariable[iRow];
4241      double value = cost[iPivot];
4242      if (value) {
4243        array[iRow] = value;
4244        index[number++] = iRow;
4245      }
4246    }
4247    arrayVector.setNumElements(number);
4248    // Extended duals before "updateTranspose"
4249    model_->clpMatrix()->dualExpanded(model_, &arrayVector, NULL, 0);
4250
4251    // Btran basic costs
4252    model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
4253
4254    // now look at dual solution
4255    for (iRow = 0; iRow < numberRows; iRow++) {
4256      // slack
4257      double value = array[iRow];
4258      if (fabs(duals[iRow] - value) > 1.0e-3)
4259        printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
4260      //duals[iRow]=value;
4261    }
4262  }
4263#endif
4264#undef CLP_DEBUG
4265  double bestDj = tolerance;
4266  int bestSequence = -1;
4267
4268  const double *cost = model_->costRegion(1);
4269
4270  model_->clpMatrix()->setOriginalWanted(numberWanted);
4271  model_->clpMatrix()->setCurrentWanted(numberWanted);
4272  int iPassR = 0, iPassC = 0;
4273  // Setup two passes
4274  // This biases towards picking row variables
4275  // This probably should be fixed
4276  int startR[4];
4277  int numberRows = model_->numberRows();
4278  startR[1] = numberColumns+numberRows;
4279  startR[2] = numberColumns;
4280  double randomR = model_->randomNumberGenerator()->randomDouble();
4281  double dstart = static_cast<double> (numberRows) * randomR;
4282  startR[0] = numberColumns + static_cast<int> (dstart);
4283  startR[3] = startR[0];
4284  double startC[4];
4285  startC[1] = 1.0;
4286  startC[2] = 0;
4287  double randomC = model_->randomNumberGenerator()->randomDouble();
4288  startC[0] = randomC;
4289  startC[3] = randomC;
4290  reducedCost = model_->djRegion(1);
4291  int sequenceOut = model_->sequenceOut();
4292  double *duals2 = duals - numberColumns;
4293  int chunk = CoinMin(1024, (numberColumns + numberRows) / 32);
4294#ifdef COIN_DETAIL
4295  if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
4296    printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
4297    int i;
4298    for (i = 0; i < 4; i++)
4299      printf("start R %d C %g ", startR[i], startC[i]);
4300    printf("\n");
4301  }
4302#endif
4303  chunk = CoinMax(chunk, 256);
4304  bool finishedR = false, finishedC = false;
4305  bool doingR = randomR > randomC;
4306  //doingR=false;
4307  int saveNumberWanted = numberWanted;
4308  while (!finishedR || !finishedC) {
4309    if (finishedR)
4310      doingR = false;
4311    if (doingR) {
4312      int saveSequence = bestSequence;
4313      int start = startR[iPassR];
4314      int end = CoinMin(startR[iPassR + 1], start + chunk / 2);
4315      int iSequence;
4316      for (iSequence = start; iSequence < end; iSequence++) {
4317        if (iSequence != sequenceOut) {
4318          double value;
4319          ClpSimplex::Status status = model_->getStatus(iSequence);
4320
4321          switch (status) {
4322
4323          case ClpSimplex::basic:
4324          case ClpSimplex::isFixed:
4325            break;
4326          case ClpSimplex::isFree:
4327          case ClpSimplex::superBasic:
4328            value = fabs(cost[iSequence] + duals2[iSequence]);
4329            if (value > FREE_ACCEPT * tolerance) {
4330              numberWanted--;
4331              // we are going to bias towards free (but only if reasonable)
4332              value *= FREE_BIAS;
4333              if (value > bestDj) {
4334                // check flagged variable and correct dj
4335                if (!model_->flagged(iSequence)) {
4336                  bestDj = value;
4337                  bestSequence = iSequence;
4338                } else {
4339                  // just to make sure we don't exit before got something
4340                  numberWanted++;
4341                }
4342              }
4343            }
4344            break;
4345          case ClpSimplex::atUpperBound:
4346            value = cost[iSequence] + duals2[iSequence];
4347            if (value > tolerance) {
4348              numberWanted--;
4349              if (value > bestDj) {
4350                // check flagged variable and correct dj
4351                if (!model_->flagged(iSequence)) {
4352                  bestDj = value;
4353                  bestSequence = iSequence;
4354                } else {
4355                  // just to make sure we don't exit before got something
4356                  numberWanted++;
4357                }
4358              }
4359            }
4360            break;
4361          case ClpSimplex::atLowerBound:
4362            value = -(cost[iSequence] + duals2[iSequence]);
4363            if (value > tolerance) {
4364              numberWanted--;
4365              if (value > bestDj) {
4366                // check flagged variable and correct dj
4367                if (!model_->flagged(iSequence)) {
4368                  bestDj = value;
4369                  bestSequence = iSequence;
4370                } else {
4371                  // just to make sure we don't exit before got something
4372                  numberWanted++;
4373                }
4374              }
4375            }
4376            break;
4377          }
4378        }
4379        if (!numberWanted)
4380          break;
4381      }
4382      numberLook -= (end - start);
4383      if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
4384        numberWanted = 0; // give up
4385      if (saveSequence != bestSequence) {
4386        // dj
4387        reducedCost[bestSequence] = cost[bestSequence] + duals[bestSequence - numberColumns];
4388        bestDj = fabs(reducedCost[bestSequence]);
4389        model_->clpMatrix()->setSavedBestSequence(bestSequence);
4390        model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
4391      }
4392      model_->clpMatrix()->setCurrentWanted(numberWanted);
4393      if (!numberWanted)
4394        break;
4395      doingR = false;
4396      // update start
4397      startR[iPassR] = iSequence;
4398      if (iSequence >= startR[iPassR + 1]) {
4399        if (iPassR)
4400          finishedR = true;
4401        else
4402          iPassR = 2;
4403      }
4404    }
4405    if (finishedC)
4406      doingR = true;
4407    if (!doingR) {
4408      int saveSequence = bestSequence;
4409      // Columns
4410      double start = startC[iPassC];
4411      // If we put this idea back then each function needs to update endFraction **
4412#if 0
4413               double dchunk = (static_cast<double> chunk) / (static_cast<double> numberColumns);
4414               double end = CoinMin(startC[iPassC+1], start + dchunk);;
4415#else
4416      double end = startC[iPassC + 1]; // force end
4417#endif
4418      model_->clpMatrix()->partialPricing(model_, start, end, bestSequence, numberWanted);
4419      numberWanted = model_->clpMatrix()->currentWanted();
4420      numberLook -= static_cast< int >((end - start) * numberColumns);
4421      if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
4422        numberWanted = 0; // give up
4423      if (saveSequence != bestSequence) {
4424        // dj
4425        bestDj = fabs(model_->clpMatrix()->reducedCost(model_, bestSequence));
4426      }
4427      if (!numberWanted)
4428        break;
4429      doingR = true;
4430      // update start
4431      startC[iPassC] = end;
4432      if (end >= startC[iPassC + 1] - 1.0e-8) {
4433        if (iPassC)
4434          finishedC = true;
4435        else
4436          iPassC = 2;
4437      }
4438    }
4439  }
4440  updates->setNumElements(0);
4441
4442  // Restore tolerance
4443  model_->setCurrentDualTolerance(saveTolerance);
4444  // Now create variable if column generation
4445  model_->clpMatrix()->createVariable(model_, bestSequence);
4446#ifndef NDEBUG
4447  if (bestSequence >= 0) {
4448    if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
4449      assert(model_->reducedCost(bestSequence) < 0.0);
4450    if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound)
4451      assert(model_->reducedCost(bestSequence) > 0.0);
4452  }
4453#endif
4454  return bestSequence;
4455}
4456
4457/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
4458*/
Note: See TracBrowser for help on using the repository browser.