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

Last change on this file since 2470 was 2385, checked in by unxusr, 10 months ago

formatting

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