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

Last change on this file since 2373 was 2373, checked in by forrest, 5 months ago

fix bugs people found

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 175.1 KB
Line 
1/* $Id: ClpPrimalColumnSteepest.cpp 2373 2018-11-19 16:24:32Z 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 (!model_->numberIterations())
3215         pivotSequence_=-1;
3216#if ABOCA_LITE
3217       int numberThreads=abcState();
3218       if (numberThreads&&mode_>1)
3219         mode_=0; // force exact
3220#endif
3221          if(weights_) {
3222               // Check if size has changed
3223               if (infeasible_->capacity() == numberRows + numberColumns &&
3224                         alternateWeights_->capacity() == numberRows +
3225                         model_->factorization()->maximumPivots()) {
3226                    //alternateWeights_->clear();
3227                    if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
3228#if ALT_UPDATE_WEIGHTS !=2
3229                         // save pivot order
3230                         CoinMemcpyN(pivotVariable,
3231                                     numberRows, alternateWeights_->getIndices());
3232#endif
3233                         // change from pivot row number to sequence number
3234                         pivotSequence_ = pivotVariable[pivotSequence_];
3235                    } else {
3236                         pivotSequence_ = -1;
3237                    }
3238                    state_ = 1;
3239               } else {
3240                    // size has changed - clear everything
3241                    delete [] weights_;
3242                    weights_ = NULL;
3243                    delete infeasible_;
3244                    infeasible_ = NULL;
3245                    delete alternateWeights_;
3246                    alternateWeights_ = NULL;
3247                    delete [] savedWeights_;
3248                    savedWeights_ = NULL;
3249                    delete [] reference_;
3250                    reference_ = NULL;
3251                    state_ = -1;
3252                    pivotSequence_ = -1;
3253               }
3254          }
3255     } else if (mode == 2 || mode == 4 || mode == 5) {
3256          // restore
3257          if (!weights_ || state_ == -1 || mode == 5) {
3258               // Partial is only allowed with certain types of matrix
3259               if ((mode_ != 4 && mode_ != 5) || numberSwitched_ || !model_->clpMatrix()->canDoPartialPricing()) {
3260                    // initialize weights
3261                    delete [] weights_;
3262                    delete alternateWeights_;
3263                    weights_ = new double[numberRows+numberColumns];
3264                    alternateWeights_ = new CoinIndexedVector();
3265                    // enough space so can use it for factorization
3266                    alternateWeights_->reserve(numberRows +
3267                                               model_->factorization()->maximumPivots());
3268                    initializeWeights();
3269                    // create saved weights
3270                    delete [] savedWeights_;
3271                    savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
3272                    // just do initialization
3273                    mode = 3;
3274               } else {
3275                    // Partial pricing
3276                    // use region as somewhere to save non-fixed slacks
3277                    // set up infeasibilities
3278                    if (!infeasible_) {
3279                         infeasible_ = new CoinIndexedVector();
3280                         infeasible_->reserve(numberColumns + numberRows);
3281                    }
3282                    infeasible_->clear();
3283                    int number = model_->numberRows() + model_->numberColumns();
3284                    int iSequence;
3285                    int numberLook = 0;
3286                    int * which = infeasible_->getIndices();
3287                    for (iSequence = model_->numberColumns(); iSequence < number; iSequence++) {
3288                         ClpSimplex::Status status = model_->getStatus(iSequence);
3289                         if (status != ClpSimplex::isFixed)
3290                              which[numberLook++] = iSequence;
3291                    }
3292                    infeasible_->setNumElements(numberLook);
3293                    doInfeasibilities = false;
3294               }
3295               savedPivotSequence_ = -2;
3296               savedSequenceOut_ = -2;
3297               if (pivotSequence_ < 0 || pivotSequence_ >= 
3298                   numberRows+numberColumns) 
3299                 pivotSequence_ = -1;
3300
3301          } else {
3302               if (mode != 4) {
3303                    // save
3304                    CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
3305                    savedPivotSequence_ = pivotSequence_;
3306                    savedSequenceOut_ = model_->sequenceOut();
3307               } else {
3308                    // restore
3309                    CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
3310                    // was - but I think should not be
3311                    //pivotSequence_= savedPivotSequence_;
3312                    //model_->setSequenceOut(savedSequenceOut_);
3313                    pivotSequence_ = -1;
3314                    model_->setSequenceOut(-1);
3315                    // indices are wrong so clear by hand
3316                    //alternateWeights_->clear();
3317                    CoinZeroN(alternateWeights_->denseVector(),
3318                              alternateWeights_->capacity());
3319                    alternateWeights_->setNumElements(0);
3320               }
3321          }
3322          state_ = 0;
3323          // set up infeasibilities
3324          if (!infeasible_) {
3325               infeasible_ = new CoinIndexedVector();
3326               infeasible_->reserve(numberColumns + numberRows);
3327          }
3328     }
3329     if (mode >= 2 && mode != 5) {
3330          if (mode != 3) {
3331               if (pivotSequence_ >= 0) {
3332#if ALT_UPDATE_WEIGHTS !=2
3333                    // restore pivot row
3334                    int iRow;
3335                    // permute alternateWeights
3336                    double * temp = model_->rowArray(3)->denseVector();;
3337                    double * work = alternateWeights_->denseVector();
3338                    int * savePivotOrder = model_->rowArray(3)->getIndices();
3339                    int * oldPivotOrder = alternateWeights_->getIndices();
3340                    for (iRow = 0; iRow < numberRows; iRow++) {
3341                         int iPivot = oldPivotOrder[iRow];
3342                         temp[iPivot] = work[iRow];
3343                         savePivotOrder[iRow] = iPivot;
3344                    }
3345                    int number = 0;
3346                    int found = -1;
3347                    int * which = oldPivotOrder;
3348                    // find pivot row and re-create alternateWeights
3349                    for (iRow = 0; iRow < numberRows; iRow++) {
3350                         int iPivot = pivotVariable[iRow];
3351                         if (iPivot == pivotSequence_)
3352                              found = iRow;
3353                         work[iRow] = temp[iPivot];
3354                         if (work[iRow])
3355                              which[number++] = iRow;
3356                    }
3357                    alternateWeights_->setNumElements(number);
3358#ifdef CLP_DEBUG
3359                    // Can happen but I should clean up
3360                    assert(found >= 0);
3361#endif
3362                    pivotSequence_ = found;
3363                    for (iRow = 0; iRow < numberRows; iRow++) {
3364                         int iPivot = savePivotOrder[iRow];
3365                         temp[iPivot] = 0.0;
3366                    }
3367#else
3368                    for (int iRow = 0; iRow < numberRows; iRow++) {
3369                         int iPivot = pivotVariable[iRow];
3370                         if (iPivot == pivotSequence_) {
3371                           pivotSequence_ = iRow;
3372                           break;
3373                         }
3374                    }
3375#endif
3376               } else {
3377                    // Just clean up
3378                    /* If this happens when alternateWeights_ is
3379                       in "save" mode then alternateWeights_->clear()
3380                       is disastrous.
3381                       As will be fairly dense anyway and this
3382                       rarely happens just zero out */
3383                    if (alternateWeights_ &&
3384                        alternateWeights_->getNumElements()) {
3385                      //alternateWeights_->clear();
3386                      CoinZeroN(alternateWeights_->denseVector(),
3387                                alternateWeights_->capacity());
3388                      alternateWeights_->setNumElements(0);
3389                    }
3390               }
3391          }
3392          // Save size of factorization
3393          if (!model_->factorization()->pivots())
3394               sizeFactorization_ = model_->factorization()->numberElements();
3395          if(!doInfeasibilities)
3396               return; // don't disturb infeasibilities
3397          double * COIN_RESTRICT infeas = infeasible_->denseVector();
3398          int * COIN_RESTRICT index = infeasible_->getIndices();
3399          int numberNonZero = 0;
3400          infeasibilitiesState_ = 0;
3401          double tolerance = model_->currentDualTolerance();
3402          int number = model_->numberRows() + model_->numberColumns();
3403          int iSequence;
3404
3405          double * reducedCost = model_->djRegion();
3406          const double * lower = model_->lowerRegion();
3407          const double * upper = model_->upperRegion();
3408          const double * solution = model_->solutionRegion();
3409          double primalTolerance = model_->currentPrimalTolerance();
3410
3411          if (!model_->nonLinearCost()->lookBothWays()) {
3412            const unsigned char * COIN_RESTRICT status =
3413              model_->statusArray();
3414#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
3415               for (iSequence = 0; iSequence < number; iSequence++) {
3416                    double value = reducedCost[iSequence];
3417                    infeas[iSequence] = 0.0;
3418                    unsigned char thisStatus=status[iSequence]&7;
3419                    if (thisStatus==3) {
3420                    } else if ((thisStatus&1)!=0) {
3421                      // basic or fixed
3422                      value=0.0;
3423                    } else if (thisStatus==2) {
3424                      value=-value;
3425                    } else {
3426                      // free or superbasic
3427                      if (fabs(value) > FREE_ACCEPT * tolerance) {
3428                        // we are going to bias towards free (but only if reasonable)
3429                        value = -fabs(value)*FREE_BIAS;
3430                      } else {
3431                        value=0.0;
3432                      }
3433                    }
3434                    if (value < -tolerance) {
3435                      infeas[iSequence] = value*value;
3436                      index[numberNonZero++] = iSequence;
3437                    }
3438               }
3439#else
3440               // Columns
3441               int numberColumns = model_->numberColumns();
3442               for (iSequence = 0; iSequence < numberColumns; iSequence++) {
3443                    infeas[iSequence] = 0.0;
3444                    double value = reducedCost[iSequence];
3445                    unsigned char thisStatus=status[iSequence]&7;
3446                    if (thisStatus==3) {
3447                    } else if ((thisStatus&1)!=0) {
3448                      // basic or fixed
3449                      value=0.0;
3450                    } else if (thisStatus==2) {
3451                      value=-value;
3452                    } else {
3453                      // free or superbasic
3454                      if (fabs(value) > FREE_ACCEPT * tolerance) {
3455                        // check hasn't slipped through
3456                        if (solution[iSequence]<lower[iSequence]+primalTolerance) {
3457                          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
3458                        } else if (solution[iSequence]>upper[iSequence]-primalTolerance) {
3459                          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
3460                          value = - value;
3461                        } else {
3462                          // we are going to bias towards free (but only if reasonable)
3463                          value = -fabs(value)*FREE_BIAS;
3464                        }
3465                      } else {
3466                        value=0.0;
3467                      }
3468                    }
3469                    if (value < -tolerance) {
3470                      infeas[iSequence] = value*value;
3471                      index[numberNonZero++] = iSequence;
3472                    }
3473               }
3474               // Rows
3475               for ( ; iSequence < number; iSequence++) {
3476                    double value = reducedCost[iSequence];
3477                    infeas[iSequence]=0.0;
3478                    unsigned char thisStatus=status[iSequence]&7;
3479                    if (thisStatus==3) {
3480                    } else if ((thisStatus&1)!=0) {
3481                      // basic or fixed
3482                      value=0.0;
3483                    } else if (thisStatus==2) {
3484                      value=-value;
3485                    } else {
3486                      // free or superbasic
3487                      if (fabs(value) > FREE_ACCEPT * tolerance) {
3488                        // we are going to bias towards free (but only if reasonable)
3489                        value = -fabs(value)*FREE_BIAS;
3490                      } else {
3491                        value=0.0;
3492                      }
3493                    }
3494                    if (value < -tolerance) {
3495                      infeas[iSequence] = value*value;
3496                      index[numberNonZero++] = iSequence;
3497                    }
3498               }
3499#endif
3500               infeasible_->setNumElements(numberNonZero);
3501          } else {
3502               ClpNonLinearCost * nonLinear = model_->nonLinearCost();
3503               infeasible_->clear();
3504               // can go both ways
3505               for (iSequence = 0; iSequence < number; iSequence++) {
3506                    double value = reducedCost[iSequence];
3507                    ClpSimplex::Status status = model_->getStatus(iSequence);
3508
3509                    switch(status) {
3510
3511                    case ClpSimplex::basic:
3512                    case ClpSimplex::isFixed:
3513                         break;
3514                    case ClpSimplex::isFree:
3515                    case ClpSimplex::superBasic:
3516                         if (fabs(value) > FREE_ACCEPT * tolerance) {
3517                              // we are going to bias towards free (but only if reasonable)
3518                              value *= FREE_BIAS;
3519                              // store square in list
3520                              infeasible_->quickAdd(iSequence, value * value);
3521                         }
3522                         break;
3523                    case ClpSimplex::atUpperBound:
3524                         if (value > tolerance) {
3525                              infeasible_->quickAdd(iSequence, value * value);
3526                         } else {
3527                              // look other way - change up should be negative
3528                              value -= nonLinear->changeUpInCost(iSequence);
3529                              if (value < -tolerance) {
3530                                   // store square in list
3531                                   infeasible_->quickAdd(iSequence, value * value);
3532                              }
3533                         }
3534                         break;
3535                    case ClpSimplex::atLowerBound:
3536                         if (value < -tolerance) {
3537                              infeasible_->quickAdd(iSequence, value * value);
3538                         } else {
3539                              // look other way - change down should be positive
3540                              value -= nonLinear->changeDownInCost(iSequence);
3541                              if (value > tolerance) {
3542                                   // store square in list
3543                                   infeasible_->quickAdd(iSequence, value * value);
3544                              }
3545                         }
3546                    }
3547               }
3548          }
3549     }
3550}
3551// Gets rid of last update
3552void
3553ClpPrimalColumnSteepest::unrollWeights()
3554{
3555     if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3556          return;
3557     double * saved = alternateWeights_->denseVector();
3558     int number = alternateWeights_->getNumElements();
3559     int * which = alternateWeights_->getIndices();
3560     int i;
3561     for (i = 0; i < number; i++) {
3562          int iRow = which[i];
3563          weights_[iRow] = saved[iRow];
3564          saved[iRow] = 0.0;
3565     }
3566     alternateWeights_->setNumElements(0);
3567}
3568
3569//-------------------------------------------------------------------
3570// Clone
3571//-------------------------------------------------------------------
3572ClpPrimalColumnPivot * ClpPrimalColumnSteepest::clone(bool CopyData) const
3573{
3574     if (CopyData) {
3575          return new ClpPrimalColumnSteepest(*this);
3576     } else {
3577          return new ClpPrimalColumnSteepest();
3578     }
3579}
3580void
3581ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
3582{
3583     // Local copy of mode so can decide what to do
3584     int switchType = mode_;
3585     if (mode_ == 4 && numberSwitched_)
3586          switchType = 3;
3587     else if (mode_ == 4 || mode_ == 5)
3588          return;
3589     int number = input->getNumElements();
3590     int * which = input->getIndices();
3591     double * work = input->denseVector();
3592#if ALT_UPDATE_WEIGHTS !=2
3593     int newNumber = 0;
3594     int * newWhich = alternateWeights_->getIndices();
3595     double * newWork = alternateWeights_->denseVector();
3596#endif
3597#ifdef ALT_UPDATE_WEIGHTSz
3598     {
3599       //int newNumber2 = 0;
3600       if (!altVector[0]) {
3601         altVector[0]=new CoinIndexedVector(2000);
3602         altVector[1]=new CoinIndexedVector(2000);
3603         altVector[2]=new CoinIndexedVector(2000);
3604       }
3605       altVector[0]->clear();
3606       // get updated pivot row
3607       int pivotRow = model_->pivotRow();
3608       // should it be - or what
3609       altVector[0]->quickAdd(pivotRow,model_->dualIn());
3610       model_->factorization()->updateColumnTranspose(altVector[2],
3611                                                      altVector[0]);
3612       double * work2 = altVector[0]->denseVector();
3613       //altVector[1]->clear();
3614       int * newWhich2 = altVector[1]->getIndices();
3615       double * newWork2 = altVector[1]->denseVector();
3616       int number2 = altVector[1]->getNumElements();
3617       int nRow=model_->numberRows();
3618       int nCol=model_->numberColumns();
3619       int nTotal=nRow+nCol;
3620       double * temp=new double[2*nTotal+nRow];
3621       memset(temp,0,(2*nTotal+nRow)*sizeof(double));
3622       double * pivRow = temp+nTotal;
3623       double * temp2 = temp+nCol;
3624       double * temp2P=pivRow+nCol;
3625       double * piU=pivRow+nTotal;
3626       double devex=0.0;
3627       double scaleFactor = 1.0/model_->dualIn();
3628       const int * pivotVariable = model_->pivotVariable();
3629       for (int i = 0; i < number; i++) {
3630         int iRow = which[i];
3631         int iPivot = pivotVariable[iRow];
3632         if (reference(iPivot)) {
3633           devex += work[iRow] * work[iRow];
3634         }
3635       }
3636       int sequenceIn = model_->sequenceIn();
3637       int sequenceOut = model_->sequenceOut();
3638       for (int i=0;i<number2;i++) {
3639         int iRow=newWhich2[i];
3640         temp2[iRow] = newWork2[iRow];
3641       }
3642       //if (!input->packedMode()) {
3643       for (int i=0;i<number;i++) {
3644         int iRow=which[i];
3645         piU[iRow] = work2[iRow];
3646         temp2P[iRow]=work2[iRow];
3647       }
3648       double alpha=model_->alpha();
3649       const int * row = model_->matrix()->getIndices();
3650       const CoinBigIndex * columnStart = model_->matrix()->getVectorStarts();
3651       const int * columnLength = model_->matrix()->getVectorLengths();
3652       const double * element = model_->matrix()->getElements();
3653       for (int i=0;i<nCol;i++) {
3654         CoinBigIndex start = columnStart[i];
3655         CoinBigIndex end = start + columnLength[i];
3656         double value=0.0;
3657         double value2=0.0;
3658         for (CoinBigIndex j = start; j < end; j++) {
3659           int iRow = row[j];
3660           value -= piU[iRow] * element[j];
3661           value2 -= newWork2[iRow] * element[j];
3662         }
3663         pivRow[i]=value;
3664         temp[i]=value2;
3665       }
3666       const unsigned char * COIN_RESTRICT statusArray = model_->statusArray();
3667       for (int i=0;i<nTotal;i++) {
3668        unsigned char thisStatus=statusArray[i]&7;
3669#if 0
3670        if (thisStatus==3) {
3671          // lower
3672        } else if ((thisStatus&1)!=0) {
3673          // basic or fixed
3674          //value=0.0;
3675        } else if (thisStatus==2) {
3676          // uppervalue=-value;
3677        } else {
3678          // free or superbasic
3679          //if (fabs(value) > FREE_ACCEPT * -dualTolerance) {
3680            // we are going to bias towards free (but only if reasonable)
3681            //value = -fabs(value)*FREE_BIAS;
3682          //} else {
3683          //value=0.0;
3684          //}
3685        }
3686#else
3687        if (thisStatus!=1&&thisStatus!=5) {
3688          double pivot = pivRow[i]*scaleFactor;
3689          double modification = temp[i];
3690          double thisWeight=weights_[i];
3691          double pivotSquared = pivot*pivot;
3692          double newWeight = thisWeight + pivotSquared * devexA - 2.0*pivot * modification;
3693          temp[i]=newWeight;
3694        } else {
3695          temp[i]=1.0;
3696        }
3697#endif
3698       }
3699       temp[sequenceOut]=devexA/(alpha*alpha);
3700       // to keep clean for debug
3701#ifndef NDEBUG
3702       {
3703         if (sequenceOut<nCol) {
3704           if (model_->columnLower()[sequenceOut]==
3705               model_->columnUpper()[sequenceOut])
3706             temp[sequenceOut]=1.0;
3707         } else {
3708           if (model_->rowLower()[sequenceOut-nCol]==
3709               model_->rowUpper()[sequenceOut-nCol])
3710             temp[sequenceOut]=1.0;
3711         }
3712       }
3713#endif
3714       temp[sequenceIn]=1.0;
3715       for (int i=0;i<nTotal;i++) {
3716         printf("%g ",temp[i]);
3717         if ((i%10)==9)
3718           printf("\n");
3719       }
3720       if (((nTotal-1)%10)!=9)
3721         printf("\n");
3722       delete [] temp;
3723     }
3724#endif
3725     int i;
3726     int sequenceIn = model_->sequenceIn();
3727     int sequenceOut = model_->sequenceOut();
3728     const int * pivotVariable = model_->pivotVariable();
3729
3730     int pivotRow = model_->pivotRow();
3731     pivotSequence_ = pivotRow;
3732
3733     devex_ = 0.0;
3734     // Can't create alternateWeights_ as packed as needed unpacked
3735     if (!input->packedMode()) {
3736          if (pivotRow >= 0) {
3737               if (switchType == 1) {
3738                    for (i = 0; i < number; i++) {
3739                         int iRow = which[i];
3740                         devex_ += work[iRow] * work[iRow];
3741#if ALT_UPDATE_WEIGHTS !=2
3742                         newWork[iRow] = -2.0 * work[iRow];
3743#endif
3744                    }
3745#if ALT_UPDATE_WEIGHTS !=2
3746                    newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3747#endif
3748                    devex_ += ADD_ONE;
3749                    weights_[sequenceOut] = 1.0 + ADD_ONE;
3750#if ALT_UPDATE_WEIGHTS !=2
3751                    CoinMemcpyN(which, number, newWhich);
3752                    alternateWeights_->setNumElements(number);
3753#endif
3754               } else {
3755                    if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3756                         for (i = 0; i < number; i++) {
3757                              int iRow = which[i];
3758                              int iPivot = pivotVariable[iRow];
3759                              if (reference(iPivot)) {
3760                                   devex_ += work[iRow] * work[iRow];
3761#if ALT_UPDATE_WEIGHTS !=2
3762                                   newWork[iRow] = -2.0 * work[iRow];
3763                                   newWhich[newNumber++] = iRow;
3764#endif
3765                              }
3766                         }
3767#if ALT_UPDATE_WEIGHTS !=2
3768                         if (!newWork[pivotRow] && devex_ > 0.0)
3769                              newWhich[newNumber++] = pivotRow; // add if not already in
3770                         newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3771#endif
3772                    } else {
3773                         for (i = 0; i < number; i++) {
3774                              int iRow = which[i];
3775                              int iPivot = pivotVariable[iRow];
3776                              if (reference(iPivot))
3777                                   devex_ += work[iRow] * work[iRow];
3778                         }
3779                    }
3780                    if (reference(sequenceIn)) {
3781                         devex_ += 1.0;
3782                    } else {
3783                    }
3784                    if (reference(sequenceOut)) {
3785                         weights_[sequenceOut] = 1.0 + 1.0;
3786                    } else {
3787                         weights_[sequenceOut] = 1.0;
3788                    }
3789#if ALT_UPDATE_WEIGHTS !=2
3790                    alternateWeights_->setNumElements(newNumber);
3791#endif
3792               }
3793          } else {
3794               if (switchType == 1) {
3795                    for (i = 0; i < number; i++) {
3796                         int iRow = which[i];
3797                         devex_ += work[iRow] * work[iRow];
3798                    }
3799                    devex_ += ADD_ONE;
3800               } else {
3801                    for (i = 0; i < number; i++) {
3802                         int iRow = which[i];
3803                         int iPivot = pivotVariable[iRow];
3804                         if (reference(iPivot)) {
3805                              devex_ += work[iRow] * work[iRow];
3806                         }
3807                    }
3808                    if (reference(sequenceIn))
3809                         devex_ += 1.0;
3810               }
3811          }
3812     } else {
3813          // packed input
3814          if (pivotRow >= 0) {
3815               if (switchType == 1) {
3816                    for (i = 0; i < number; i++) {
3817#if ALT_UPDATE_WEIGHTS !=2
3818                         int iRow = which[i];
3819#endif
3820                         devex_ += work[i] * work[i];
3821#if ALT_UPDATE_WEIGHTS !=2
3822                         newWork[iRow] = -2.0 * work[i];
3823#endif
3824                    }
3825#if ALT_UPDATE_WEIGHTS !=2
3826                    newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3827#endif
3828                    devex_ += ADD_ONE;
3829                    weights_[sequenceOut] = 1.0 + ADD_ONE;
3830#if ALT_UPDATE_WEIGHTS !=2
3831                    CoinMemcpyN(which, number, newWhich);
3832                    alternateWeights_->setNumElements(number);
3833#endif
3834               } else {
3835                    if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
3836                         for (i = 0; i < number; i++) {
3837                              int iRow = which[i];
3838                              int iPivot = pivotVariable[iRow];
3839                              if (reference(iPivot)) {
3840                                   devex_ += work[i] * work[i];
3841#if ALT_UPDATE_WEIGHTS !=2
3842                                   newWork[iRow] = -2.0 * work[i];
3843                                   newWhich[newNumber++] = iRow;
3844#endif
3845                              }
3846                         }
3847#if ALT_UPDATE_WEIGHTS !=2
3848                         if (!newWork[pivotRow] && devex_ > 0.0)
3849                              newWhich[newNumber++] = pivotRow; // add if not already in
3850                         newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
3851#endif
3852                    } else {
3853                         for (i = 0; i < number; i++) {
3854                              int iRow = which[i];
3855                              int iPivot = pivotVariable[iRow];
3856                              if (reference(iPivot))
3857                                   devex_ += work[i] * work[i];
3858                         }
3859                    }
3860                    if (reference(sequenceIn)) {
3861                         devex_ += 1.0;
3862                    } else {
3863                    }
3864                    if (reference(sequenceOut)) {
3865                         weights_[sequenceOut] = 1.0 + 1.0;
3866                    } else {
3867                         weights_[sequenceOut] = 1.0;
3868                    }
3869#if ALT_UPDATE_WEIGHTS !=2
3870                    alternateWeights_->setNumElements(newNumber);
3871#endif
3872               }
3873          } else {
3874               if (switchType == 1) {
3875                    for (i = 0; i < number; i++) {
3876                         devex_ += work[i] * work[i];
3877                    }
3878                    devex_ += ADD_ONE;
3879               } else {
3880                    for (i = 0; i < number; i++) {
3881                         int iRow = which[i];
3882                         int iPivot = pivotVariable[iRow];
3883                         if (reference(iPivot)) {
3884                              devex_ += work[i] * work[i];
3885                         }
3886                    }
3887                    if (reference(sequenceIn))
3888                         devex_ += 1.0;
3889               }
3890          }
3891     }
3892     if (devex_<1.001e-30) {
3893       COIN_DETAIL_PRINT(printf("devex of incoming tiny %d %g\n",sequenceIn,devex_));
3894       devex_=1.0e-30;
3895     }
3896     double oldDevex = weights_[sequenceIn];
3897#ifdef CLP_DEBUG
3898     if ((model_->messageHandler()->logLevel() & 32))
3899          printf("old weight %g, new %g\n", oldDevex, devex_);
3900#endif
3901     double check = CoinMax(devex_, oldDevex) + 0.1;
3902     weights_[sequenceIn] = devex_;
3903     double testValue = 0.1;
3904     if (mode_ == 4 && numberSwitched_ == 1)
3905          testValue = 0.5;
3906     if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
3907#ifdef CLP_DEBUG
3908          if ((model_->messageHandler()->logLevel() & 48) == 16)
3909               printf("old weight %g, new %g\n", oldDevex, devex_);
3910#endif
3911          //printf("old weight %g, new %g\n",oldDevex,devex_);
3912          testValue = 0.99;
3913          if (mode_ == 1)
3914               testValue = 1.01e1; // make unlikely to do if steepest
3915          else if (mode_ == 4 && numberSwitched_ == 1)
3916               testValue = 0.9;
3917          double difference = fabs(devex_ - oldDevex);
3918          if ( difference > testValue * check ) {
3919               // need to redo
3920               model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
3921                                                 *model_->messagesPointer())
3922                         << oldDevex << devex_
3923                         << CoinMessageEol;
3924               initializeWeights();
3925               // redo devex_
3926               if (pivotRow >= 0) 
3927                 devex_=1.0;
3928          }
3929     }
3930     if (pivotRow >= 0) {
3931          // set outgoing weight here
3932       double alpha = model_->alpha();
3933       if (fabs(alpha)>1.0e15) {
3934         COIN_DETAIL_PRINT(printf("alpha %g for %d !!\n",alpha,model_->sequenceOut()));
3935         alpha = 1.0e15;
3936       }
3937       weights_[model_->sequenceOut()] = devex_ / (alpha*alpha);
3938     }
3939}
3940// Checks accuracy - just for debug
3941void
3942ClpPrimalColumnSteepest::checkAccuracy(int sequence,
3943                                       double relativeTolerance,
3944                                       CoinIndexedVector * rowArray1,
3945                                       CoinIndexedVector * rowArray2)
3946{
3947     if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
3948          return;
3949     model_->unpack(rowArray1, sequence);
3950     model_->factorization()->updateColumn(rowArray2, rowArray1);
3951     int number = rowArray1->getNumElements();
3952     int * which = rowArray1->getIndices();
3953     double * work = rowArray1->denseVector();
3954     const int * pivotVariable = model_->pivotVariable();
3955
3956     double devex = 0.0;
3957     int i;
3958
3959     if (mode_ == 1) {
3960          for (i = 0; i < number; i++) {
3961               int iRow = which[i];
3962               devex += work[iRow] * work[iRow];
3963               work[iRow] = 0.0;
3964          }
3965          devex += ADD_ONE;
3966     } else {
3967          for (i = 0; i < number; i++) {
3968               int iRow = which[i];
3969               int iPivot = pivotVariable[iRow];
3970               if (reference(iPivot)) {
3971                    devex += work[iRow] * work[iRow];
3972               }
3973               work[iRow] = 0.0;
3974          }
3975          if (reference(sequence))
3976               devex += 1.0;
3977     }
3978
3979     double oldDevex = CoinMax(weights_[sequence],1.0e-4);
3980     devex = CoinMax(devex,1.0e-4);
3981     double check = CoinMax(devex, oldDevex);;
3982     rowArray1->setNumElements(0);
3983     if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
3984       //COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
3985       printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex);
3986       if (mode_==0) {
3987         rowArray1->setNumElements(0);
3988         model_->unpack(rowArray1, sequence);
3989         number = rowArray1->getNumElements();
3990         for (i = 0; i< number;i++)
3991           printf("(%d,%g) ",which[i],work[which[i]]);
3992         printf("\n");
3993         model_->factorization()->updateColumn(rowArray2, rowArray1);
3994         number = rowArray1->getNumElements();
3995         for (i = 0; i< number;i++)
3996           printf("(%d,%g) ",which[i],work[which[i]]);
3997         printf("\n");
3998         devex=0.0;
3999         for (i = 0; i < number; i++) {
4000           int iRow = which[i];
4001           int iPivot = pivotVariable[iRow];
4002           if (reference(iPivot)) {
4003             devex += work[iRow] * work[iRow];
4004           }
4005           work[iRow] = 0.0;
4006         }
4007         if (reference(sequence))
4008           devex += 1.0;
4009       }
4010          // update so won't print again
4011          weights_[sequence] = devex;
4012     }
4013}
4014
4015// Initialize weights
4016void
4017ClpPrimalColumnSteepest::initializeWeights()
4018{
4019     int numberRows = model_->numberRows();
4020     int numberColumns = model_->numberColumns();
4021     int number = numberRows + numberColumns;
4022     int iSequence;
4023     if (mode_ != 1) {
4024          // initialize to 1.0
4025          // and set reference framework
4026          if (!reference_) {
4027               int nWords = (number + 31) >> 5;
4028               reference_ = new unsigned int[nWords];
4029               CoinZeroN(reference_, nWords);
4030          }
4031
4032          for (iSequence = 0; iSequence < number; iSequence++) {
4033               weights_[iSequence] = 1.0;
4034               if (model_->getStatus(iSequence) == ClpSimplex::basic) {
4035                    setReference(iSequence, false);
4036               } else {
4037                    setReference(iSequence, true);
4038               }
4039          }
4040     } else {
4041          CoinIndexedVector * temp = new CoinIndexedVector();
4042          temp->reserve(numberRows +
4043                        model_->factorization()->maximumPivots());
4044          double * array = alternateWeights_->denseVector();
4045          int * which = alternateWeights_->getIndices();
4046
4047          for (iSequence = 0; iSequence < number; iSequence++) {
4048               weights_[iSequence] = 1.0 + ADD_ONE;
4049               if (model_->getStatus(iSequence) != ClpSimplex::basic &&
4050                         model_->getStatus(iSequence) != ClpSimplex::isFixed) {
4051                    model_->unpack(alternateWeights_, iSequence);
4052                    double value = ADD_ONE;
4053                    model_->factorization()->updateColumn(temp, alternateWeights_);
4054                    int number = alternateWeights_->getNumElements();
4055                    int j;
4056                    for (j = 0; j < number; j++) {
4057                         int iRow = which[j];
4058                         value += array[iRow] * array[iRow];
4059                         array[iRow] = 0.0;
4060                    }
4061                    alternateWeights_->setNumElements(0);
4062                    weights_[iSequence] = value;
4063               }
4064          }
4065          delete temp;
4066     }
4067}
4068// Gets rid of all arrays
4069void
4070ClpPrimalColumnSteepest::clearArrays()
4071{
4072     if (persistence_ == normal) {
4073          delete [] weights_;
4074          weights_ = NULL;
4075          delete infeasible_;
4076          infeasible_ = NULL;
4077          delete alternateWeights_;
4078          alternateWeights_ = NULL;
4079          delete [] savedWeights_;
4080          savedWeights_ = NULL;
4081          delete [] reference_;
4082          reference_ = NULL;
4083     }
4084     pivotSequence_ = -1;
4085     state_ = -1;
4086     savedPivotSequence_ = -1;
4087     savedSequenceOut_ = -1;
4088     devex_ = 0.0;
4089}
4090// Returns true if would not find any column
4091bool
4092ClpPrimalColumnSteepest::looksOptimal() const
4093{
4094     if (looksOptimal_)
4095          return true; // user overrode
4096     //**** THIS MUST MATCH the action coding above
4097     double tolerance = model_->currentDualTolerance();
4098     // we can't really trust infeasibilities if there is dual error
4099     // this coding has to mimic coding in checkDualSolution
4100     double error = CoinMin(1.0e-2, model_->largestDualError());
4101     // allow tolerance at least slightly bigger than standard
4102     tolerance = tolerance +  error;
4103     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
4104          // we can't really trust infeasibilities if there is dual error
4105          double checkTolerance = 1.0e-8;
4106          if (!model_->factorization()->pivots())
4107               checkTolerance = 1.0e-6;
4108          if (model_->largestDualError() > checkTolerance)
4109               tolerance *= model_->largestDualError() / checkTolerance;
4110          // But cap
4111          tolerance = CoinMin(1000.0, tolerance);
4112     }
4113     int number = model_->numberRows() + model_->numberColumns();
4114     int iSequence;
4115
4116     double * reducedCost = model_->djRegion();
4117     int numberInfeasible = 0;
4118     if (!model_->nonLinearCost()->lookBothWays()) {
4119          for (iSequence = 0; iSequence < number; iSequence++) {
4120               double value = reducedCost[iSequence];
4121               ClpSimplex::Status status = model_->getStatus(iSequence);
4122
4123               switch(status) {
4124
4125               case ClpSimplex::basic:
4126               case ClpSimplex::isFixed:
4127                    break;
4128               case ClpSimplex::isFree:
4129               case ClpSimplex::superBasic:
4130                    if (fabs(value) > FREE_ACCEPT * tolerance)
4131                         numberInfeasible++;
4132                    break;
4133               case ClpSimplex::atUpperBound:
4134                    if (value > tolerance)
4135                         numberInfeasible++;
4136                    break;
4137               case ClpSimplex::atLowerBound:
4138                    if (value < -tolerance)
4139                         numberInfeasible++;
4140               }
4141          }
4142     } else {
4143          ClpNonLinearCost * nonLinear = model_->nonLinearCost();
4144          // can go both ways
4145          for (iSequence = 0; iSequence < number; iSequence++) {
4146               double value = reducedCost[iSequence];
4147               ClpSimplex::Status status = model_->getStatus(iSequence);
4148
4149               switch(status) {
4150
4151               case ClpSimplex::basic:
4152               case ClpSimplex::isFixed:
4153                    break;
4154               case ClpSimplex::isFree:
4155               case ClpSimplex::superBasic:
4156                    if (fabs(value) > FREE_ACCEPT * tolerance)
4157                         numberInfeasible++;
4158                    break;
4159               case ClpSimplex::atUpperBound:
4160                    if (value > tolerance) {
4161                         numberInfeasible++;
4162                    } else {
4163                         // look other way - change up should be negative
4164                         value -= nonLinear->changeUpInCost(iSequence);
4165                         if (value < -tolerance)
4166                              numberInfeasible++;
4167                    }
4168                    break;
4169               case ClpSimplex::atLowerBound:
4170                    if (value < -tolerance) {
4171                         numberInfeasible++;
4172                    } else {
4173                         // look other way - change down should be positive
4174                         value -= nonLinear->changeDownInCost(iSequence);
4175                         if (value > tolerance)
4176                              numberInfeasible++;
4177                    }
4178               }
4179          }
4180     }
4181     return numberInfeasible == 0;
4182}
4183/* Returns number of extra columns for sprint algorithm - 0 means off.
4184   Also number of iterations before recompute
4185*/
4186int
4187ClpPrimalColumnSteepest::numberSprintColumns(int & numberIterations) const
4188{
4189     numberIterations = 0;
4190     int numberAdd = 0;
4191     if (!numberSwitched_ && mode_ >= 10) {
4192          numberIterations = CoinMin(2000, model_->numberRows() / 5);
4193          numberIterations = CoinMax(numberIterations, model_->factorizationFrequency());
4194          numberIterations = CoinMax(numberIterations, 500);
4195          if (mode_ == 10) {
4196               numberAdd = CoinMax(300, model_->numberColumns() / 10);
4197               numberAdd = CoinMax(numberAdd, model_->numberRows() / 5);
4198               // fake all
4199               //numberAdd=1000000;
4200               numberAdd = CoinMin(numberAdd, model_->numberColumns());
4201          } else {
4202               abort();
4203          }
4204     }
4205     return numberAdd;
4206}
4207// Switch off sprint idea
4208void
4209ClpPrimalColumnSteepest::switchOffSprint()
4210{
4211     numberSwitched_ = 10;
4212}
4213// Update djs doing partial pricing (dantzig)
4214int
4215ClpPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
4216                                        CoinIndexedVector * spareRow2,
4217                                        int numberWanted,
4218                                        int numberLook)
4219{
4220     int number = 0;
4221     int * index;
4222     double * updateBy;
4223     double * reducedCost;
4224     double saveTolerance = model_->currentDualTolerance();
4225     double tolerance = model_->currentDualTolerance();
4226     // we can't really trust infeasibilities if there is dual error
4227     // this coding has to mimic coding in checkDualSolution
4228     double error = CoinMin(1.0e-2, model_->largestDualError());
4229     // allow tolerance at least slightly bigger than standard
4230     tolerance = tolerance +  error;
4231     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
4232          // we can't really trust infeasibilities if there is dual error
4233          double checkTolerance = 1.0e-8;
4234          if (!model_->factorization()->pivots())
4235               checkTolerance = 1.0e-6;
4236          if (model_->largestDualError() > checkTolerance)
4237               tolerance *= model_->largestDualError() / checkTolerance;
4238          // But cap
4239          tolerance = CoinMin(1000.0, tolerance);
4240     }
4241     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
4242          tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
4243     // So partial pricing can use
4244     model_->setCurrentDualTolerance(tolerance);
4245     model_->factorization()->updateColumnTranspose(spareRow2, updates);
4246     int numberColumns = model_->numberColumns();
4247
4248     // Rows
4249     reducedCost = model_->djRegion(0);
4250
4251     number = updates->getNumElements();
4252     index = updates->getIndices();
4253     updateBy = updates->denseVector();
4254     int j;
4255     double * duals = model_->dualRowSolution();
4256     for (j = 0; j < number; j++) {
4257          int iSequence = index[j];
4258          double value = duals[iSequence];
4259          value -= updateBy[j];
4260          updateBy[j] = 0.0;
4261          duals[iSequence] = value;
4262     }
4263     //#define CLP_DEBUG
4264#ifdef CLP_DEBUG
4265     // check duals
4266     {
4267          int numberRows = model_->numberRows();
4268          //work space
4269          CoinIndexedVector arrayVector;
4270          arrayVector.reserve(numberRows + 1000);
4271          CoinIndexedVector workSpace;
4272          workSpace.reserve(numberRows + 1000);
4273
4274
4275          int iRow;
4276          double * array = arrayVector.denseVector();
4277          int * index = arrayVector.getIndices();
4278          int number = 0;
4279          int * pivotVariable = model_->pivotVariable();
4280          double * cost = model_->costRegion();
4281          for (iRow = 0; iRow < numberRows; iRow++) {
4282               int iPivot = pivotVariable[iRow];
4283               double value = cost[iPivot];
4284               if (value) {
4285                    array[iRow] = value;
4286                    index[number++] = iRow;
4287               }
4288          }
4289          arrayVector.setNumElements(number);
4290          // Extended duals before "updateTranspose"
4291          model_->clpMatrix()->dualExpanded(model_, &arrayVector, NULL, 0);
4292
4293          // Btran basic costs
4294          model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
4295
4296          // now look at dual solution
4297          for (iRow = 0; iRow < numberRows; iRow++) {
4298               // slack
4299               double value = array[iRow];
4300               if (fabs(duals[iRow] - value) > 1.0e-3)
4301                    printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
4302               //duals[iRow]=value;
4303          }
4304     }
4305#endif
4306#undef CLP_DEBUG
4307     double bestDj = tolerance;
4308     int bestSequence = -1;
4309
4310     const double * cost = model_->costRegion(1);
4311
4312     model_->clpMatrix()->setOriginalWanted(numberWanted);
4313     model_->clpMatrix()->setCurrentWanted(numberWanted);
4314     int iPassR = 0, iPassC = 0;
4315     // Setup two passes
4316     // This biases towards picking row variables
4317     // This probably should be fixed
4318     int startR[4];
4319     const int * which = infeasible_->getIndices();
4320     int nSlacks = infeasible_->getNumElements();
4321     startR[1] = nSlacks;
4322     startR[2] = 0;
4323     double randomR = model_->randomNumberGenerator()->randomDouble();
4324     double dstart = static_cast<double> (nSlacks) * randomR;
4325     startR[0] = static_cast<int> (dstart);
4326     startR[3] = startR[0];
4327     double startC[4];
4328     startC[1] = 1.0;
4329     startC[2] = 0;
4330     double randomC = model_->randomNumberGenerator()->randomDouble();
4331     startC[0] = randomC;
4332     startC[3] = randomC;
4333     reducedCost = model_->djRegion(1);
4334     int sequenceOut = model_->sequenceOut();
4335     double * duals2 = duals - numberColumns;
4336     int chunk = CoinMin(1024, (numberColumns + nSlacks) / 32);
4337#ifdef COIN_DETAIL
4338     if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
4339          printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
4340          int i;
4341          for (i = 0; i < 4; i++)
4342               printf("start R %d C %g ", startR[i], startC[i]);
4343          printf("\n");
4344     }
4345#endif
4346     chunk = CoinMax(chunk, 256);
4347     bool finishedR = false, finishedC = false;
4348     bool doingR = randomR > randomC;
4349     //doingR=false;
4350     int saveNumberWanted = numberWanted;
4351     while (!finishedR || !finishedC) {
4352          if (finishedR)
4353               doingR = false;
4354          if (doingR) {
4355               int saveSequence = bestSequence;
4356               int start = startR[iPassR];
4357               int end = CoinMin(startR[iPassR+1], start + chunk / 2);
4358               int jSequence;
4359               for (jSequence = start; jSequence < end; jSequence++) {
4360                    int iSequence = which[jSequence];
4361                    if (iSequence != sequenceOut) {
4362                         double value;
4363                         ClpSimplex::Status status = model_->getStatus(iSequence);
4364
4365                         switch(status) {
4366
4367                         case ClpSimplex::basic:
4368                         case ClpSimplex::isFixed:
4369                              break;
4370                         case ClpSimplex::isFree:
4371                         case ClpSimplex::superBasic:
4372                              value = fabs(cost[iSequence] + duals2[iSequence]);
4373                              if (value > FREE_ACCEPT * tolerance) {
4374                                   numberWanted--;
4375                                   // we are going to bias towards free (but only if reasonable)
4376                                   value *= FREE_BIAS;
4377                                   if (value > bestDj) {
4378                                        // check flagged variable and correct dj
4379                                        if (!model_->flagged(iSequence)) {
4380                                             bestDj = value;
4381                                             bestSequence = iSequence;
4382                                        } else {
4383                                             // just to make sure we don't exit before got something
4384                                             numberWanted++;
4385                                        }
4386                                   }
4387                              }
4388                              break;
4389                         case ClpSimplex::atUpperBound:
4390                              value = cost[iSequence] + duals2[iSequence];
4391                              if (value > tolerance) {
4392                                   numberWanted--;
4393                                   if (value > bestDj) {
4394                                        // check flagged variable and correct dj
4395                                        if (!model_->flagged(iSequence)) {
4396                                             bestDj = value;
4397                                             bestSequence = iSequence;
4398                                        } else {
4399                                             // just to make sure we don't exit before got something
4400                                             numberWanted++;
4401                                        }
4402                                   }
4403                              }
4404                              break;
4405                         case ClpSimplex::atLowerBound:
4406                              value = -(cost[iSequence] + duals2[iSequence]);
4407                              if (value > tolerance) {
4408                                   numberWanted--;
4409                                   if (value > bestDj) {
4410                                        // check flagged variable and correct dj
4411                                        if (!model_->flagged(iSequence)) {
4412                                             bestDj = value;
4413                                             bestSequence = iSequence;
4414                                        } else {
4415                                             // just to make sure we don't exit before got something
4416                                             numberWanted++;
4417                                        }
4418                                   }
4419                              }
4420                              break;
4421                         }
4422                    }
4423                    if (!numberWanted)
4424                         break;
4425               }
4426               numberLook -= (end - start);
4427               if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
4428                    numberWanted = 0; // give up
4429               if (saveSequence != bestSequence) {
4430                    // dj
4431                    reducedCost[bestSequence] = cost[bestSequence] + duals[bestSequence-numberColumns];
4432                    bestDj = fabs(reducedCost[bestSequence]);
4433                    model_->clpMatrix()->setSavedBestSequence(bestSequence);
4434                    model_->clpMatrix()->setSavedBestDj(reducedCost[bestSequence]);
4435               }
4436               model_->clpMatrix()->setCurrentWanted(numberWanted);
4437               if (!numberWanted)
4438                    break;
4439               doingR = false;
4440               // update start
4441               startR[iPassR] = jSequence;
4442               if (jSequence >= startR[iPassR+1]) {
4443                    if (iPassR)
4444                         finishedR = true;
4445                    else
4446                         iPassR = 2;
4447               }
4448          }
4449          if (finishedC)
4450               doingR = true;
4451          if (!doingR) {
4452               int saveSequence = bestSequence;
4453               // Columns
4454               double start = startC[iPassC];
4455               // If we put this idea back then each function needs to update endFraction **
4456#if 0
4457               double dchunk = (static_cast<double> chunk) / (static_cast<double> numberColumns);
4458               double end = CoinMin(startC[iPassC+1], start + dchunk);;
4459#else
4460               double end = startC[iPassC+1]; // force end
4461#endif
4462               model_->clpMatrix()->partialPricing(model_, start, end, bestSequence, numberWanted);
4463               numberWanted = model_->clpMatrix()->currentWanted();
4464               numberLook -= static_cast<int> ((end - start) * numberColumns);
4465               if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
4466                    numberWanted = 0; // give up
4467               if (saveSequence != bestSequence) {
4468                    // dj
4469                    bestDj = fabs(model_->clpMatrix()->reducedCost(model_, bestSequence));
4470               }
4471               if (!numberWanted)
4472                    break;
4473               doingR = true;
4474               // update start
4475               startC[iPassC] = end;
4476               if (end >= startC[iPassC+1] - 1.0e-8) {
4477                    if (iPassC)
4478                         finishedC = true;
4479                    else
4480                         iPassC = 2;
4481               }
4482          }
4483     }
4484     updates->setNumElements(0);
4485
4486     // Restore tolerance
4487     model_->setCurrentDualTolerance(saveTolerance);
4488     // Now create variable if column generation
4489     model_->clpMatrix()->createVariable(model_, bestSequence);
4490#ifndef NDEBUG
4491     if (bestSequence >= 0) {
4492          if (model_->getStatus(bestSequence) == ClpSimplex::atLowerBound)
4493               assert(model_->reducedCost(bestSequence) < 0.0);
4494          if (model_->getStatus(bestSequence) == ClpSimplex::atUpperBound)
4495               assert(model_->reducedCost(bestSequence) > 0.0);
4496     }
4497#endif
4498     return bestSequence;
4499}
Note: See TracBrowser for help on using the repository browser.