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

Last change on this file since 2235 was 2235, checked in by forrest, 4 years ago

stuff for vector matrix

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