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

Last change on this file since 2289 was 2289, checked in by forrest, 8 months ago

zero divide and last change to ClpSimplexDual? not such a good idea - no idea why

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