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

Last change on this file since 2030 was 2025, checked in by forrest, 6 years ago

mainly dantzig-wolfe - also Clp strong branching (rather than OsiClp?)

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