source: trunk/Clp/src/AbcPrimalColumnSteepest.cpp @ 1878

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

File size: 84.4 KB
Line 
1/* $Id: AbcPrimalColumnSteepest.cpp 1732 2011-05-31 08:09:41Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#if CLP_HAS_ABC
7#include "CoinPragma.hpp"
8
9#include "AbcSimplex.hpp"
10#include "AbcPrimalColumnSteepest.hpp"
11#include "CoinIndexedVector.hpp"
12#include "AbcSimplexFactorization.hpp"
13#include "AbcNonLinearCost.hpp"
14#include "ClpMessage.hpp"
15#include "CoinHelperFunctions.hpp"
16#undef COIN_DETAIL_PRINT
17#define COIN_DETAIL_PRINT(s) s
18#include <stdio.h>
19//#define CLP_DEBUG
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27AbcPrimalColumnSteepest::AbcPrimalColumnSteepest (int mode)
28     : AbcPrimalColumnPivot(),
29       devex_(0.0),
30       weights_(NULL),
31       infeasible_(NULL),
32       alternateWeights_(NULL),
33       savedWeights_(NULL),
34       reference_(NULL),
35       state_(-1),
36       mode_(mode),
37       persistence_(normal),
38       numberSwitched_(0),
39       pivotSequence_(-1),
40       savedPivotSequence_(-1),
41       savedSequenceOut_(-1),
42       sizeFactorization_(0)
43{
44     type_ = 2 + 64 * mode;
45}
46//-------------------------------------------------------------------
47// Copy constructor
48//-------------------------------------------------------------------
49AbcPrimalColumnSteepest::AbcPrimalColumnSteepest (const AbcPrimalColumnSteepest & rhs)
50     : AbcPrimalColumnPivot(rhs)
51{
52     state_ = rhs.state_;
53     mode_ = rhs.mode_;
54     persistence_ = rhs.persistence_;
55     numberSwitched_ = rhs.numberSwitched_;
56     model_ = rhs.model_;
57     pivotSequence_ = rhs.pivotSequence_;
58     savedPivotSequence_ = rhs.savedPivotSequence_;
59     savedSequenceOut_ = rhs.savedSequenceOut_;
60     sizeFactorization_ = rhs.sizeFactorization_;
61     devex_ = rhs.devex_;
62     if ((model_ && model_->whatsChanged() & 1) != 0) {
63          if (rhs.infeasible_) {
64               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
65          } else {
66               infeasible_ = NULL;
67          }
68          reference_ = NULL;
69          if (rhs.weights_) {
70               assert(model_);
71               int number = model_->numberRows() + model_->numberColumns();
72               assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
73               weights_ = new double[number];
74               CoinMemcpyN(rhs.weights_, number, weights_);
75               savedWeights_ = new double[number];
76               CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
77               if (mode_ != 1) {
78                    reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
79               }
80          } else {
81               weights_ = NULL;
82               savedWeights_ = NULL;
83          }
84          if (rhs.alternateWeights_) {
85               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
86          } else {
87               alternateWeights_ = NULL;
88          }
89     } else {
90          infeasible_ = NULL;
91          reference_ = NULL;
92          weights_ = NULL;
93          savedWeights_ = NULL;
94          alternateWeights_ = NULL;
95     }
96}
97
98//-------------------------------------------------------------------
99// Destructor
100//-------------------------------------------------------------------
101AbcPrimalColumnSteepest::~AbcPrimalColumnSteepest ()
102{
103     delete [] weights_;
104     delete infeasible_;
105     delete alternateWeights_;
106     delete [] savedWeights_;
107     delete [] reference_;
108}
109
110//----------------------------------------------------------------
111// Assignment operator
112//-------------------------------------------------------------------
113AbcPrimalColumnSteepest &
114AbcPrimalColumnSteepest::operator=(const AbcPrimalColumnSteepest& rhs)
115{
116     if (this != &rhs) {
117          AbcPrimalColumnPivot::operator=(rhs);
118          state_ = rhs.state_;
119          mode_ = rhs.mode_;
120          persistence_ = rhs.persistence_;
121          numberSwitched_ = rhs.numberSwitched_;
122          model_ = rhs.model_;
123          pivotSequence_ = rhs.pivotSequence_;
124          savedPivotSequence_ = rhs.savedPivotSequence_;
125          savedSequenceOut_ = rhs.savedSequenceOut_;
126          sizeFactorization_ = rhs.sizeFactorization_;
127          devex_ = rhs.devex_;
128          delete [] weights_;
129          delete [] reference_;
130          reference_ = NULL;
131          delete infeasible_;
132          delete alternateWeights_;
133          delete [] savedWeights_;
134          savedWeights_ = NULL;
135          if (rhs.infeasible_ != NULL) {
136               infeasible_ = new CoinIndexedVector(rhs.infeasible_);
137          } else {
138               infeasible_ = NULL;
139          }
140          if (rhs.weights_ != NULL) {
141               assert(model_);
142               int number = model_->numberRows() + model_->numberColumns();
143               assert (number == rhs.model_->numberRows() + rhs.model_->numberColumns());
144               weights_ = new double[number];
145               CoinMemcpyN(rhs.weights_, number, weights_);
146               savedWeights_ = new double[number];
147               CoinMemcpyN(rhs.savedWeights_, number, savedWeights_);
148               if (mode_ != 1) {
149                    reference_ = CoinCopyOfArray(rhs.reference_, (number + 31) >> 5);
150               }
151          } else {
152               weights_ = NULL;
153          }
154          if (rhs.alternateWeights_ != NULL) {
155               alternateWeights_ = new CoinIndexedVector(rhs.alternateWeights_);
156          } else {
157               alternateWeights_ = NULL;
158          }
159     }
160     return *this;
161}
162// These have to match ClpPackedMatrix version
163#define TRY_NORM 1.0e-4
164#define ADD_ONE 1.0
165// Returns pivot column, -1 if none
166/*      The Packed CoinIndexedVector updates has cost updates - for normal LP
167        that is just +-weight where a feasibility changed.  It also has
168        reduced cost from last iteration in pivot row*/
169int
170AbcPrimalColumnSteepest::pivotColumn(CoinPartitionedVector * updates,
171                                     CoinPartitionedVector * spareRow2,
172                                     CoinPartitionedVector * spareColumn1)
173{
174     assert(model_);
175     int number = 0;
176     int * index;
177     double tolerance = model_->currentDualTolerance();
178     // we can't really trust infeasibilities if there is dual error
179     // this coding has to mimic coding in checkDualSolution
180     double error = CoinMin(1.0e-2, model_->largestDualError());
181     // allow tolerance at least slightly bigger than standard
182     tolerance = tolerance +  error;
183     int pivotRow = model_->pivotRow();
184     int anyUpdates;
185     double * infeas = infeasible_->denseVector();
186
187     // Local copy of mode so can decide what to do
188     int switchType;
189     if (mode_ == 4)
190          switchType = 5 - numberSwitched_;
191     else if (mode_ >= 10)
192          switchType = 3;
193     else
194          switchType = mode_;
195     /* switchType -
196        0 - all exact devex
197        1 - all steepest
198        2 - some exact devex
199        3 - auto some exact devex
200        4 - devex
201        5 - dantzig
202        10 - can go to mini-sprint
203     */
204     if (updates->getNumElements() > 1) {
205          // would have to have two goes for devex, three for steepest
206          anyUpdates = 2;
207     } else if (updates->getNumElements()) {
208          if (updates->getIndices()[0] == pivotRow
209              && fabs(updates->denseVector()[pivotRow]) > 1.0e-6) {
210            //&& fabs(updates->denseVector()[pivotRow]) < 1.0e6) {
211               // reasonable size
212               anyUpdates = 1;
213               //if (fabs(model_->dualIn())<1.0e-4||fabs(fabs(model_->dualIn())-fabs(updates->denseVector()[0]))>1.0e-5)
214               //printf("dualin %g pivot %g\n",model_->dualIn(),updates->denseVector()[0]);
215          } else {
216               // too small
217               anyUpdates = 2;
218          }
219     } else if (pivotSequence_ >= 0) {
220          // just after re-factorization
221          anyUpdates = -1;
222     } else {
223          // sub flip - nothing to do
224          anyUpdates = 0;
225     }
226     int sequenceOut = model_->sequenceOut();
227     if (switchType == 5) {
228               pivotSequence_ = -1;
229               pivotRow = -1;
230               // See if to switch
231               int numberRows = model_->numberRows();
232               int numberWanted = 10;
233               int numberColumns = model_->numberColumns();
234               double ratio = static_cast<double> (sizeFactorization_) /
235                              static_cast<double> (numberRows);
236               // Number of dual infeasibilities at last invert
237               int numberDual = model_->numberDualInfeasibilities();
238               int numberLook = CoinMin(numberDual, numberColumns / 10);
239               if (ratio < 1.0) {
240                    numberWanted = 100;
241                    numberLook /= 20;
242                    numberWanted = CoinMax(numberWanted, numberLook);
243               } else if (ratio < 3.0) {
244                    numberWanted = 500;
245                    numberLook /= 15;
246                    numberWanted = CoinMax(numberWanted, numberLook);
247               } else if (ratio < 4.0 || mode_ == 5) {
248                    numberWanted = 1000;
249                    numberLook /= 10;
250                    numberWanted = CoinMax(numberWanted, numberLook);
251               } else if (mode_ != 5) {
252                    switchType = 4;
253                    // initialize
254                    numberSwitched_++;
255                    // Make sure will re-do
256                    delete [] weights_;
257                    weights_ = NULL;
258                    //work space
259                    int whichArray[2];
260                    for (int i=0;i<2;i++)
261                      whichArray[i]=model_->getAvailableArray();
262                    CoinIndexedVector * array1 = model_->usefulArray(whichArray[0]);
263                    CoinIndexedVector * array2 = model_->usefulArray(whichArray[1]);
264                    model_->computeDuals(NULL,array1,array2);
265                    for (int i=0;i<2;i++)
266                      model_->setAvailableArray(whichArray[i]);
267                    saveWeights(model_, 4);
268                    anyUpdates = 0;
269                    COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
270               }
271               if (switchType == 5) {
272                    numberLook *= 5; // needs tuning for gub
273                    if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
274                         COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d look %d\n",
275                                                  sizeFactorization_, ratio, numberWanted, numberLook));
276                    }
277                    // Update duals and row djs
278                    // Do partial pricing
279                    return partialPricing(updates, numberWanted, numberLook);
280               }
281     }
282     if (switchType == 5) {
283          if (anyUpdates > 0) {
284               justDjs(updates,spareColumn1);
285          }
286     } else if (anyUpdates == 1) {
287          if (switchType < 4) {
288               // exact etc when can use dj
289            doSteepestWork(updates,spareRow2,spareColumn1,2);
290          } else {
291               // devex etc when can use dj
292               djsAndDevex(updates, spareRow2,
293                           spareColumn1);
294          }
295     } else if (anyUpdates == -1) {
296          if (switchType < 4) {
297               // exact etc when djs okay
298            //if (model_->maximumIterations()!=1000)
299            doSteepestWork(updates,spareRow2,spareColumn1,1);
300          } else {
301               // devex etc when djs okay
302               justDevex(updates, 
303                         spareColumn1);
304          }
305     } else if (anyUpdates == 2) {
306          if (switchType < 4) {
307               // exact etc when have to use pivot
308            doSteepestWork(updates,spareRow2,spareColumn1,3);
309          } else {
310               // devex etc when have to use pivot
311               djsAndDevex2(updates,
312                            spareColumn1);
313          }
314     }
315#ifdef CLP_DEBUG
316     alternateWeights_->checkClear();
317#endif
318     // make sure outgoing from last iteration okay
319     if (sequenceOut >= 0) {
320          AbcSimplex::Status status = model_->getInternalStatus(sequenceOut);
321          double value = model_->reducedCost(sequenceOut);
322
323          switch(status) {
324
325          case AbcSimplex::basic:
326          case AbcSimplex::isFixed:
327               break;
328          case AbcSimplex::isFree:
329          case AbcSimplex::superBasic:
330               if (fabs(value) > FREE_ACCEPT * tolerance) {
331                    // we are going to bias towards free (but only if reasonable)
332                    value *= FREE_BIAS;
333                    // store square in list
334                    if (infeas[sequenceOut])
335                         infeas[sequenceOut] = value * value; // already there
336                    else
337                         infeasible_->quickAdd(sequenceOut, value * value);
338               } else {
339                    infeasible_->zero(sequenceOut);
340               }
341               break;
342          case AbcSimplex::atUpperBound:
343               if (value > tolerance) {
344                    // store square in list
345                    if (infeas[sequenceOut])
346                         infeas[sequenceOut] = value * value; // already there
347                    else
348                         infeasible_->quickAdd(sequenceOut, value * value);
349               } else {
350                    infeasible_->zero(sequenceOut);
351               }
352               break;
353          case AbcSimplex::atLowerBound:
354               if (value < -tolerance) {
355                    // store square in list
356                    if (infeas[sequenceOut])
357                         infeas[sequenceOut] = value * value; // already there
358                    else
359                         infeasible_->quickAdd(sequenceOut, value * value);
360               } else {
361                    infeasible_->zero(sequenceOut);
362               }
363          }
364     }
365     // update of duals finished - now do pricing
366     // See what sort of pricing
367     int numberWanted = 10;
368     number = infeasible_->getNumElements();
369     int numberColumns = model_->numberColumns();
370     if (switchType == 5) {
371          pivotSequence_ = -1;
372          pivotRow = -1;
373          // See if to switch
374          int numberRows = model_->numberRows();
375          // ratio is done on number of columns here
376          //double ratio = static_cast<double> sizeFactorization_/static_cast<double> numberColumns;
377          double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
378          if (ratio < 1.0) {
379               numberWanted = CoinMax(100, number / 200);
380          } else if (ratio < 2.0 - 1.0) {
381               numberWanted = CoinMax(500, number / 40);
382          } else if (ratio < 4.0 - 3.0 || mode_ == 5) {
383               numberWanted = CoinMax(2000, number / 10);
384               numberWanted = CoinMax(numberWanted, numberColumns / 30);
385          } else if (mode_ != 5) {
386               switchType = 4;
387               // initialize
388               numberSwitched_++;
389               // Make sure will re-do
390               delete [] weights_;
391               weights_ = NULL;
392               saveWeights(model_, 4);
393               COIN_DETAIL_PRINT(printf("switching to devex %d nel ratio %g\n", sizeFactorization_, ratio));
394          }
395          //if (model_->numberIterations()%1000==0)
396          //printf("numels %d ratio %g wanted %d\n",sizeFactorization_,ratio,numberWanted);
397     }
398     int numberRows = model_->numberRows();
399     // ratio is done on number of rows here
400     double ratio = static_cast<double> (sizeFactorization_) / static_cast<double> (numberRows);
401     if(switchType == 4) {
402          // Still in devex mode
403          // Go to steepest if lot of iterations?
404          if (ratio < 5.0) {
405               numberWanted = CoinMax(2000, number / 10);
406               numberWanted = CoinMax(numberWanted, numberColumns / 20);
407          } else if (ratio < 7.0) {
408               numberWanted = CoinMax(2000, number / 5);
409               numberWanted = CoinMax(numberWanted, numberColumns / 10);
410          } else {
411               // we can zero out
412               updates->clear();
413               spareColumn1->clear();
414               switchType = 3;
415               // initialize
416               pivotSequence_ = -1;
417               pivotRow = -1;
418               numberSwitched_++;
419               // Make sure will re-do
420               delete [] weights_;
421               weights_ = NULL;
422               saveWeights(model_, 4);
423               COIN_DETAIL_PRINT(printf("switching to exact %d nel ratio %g\n", sizeFactorization_, ratio));
424               updates->clear();
425          }
426          if (model_->numberIterations() % 1000 == 0)
427            COIN_DETAIL_PRINT(printf("numels %d ratio %g wanted %d type x\n", sizeFactorization_, ratio, numberWanted));
428     }
429     if (switchType < 4) {
430          if (switchType < 2 ) {
431               numberWanted = number + 1;
432          } else if (switchType == 2) {
433               numberWanted = CoinMax(2000, number / 8);
434          } else {
435               if (ratio < 1.0) {
436                    numberWanted = CoinMax(2000, number / 20);
437               } else if (ratio < 5.0) {
438                    numberWanted = CoinMax(2000, number / 10);
439                    numberWanted = CoinMax(numberWanted, numberColumns / 40);
440               } else if (ratio < 10.0) {
441                    numberWanted = CoinMax(2000, number / 8);
442                    numberWanted = CoinMax(numberWanted, numberColumns / 20);
443               } else {
444                    ratio = number * (ratio / 80.0);
445                    if (ratio > number) {
446                         numberWanted = number + 1;
447                    } else {
448                         numberWanted = CoinMax(2000, static_cast<int> (ratio));
449                         numberWanted = CoinMax(numberWanted, numberColumns / 10);
450                    }
451               }
452          }
453          //if (model_->numberIterations()%1000==0)
454          //printf("numels %d ratio %g wanted %d type %d\n",sizeFactorization_,ratio,numberWanted,
455          //switchType);
456     }
457
458
459     double bestDj = 1.0e-30;
460     int bestSequence = -1;
461
462     int i, iSequence;
463     index = infeasible_->getIndices();
464     number = infeasible_->getNumElements();
465     if(model_->numberIterations() < model_->lastBadIteration() + 200 &&
466               model_->factorization()->pivots() > 10) {
467          // we can't really trust infeasibilities if there is dual error
468          double checkTolerance = 1.0e-8;
469          if (model_->largestDualError() > checkTolerance)
470               tolerance *= model_->largestDualError() / checkTolerance;
471          // But cap
472          tolerance = CoinMin(1000.0, tolerance);
473     }
474#ifdef CLP_DEBUG
475     if (model_->numberDualInfeasibilities() == 1)
476          printf("** %g %g %g %x %x %d\n", tolerance, model_->dualTolerance(),
477                 model_->largestDualError(), model_, model_->messageHandler(),
478                 number);
479#endif
480     // stop last one coming immediately
481     double saveOutInfeasibility = 0.0;
482     if (sequenceOut >= 0) {
483          saveOutInfeasibility = infeas[sequenceOut];
484          infeas[sequenceOut] = 0.0;
485     }
486     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
487          tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
488     tolerance *= tolerance; // as we are using squares
489
490     int iPass;
491     // Setup two passes
492     int start[4];
493     start[1] = number;
494     start[2] = 0;
495     double dstart = static_cast<double> (number) * model_->randomNumberGenerator()->randomDouble();
496     start[0] = static_cast<int> (dstart);
497     start[3] = start[0];
498     //double largestWeight=0.0;
499     //double smallestWeight=1.0e100;
500     for (iPass = 0; iPass < 2; iPass++) {
501          int end = start[2*iPass+1];
502          if (switchType < 5) {
503               for (i = start[2*iPass]; i < end; i++) {
504                    iSequence = index[i];
505                    double value = infeas[iSequence];
506                    double weight = weights_[iSequence];
507                    if (value > tolerance) {
508                         //weight=1.0;
509                         if (value > bestDj * weight) {
510                              // check flagged variable and correct dj
511                              if (!model_->flagged(iSequence)) {
512                                   bestDj = value / weight;
513                                   bestSequence = iSequence;
514                              } else {
515                                   // just to make sure we don't exit before got something
516                                   numberWanted++;
517                              }
518                         }
519                         numberWanted--;
520                    }
521                    if (!numberWanted)
522                         break;
523               }
524          } else {
525               // Dantzig
526               for (i = start[2*iPass]; i < end; i++) {
527                    iSequence = index[i];
528                    double value = infeas[iSequence];
529                    if (value > tolerance) {
530                         if (value > bestDj) {
531                              // check flagged variable and correct dj
532                              if (!model_->flagged(iSequence)) {
533                                   bestDj = value;
534                                   bestSequence = iSequence;
535                              } else {
536                                   // just to make sure we don't exit before got something
537                                   numberWanted++;
538                              }
539                         }
540                         numberWanted--;
541                    }
542                    if (!numberWanted)
543                         break;
544               }
545          }
546          if (!numberWanted)
547               break;
548     }
549     if (sequenceOut >= 0) {
550          infeas[sequenceOut] = saveOutInfeasibility;
551     }
552     /*if (model_->numberIterations()%100==0)
553       printf("%d best %g\n",bestSequence,bestDj);*/
554
555#ifndef NDEBUG
556     if (bestSequence >= 0) {
557          if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
558               assert(model_->reducedCost(bestSequence) < 0.0);
559          if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound) {
560               assert(model_->reducedCost(bestSequence) > 0.0);
561          }
562     }
563#endif
564#if 1
565     if (model_->logLevel()==127) {
566       double * reducedCost=model_->djRegion();
567       for (int i=0;i<numberRows;i++)
568         printf("row %d weight %g infeas %g dj %g\n",i,weights_[i],infeas[i]>1.0e-13 ? infeas[i]: 0.0,reducedCost[i]);
569       for (int i=numberRows;i<numberRows+numberColumns;i++)
570         printf("column %d weight %g infeas %g dj %g\n",i-numberRows,weights_[i],infeas[i]>1.0e-13 ? infeas[i]: 0.0,reducedCost[i]);
571     }
572#endif
573#if SOME_DEBUG_1
574     if (switchType<4) {
575       // check for accuracy
576       int iCheck = 892;
577       //printf("weight for iCheck is %g\n",weights_[iCheck]);
578       int numberRows = model_->numberRows();
579       int numberColumns = model_->numberColumns();
580       for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
581         if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
582             model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
583           checkAccuracy(iCheck, 1.0e-1, updates);
584       }
585     }
586#endif
587#if 1
588     for (int iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
589       if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
590           model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
591         assert (fabs(model_->costRegion()[iCheck])<1.0e9);
592     }
593#endif
594     return bestSequence;
595}
596/* Does steepest work
597   type -
598   0 - just djs
599   1 - just steepest
600   2 - both using scaleFactor
601   3 - both using extra array
602*/
603int 
604AbcPrimalColumnSteepest::doSteepestWork(CoinPartitionedVector * updates,
605                                        CoinPartitionedVector * spareRow2,
606                                        CoinPartitionedVector * spareColumn1,
607                                        int type)
608{
609  int pivotRow = model_->pivotRow(); 
610  if (type==1) {
611    pivotRow=pivotSequence_;
612    if (pivotRow<0)
613      type=-1;
614  } else if (pivotSequence_<0) {
615    assert (model_->sequenceIn()==model_->sequenceOut());
616    type=0;
617  }
618  // see if reference
619  int sequenceIn = pivotRow>=0 ? model_->sequenceIn() : -1;
620  // save outgoing weight round update
621  double outgoingWeight = 0.0;
622  int sequenceOut = model_->sequenceOut();
623  if (sequenceOut >= 0)
624    outgoingWeight = weights_[sequenceOut];
625  double referenceIn;
626  if (mode_ != 1) {
627    if(sequenceIn>=0&&reference(sequenceIn))
628      referenceIn = 1.0;
629    else
630      referenceIn = 0.0;
631  } else {
632    referenceIn = -1.0;
633  }
634  double * infeasibilities=infeasible_->denseVector();
635  if (sequenceIn>=0)
636    infeasibilities[sequenceIn]=0.0;
637  double * array=updates->denseVector();
638  double scaleFactor;
639  double devex=devex_;
640  if (type==0) {
641    // just djs - to keep clean swap
642    scaleFactor=1.0;
643    CoinPartitionedVector * temp = updates;
644    updates=spareRow2;
645    spareRow2=temp;
646    assert (!alternateWeights_->getNumElements());
647    devex=0.0;
648  } else if (type==1) {
649    // just steepest
650    updates->clear();
651    scaleFactor=COIN_DBL_MAX;
652    // might as well set dj to 1
653    double dj = -1.0;
654    assert (pivotRow>=0);
655    spareRow2->createUnpacked(1, &pivotRow, &dj);
656  } else if (type==2) {
657    // using scaleFactor - swap
658    assert (pivotRow>=0);
659    scaleFactor=-array[pivotRow];
660    array[pivotRow]=-1.0;
661    CoinPartitionedVector * temp = updates;
662    updates=spareRow2;
663    spareRow2=temp;
664  } else if (type==3) {
665    // need two arrays
666    scaleFactor=0.0;
667    // might as well set dj to 1
668    double dj = -1.0;
669    assert (pivotRow>=0);
670    spareRow2->createUnpacked(1, &pivotRow, &dj);
671  }
672  if (type>=0) {
673    // parallelize this
674    if (type==0) {
675      model_->factorization()->updateColumnTranspose(*spareRow2);
676    } else if (type<3) {
677      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2,0);
678      model_->factorization()->updateColumnTransposeCpu(*alternateWeights_,1);
679      cilk_sync;
680    } else {
681      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*updates,0);
682      cilk_spawn model_->factorization()->updateColumnTransposeCpu(*spareRow2,1);
683      model_->factorization()->updateColumnTransposeCpu(*alternateWeights_,2);
684      cilk_sync;
685    }
686    model_->abcMatrix()->primalColumnDouble(*spareRow2,
687                                            *updates,
688                                            *alternateWeights_,
689                                            *spareColumn1,
690                                            *infeasible_,referenceIn,devex,
691                                            reference_,weights_,scaleFactor);
692  }
693  pivotSequence_=-1;
694  // later do pricing here
695  // later move pricing into abcmatrix
696  // restore outgoing weight
697  if (sequenceOut >= 0)
698    weights_[sequenceOut] = outgoingWeight;
699  alternateWeights_->clear();
700  updates->clear();
701  spareRow2->clear();
702  return 0;
703}
704// Just update djs
705void
706AbcPrimalColumnSteepest::justDjs(CoinIndexedVector * updates,
707                                 CoinIndexedVector * spareColumn1)
708{
709     int iSection, j;
710     int number = 0;
711     double multiplier;
712     int * index;
713     double * updateBy;
714     double * reducedCost;
715     double tolerance = model_->currentDualTolerance();
716     // we can't really trust infeasibilities if there is dual error
717     // this coding has to mimic coding in checkDualSolution
718     double error = CoinMin(1.0e-2, model_->largestDualError());
719     // allow tolerance at least slightly bigger than standard
720     tolerance = tolerance +  error;
721     int pivotRow = model_->pivotRow();
722     double * infeas = infeasible_->denseVector();
723     //updates->scanAndPack();
724     model_->factorization()->updateColumnTranspose(*updates);
725
726     // put row of tableau in rowArray and columnArray
727     model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
728     // normal
729     reducedCost = model_->djRegion();
730     for (iSection = 0; iSection < 2; iSection++) {
731
732          int addSequence;
733#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
734          double slack_multiplier;
735#endif
736
737          if (!iSection) {
738               number = updates->getNumElements();
739               index = updates->getIndices();
740               updateBy = updates->denseVector();
741               addSequence = 0;
742               multiplier=-1;
743#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
744               slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
745#endif
746          } else {
747               number = spareColumn1->getNumElements();
748               index = spareColumn1->getIndices();
749               updateBy = spareColumn1->denseVector();
750               addSequence = model_->maximumAbcNumberRows();
751               multiplier=1;
752#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
753               slack_multiplier = 1.0;
754#endif
755          }
756
757          for (j = 0; j < number; j++) {
758               int iSequence = index[j];
759               double tableauValue=updateBy[iSequence];
760               updateBy[iSequence] = 0.0;
761               iSequence += addSequence;
762               double value = reducedCost[iSequence];
763               value -= multiplier*tableauValue;
764               reducedCost[iSequence] = value;
765               AbcSimplex::Status status = model_->getInternalStatus(iSequence);
766
767               switch(status) {
768
769               case AbcSimplex::basic:
770                    infeasible_->zero(iSequence);
771               case AbcSimplex::isFixed:
772                    break;
773               case AbcSimplex::isFree:
774               case AbcSimplex::superBasic:
775                    if (fabs(value) > FREE_ACCEPT * tolerance) {
776                         // we are going to bias towards free (but only if reasonable)
777                         value *= FREE_BIAS;
778                         // store square in list
779                         if (infeas[iSequence])
780                              infeas[iSequence] = value * value; // already there
781                         else
782                              infeasible_->quickAdd(iSequence , value * value);
783                    } else {
784                         infeasible_->zero(iSequence);
785                    }
786                    break;
787               case AbcSimplex::atUpperBound:
788                    if (value > tolerance) {
789#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
790                        value *= value*slack_multiplier;
791#else
792                        value *= value;
793#endif
794                         // store square in list
795                         if (infeas[iSequence])
796                              infeas[iSequence] = value; // already there
797                         else
798                              infeasible_->quickAdd(iSequence, value);
799                    } else {
800                         infeasible_->zero(iSequence);
801                    }
802                    break;
803               case AbcSimplex::atLowerBound:
804                    if (value < -tolerance) {
805#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
806                        value *= value*slack_multiplier;
807#else
808                        value *= value;
809#endif
810                         // store square in list
811                         if (infeas[iSequence])
812                              infeas[iSequence] = value; // already there
813                         else
814                              infeasible_->quickAdd(iSequence, value);
815                    } else {
816                         infeasible_->zero(iSequence);
817                    }
818               }
819          }
820     }
821     updates->setNumElements(0);
822     spareColumn1->setNumElements(0);
823     if (pivotRow >= 0) {
824          // make sure infeasibility on incoming is 0.0
825          int sequenceIn = model_->sequenceIn();
826          infeasible_->zero(sequenceIn);
827     }
828}
829// Update djs, weights for Devex
830void
831AbcPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates,
832                                     CoinIndexedVector * spareRow2,
833                                     CoinIndexedVector * spareColumn1)
834{
835     int j;
836     int number = 0;
837     int * index;
838     double * updateBy;
839     double * reducedCost;
840     double tolerance = model_->currentDualTolerance();
841     // we can't really trust infeasibilities if there is dual error
842     // this coding has to mimic coding in checkDualSolution
843     double error = CoinMin(1.0e-2, model_->largestDualError());
844     // allow tolerance at least slightly bigger than standard
845     tolerance = tolerance +  error;
846     // for weights update we use pivotSequence
847     // unset in case sub flip
848     assert (pivotSequence_ >= 0);
849     assert (model_->pivotVariable()[pivotSequence_] == model_->sequenceIn());
850     double scaleFactor = 1.0 / updates->denseVector()[pivotSequence_]; // as formula is with 1.0
851     pivotSequence_ = -1;
852     double * infeas = infeasible_->denseVector();
853     //updates->scanAndPack();
854     model_->factorization()->updateColumnTranspose(*updates);
855     // and we can see if reference
856     //double referenceIn = 0.0;
857     int sequenceIn = model_->sequenceIn();
858     //if (mode_ != 1 && reference(sequenceIn))
859     //   referenceIn = 1.0;
860     // save outgoing weight round update
861     double outgoingWeight = 0.0;
862     int sequenceOut = model_->sequenceOut();
863     if (sequenceOut >= 0)
864          outgoingWeight = weights_[sequenceOut];
865
866     // put row of tableau in rowArray and columnArray
867     model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
868     // update weights
869     double * weight;
870     // rows
871     reducedCost = model_->djRegion();
872
873     number = updates->getNumElements();
874     index = updates->getIndices();
875     updateBy = updates->denseVector();
876     weight = weights_;
877     // Devex
878     for (j = 0; j < number; j++) {
879          double thisWeight;
880          double pivot;
881          double value3;
882          int iSequence = index[j];
883          double value = reducedCost[iSequence];
884          double value2 = updateBy[iSequence];
885          updateBy[iSequence] = 0.0;
886          value -= -value2;
887          reducedCost[iSequence] = value;
888          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
889
890          switch(status) {
891
892          case AbcSimplex::basic:
893               infeasible_->zero(iSequence);
894          case AbcSimplex::isFixed:
895               break;
896          case AbcSimplex::isFree:
897          case AbcSimplex::superBasic:
898               thisWeight = weight[iSequence];
899               // row has -1
900               pivot = value2 * scaleFactor;
901               value3 = pivot * pivot * devex_;
902               if (reference(iSequence))
903                    value3 += 1.0;
904               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
905               if (fabs(value) > FREE_ACCEPT * tolerance) {
906                    // we are going to bias towards free (but only if reasonable)
907                    value *= FREE_BIAS;
908                    // store square in list
909                    if (infeas[iSequence])
910                         infeas[iSequence] = value * value; // already there
911                    else
912                         infeasible_->quickAdd(iSequence , value * value);
913               } else {
914                    infeasible_->zero(iSequence);
915               }
916               break;
917          case AbcSimplex::atUpperBound:
918               thisWeight = weight[iSequence];
919               // row has -1
920               pivot = value2 * scaleFactor;
921               value3 = pivot * pivot * devex_;
922               if (reference(iSequence))
923                    value3 += 1.0;
924               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
925               if (value > tolerance) {
926                    // store square in list
927#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
928                    value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
929#else
930                    value *= value;
931#endif
932                    if (infeas[iSequence])
933                         infeas[iSequence] = value; // already there
934                    else
935                         infeasible_->quickAdd(iSequence , value);
936               } else {
937                    infeasible_->zero(iSequence);
938               }
939               break;
940          case AbcSimplex::atLowerBound:
941               thisWeight = weight[iSequence];
942               // row has -1
943               pivot = value2 * scaleFactor;
944               value3 = pivot * pivot * devex_;
945               if (reference(iSequence))
946                    value3 += 1.0;
947               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
948               if (value < -tolerance) {
949                    // store square in list
950#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
951                    value *= value*CLP_PRIMAL_SLACK_MULTIPLIER;
952#else
953                    value *= value;
954#endif
955                    if (infeas[iSequence])
956                         infeas[iSequence] = value; // already there
957                    else
958                         infeasible_->quickAdd(iSequence , value);
959               } else {
960                    infeasible_->zero(iSequence);
961               }
962          }
963     }
964
965     // columns
966     int addSequence = model_->maximumAbcNumberRows();
967
968     scaleFactor = -scaleFactor;
969     number = spareColumn1->getNumElements();
970     index = spareColumn1->getIndices();
971     updateBy = spareColumn1->denseVector();
972
973     // Devex
974
975     for (j = 0; j < number; j++) {
976          double thisWeight;
977          double pivot;
978          double value3;
979          int iSequence = index[j];
980          double value2=updateBy[iSequence];
981          updateBy[iSequence] = 0.0;
982          iSequence += addSequence;
983          double value = reducedCost[iSequence];
984          value -= value2;
985          reducedCost[iSequence] = value;
986          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
987
988          switch(status) {
989
990          case AbcSimplex::basic:
991               infeasible_->zero(iSequence);
992          case AbcSimplex::isFixed:
993               break;
994          case AbcSimplex::isFree:
995          case AbcSimplex::superBasic:
996               thisWeight = weight[iSequence];
997               // row has -1
998               pivot = value2 * scaleFactor;
999               value3 = pivot * pivot * devex_;
1000               if (reference(iSequence))
1001                    value3 += 1.0;
1002               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1003               if (fabs(value) > FREE_ACCEPT * tolerance) {
1004                    // we are going to bias towards free (but only if reasonable)
1005                    value *= FREE_BIAS;
1006                    // store square in list
1007                    if (infeas[iSequence])
1008                         infeas[iSequence] = value * value; // already there
1009                    else
1010                         infeasible_->quickAdd(iSequence, value * value);
1011               } else {
1012                    infeasible_->zero(iSequence);
1013               }
1014               break;
1015          case AbcSimplex::atUpperBound:
1016               thisWeight = weight[iSequence];
1017               // row has -1
1018               pivot = value2 * scaleFactor;
1019               value3 = pivot * pivot * devex_;
1020               if (reference(iSequence))
1021                    value3 += 1.0;
1022               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1023               if (value > tolerance) {
1024                    // store square in list
1025                    if (infeas[iSequence])
1026                         infeas[iSequence] = value * value; // already there
1027                    else
1028                         infeasible_->quickAdd(iSequence, value * value);
1029               } else {
1030                    infeasible_->zero(iSequence);
1031               }
1032               break;
1033          case AbcSimplex::atLowerBound:
1034               thisWeight = weight[iSequence];
1035               // row has -1
1036               pivot = value2 * scaleFactor;
1037               value3 = pivot * pivot * devex_;
1038               if (reference(iSequence))
1039                    value3 += 1.0;
1040               weight[iSequence] = CoinMax(0.99 * thisWeight, value3);
1041               if (value < -tolerance) {
1042                    // store square in list
1043                    if (infeas[iSequence])
1044                         infeas[iSequence] = value * value; // already there
1045                    else
1046                         infeasible_->quickAdd(iSequence, value * value);
1047               } else {
1048                    infeasible_->zero(iSequence);
1049               }
1050          }
1051     }
1052     // restore outgoing weight
1053     if (sequenceOut >= 0)
1054          weights_[sequenceOut] = outgoingWeight;
1055     // make sure infeasibility on incoming is 0.0
1056     infeasible_->zero(sequenceIn);
1057     spareRow2->setNumElements(0);
1058     //#define SOME_DEBUG_1
1059#ifdef SOME_DEBUG_1
1060     // check for accuracy
1061     int iCheck = 892;
1062     //printf("weight for iCheck is %g\n",weights_[iCheck]);
1063     int numberRows = model_->numberRows();
1064     //int numberColumns = model_->numberColumns();
1065     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1066          if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
1067                    !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
1068               checkAccuracy(iCheck, 1.0e-1, updates);
1069     }
1070#endif
1071     updates->setNumElements(0);
1072     spareColumn1->setNumElements(0);
1073}
1074// Update djs, weights for Devex
1075void
1076AbcPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates,
1077                                      CoinIndexedVector * spareColumn1)
1078{
1079     int iSection, j;
1080     int number = 0;
1081     double multiplier;
1082     int * index;
1083     double * updateBy;
1084     double * reducedCost;
1085     // dj could be very small (or even zero - take care)
1086     double dj = model_->dualIn();
1087     double tolerance = model_->currentDualTolerance();
1088     // we can't really trust infeasibilities if there is dual error
1089     // this coding has to mimic coding in checkDualSolution
1090     double error = CoinMin(1.0e-2, model_->largestDualError());
1091     // allow tolerance at least slightly bigger than standard
1092     tolerance = tolerance +  error;
1093     int pivotRow = model_->pivotRow();
1094     double * infeas = infeasible_->denseVector();
1095     //updates->scanAndPack();
1096     model_->factorization()->updateColumnTranspose(*updates);
1097
1098     // put row of tableau in rowArray and columnArray
1099     model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
1100     // normal
1101     reducedCost = model_->djRegion();
1102     for (iSection = 0; iSection < 2; iSection++) {
1103
1104          int addSequence;
1105#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1106          double slack_multiplier;
1107#endif
1108
1109          if (!iSection) {
1110               number = updates->getNumElements();
1111               index = updates->getIndices();
1112               updateBy = updates->denseVector();
1113               addSequence = 0;
1114               multiplier=-1;
1115#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1116               slack_multiplier = CLP_PRIMAL_SLACK_MULTIPLIER;
1117#endif
1118          } else {
1119               number = spareColumn1->getNumElements();
1120               index = spareColumn1->getIndices();
1121               updateBy = spareColumn1->denseVector();
1122               addSequence = model_->maximumAbcNumberRows();
1123               multiplier=1;
1124#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1125               slack_multiplier = 1.0;
1126#endif
1127          }
1128
1129          for (j = 0; j < number; j++) {
1130               int iSequence = index[j];
1131               double tableauValue=updateBy[iSequence];
1132               updateBy[iSequence] = 0.0;
1133               iSequence += addSequence;
1134               double value = reducedCost[iSequence];
1135               value -= multiplier*tableauValue;
1136               AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1137
1138               switch(status) {
1139
1140               case AbcSimplex::basic:
1141                    infeasible_->zero(iSequence);
1142               case AbcSimplex::isFixed:
1143                    break;
1144               case AbcSimplex::isFree:
1145               case AbcSimplex::superBasic:
1146                 reducedCost[iSequence] = value;
1147                    if (fabs(value) > FREE_ACCEPT * tolerance) {
1148                         // we are going to bias towards free (but only if reasonable)
1149                         value *= FREE_BIAS;
1150                         // store square in list
1151                         if (infeas[iSequence])
1152                              infeas[iSequence] = value * value; // already there
1153                         else
1154                              infeasible_->quickAdd(iSequence , value * value);
1155                    } else {
1156                         infeasible_->zero(iSequence );
1157                    }
1158                    break;
1159               case AbcSimplex::atUpperBound:
1160                 reducedCost[iSequence] = value;
1161                    if (value > tolerance) {
1162#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1163                        value *= value*slack_multiplier;
1164#else
1165                        value *= value;
1166#endif
1167                         // store square in list
1168                         if (infeas[iSequence])
1169                              infeas[iSequence] = value; // already there
1170                         else
1171                              infeasible_->quickAdd(iSequence, value);
1172                    } else {
1173                         infeasible_->zero(iSequence);
1174                    }
1175                    break;
1176               case AbcSimplex::atLowerBound:
1177                 reducedCost[iSequence] = value;
1178                    if (value < -tolerance) {
1179#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
1180                        value *= value*slack_multiplier;
1181#else
1182                        value *= value;
1183#endif
1184                         // store square in list
1185                         if (infeas[iSequence])
1186                              infeas[iSequence] = value; // already there
1187                         else
1188                              infeasible_->quickAdd(iSequence, value);
1189                    } else {
1190                         infeasible_->zero(iSequence);
1191                    }
1192               }
1193          }
1194     }
1195     // They are empty
1196     updates->setNumElements(0);
1197     spareColumn1->setNumElements(0);
1198     // make sure infeasibility on incoming is 0.0
1199     int sequenceIn = model_->sequenceIn();
1200     infeasible_->zero(sequenceIn);
1201     // for weights update we use pivotSequence
1202     if (pivotSequence_ >= 0) {
1203          pivotRow = pivotSequence_;
1204          // unset in case sub flip
1205          pivotSequence_ = -1;
1206          // make sure infeasibility on incoming is 0.0
1207          const int * pivotVariable = model_->pivotVariable();
1208          sequenceIn = pivotVariable[pivotRow];
1209          infeasible_->zero(sequenceIn);
1210          // and we can see if reference
1211          //double referenceIn = 0.0;
1212          //if (mode_ != 1 && reference(sequenceIn))
1213          //   referenceIn = 1.0;
1214          // save outgoing weight round update
1215          double outgoingWeight = 0.0;
1216          int sequenceOut = model_->sequenceOut();
1217          if (sequenceOut >= 0)
1218               outgoingWeight = weights_[sequenceOut];
1219          // update weights
1220          updates->setNumElements(0);
1221          spareColumn1->setNumElements(0);
1222          // might as well set dj to 1
1223          dj = 1.0;
1224          updates->insert(pivotRow, -dj);
1225          model_->factorization()->updateColumnTranspose(*updates);
1226          // put row of tableau in rowArray and columnArray
1227          model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
1228          double * weight;
1229          int numberColumns = model_->numberColumns();
1230          // rows
1231          number = updates->getNumElements();
1232          index = updates->getIndices();
1233          updateBy = updates->denseVector();
1234          weight = weights_;
1235
1236          assert (devex_ > 0.0);
1237          for (j = 0; j < number; j++) {
1238               int iSequence = index[j];
1239               double thisWeight = weight[iSequence];
1240               // row has -1
1241               double pivot = - updateBy[iSequence];
1242               updateBy[iSequence] = 0.0;
1243               double value = pivot * pivot * devex_;
1244               if (reference(iSequence + numberColumns))
1245                    value += 1.0;
1246               weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1247          }
1248
1249          // columns
1250          number = spareColumn1->getNumElements();
1251          index = spareColumn1->getIndices();
1252          updateBy = spareColumn1->denseVector();
1253          int addSequence = model_->maximumAbcNumberRows();
1254          for (j = 0; j < number; j++) {
1255               int iSequence = index[j];
1256               double pivot=updateBy[iSequence];
1257               updateBy[iSequence] = 0.0;
1258               iSequence += addSequence;
1259               double thisWeight = weight[iSequence];
1260               // row has -1
1261               double value = pivot * pivot * devex_;
1262               if (reference(iSequence))
1263                    value += 1.0;
1264               weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1265          }
1266          // restore outgoing weight
1267          if (sequenceOut >= 0)
1268               weights_[sequenceOut] = outgoingWeight;
1269          //#define SOME_DEBUG_1
1270#ifdef SOME_DEBUG_1
1271          // check for accuracy
1272          int iCheck = 892;
1273          //printf("weight for iCheck is %g\n",weights_[iCheck]);
1274          int numberRows = model_->numberRows();
1275          //int numberColumns = model_->numberColumns();
1276          for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1277               if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
1278                         !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
1279                    checkAccuracy(iCheck, 1.0e-1, updates);
1280          }
1281#endif
1282          updates->setNumElements(0);
1283          spareColumn1->setNumElements(0);
1284     }
1285}
1286// Update weights for Devex
1287void
1288AbcPrimalColumnSteepest::justDevex(CoinIndexedVector * updates,
1289                                   CoinIndexedVector * spareColumn1)
1290{
1291     int j;
1292     int number = 0;
1293     int * index;
1294     double * updateBy;
1295     // dj could be very small (or even zero - take care)
1296     double dj = model_->dualIn();
1297     double tolerance = model_->currentDualTolerance();
1298     // we can't really trust infeasibilities if there is dual error
1299     // this coding has to mimic coding in checkDualSolution
1300     double error = CoinMin(1.0e-2, model_->largestDualError());
1301     // allow tolerance at least slightly bigger than standard
1302     tolerance = tolerance +  error;
1303     int pivotRow = model_->pivotRow();
1304     // for weights update we use pivotSequence
1305     pivotRow = pivotSequence_;
1306     assert (pivotRow >= 0);
1307     // make sure infeasibility on incoming is 0.0
1308     const int * pivotVariable = model_->pivotVariable();
1309     int sequenceIn = pivotVariable[pivotRow];
1310     infeasible_->zero(sequenceIn);
1311     // and we can see if reference
1312     //double referenceIn = 0.0;
1313     //if (mode_ != 1 && reference(sequenceIn))
1314     //   referenceIn = 1.0;
1315     // save outgoing weight round update
1316     double outgoingWeight = 0.0;
1317     int sequenceOut = model_->sequenceOut();
1318     if (sequenceOut >= 0)
1319          outgoingWeight = weights_[sequenceOut];
1320     assert (!updates->getNumElements()); 
1321     assert (!spareColumn1->getNumElements());
1322     // unset in case sub flip
1323     pivotSequence_ = -1;
1324     // might as well set dj to 1
1325     dj = -1.0;
1326     updates->createUnpacked(1, &pivotRow, &dj);
1327     model_->factorization()->updateColumnTranspose(*updates);
1328     // put row of tableau in rowArray and columnArray
1329     model_->abcMatrix()->transposeTimes(*updates, *spareColumn1);
1330     double * weight;
1331     // rows
1332     number = updates->getNumElements();
1333     index = updates->getIndices();
1334     updateBy = updates->denseVector();
1335     weight = weights_;
1336
1337     // Devex
1338     assert (devex_ > 0.0);
1339     for (j = 0; j < number; j++) {
1340          int iSequence = index[j];
1341          double thisWeight = weight[iSequence];
1342          // row has -1
1343          double pivot = - updateBy[iSequence];
1344          updateBy[iSequence] = 0.0;
1345          double value = pivot * pivot * devex_;
1346          if (reference(iSequence))
1347               value += 1.0;
1348          weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1349     }
1350
1351     // columns
1352     int addSequence = model_->maximumAbcNumberRows();
1353     number = spareColumn1->getNumElements();
1354     index = spareColumn1->getIndices();
1355     updateBy = spareColumn1->denseVector();
1356     // Devex
1357     for (j = 0; j < number; j++) {
1358          int iSequence = index[j];
1359          double pivot=updateBy[iSequence];
1360          updateBy[iSequence] = 0.0;
1361          iSequence += addSequence;
1362          double thisWeight = weight[iSequence];
1363          // row has -1
1364          double value = pivot * pivot * devex_;
1365          if (reference(iSequence))
1366               value += 1.0;
1367          weight[iSequence] = CoinMax(0.99 * thisWeight, value);
1368     }
1369     // restore outgoing weight
1370     if (sequenceOut >= 0)
1371          weights_[sequenceOut] = outgoingWeight;
1372     //#define SOME_DEBUG_1
1373#ifdef SOME_DEBUG_1
1374     // check for accuracy
1375     int iCheck = 892;
1376     //printf("weight for iCheck is %g\n",weights_[iCheck]);
1377     int numberRows = model_->numberRows();
1378     //int numberColumns = model_->numberColumns();
1379     for (iCheck = 0; iCheck < numberRows + numberColumns; iCheck++) {
1380          if (model_->getInternalStatus(iCheck) != AbcSimplex::basic &&
1381                    !model_->getInternalStatus(iCheck) != AbcSimplex::isFixed)
1382               checkAccuracy(iCheck, 1.0e-1, updates);
1383     }
1384#endif
1385     updates->setNumElements(0);
1386     spareColumn1->setNumElements(0);
1387}
1388// Called when maximum pivots changes
1389void
1390AbcPrimalColumnSteepest::maximumPivotsChanged()
1391{
1392     if (alternateWeights_ &&
1393               alternateWeights_->capacity() != model_->numberRows() +
1394               model_->factorization()->maximumPivots()) {
1395          delete alternateWeights_;
1396          alternateWeights_ = new CoinIndexedVector();
1397          // enough space so can use it for factorization
1398          // enoughfor ordered
1399          alternateWeights_->reserve(2*model_->numberRows() +
1400                                     model_->factorization()->maximumPivots());
1401     }
1402}
1403/*
1404   1) before factorization
1405   2) after factorization
1406   3) just redo infeasibilities
1407   4) restore weights
1408   5) at end of values pass (so need initialization)
1409*/
1410void
1411AbcPrimalColumnSteepest::saveWeights(AbcSimplex * model, int mode)
1412{
1413  model_ = model;
1414  if (mode_ == 4 || mode_ == 5) {
1415    if (mode == 1 && !weights_)
1416      numberSwitched_ = 0; // Reset
1417  }
1418  // alternateWeights_ is defined as indexed but is treated oddly
1419  // at times
1420  int numberRows = model_->numberRows();
1421  int numberColumns = model_->numberColumns();
1422  int maximumRows=model_->maximumAbcNumberRows();
1423  const int * pivotVariable = model_->pivotVariable();
1424  bool doInfeasibilities = true;
1425  if (mode == 1) {
1426    if(weights_) {
1427      // Check if size has changed
1428      if (infeasible_->capacity() == numberRows + numberColumns &&
1429          alternateWeights_->capacity() == 2*numberRows +
1430          model_->factorization()->maximumPivots()) {
1431        //alternateWeights_->clear();
1432        if (pivotSequence_ >= 0 && pivotSequence_ < numberRows) {
1433          // save pivot order
1434          CoinMemcpyN(pivotVariable,
1435                      numberRows, alternateWeights_->getIndices());
1436          // change from pivot row number to sequence number
1437          pivotSequence_ = pivotVariable[pivotSequence_];
1438        } else {
1439          pivotSequence_ = -1;
1440        }
1441        state_ = 1;
1442      } else {
1443        // size has changed - clear everything
1444        delete [] weights_;
1445        weights_ = NULL;
1446        delete infeasible_;
1447        infeasible_ = NULL;
1448        delete alternateWeights_;
1449        alternateWeights_ = NULL;
1450        delete [] savedWeights_;
1451        savedWeights_ = NULL;
1452        delete [] reference_;
1453        reference_ = NULL;
1454        state_ = -1;
1455        pivotSequence_ = -1;
1456      }
1457    }
1458  } else if (mode == 2 || mode == 4 || mode == 5) {
1459    // restore
1460    if (!weights_ || state_ == -1 || mode == 5) {
1461      // Partial is only allowed with certain types of matrix
1462      if ((mode_ != 4 && mode_ != 5) || numberSwitched_ ) {
1463        // initialize weights
1464        delete [] weights_;
1465        delete alternateWeights_;
1466        weights_ = new double[numberRows+numberColumns];
1467        alternateWeights_ = new CoinIndexedVector();
1468        // enough space so can use it for factorization
1469        alternateWeights_->reserve(2*numberRows +
1470                                   model_->factorization()->maximumPivots());
1471        initializeWeights();
1472        // create saved weights
1473        delete [] savedWeights_;
1474        savedWeights_ = CoinCopyOfArray(weights_, numberRows + numberColumns);
1475        // just do initialization
1476        mode = 3;
1477      } else {
1478        // Partial pricing
1479        // use region as somewhere to save non-fixed slacks
1480        // set up infeasibilities
1481        if (!infeasible_) {
1482          infeasible_ = new CoinIndexedVector();
1483          infeasible_->reserve(numberColumns + numberRows);
1484        }
1485        infeasible_->clear();
1486        int number = model_->numberRows();
1487        int iSequence;
1488        int numberLook = 0;
1489        int * which = infeasible_->getIndices();
1490        for (iSequence = 0; iSequence < number; iSequence++) {
1491          AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1492          if (status != AbcSimplex::isFixed)
1493            which[numberLook++] = iSequence;
1494        }
1495        infeasible_->setNumElements(numberLook);
1496        doInfeasibilities = false;
1497      }
1498      savedPivotSequence_ = -2;
1499      savedSequenceOut_ = -2;
1500    } else {
1501      if (mode != 4) {
1502        // save
1503        CoinMemcpyN(weights_, (numberRows + numberColumns), savedWeights_);
1504        savedPivotSequence_ = pivotSequence_;
1505        savedSequenceOut_ = model_->sequenceOut();
1506      } else {
1507        // restore
1508        CoinMemcpyN(savedWeights_, (numberRows + numberColumns), weights_);
1509        // was - but I think should not be
1510        //pivotSequence_= savedPivotSequence_;
1511        //model_->setSequenceOut(savedSequenceOut_);
1512        pivotSequence_ = -1;
1513        model_->setSequenceOut(-1);
1514        // indices are wrong so clear by hand
1515        //alternateWeights_->clear();
1516        CoinZeroN(alternateWeights_->denseVector(),
1517                  alternateWeights_->capacity());
1518        alternateWeights_->setNumElements(0);
1519      }
1520    }
1521    state_ = 0;
1522    // set up infeasibilities
1523    if (!infeasible_) {
1524      infeasible_ = new CoinIndexedVector();
1525      infeasible_->reserve(numberColumns + numberRows);
1526    }
1527  }
1528  if (mode >= 2 && mode != 5) {
1529    if (mode != 3) {
1530      if (pivotSequence_ >= 0) {
1531        // restore pivot row
1532        int iRow;
1533        // permute alternateWeights
1534        int iVector=model_->getAvailableArray();
1535        double * temp = model_->usefulArray(iVector)->denseVector();;
1536        double * work = alternateWeights_->denseVector();
1537        int * savePivotOrder = model_->usefulArray(iVector)->getIndices();
1538        int * oldPivotOrder = alternateWeights_->getIndices();
1539        for (iRow = 0; iRow < numberRows; iRow++) {
1540          int iPivot = oldPivotOrder[iRow];
1541          temp[iPivot] = work[iRow];
1542          savePivotOrder[iRow] = iPivot;
1543        }
1544        int number = 0;
1545        int found = -1;
1546        int * which = oldPivotOrder;
1547        // find pivot row and re-create alternateWeights
1548        for (iRow = 0; iRow < numberRows; iRow++) {
1549          int iPivot = pivotVariable[iRow];
1550          if (iPivot == pivotSequence_)
1551            found = iRow;
1552          work[iRow] = temp[iPivot];
1553          if (work[iRow])
1554            which[number++] = iRow;
1555        }
1556        alternateWeights_->setNumElements(number);
1557#ifdef CLP_DEBUG
1558        // Can happen but I should clean up
1559        assert(found >= 0);
1560#endif
1561        pivotSequence_ = found;
1562        for (iRow = 0; iRow < numberRows; iRow++) {
1563          int iPivot = savePivotOrder[iRow];
1564          temp[iPivot] = 0.0;
1565        }
1566        model_->setAvailableArray(iVector);
1567      } else {
1568        // Just clean up
1569        if (alternateWeights_)
1570          alternateWeights_->clear();
1571      }
1572    }
1573    // Save size of factorization
1574    if (!model->factorization()->pivots())
1575      sizeFactorization_ = model_->factorization()->numberElements();
1576    if(!doInfeasibilities)
1577      return; // don't disturb infeasibilities
1578    infeasible_->clear();
1579    double tolerance = model_->currentDualTolerance();
1580    int number = model_->numberRows() + model_->numberColumns();
1581    int iSequence;
1582   
1583    double * reducedCost = model_->djRegion();
1584#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
1585    for (iSequence = 0; iSequence < number; iSequence++) {
1586      double value = reducedCost[iSequence];
1587      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1588     
1589      switch(status) {
1590       
1591      case AbcSimplex::basic:
1592      case AbcSimplex::isFixed:
1593        break;
1594      case AbcSimplex::isFree:
1595      case AbcSimplex::superBasic:
1596        if (fabs(value) > FREE_ACCEPT * tolerance) {
1597          // we are going to bias towards free (but only if reasonable)
1598          value *= FREE_BIAS;
1599          // store square in list
1600          infeasible_->quickAdd(iSequence, value * value);
1601        }
1602        break;
1603      case AbcSimplex::atUpperBound:
1604        if (value > tolerance) {
1605          infeasible_->quickAdd(iSequence, value * value);
1606        }
1607        break;
1608      case AbcSimplex::atLowerBound:
1609        if (value < -tolerance) {
1610          infeasible_->quickAdd(iSequence, value * value);
1611        }
1612      }
1613    }
1614#else
1615    // Columns
1616    for (iSequence = maximumRows; iSequence < number; iSequence++) {
1617      double value = reducedCost[iSequence];
1618      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1619     
1620      switch(status) {
1621       
1622      case AbcSimplex::basic:
1623      case AbcSimplex::isFixed:
1624        break;
1625      case AbcSimplex::isFree:
1626      case AbcSimplex::superBasic:
1627        if (fabs(value) > FREE_ACCEPT * tolerance) {
1628          // we are going to bias towards free (but only if reasonable)
1629          value *= FREE_BIAS;
1630          // store square in list
1631          infeasible_->quickAdd(iSequence, value * value);
1632        }
1633        break;
1634      case AbcSimplex::atUpperBound:
1635        if (value > tolerance) {
1636          infeasible_->quickAdd(iSequence, value * value);
1637        }
1638        break;
1639      case AbcSimplex::atLowerBound:
1640        if (value < -tolerance) {
1641          infeasible_->quickAdd(iSequence, value * value);
1642        }
1643      }
1644    }
1645    // Rows
1646    for (iSequence=0 ; iSequence < numberRows; iSequence++) {
1647      double value = reducedCost[iSequence];
1648      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1649     
1650      switch(status) {
1651       
1652      case AbcSimplex::basic:
1653      case AbcSimplex::isFixed:
1654        break;
1655      case AbcSimplex::isFree:
1656      case AbcSimplex::superBasic:
1657        if (fabs(value) > FREE_ACCEPT * tolerance) {
1658          // we are going to bias towards free (but only if reasonable)
1659          value *= FREE_BIAS;
1660          // store square in list
1661          infeasible_->quickAdd(iSequence, value * value);
1662        }
1663        break;
1664      case AbcSimplex::atUpperBound:
1665        if (value > tolerance) {
1666          infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
1667        }
1668        break;
1669      case AbcSimplex::atLowerBound:
1670        if (value < -tolerance) {
1671          infeasible_->quickAdd(iSequence, value * value * CLP_PRIMAL_SLACK_MULTIPLIER);
1672        }
1673      }
1674    }
1675#endif
1676  } 
1677}
1678// Gets rid of last update
1679void
1680AbcPrimalColumnSteepest::unrollWeights()
1681{
1682     if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
1683          return;
1684     double * saved = alternateWeights_->denseVector();
1685     int number = alternateWeights_->getNumElements();
1686     int * which = alternateWeights_->getIndices();
1687     int i;
1688     for (i = 0; i < number; i++) {
1689          int iRow = which[i];
1690          weights_[iRow] = saved[iRow];
1691          saved[iRow] = 0.0;
1692     }
1693     alternateWeights_->setNumElements(0);
1694}
1695
1696//-------------------------------------------------------------------
1697// Clone
1698//-------------------------------------------------------------------
1699AbcPrimalColumnPivot * AbcPrimalColumnSteepest::clone(bool CopyData) const
1700{
1701     if (CopyData) {
1702          return new AbcPrimalColumnSteepest(*this);
1703     } else {
1704          return new AbcPrimalColumnSteepest();
1705     }
1706}
1707void
1708AbcPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
1709{
1710     // Local copy of mode so can decide what to do
1711     int switchType = mode_;
1712     if (mode_ == 4 && numberSwitched_)
1713          switchType = 3;
1714     else if (mode_ == 4 || mode_ == 5)
1715          return;
1716     int number = input->getNumElements();
1717     int * which = input->getIndices();
1718     double * work = input->denseVector();
1719     int newNumber = 0;
1720     int * newWhich = alternateWeights_->getIndices();
1721     double * newWork = alternateWeights_->denseVector();
1722     int i;
1723     int sequenceIn = model_->sequenceIn();
1724     int sequenceOut = model_->sequenceOut();
1725     const int * pivotVariable = model_->pivotVariable();
1726
1727     int pivotRow = model_->pivotRow();
1728     pivotSequence_ = pivotRow;
1729
1730     devex_ = 0.0;
1731          if (pivotRow >= 0) {
1732               if (switchType == 1) {
1733                    for (i = 0; i < number; i++) {
1734                         int iRow = which[i];
1735                         devex_ += work[iRow] * work[iRow];
1736                         newWork[iRow] = -2.0 * work[iRow];
1737                    }
1738                    newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
1739                    devex_ += ADD_ONE;
1740                    weights_[sequenceOut] = 1.0 + ADD_ONE;
1741                    CoinMemcpyN(which, number, newWhich);
1742                    alternateWeights_->setNumElements(number);
1743               } else {
1744                    if ((mode_ != 4 && mode_ != 5) || numberSwitched_ > 1) {
1745                         for (i = 0; i < number; i++) {
1746                              int iRow = which[i];
1747                              int iPivot = pivotVariable[iRow];
1748                              if (reference(iPivot)) {
1749                                   devex_ += work[iRow] * work[iRow];
1750                                   newWork[iRow] = -2.0 * work[iRow];
1751                                   newWhich[newNumber++] = iRow;
1752                              }
1753                         }
1754                         if (!newWork[pivotRow] && devex_ > 0.0)
1755                              newWhich[newNumber++] = pivotRow; // add if not already in
1756                         newWork[pivotRow] = -2.0 * CoinMax(devex_, 0.0);
1757                    } else {
1758                         for (i = 0; i < number; i++) {
1759                              int iRow = which[i];
1760                              int iPivot = pivotVariable[iRow];
1761                              if (reference(iPivot))
1762                                   devex_ += work[iRow] * work[iRow];
1763                         }
1764                    }
1765                    if (reference(sequenceIn)) {
1766                         devex_ += 1.0;
1767                    } else {
1768                    }
1769                    if (reference(sequenceOut)) {
1770                         weights_[sequenceOut] = 1.0 + 1.0;
1771                    } else {
1772                         weights_[sequenceOut] = 1.0;
1773                    }
1774                    alternateWeights_->setNumElements(newNumber);
1775               }
1776          } else {
1777               if (switchType == 1) {
1778                    for (i = 0; i < number; i++) {
1779                         int iRow = which[i];
1780                         devex_ += work[iRow] * work[iRow];
1781                    }
1782                    devex_ += ADD_ONE;
1783               } else {
1784                    for (i = 0; i < number; i++) {
1785                         int iRow = which[i];
1786                         int iPivot = pivotVariable[iRow];
1787                         if (reference(iPivot)) {
1788                              devex_ += work[iRow] * work[iRow];
1789                         }
1790                    }
1791                    if (reference(sequenceIn))
1792                         devex_ += 1.0;
1793               }
1794          }
1795     double oldDevex = weights_[sequenceIn];
1796#ifdef CLP_DEBUG
1797     if ((model_->messageHandler()->logLevel() & 32))
1798          printf("old weight %g, new %g\n", oldDevex, devex_);
1799#endif
1800     double check = CoinMax(devex_, oldDevex) + 0.1;
1801     weights_[sequenceIn] = devex_;
1802     double testValue = 0.1;
1803     if (mode_ == 4 && numberSwitched_ == 1)
1804          testValue = 0.5;
1805     if ( fabs ( devex_ - oldDevex ) > testValue * check ) {
1806#ifdef CLP_DEBUG
1807          if ((model_->messageHandler()->logLevel() & 48) == 16)
1808               printf("old weight %g, new %g\n", oldDevex, devex_);
1809#endif
1810          //printf("old weight %g, new %g\n",oldDevex,devex_);
1811          testValue = 0.99;
1812          if (mode_ == 1)
1813               testValue = 1.01e1; // make unlikely to do if steepest
1814          else if (mode_ == 4 && numberSwitched_ == 1)
1815               testValue = 0.9;
1816          double difference = fabs(devex_ - oldDevex);
1817          if ( difference > testValue * check ) {
1818               // need to redo
1819               model_->messageHandler()->message(CLP_INITIALIZE_STEEP,
1820                                                 *model_->messagesPointer())
1821                         << oldDevex << devex_
1822                         << CoinMessageEol;
1823               initializeWeights();
1824          }
1825     }
1826     if (pivotRow >= 0) {
1827          // set outgoing weight here
1828          weights_[model_->sequenceOut()] = devex_ / (model_->alpha() * model_->alpha());
1829     }
1830}
1831// Checks accuracy - just for debug
1832void
1833AbcPrimalColumnSteepest::checkAccuracy(int sequence,
1834                                       double relativeTolerance,
1835                                       CoinIndexedVector * rowArray1)
1836{
1837     if ((mode_ == 4 || mode_ == 5) && !numberSwitched_)
1838          return;
1839     model_->unpack(*rowArray1, sequence);
1840     model_->factorization()->updateColumn(*rowArray1);
1841     int number = rowArray1->getNumElements();
1842     int * which = rowArray1->getIndices();
1843     double * work = rowArray1->denseVector();
1844     const int * pivotVariable = model_->pivotVariable();
1845
1846     double devex = 0.0;
1847     int i;
1848
1849     if (mode_ == 1) {
1850          for (i = 0; i < number; i++) {
1851               int iRow = which[i];
1852               devex += work[iRow] * work[iRow];
1853               work[iRow] = 0.0;
1854          }
1855          devex += ADD_ONE;
1856     } else {
1857          for (i = 0; i < number; i++) {
1858               int iRow = which[i];
1859               int iPivot = pivotVariable[iRow];
1860               if (reference(iPivot)) {
1861                    devex += work[iRow] * work[iRow];
1862               }
1863               work[iRow] = 0.0;
1864          }
1865          if (reference(sequence))
1866               devex += 1.0;
1867     }
1868
1869     double oldDevex = weights_[sequence];
1870     double check = CoinMax(devex, oldDevex);;
1871     if ( fabs ( devex - oldDevex ) > relativeTolerance * check ) {
1872       COIN_DETAIL_PRINT(printf("check %d old weight %g, new %g\n", sequence, oldDevex, devex));
1873          // update so won't print again
1874          weights_[sequence] = devex;
1875     }
1876     rowArray1->setNumElements(0);
1877}
1878
1879// Initialize weights
1880void
1881AbcPrimalColumnSteepest::initializeWeights()
1882{
1883     int numberRows = model_->numberRows();
1884     int numberColumns = model_->numberColumns();
1885     int number = numberRows + numberColumns;
1886     int iSequence;
1887     if (mode_ != 1) {
1888          // initialize to 1.0
1889          // and set reference framework
1890          if (!reference_) {
1891               int nWords = (number + 31) >> 5;
1892               reference_ = new unsigned int[nWords];
1893               CoinZeroN(reference_, nWords);
1894          }
1895
1896          for (iSequence = 0; iSequence < number; iSequence++) {
1897               weights_[iSequence] = 1.0;
1898               if (model_->getInternalStatus(iSequence) == AbcSimplex::basic) {
1899                    setReference(iSequence, false);
1900               } else {
1901                    setReference(iSequence, true);
1902               }
1903          }
1904     } else {
1905          CoinIndexedVector * temp = new CoinIndexedVector();
1906          temp->reserve(numberRows +
1907                        model_->factorization()->maximumPivots());
1908          double * array = alternateWeights_->denseVector();
1909          int * which = alternateWeights_->getIndices();
1910
1911          for (iSequence = 0; iSequence < number; iSequence++) {
1912               weights_[iSequence] = 1.0 + ADD_ONE;
1913               if (model_->getInternalStatus(iSequence) != AbcSimplex::basic &&
1914                         model_->getInternalStatus(iSequence) != AbcSimplex::isFixed) {
1915                    model_->unpack(*alternateWeights_, iSequence);
1916                    double value = ADD_ONE;
1917                    model_->factorization()->updateColumn(*alternateWeights_);
1918                    int number = alternateWeights_->getNumElements();
1919                    int j;
1920                    for (j = 0; j < number; j++) {
1921                         int iRow = which[j];
1922                         value += array[iRow] * array[iRow];
1923                         array[iRow] = 0.0;
1924                    }
1925                    alternateWeights_->setNumElements(0);
1926                    weights_[iSequence] = value;
1927               }
1928          }
1929          delete temp;
1930     }
1931}
1932// Gets rid of all arrays
1933void
1934AbcPrimalColumnSteepest::clearArrays()
1935{
1936     if (persistence_ == normal) {
1937          delete [] weights_;
1938          weights_ = NULL;
1939          delete infeasible_;
1940          infeasible_ = NULL;
1941          delete alternateWeights_;
1942          alternateWeights_ = NULL;
1943          delete [] savedWeights_;
1944          savedWeights_ = NULL;
1945          delete [] reference_;
1946          reference_ = NULL;
1947     }
1948     pivotSequence_ = -1;
1949     state_ = -1;
1950     savedPivotSequence_ = -1;
1951     savedSequenceOut_ = -1;
1952     devex_ = 0.0;
1953}
1954// Returns true if would not find any column
1955bool
1956AbcPrimalColumnSteepest::looksOptimal() const
1957{
1958     if (looksOptimal_)
1959          return true; // user overrode
1960     //**** THIS MUST MATCH the action coding above
1961     double tolerance = model_->currentDualTolerance();
1962     // we can't really trust infeasibilities if there is dual error
1963     // this coding has to mimic coding in checkDualSolution
1964     double error = CoinMin(1.0e-2, model_->largestDualError());
1965     // allow tolerance at least slightly bigger than standard
1966     tolerance = tolerance +  error;
1967     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
1968          // we can't really trust infeasibilities if there is dual error
1969          double checkTolerance = 1.0e-8;
1970          if (!model_->factorization()->pivots())
1971               checkTolerance = 1.0e-6;
1972          if (model_->largestDualError() > checkTolerance)
1973               tolerance *= model_->largestDualError() / checkTolerance;
1974          // But cap
1975          tolerance = CoinMin(1000.0, tolerance);
1976     }
1977     int number = model_->numberRows() + model_->numberColumns();
1978     int iSequence;
1979
1980     double * reducedCost = model_->djRegion();
1981     int numberInfeasible = 0;
1982          for (iSequence = 0; iSequence < number; iSequence++) {
1983               double value = reducedCost[iSequence];
1984               AbcSimplex::Status status = model_->getInternalStatus(iSequence);
1985
1986               switch(status) {
1987
1988               case AbcSimplex::basic:
1989               case AbcSimplex::isFixed:
1990                    break;
1991               case AbcSimplex::isFree:
1992               case AbcSimplex::superBasic:
1993                    if (fabs(value) > FREE_ACCEPT * tolerance)
1994                         numberInfeasible++;
1995                    break;
1996               case AbcSimplex::atUpperBound:
1997                    if (value > tolerance)
1998                         numberInfeasible++;
1999                    break;
2000               case AbcSimplex::atLowerBound:
2001                    if (value < -tolerance)
2002                         numberInfeasible++;
2003               }
2004          }
2005     return numberInfeasible == 0;
2006}
2007// Update djs doing partial pricing (dantzig)
2008int
2009AbcPrimalColumnSteepest::partialPricing(CoinIndexedVector * updates,
2010                                        int numberWanted,
2011                                        int numberLook)
2012{
2013     int number = 0;
2014     int * index;
2015     double * updateBy;
2016     double * reducedCost;
2017     double saveTolerance = model_->currentDualTolerance();
2018     double tolerance = model_->currentDualTolerance();
2019     // we can't really trust infeasibilities if there is dual error
2020     // this coding has to mimic coding in checkDualSolution
2021     double error = CoinMin(1.0e-2, model_->largestDualError());
2022     // allow tolerance at least slightly bigger than standard
2023     tolerance = tolerance +  error;
2024     if(model_->numberIterations() < model_->lastBadIteration() + 200) {
2025          // we can't really trust infeasibilities if there is dual error
2026          double checkTolerance = 1.0e-8;
2027          if (!model_->factorization()->pivots())
2028               checkTolerance = 1.0e-6;
2029          if (model_->largestDualError() > checkTolerance)
2030               tolerance *= model_->largestDualError() / checkTolerance;
2031          // But cap
2032          tolerance = CoinMin(1000.0, tolerance);
2033     }
2034     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
2035       tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
2036     // So partial pricing can use
2037     model_->setCurrentDualTolerance(tolerance);
2038     model_->factorization()->updateColumnTranspose(*updates);
2039     int numberColumns = model_->numberColumns();
2040
2041     // Rows
2042     reducedCost = model_->djRegion();
2043
2044     number = updates->getNumElements();
2045     index = updates->getIndices();
2046     updateBy = updates->denseVector();
2047     int j;
2048     double * duals = model_->dualRowSolution();
2049     for (j = 0; j < number; j++) {
2050          int iSequence = index[j];
2051          double value = duals[iSequence];
2052          value -= updateBy[iSequence];
2053          updateBy[iSequence] = 0.0;
2054          duals[iSequence] = value;
2055     }
2056     //#define CLP_DEBUG
2057#ifdef CLP_DEBUG
2058     // check duals
2059     {
2060          //work space
2061          CoinIndexedVector arrayVector;
2062          arrayVector.reserve(numberRows + 1000);
2063          CoinIndexedVector workSpace;
2064          workSpace.reserve(numberRows + 1000);
2065
2066
2067          int iRow;
2068          double * array = arrayVector.denseVector();
2069          int * index = arrayVector.getIndices();
2070          int number = 0;
2071          int * pivotVariable = model_->pivotVariable();
2072          double * cost = model_->costRegion();
2073          for (iRow = 0; iRow < numberRows; iRow++) {
2074               int iPivot = pivotVariable[iRow];
2075               double value = cost[iPivot];
2076               if (value) {
2077                    array[iRow] = value;
2078                    index[number++] = iRow;
2079               }
2080          }
2081          arrayVector.setNumElements(number);
2082
2083          // Btran basic costs
2084          model_->factorization()->updateColumnTranspose(&workSpace, &arrayVector);
2085
2086          // now look at dual solution
2087          for (iRow = 0; iRow < numberRows; iRow++) {
2088               // slack
2089               double value = array[iRow];
2090               if (fabs(duals[iRow] - value) > 1.0e-3)
2091                    printf("bad row %d old dual %g new %g\n", iRow, duals[iRow], value);
2092               //duals[iRow]=value;
2093          }
2094     }
2095#endif
2096#undef CLP_DEBUG
2097     double bestDj = tolerance;
2098     int bestSequence = -1;
2099
2100     const double * cost = model_->costRegion();
2101
2102     model_->abcMatrix()->setOriginalWanted(numberWanted);
2103     model_->abcMatrix()->setCurrentWanted(numberWanted);
2104     int iPassR = 0, iPassC = 0;
2105     // Setup two passes
2106     // This biases towards picking row variables
2107     // This probably should be fixed
2108     int startR[4];
2109     const int * which = infeasible_->getIndices();
2110     int nSlacks = infeasible_->getNumElements();
2111     startR[1] = nSlacks;
2112     startR[2] = 0;
2113     double randomR = model_->randomNumberGenerator()->randomDouble();
2114     double dstart = static_cast<double> (nSlacks) * randomR;
2115     startR[0] = static_cast<int> (dstart);
2116     startR[3] = startR[0];
2117     double startC[4];
2118     startC[1] = 1.0;
2119     startC[2] = 0;
2120     double randomC = model_->randomNumberGenerator()->randomDouble();
2121     startC[0] = randomC;
2122     startC[3] = randomC;
2123     reducedCost = model_->djRegion();
2124     int sequenceOut = model_->sequenceOut();
2125     int chunk = CoinMin(1024, (numberColumns + nSlacks) / 32);
2126#ifdef COIN_DETAIL
2127     if (model_->numberIterations() % 1000 == 0 && model_->logLevel() > 1) {
2128          printf("%d wanted, nSlacks %d, chunk %d\n", numberWanted, nSlacks, chunk);
2129          int i;
2130          for (i = 0; i < 4; i++)
2131               printf("start R %d C %g ", startR[i], startC[i]);
2132          printf("\n");
2133     }
2134#endif
2135     chunk = CoinMax(chunk, 256);
2136     bool finishedR = false, finishedC = false;
2137     bool doingR = randomR > randomC;
2138     //doingR=false;
2139     int saveNumberWanted = numberWanted;
2140     while (!finishedR || !finishedC) {
2141          if (finishedR)
2142               doingR = false;
2143          if (doingR) {
2144               int saveSequence = bestSequence;
2145               int start = startR[iPassR];
2146               int end = CoinMin(startR[iPassR+1], start + chunk / 2);
2147               int jSequence;
2148               for (jSequence = start; jSequence < end; jSequence++) {
2149                    int iSequence = which[jSequence];
2150                    if (iSequence != sequenceOut) {
2151                         double value;
2152                         AbcSimplex::Status status = model_->getInternalStatus(iSequence);
2153
2154                         switch(status) {
2155
2156                         case AbcSimplex::basic:
2157                         case AbcSimplex::isFixed:
2158                              break;
2159                         case AbcSimplex::isFree:
2160                         case AbcSimplex::superBasic:
2161                              value = fabs(cost[iSequence] - duals[iSequence]);
2162                              if (value > FREE_ACCEPT * tolerance) {
2163                                   numberWanted--;
2164                                   // we are going to bias towards free (but only if reasonable)
2165                                   value *= FREE_BIAS;
2166                                   if (value > bestDj) {
2167                                        // check flagged variable and correct dj
2168                                        if (!model_->flagged(iSequence)) {
2169                                             bestDj = value;
2170                                             bestSequence = iSequence;
2171                                        } else {
2172                                             // just to make sure we don't exit before got something
2173                                             numberWanted++;
2174                                        }
2175                                   }
2176                              }
2177                              break;
2178                         case AbcSimplex::atUpperBound:
2179                              value = cost[iSequence] - duals[iSequence];
2180                              if (value > tolerance) {
2181                                   numberWanted--;
2182                                   if (value > bestDj) {
2183                                        // check flagged variable and correct dj
2184                                        if (!model_->flagged(iSequence)) {
2185                                             bestDj = value;
2186                                             bestSequence = iSequence;
2187                                        } else {
2188                                             // just to make sure we don't exit before got something
2189                                             numberWanted++;
2190                                        }
2191                                   }
2192                              }
2193                              break;
2194                         case AbcSimplex::atLowerBound:
2195                              value = -(cost[iSequence] - duals[iSequence]);
2196                              if (value > tolerance) {
2197                                   numberWanted--;
2198                                   if (value > bestDj) {
2199                                        // check flagged variable and correct dj
2200                                        if (!model_->flagged(iSequence)) {
2201                                             bestDj = value;
2202                                             bestSequence = iSequence;
2203                                        } else {
2204                                             // just to make sure we don't exit before got something
2205                                             numberWanted++;
2206                                        }
2207                                   }
2208                              }
2209                              break;
2210                         }
2211                    }
2212                    if (!numberWanted)
2213                         break;
2214               }
2215               numberLook -= (end - start);
2216               if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
2217                    numberWanted = 0; // give up
2218               if (saveSequence != bestSequence) {
2219                    // dj
2220                    reducedCost[bestSequence] = cost[bestSequence] - duals[bestSequence];
2221                    bestDj = fabs(reducedCost[bestSequence]);
2222                    model_->abcMatrix()->setSavedBestSequence(bestSequence);
2223                    model_->abcMatrix()->setSavedBestDj(reducedCost[bestSequence]);
2224               }
2225               model_->abcMatrix()->setCurrentWanted(numberWanted);
2226               if (!numberWanted)
2227                    break;
2228               doingR = false;
2229               // update start
2230               startR[iPassR] = jSequence;
2231               if (jSequence >= startR[iPassR+1]) {
2232                    if (iPassR)
2233                         finishedR = true;
2234                    else
2235                         iPassR = 2;
2236               }
2237          }
2238          if (finishedC)
2239               doingR = true;
2240          if (!doingR) {
2241            // temp
2242            int saveSequence = bestSequence;
2243            // Columns
2244            double start = startC[iPassC];
2245            // If we put this idea back then each function needs to update endFraction **
2246#if 0
2247            double dchunk = (static_cast<double> chunk) / (static_cast<double> numberColumns);
2248            double end = CoinMin(startC[iPassC+1], start + dchunk);;
2249#else
2250            double end = startC[iPassC+1]; // force end
2251#endif
2252            model_->abcMatrix()->partialPricing(start, end, bestSequence, numberWanted);
2253            numberWanted = model_->abcMatrix()->currentWanted();
2254            numberLook -= static_cast<int> ((end - start) * numberColumns);
2255            if (numberLook < 0 && (10 * (saveNumberWanted - numberWanted) > saveNumberWanted))
2256              numberWanted = 0; // give up
2257            if (bestSequence!=saveSequence) {
2258              // dj
2259              bestDj = fabs(reducedCost[bestSequence]);
2260            }
2261            if (!numberWanted)
2262              break;
2263            doingR = true;
2264            // update start
2265            startC[iPassC] = end;
2266            if (end >= startC[iPassC+1] - 1.0e-8) {
2267              if (iPassC)
2268                finishedC = true;
2269              else
2270                iPassC = 2;
2271            }
2272          }
2273     }
2274     updates->setNumElements(0);
2275
2276     // Restore tolerance
2277     model_->setCurrentDualTolerance(saveTolerance);
2278#ifndef NDEBUG
2279     if (bestSequence >= 0) {
2280          if (model_->getInternalStatus(bestSequence) == AbcSimplex::atLowerBound)
2281               assert(model_->reducedCost(bestSequence) < 0.0);
2282          if (model_->getInternalStatus(bestSequence) == AbcSimplex::atUpperBound)
2283               assert(model_->reducedCost(bestSequence) > 0.0);
2284     }
2285#endif
2286     return bestSequence;
2287}
2288#endif
Note: See TracBrowser for help on using the repository browser.