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

Last change on this file since 2470 was 2429, checked in by stefan, 8 months ago

change more reinterpret_cast from NULL to C-cast, see also #93

  • Property svn:keywords set to Id
File size: 226.0 KB
RevLine 
[1910]1/* $Id: AbcSimplex.cpp 2429 2019-03-11 16:34:25Z stefan $ */
[1878]2// Copyright (C) 2000, 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 "ClpConfig.h"
7
8#include "CoinPragma.hpp"
9#include <math.h>
10//#define ABC_DEBUG 2
11
[2385]12#if SLIM_CLP == 2
[1878]13#define SLIM_NOIO
14#endif
15#include "CoinHelperFunctions.hpp"
16#include "CoinFloatEqual.hpp"
17#include "ClpSimplex.hpp"
18#include "AbcSimplex.hpp"
19#include "AbcSimplexDual.hpp"
20#include "AbcSimplexFactorization.hpp"
21#include "AbcNonLinearCost.hpp"
22#include "CoinAbcCommon.hpp"
23#include "AbcMatrix.hpp"
24#include "CoinIndexedVector.hpp"
25#include "AbcDualRowDantzig.hpp"
26#include "AbcDualRowSteepest.hpp"
27#include "AbcPrimalColumnDantzig.hpp"
28#include "AbcPrimalColumnSteepest.hpp"
29#include "ClpMessage.hpp"
30#include "ClpEventHandler.hpp"
31#include "ClpLinearObjective.hpp"
32#include "CoinAbcHelperFunctions.hpp"
33#include "CoinModel.hpp"
34#include "CoinLpIO.hpp"
35#include <cfloat>
36
37#include <string>
38#include <stdio.h>
39#include <iostream>
40//#############################################################################
[2385]41AbcSimplex::AbcSimplex(bool emptyMessages)
42  :
43
[1878]44  ClpSimplex(emptyMessages)
45{
46  gutsOfDelete(0);
[2385]47  gutsOfInitialize(0, 0, true);
48  assert(maximumAbcNumberRows_ >= 0);
[1878]49  //printf("XX %x AbcSimplex constructor\n",this);
50}
51
52//-----------------------------------------------------------------------------
53
[2385]54AbcSimplex::~AbcSimplex()
[1878]55{
56  //printf("XX %x AbcSimplex destructor\n",this);
57  gutsOfDelete(1);
58}
59// Copy constructor.
[2385]60AbcSimplex::AbcSimplex(const AbcSimplex &rhs)
61  : ClpSimplex(rhs)
[1878]62{
63  //printf("XX %x AbcSimplex constructor from %x\n",this,&rhs);
64  gutsOfDelete(0);
[2385]65  gutsOfInitialize(numberRows_, numberColumns_, false);
[1878]66  gutsOfCopy(rhs);
[2385]67  assert(maximumAbcNumberRows_ >= 0);
[1878]68}
69#include "ClpDualRowSteepest.hpp"
70#include "ClpPrimalColumnSteepest.hpp"
71#include "ClpFactorization.hpp"
72// Copy constructor from model
[2385]73AbcSimplex::AbcSimplex(const ClpSimplex &rhs)
74  : ClpSimplex(rhs)
[1878]75{
76  //printf("XX %x AbcSimplex constructor from ClpSimplex\n",this);
77  gutsOfDelete(0);
[2385]78  gutsOfInitialize(numberRows_, numberColumns_, true);
79  gutsOfResize(numberRows_, numberColumns_);
[1878]80  // Set up row/column selection and factorization type
[2385]81  ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(rhs.dualRowPivot());
[1878]82  if (pivot) {
83    AbcDualRowSteepest steep(pivot->mode());
84    setDualRowPivotAlgorithm(steep);
85  } else {
86    AbcDualRowDantzig dantzig;
87    setDualRowPivotAlgorithm(dantzig);
88  }
[2385]89  ClpPrimalColumnSteepest *pivotColumn = dynamic_cast< ClpPrimalColumnSteepest * >(rhs.primalColumnPivot());
[1878]90  if (pivotColumn) {
91    AbcPrimalColumnSteepest steep(pivotColumn->mode());
92    setPrimalColumnPivotAlgorithm(steep);
93  } else {
94    AbcPrimalColumnDantzig dantzig;
95    setPrimalColumnPivotAlgorithm(dantzig);
96  }
97  //if (rhs.factorization()->)
98  //factorization_->forceOtherFactorization();
[2385]99  factorization()->synchronize(rhs.factorization(), this);
[1878]100  //factorization_->setGoDenseThreshold(rhs.factorization()->goDenseThreshold());
101  //factorization_->setGoSmallThreshold(rhs.factorization()->goSmallThreshold());
[2385]102  translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
103  assert(maximumAbcNumberRows_ >= 0);
[1878]104}
105// Assignment operator. This copies the data
106AbcSimplex &
[2385]107AbcSimplex::operator=(const AbcSimplex &rhs)
[1878]108{
109  if (this != &rhs) {
110    gutsOfDelete(1);
111    ClpSimplex::operator=(rhs);
112    gutsOfCopy(rhs);
[2385]113    assert(maximumAbcNumberRows_ >= 0);
[1878]114  }
115  return *this;
116}
117// fills in perturbationSaved_ from start with 0.5+random
[2385]118void AbcSimplex::fillPerturbation(int start, int number)
[1878]119{
[2385]120  double *array = perturbationSaved_ + start;
121  for (int i = 0; i < number; i++)
122    array[i] = 0.5 + 0.5 * randomNumberGenerator_.randomDouble();
[1878]123}
124// Sets up all extra pointers
[2385]125void AbcSimplex::setupPointers(int maxRows, int maxColumns)
[1878]126{
[2385]127  int numberTotal = maxRows + maxColumns;
128  scaleToExternal_ = scaleFromExternal_ + numberTotal;
129  tempArray_ = offset_ + numberTotal;
130  lowerSaved_ = abcLower_ + numberTotal;
131  upperSaved_ = abcUpper_ + numberTotal;
132  costSaved_ = abcCost_ + numberTotal;
133  solutionSaved_ = abcSolution_ + numberTotal;
134  djSaved_ = abcDj_ + numberTotal;
135  internalStatusSaved_ = internalStatus_ + numberTotal;
136  perturbationSaved_ = abcPerturbation_ + numberTotal;
137  offsetRhs_ = tempArray_ + numberTotal;
138  lowerBasic_ = lowerSaved_ + numberTotal;
139  upperBasic_ = upperSaved_ + numberTotal;
140  costBasic_ = costSaved_ + 2 * numberTotal;
141  solutionBasic_ = solutionSaved_ + numberTotal;
142  djBasic_ = djSaved_ + numberTotal;
143  perturbationBasic_ = perturbationSaved_ + numberTotal;
144  columnUseScale_ = scaleToExternal_ + maxRows;
145  inverseColumnUseScale_ = scaleFromExternal_ + maxRows;
[1878]146}
147// Copies all saved versions to working versions and may do something for perturbation
[2385]148void AbcSimplex::copyFromSaved(int which)
[1878]149{
[2385]150  if ((which & 1) != 0)
151    CoinAbcMemcpy(abcSolution_, solutionSaved_, maximumNumberTotal_);
152  if ((which & 2) != 0)
153    CoinAbcMemcpy(abcCost_, costSaved_, maximumNumberTotal_);
154  if ((which & 4) != 0)
155    CoinAbcMemcpy(abcLower_, lowerSaved_, maximumNumberTotal_);
156  if ((which & 8) != 0)
157    CoinAbcMemcpy(abcUpper_, upperSaved_, maximumNumberTotal_);
158  if ((which & 16) != 0)
159    CoinAbcMemcpy(abcDj_, djSaved_, maximumNumberTotal_);
160  if ((which & 32) != 0) {
161    CoinAbcMemcpy(abcLower_, abcPerturbation_, numberTotal_);
162    CoinAbcMemcpy(abcUpper_, abcPerturbation_ + numberTotal_, numberTotal_);
[1878]163  }
164}
[2385]165void AbcSimplex::gutsOfCopy(const AbcSimplex &rhs)
[1878]166{
[1910]167#ifdef ABC_INHERIT
[2385]168  abcSimplex_ = this;
[1910]169#endif
[2385]170  assert(numberRows_ == rhs.numberRows_);
171  assert(numberColumns_ == rhs.numberColumns_);
[1878]172  numberTotal_ = rhs.numberTotal_;
173  maximumNumberTotal_ = rhs.maximumNumberTotal_;
174  // special options here to be safe
[2385]175  specialOptions_ = rhs.specialOptions_;
[1878]176  //assert (maximumInternalRows_ >= numberRows_);
177  //assert (maximumInternalColumns_ >= numberColumns_);
178  // Copy all scalars
[2385]179  CoinAbcMemcpy(reinterpret_cast< int * >(&sumNonBasicCosts_),
180    reinterpret_cast< const int * >(&rhs.sumNonBasicCosts_),
181    (&swappedAlgorithm_ - reinterpret_cast< int * >(&sumNonBasicCosts_)) + 1);
[1878]182  // could add 2 so can go off end
[2385]183  int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_;
184  internalStatus_ = ClpCopyOfArray(rhs.internalStatus_, sizeArray + maximumNumberTotal_);
[1878]185  abcLower_ = ClpCopyOfArray(rhs.abcLower_, sizeArray);
186  abcUpper_ = ClpCopyOfArray(rhs.abcUpper_, sizeArray);
[2385]187  abcCost_ = ClpCopyOfArray(rhs.abcCost_, sizeArray + maximumNumberTotal_);
[1878]188  abcDj_ = ClpCopyOfArray(rhs.abcDj_, sizeArray);
[1990]189
[2385]190  abcSolution_ = ClpCopyOfArray(rhs.abcSolution_, sizeArray + maximumNumberTotal_);
191  abcPerturbation_ = ClpCopyOfArray(rhs.abcPerturbation_, sizeArray);
192  abcPivotVariable_ = ClpCopyOfArray(rhs.abcPivotVariable_, maximumAbcNumberRows_);
[1878]193  //fromExternal_ = ClpCopyOfArray(rhs.fromExternal_,sizeArray);
194  //toExternal_ = ClpCopyOfArray(rhs.toExternal_,sizeArray);
[2385]195  scaleFromExternal_ = ClpCopyOfArray(rhs.scaleFromExternal_, sizeArray);
196  offset_ = ClpCopyOfArray(rhs.offset_, sizeArray);
197  setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
[1878]198  if (rhs.abcMatrix_)
199    abcMatrix_ = new AbcMatrix(*rhs.abcMatrix_);
200  else
201    abcMatrix_ = NULL;
202  for (int i = 0; i < ABC_NUMBER_USEFUL; i++) {
203    usefulArray_[i] = rhs.usefulArray_[i];
204  }
205  if (rhs.abcFactorization_) {
206    setFactorization(*rhs.abcFactorization_);
207  } else {
208    delete abcFactorization_;
209    abcFactorization_ = NULL;
210  }
211#ifdef EARLY_FACTORIZE
212  delete abcEarlyFactorization_;
213  if (rhs.abcEarlyFactorization_) {
214    abcEarlyFactorization_ = new AbcSimplexFactorization(*rhs.abcEarlyFactorization_);
215  } else {
216    abcEarlyFactorization_ = NULL;
217  }
218#endif
219  abcDualRowPivot_ = rhs.abcDualRowPivot_->clone(true);
220  abcDualRowPivot_->setModel(this);
221  abcPrimalColumnPivot_ = rhs.abcPrimalColumnPivot_->clone(true);
222  abcPrimalColumnPivot_->setModel(this);
223  if (rhs.abcBaseModel_) {
224    abcBaseModel_ = new AbcSimplex(*rhs.abcBaseModel_);
225  } else {
226    abcBaseModel_ = NULL;
227  }
228  if (rhs.clpModel_) {
229    clpModel_ = new ClpSimplex(*rhs.clpModel_);
230  } else {
231    clpModel_ = NULL;
232  }
233  abcProgress_ = rhs.abcProgress_;
234  solveType_ = rhs.solveType_;
235}
236// type == 0 nullify, 1 also delete
[2385]237void AbcSimplex::gutsOfDelete(int type)
[1878]238{
239  if (type) {
[2385]240    delete[] abcLower_;
241    delete[] abcUpper_;
242    delete[] abcCost_;
243    delete[] abcDj_;
244    delete[] abcSolution_;
[1878]245    //delete [] fromExternal_ ;
246    //delete [] toExternal_ ;
[2385]247    delete[] scaleFromExternal_;
[1878]248    //delete [] scaleToExternal_ ;
[2385]249    delete[] offset_;
250    delete[] internalStatus_;
251    delete[] abcPerturbation_;
252    delete abcMatrix_;
[1878]253    delete abcFactorization_;
254#ifdef EARLY_FACTORIZE
255    delete abcEarlyFactorization_;
256#endif
[2385]257    delete[] abcPivotVariable_;
[1878]258    delete abcDualRowPivot_;
259    delete abcPrimalColumnPivot_;
260    delete abcBaseModel_;
261    delete clpModel_;
262    delete abcNonLinearCost_;
263  }
[2385]264  CoinAbcMemset0(reinterpret_cast< char * >(&scaleToExternal_),
265    reinterpret_cast< char * >(&usefulArray_[0]) - reinterpret_cast< char * >(&scaleToExternal_));
[1878]266}
[2385]267template < class T >
268T *newArray(T * /*nullArray*/, int size)
[1878]269{
270  if (size) {
[2385]271    T *arrayNew = new T[size];
[1878]272#ifndef NDEBUG
[2385]273    memset(arrayNew, 'A', (size) * sizeof(T));
[1878]274#endif
275    return arrayNew;
276  } else {
277    return NULL;
278  }
279}
[2385]280template < class T >
281T *resizeArray(T *array, int oldSize1, int oldSize2, int newSize1, int newSize2, int extra)
[1878]282{
[2385]283  if ((array || !oldSize1) && (newSize1 != oldSize1 || newSize2 != oldSize2)) {
284    int newTotal = newSize1 + newSize2;
285    int oldTotal = oldSize1 + oldSize2;
286    T *arrayNew;
287    if (newSize1 > oldSize1 || newSize2 > oldSize2) {
288      arrayNew = new T[2 * newTotal + newSize1 + extra];
[1878]289#ifndef NDEBUG
[2385]290      memset(arrayNew, 'A', (2 * newTotal + newSize1 + extra) * sizeof(T));
[1878]291#endif
[2385]292      CoinAbcMemcpy(arrayNew, array, oldSize1);
293      CoinAbcMemcpy(arrayNew + newSize1, array + oldSize1, oldSize2);
[1878]294      // and second half
[2385]295      CoinAbcMemcpy(arrayNew + newTotal, array, oldSize1 + oldTotal);
296      CoinAbcMemcpy(arrayNew + newSize1 + newTotal, array + oldSize1 + oldTotal, oldSize2);
297      delete[] array;
[1878]298    } else {
[2385]299      arrayNew = array;
300      for (int i = 0; i < newSize2; i++)
301        array[newSize1 + i] = array[oldSize1 + i];
302      for (int i = 0; i < newSize1; i++)
303        array[newTotal + i] = array[oldTotal + i];
304      for (int i = 0; i < newSize2; i++)
305        array[newSize1 + newTotal + i] = array[oldSize1 + oldTotal + i];
[1878]306    }
307    return arrayNew;
308  } else {
309    return array;
310  }
311}
312// Initializes arrays
[2385]313void AbcSimplex::gutsOfInitialize(int numberRows, int numberColumns, bool doMore)
[1878]314{
[1910]315#ifdef ABC_INHERIT
[2385]316  abcSimplex_ = this;
[1910]317#endif
[1878]318  // Zero all
[2385]319  CoinAbcMemset0(&sumNonBasicCosts_, (reinterpret_cast< double * >(&usefulArray_[0]) - &sumNonBasicCosts_));
[1878]320  zeroTolerance_ = 1.0e-13;
321  bestObjectiveValue_ = -COIN_DBL_MAX;
322  primalToleranceToGetOptimal_ = -1.0;
323  primalTolerance_ = 1.0e-6;
324  //dualTolerance_ = 1.0e-6;
325  alphaAccuracy_ = -1.0;
326  upperIn_ = -COIN_DBL_MAX;
327  lowerOut_ = -1;
328  valueOut_ = -1;
329  upperOut_ = -1;
330  dualOut_ = -1;
331  acceptablePivot_ = 1.0e-8;
332  //dualBound_=1.0e9;
333  sequenceIn_ = -1;
334  directionIn_ = -1;
335  sequenceOut_ = -1;
336  directionOut_ = -1;
337  pivotRow_ = -1;
338  lastGoodIteration_ = -100;
339  numberPrimalInfeasibilities_ = 100;
[2385]340  numberRefinements_ = 1;
[1878]341  changeMade_ = 1;
342  forceFactorization_ = -1;
[2385]343  if (perturbation_ < 50 || (perturbation_ > 60 && perturbation_ < 100))
[1878]344    perturbation_ = 50;
345  lastBadIteration_ = -999999;
346  lastFlaggedIteration_ = -999999; // doesn't seem to be used
347  firstFree_ = -1;
348  incomingInfeasibility_ = 1.0;
349  allowedInfeasibility_ = 10.0;
350  solveType_ = 1; // say simplex based life form
351  //specialOptions_|=65536;
352  //ClpModel::startPermanentArrays();
[2385]353  maximumInternalRows_ = 0;
354  maximumInternalColumns_ = 0;
355  numberRows_ = numberRows;
356  numberColumns_ = numberColumns;
357  numberTotal_ = numberRows_ + numberColumns_;
358  maximumAbcNumberRows_ = numberRows;
359  maximumAbcNumberColumns_ = numberColumns;
360  maximumNumberTotal_ = numberTotal_;
361  int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_;
[1878]362  if (doMore) {
363    // say Steepest pricing
364    abcDualRowPivot_ = new AbcDualRowSteepest();
365    abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
[2429]366    internalStatus_ = newArray((unsigned char *)NULL,
[2385]367      sizeArray + maximumNumberTotal_);
[2429]368    abcLower_ = newArray((double *)NULL, sizeArray);
369    abcUpper_ = newArray((double *)NULL, sizeArray);
370    abcCost_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_);
371    abcDj_ = newArray((double *)NULL, sizeArray);
372    abcSolution_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_);
373    //fromExternal_ = newArray((int *)NULL,sizeArray);
374    //toExternal_ = newArray((int *)NULL,sizeArray);
375    scaleFromExternal_ = newArray((double *)NULL, sizeArray);
376    offset_ = newArray((double *)NULL, sizeArray);
377    abcPerturbation_ = newArray((double *)NULL, sizeArray);
378    abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_);
[1878]379    // Fill perturbation array
[2385]380    setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
381    fillPerturbation(0, maximumNumberTotal_);
[1878]382  }
383  // get an empty factorization so we can set tolerances etc
384  getEmptyFactorization();
[2385]385  for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
386    usefulArray_[i].reserve(CoinMax(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200), 2 * numberRows_));
[1878]387  //savedStatus_ = internalStatus_+maximumNumberTotal_;
388  //startPermanentArrays();
389}
390// resizes arrays
[2385]391void AbcSimplex::gutsOfResize(int numberRows, int numberColumns)
[1878]392{
393  if (!abcSolution_) {
[2385]394    numberRows_ = 0;
395    numberColumns_ = 0;
396    numberTotal_ = 0;
397    maximumAbcNumberRows_ = 0;
398    maximumAbcNumberColumns_ = 0;
399    maximumNumberTotal_ = 0;
[1878]400  }
[2385]401  if (numberRows == numberRows_ && numberColumns == numberColumns_)
[1878]402    return;
403  // can do on state bit
[2385]404  int newSize1 = CoinMax(numberRows, maximumAbcNumberRows_);
405  if ((stateOfProblem_ & ADD_A_BIT) != 0 && numberRows > maximumAbcNumberRows_)
406    newSize1 = CoinMax(numberRows, maximumAbcNumberRows_ + CoinMin(100, numberRows_ / 10));
407  int newSize2 = CoinMax(numberColumns, maximumAbcNumberColumns_);
408  numberRows_ = numberRows;
409  numberColumns_ = numberColumns;
410  numberTotal_ = numberRows_ + numberColumns_;
[1878]411  //fromExternal_ = resizeArray(fromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1);
412  //toExternal_ = resizeArray(toExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
[2385]413  scaleFromExternal_ = resizeArray(scaleFromExternal_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
[1878]414  //scaleToExternal_ = resizeArray(scaleToExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
[2385]415  internalStatus_ = resizeArray(internalStatus_, maximumAbcNumberRows_,
416    maximumAbcNumberColumns_,
417    newSize1, newSize2, numberTotal_);
418  abcLower_ = resizeArray(abcLower_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
419  abcUpper_ = resizeArray(abcUpper_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
420  abcCost_ = resizeArray(abcCost_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_);
421  abcDj_ = resizeArray(abcDj_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
422  abcSolution_ = resizeArray(abcSolution_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_);
423  abcPerturbation_ = resizeArray(abcPerturbation_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
424  offset_ = resizeArray(offset_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
425  setupPointers(newSize1, newSize2);
[1878]426  // Fill gaps in perturbation array
[2385]427  fillPerturbation(maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_);
428  fillPerturbation(newSize1 + maximumAbcNumberColumns_, newSize2 - maximumAbcNumberColumns_);
[1878]429  // Clean gap
430  //CoinIotaN(fromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
431  //CoinIotaN(toExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
[2385]432  CoinFillN(scaleFromExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0);
433  CoinFillN(scaleToExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0);
434  CoinFillN(internalStatusSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, static_cast< unsigned char >(1));
435  CoinFillN(lowerSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, -COIN_DBL_MAX);
436  CoinFillN(upperSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, COIN_DBL_MAX);
437  CoinFillN(costSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
438  CoinFillN(djSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
439  CoinFillN(solutionSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
440  CoinFillN(offset_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
441  maximumAbcNumberRows_ = newSize1;
442  maximumAbcNumberColumns_ = newSize2;
443  maximumNumberTotal_ = newSize1 + newSize2;
444  delete[] abcPivotVariable_;
[1878]445  abcPivotVariable_ = new int[maximumAbcNumberRows_];
[2385]446  for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
447    usefulArray_[i].reserve(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200));
[1878]448}
[2385]449void AbcSimplex::translate(int type)
[1878]450{
451  // clear bottom bits
[2385]452  stateOfProblem_ &= ~(VALUES_PASS - 1);
453  if ((type & DO_SCALE_AND_MATRIX) != 0) {
[1878]454    //stateOfProblem_ |= DO_SCALE_AND_MATRIX;
455    stateOfProblem_ |= DO_SCALE_AND_MATRIX;
456    delete abcMatrix_;
[2385]457    abcMatrix_ = new AbcMatrix(*matrix());
[1878]458    abcMatrix_->setModel(this);
459    abcMatrix_->scale(scalingFlag_ ? 0 : -1);
460  }
[2385]461  if ((type & DO_STATUS) != 0 && (type & DO_BASIS_AND_ORDER) == 0) {
[1878]462    // from Clp enum to Abc enum (and bound flip)
[2385]463    unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
464    unsigned char *COIN_RESTRICT statusAbc = internalStatus_;
465    const unsigned char *COIN_RESTRICT statusClp = status_;
[1878]466    int i;
[2385]467    for (i = numberRows_ - 1; i >= 0; i--) {
468      unsigned char status = statusClp[i] & 7;
469      bool basicClp = status == 1;
470      bool basicAbc = (statusAbc[i] & 7) == 4;
471      if (basicClp == basicAbc)
472        statusAbc[i] = lookupToAbcSlack[status];
[1878]473      else
[2385]474        break;
[1878]475    }
476    if (!i) {
477      // from Clp enum to Abc enum
[2385]478      unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
479      statusAbc += maximumAbcNumberRows_;
480      statusClp += maximumAbcNumberRows_;
481      for (i = numberColumns_ - 1; i >= 0; i--) {
482        unsigned char status = statusClp[i] & 7;
483        bool basicClp = status == 1;
484        bool basicAbc = (statusAbc[i] & 7) == 4;
485        if (basicClp == basicAbc)
486          statusAbc[i] = lookupToAbc[status];
487        else
488          break;
[1878]489      }
490      if (i)
[2385]491        type |= DO_BASIS_AND_ORDER;
[1878]492    } else {
[2385]493      type |= DO_BASIS_AND_ORDER;
[1878]494    }
495    stateOfProblem_ |= DO_STATUS;
496  }
[2385]497  if ((type & DO_BASIS_AND_ORDER) != 0) {
[1878]498    permuteIn();
499    permuteBasis();
500    stateOfProblem_ |= DO_BASIS_AND_ORDER;
501    type &= ~DO_SOLUTION;
502  }
[2385]503  if ((type & DO_SOLUTION) != 0) {
[1878]504    permuteBasis();
505    stateOfProblem_ |= DO_SOLUTION;
506  }
[2385]507  if ((type & DO_JUST_BOUNDS) != 0) {
[1878]508    stateOfProblem_ |= DO_JUST_BOUNDS;
509  }
510  if (!type) {
511    // just move stuff down
[2385]512    CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_);
513    CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_);
514    CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
[1878]515  }
516}
[2024]517#ifdef ABC_SPRINT
518// Overwrite to create sub problem (just internal arrays) - save full stuff
[2385]519AbcSimplex *
520AbcSimplex::createSubProblem(int numberColumns, const int *whichColumn)
[2024]521{
522  int numberFullColumns = numberColumns_;
523  numberColumns_ = numberColumns;
[2385]524  AbcSimplex *subProblem = this;
525  AbcSimplex *fullProblem = reinterpret_cast< AbcSimplex * >(new char[sizeof(AbcSimplex)]);
526  memset(fullProblem, 0, sizeof(AbcSimplex));
[2024]527  fullProblem->maximumAbcNumberRows_ = maximumAbcNumberRows_;
528  fullProblem->numberColumns_ = numberFullColumns;
529  fullProblem->numberTotal_ = numberTotal_;
530  fullProblem->maximumNumberTotal_ = maximumNumberTotal_;
531  fullProblem->numberTotalWithoutFixed_ = numberTotalWithoutFixed_;
[2385]532  fullProblem->abcPrimalColumnPivot_ = abcPrimalColumnPivot_;
533  fullProblem->internalStatus_ = internalStatus_;
534  fullProblem->abcLower_ = abcLower_;
535  fullProblem->abcUpper_ = abcUpper_;
536  fullProblem->abcCost_ = abcCost_;
537  fullProblem->abcDj_ = abcDj_;
538  fullProblem->abcSolution_ = abcSolution_;
539  fullProblem->scaleFromExternal_ = scaleFromExternal_;
540  fullProblem->offset_ = offset_;
541  fullProblem->abcPerturbation_ = abcPerturbation_;
542  fullProblem->abcPivotVariable_ = abcPivotVariable_;
543  fullProblem->abcMatrix_ = abcMatrix_;
[2024]544  fullProblem->abcNonLinearCost_ = abcNonLinearCost_;
[2385]545  fullProblem->setupPointers(maximumAbcNumberRows_, numberFullColumns);
546  subProblem->numberTotal_ = maximumAbcNumberRows_ + numberColumns;
547  subProblem->maximumNumberTotal_ = maximumAbcNumberRows_ + numberColumns;
548  subProblem->numberTotalWithoutFixed_ = subProblem->numberTotal_;
549  int sizeArray = 2 * subProblem->maximumNumberTotal_ + maximumAbcNumberRows_;
[2429]550  subProblem->internalStatus_ = newArray((unsigned char *)NULL,
[2385]551    sizeArray + subProblem->maximumNumberTotal_);
[2429]552  subProblem->abcLower_ = newArray((double *)NULL, sizeArray);
553  subProblem->abcUpper_ = newArray((double *)NULL, sizeArray);
554  subProblem->abcCost_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_);
555  subProblem->abcDj_ = newArray((double *)NULL, sizeArray);
556  subProblem->abcSolution_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_);
557  //fromExternal_ = newArray((int *)NULL,sizeArray);
558  //toExternal_ = newArray((int *)NULL,sizeArray);
559  subProblem->scaleFromExternal_ = newArray((double *)NULL, sizeArray);
560  subProblem->offset_ = newArray((double *)NULL, sizeArray);
561  subProblem->abcPerturbation_ = newArray((double *)NULL, sizeArray);
562  subProblem->abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_);
[2385]563  subProblem->setupPointers(maximumAbcNumberRows_, numberColumns);
[2024]564  // could use arrays - but for now be safe
[2385]565  int *backward = new int[numberFullColumns + numberRows_];
566  int *whichRow = backward + numberFullColumns;
567  for (int i = 0; i < numberFullColumns; i++)
568    backward[i] = -1;
569  for (int i = 0; i < numberColumns; i++)
570    backward[whichColumn[i]] = i + numberRows_;
571  for (int i = 0; i < numberRows_; i++) {
572    whichRow[i] = i;
[2024]573    int iPivot = fullProblem->abcPivotVariable_[i];
[2385]574    if (iPivot < numberRows_) {
575      subProblem->abcPivotVariable_[i] = iPivot;
[2024]576    } else {
[2385]577      assert(backward[iPivot - numberRows_] >= 0);
578      subProblem->abcPivotVariable_[i] = backward[iPivot - numberRows_];
[2024]579    }
580    subProblem->internalStatus_[i] = fullProblem->internalStatus_[i];
581    subProblem->abcLower_[i] = fullProblem->abcLower_[i];
582    subProblem->abcUpper_[i] = fullProblem->abcUpper_[i];
583    subProblem->abcCost_[i] = fullProblem->abcCost_[i];
584    subProblem->abcDj_[i] = fullProblem->abcDj_[i];
585    subProblem->abcSolution_[i] = fullProblem->abcSolution_[i];
586    subProblem->abcPerturbation_[i] = fullProblem->abcPerturbation_[i];
587    subProblem->internalStatusSaved_[i] = fullProblem->internalStatusSaved_[i];
588    subProblem->perturbationSaved_[i] = fullProblem->perturbationSaved_[i];
589    subProblem->lowerSaved_[i] = fullProblem->lowerSaved_[i];
590    subProblem->upperSaved_[i] = fullProblem->upperSaved_[i];
591    subProblem->costSaved_[i] = fullProblem->costSaved_[i];
592    subProblem->djSaved_[i] = fullProblem->djSaved_[i];
593    subProblem->solutionSaved_[i] = fullProblem->solutionSaved_[i];
594    subProblem->offset_[i] = fullProblem->offset_[i];
595    subProblem->lowerBasic_[i] = fullProblem->lowerBasic_[i];
596    subProblem->upperBasic_[i] = fullProblem->upperBasic_[i];
597    subProblem->costBasic_[i] = fullProblem->costBasic_[i];
598    subProblem->solutionBasic_[i] = fullProblem->solutionBasic_[i];
599    subProblem->djBasic_[i] = fullProblem->djBasic_[i];
600  }
[2385]601  for (int i = 0; i < numberColumns; i++) {
602    int k = whichColumn[i];
603    subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k];
604    subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k];
605    subProblem->abcLower_[maximumAbcNumberRows_ + i] = fullProblem->abcLower_[maximumAbcNumberRows_ + k];
606    subProblem->abcUpper_[maximumAbcNumberRows_ + i] = fullProblem->abcUpper_[maximumAbcNumberRows_ + k];
607    subProblem->abcCost_[maximumAbcNumberRows_ + i] = fullProblem->abcCost_[maximumAbcNumberRows_ + k];
608    subProblem->abcDj_[maximumAbcNumberRows_ + i] = fullProblem->abcDj_[maximumAbcNumberRows_ + k];
609    subProblem->abcSolution_[maximumAbcNumberRows_ + i] = fullProblem->abcSolution_[maximumAbcNumberRows_ + k];
610    subProblem->abcPerturbation_[maximumAbcNumberRows_ + i] = fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k];
611    subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i] = fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k];
612    subProblem->perturbationSaved_[maximumAbcNumberRows_ + i] = fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k];
613    subProblem->lowerSaved_[maximumAbcNumberRows_ + i] = fullProblem->lowerSaved_[maximumAbcNumberRows_ + k];
614    subProblem->upperSaved_[maximumAbcNumberRows_ + i] = fullProblem->upperSaved_[maximumAbcNumberRows_ + k];
615    subProblem->costSaved_[maximumAbcNumberRows_ + i] = fullProblem->costSaved_[maximumAbcNumberRows_ + k];
616    subProblem->djSaved_[maximumAbcNumberRows_ + i] = fullProblem->djSaved_[maximumAbcNumberRows_ + k];
617    subProblem->solutionSaved_[maximumAbcNumberRows_ + i] = fullProblem->solutionSaved_[maximumAbcNumberRows_ + k];
618    subProblem->offset_[maximumAbcNumberRows_ + i] = fullProblem->offset_[maximumAbcNumberRows_ + k];
[2024]619  }
[2385]620  subProblem->abcNonLinearCost_ = new AbcNonLinearCost(subProblem);
[2024]621  subProblem->abcNonLinearCost_->checkInfeasibilities(0.0);
[2385]622  subProblem->abcMatrix_ = new AbcMatrix(*fullProblem->abcMatrix_, numberRows_, whichRow,
623    numberColumns, whichColumn);
[2024]624  subProblem->abcMatrix_->setModel(subProblem);
625  subProblem->abcMatrix_->rebalance();
626  subProblem->abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
[2385]627  subProblem->abcPrimalColumnPivot_->saveWeights(subProblem, 2);
628  delete[] backward;
[2024]629  // swap
630  return fullProblem;
631}
632// Restore stuff from sub problem (and delete sub problem)
[2385]633void AbcSimplex::restoreFromSubProblem(AbcSimplex *fullProblem, const int *whichColumn)
[2024]634{
[2385]635  AbcSimplex *subProblem = this;
636  for (int i = 0; i < numberRows_; i++) {
[2024]637    int iPivot = subProblem->abcPivotVariable_[i];
[2385]638    if (iPivot < numberRows_) {
639      fullProblem->abcPivotVariable_[i] = iPivot;
[2024]640    } else {
[2385]641      fullProblem->abcPivotVariable_[i] = whichColumn[iPivot - numberRows_] + numberRows_;
[2024]642    }
643    fullProblem->internalStatus_[i] = subProblem->internalStatus_[i];
644    fullProblem->abcLower_[i] = subProblem->abcLower_[i];
645    fullProblem->abcUpper_[i] = subProblem->abcUpper_[i];
646    fullProblem->abcCost_[i] = subProblem->abcCost_[i];
647    fullProblem->abcDj_[i] = subProblem->abcDj_[i];
648    fullProblem->abcSolution_[i] = subProblem->abcSolution_[i];
649    fullProblem->abcPerturbation_[i] = subProblem->abcPerturbation_[i];
650    fullProblem->internalStatusSaved_[i] = subProblem->internalStatusSaved_[i];
651    fullProblem->perturbationSaved_[i] = subProblem->perturbationSaved_[i];
652    fullProblem->lowerSaved_[i] = subProblem->lowerSaved_[i];
653    fullProblem->upperSaved_[i] = subProblem->upperSaved_[i];
654    fullProblem->costSaved_[i] = subProblem->costSaved_[i];
655    fullProblem->djSaved_[i] = subProblem->djSaved_[i];
656    fullProblem->solutionSaved_[i] = subProblem->solutionSaved_[i];
657    fullProblem->offset_[i] = subProblem->offset_[i];
658    fullProblem->lowerBasic_[i] = subProblem->lowerBasic_[i];
659    fullProblem->upperBasic_[i] = subProblem->upperBasic_[i];
660    fullProblem->costBasic_[i] = subProblem->costBasic_[i];
661    fullProblem->solutionBasic_[i] = subProblem->solutionBasic_[i];
662    fullProblem->djBasic_[i] = subProblem->djBasic_[i];
663  }
664  int numberColumns = subProblem->numberColumns_;
[2385]665  for (int i = 0; i < numberColumns; i++) {
666    int k = whichColumn[i];
667    fullProblem->internalStatus_[maximumAbcNumberRows_ + k] = subProblem->internalStatus_[maximumAbcNumberRows_ + i];
668    fullProblem->abcLower_[maximumAbcNumberRows_ + k] = subProblem->abcLower_[maximumAbcNumberRows_ + i];
669    fullProblem->abcUpper_[maximumAbcNumberRows_ + k] = subProblem->abcUpper_[maximumAbcNumberRows_ + i];
670    fullProblem->abcCost_[maximumAbcNumberRows_ + k] = subProblem->abcCost_[maximumAbcNumberRows_ + i];
671    fullProblem->abcDj_[maximumAbcNumberRows_ + k] = subProblem->abcDj_[maximumAbcNumberRows_ + i];
672    fullProblem->abcSolution_[maximumAbcNumberRows_ + k] = subProblem->abcSolution_[maximumAbcNumberRows_ + i];
673    fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k] = subProblem->abcPerturbation_[maximumAbcNumberRows_ + i];
674    fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k] = subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i];
675    fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k] = subProblem->perturbationSaved_[maximumAbcNumberRows_ + i];
676    fullProblem->lowerSaved_[maximumAbcNumberRows_ + k] = subProblem->lowerSaved_[maximumAbcNumberRows_ + i];
677    fullProblem->upperSaved_[maximumAbcNumberRows_ + k] = subProblem->upperSaved_[maximumAbcNumberRows_ + i];
678    fullProblem->costSaved_[maximumAbcNumberRows_ + k] = subProblem->costSaved_[maximumAbcNumberRows_ + i];
679    fullProblem->djSaved_[maximumAbcNumberRows_ + k] = subProblem->djSaved_[maximumAbcNumberRows_ + i];
680    fullProblem->solutionSaved_[maximumAbcNumberRows_ + k] = subProblem->solutionSaved_[maximumAbcNumberRows_ + i];
681    fullProblem->offset_[maximumAbcNumberRows_ + k] = subProblem->offset_[maximumAbcNumberRows_ + i];
[2024]682  }
[2385]683  delete[] subProblem->internalStatus_;
684  delete[] subProblem->abcPerturbation_;
[2024]685  delete subProblem->abcMatrix_;
[2385]686  delete[] subProblem->abcLower_;
687  delete[] subProblem->abcUpper_;
688  delete[] subProblem->abcCost_;
689  delete[] subProblem->abcSolution_;
690  delete[] subProblem->abcDj_;
[2024]691  delete subProblem->abcPrimalColumnPivot_;
[2385]692  delete[] subProblem->scaleFromExternal_;
693  delete[] subProblem->offset_;
694  delete[] subProblem->abcPivotVariable_;
695  delete[] subProblem->reversePivotVariable_;
[2024]696  delete subProblem->abcNonLinearCost_;
697  numberColumns_ = fullProblem->numberColumns_;
[2385]698  numberTotal_ = fullProblem->numberTotal_;
699  maximumNumberTotal_ = fullProblem->maximumNumberTotal_;
700  numberTotalWithoutFixed_ = fullProblem->numberTotalWithoutFixed_;
701  abcPrimalColumnPivot_ = fullProblem->abcPrimalColumnPivot_;
702  internalStatus_ = fullProblem->internalStatus_;
703  abcLower_ = fullProblem->abcLower_;
704  abcUpper_ = fullProblem->abcUpper_;
705  abcCost_ = fullProblem->abcCost_;
706  abcDj_ = fullProblem->abcDj_;
707  abcSolution_ = fullProblem->abcSolution_;
708  scaleFromExternal_ = fullProblem->scaleFromExternal_;
709  offset_ = fullProblem->offset_;
710  abcPerturbation_ = fullProblem->abcPerturbation_;
711  abcPivotVariable_ = fullProblem->abcPivotVariable_;
712  abcMatrix_ = fullProblem->abcMatrix_;
713  setupPointers(maximumAbcNumberRows_, numberColumns);
[2024]714  // ? redo nonlinearcost
715  abcNonLinearCost_ = fullProblem->abcNonLinearCost_;
716  abcNonLinearCost_->refresh();
[2385]717  delete[] reinterpret_cast< char * >(fullProblem);
[2024]718}
719#endif
[1878]720/* Sets dual values pass djs using unscaled duals
721   type 1 - values pass
722   type 2 - just use as infeasibility weights
723   type 3 - as 2 but crash
724*/
[2385]725void AbcSimplex::setupDualValuesPass(const double *fakeDuals,
726  const double *fakePrimals,
727  int type)
[1878]728{
729  // allslack basis
[2385]730  memset(internalStatus_, 6, numberRows_);
[1878]731  // temp
[2385]732  if (type == 1) {
733    bool allEqual = true;
734    for (int i = 0; i < numberRows_; i++) {
735      if (rowUpper_[i] > rowLower_[i]) {
736        allEqual = false;
737        break;
[1878]738      }
739    }
740    if (allEqual) {
741      // just modify costs
[2385]742      transposeTimes(-1.0, fakeDuals, objective());
[1878]743      return;
744    }
745  }
746  // compute unscaled djs
[2385]747  CoinAbcMemcpy(djSaved_ + maximumAbcNumberRows_, objective(), numberColumns_);
748  matrix_->transposeTimes(-1.0, fakeDuals, djSaved_ + maximumAbcNumberRows_);
[1878]749  // save fake solution
[2385]750  assert(solution_);
[1878]751  //solution_ = new double[numberTotal_];
[2385]752  CoinAbcMemset0(solution_, numberRows_);
753  CoinAbcMemcpy(solution_ + maximumAbcNumberRows_, fakePrimals, numberColumns_);
[1878]754  // adjust
[2385]755  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++)
756    solution_[iSequence] -= offset_[iSequence];
757  matrix_->times(-1.0, solution_ + maximumAbcNumberRows_, solution_);
[1878]758  double direction = optimizationDirection_;
[2385]759  const double *COIN_RESTRICT rowScale = scaleFromExternal_;
760  const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
761  for (int iRow = 0; iRow < numberRows_; iRow++) {
762    djSaved_[iRow] = direction * fakeDuals[iRow] * inverseRowScale[iRow];
[1878]763    solution_[iRow] *= rowScale[iRow];
764  }
[2385]765  const double *COIN_RESTRICT columnScale = scaleToExternal_;
766  for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++)
767    djSaved_[iColumn] *= direction * columnScale[iColumn];
[1878]768  // Set solution values
[2385]769  const double *lower = abcLower_ + maximumAbcNumberRows_;
770  const double *upper = abcUpper_ + maximumAbcNumberRows_;
771  double *solution = abcSolution_ + maximumAbcNumberRows_;
772  const double *djs = djSaved_ + maximumAbcNumberRows_;
773  const double *inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
774  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
775    double thisValue = fakePrimals[iColumn];
776    double thisLower = columnLower_[iColumn];
777    double thisUpper = columnUpper_[iColumn];
778    double thisDj = djs[iColumn];
779    solution_[iColumn + maximumAbcNumberRows_] = solution_[iColumn + maximumAbcNumberRows_] * inverseColumnScale[iColumn];
780    if (thisLower > -1.0e30) {
781      if (thisUpper < 1.0e30) {
782        double gapUp = thisUpper - thisValue;
783        double gapDown = thisValue - thisLower;
784        bool goUp;
785        if (gapUp > gapDown && thisDj > -dualTolerance_) {
786          goUp = false;
787        } else if (gapUp < gapDown && thisDj < dualTolerance_) {
788          goUp = true;
789        } else {
790          // wild guess
791          double badUp;
792          double badDown;
793          if (gapUp > gapDown) {
794            badUp = gapUp * dualTolerance_;
795            badDown = -gapDown * thisDj;
796          } else {
797            badUp = gapUp * thisDj;
798            badDown = gapDown * dualTolerance_;
799          }
800          goUp = badDown > badUp;
801        }
802        if (goUp) {
803          solution[iColumn] = upper[iColumn];
804          setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound);
805          setStatus(iColumn, ClpSimplex::atUpperBound);
806        } else {
807          solution[iColumn] = lower[iColumn];
808          setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound);
809          setStatus(iColumn, ClpSimplex::atLowerBound);
810        }
[1878]811      } else {
[2385]812        solution[iColumn] = lower[iColumn];
813        setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound);
814        setStatus(iColumn, ClpSimplex::atLowerBound);
[1878]815      }
[2385]816    } else if (thisUpper < 1.0e30) {
817      solution[iColumn] = upper[iColumn];
818      setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound);
819      setStatus(iColumn, ClpSimplex::atUpperBound);
[1878]820    } else {
821      // free
[2385]822      solution[iColumn] = thisValue * inverseColumnScale[iColumn];
823      setInternalStatus(iColumn + maximumAbcNumberRows_, isFree);
824      setStatus(iColumn, ClpSimplex::isFree);
[1878]825    }
826  }
827  switch (type) {
828  case 1:
829    stateOfProblem_ |= VALUES_PASS;
830    break;
831  case 2:
832    stateOfProblem_ |= VALUES_PASS2;
[2385]833    delete[] solution_;
834    solution_ = NULL;
[1878]835    break;
836  case 3:
837    stateOfProblem_ |= VALUES_PASS2;
838    break;
839  }
840}
841//#############################################################################
[2385]842int AbcSimplex::computePrimals(CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
[1878]843{
[2385]844
[1878]845  arrayVector->clear();
846  previousVector->clear();
847  // accumulate non basic stuff
[2385]848  double *COIN_RESTRICT array = arrayVector->denseVector();
849  CoinAbcScatterZeroTo(abcSolution_, abcPivotVariable_, numberRows_);
[1878]850  abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, array);
851#if 0
852  static int xxxxxx=0;
853  if (xxxxxx==0)
854  for (int i=0;i<numberRows_;i++)
855    printf("%d %.18g\n",i,array[i]);
856#endif
857  //arrayVector->scan(0,numberRows_,zeroTolerance_);
858  // Ftran adjusted RHS and iterate to improve accuracy
859  double lastError = COIN_DBL_MAX;
[2385]860  CoinIndexedVector *thisVector = arrayVector;
861  CoinIndexedVector *lastVector = previousVector;
[1878]862  //if (arrayVector->getNumElements())
863#if 0
864  double largest=0.0;
865  int iLargest=-1;
866  for (int i=0;i<numberRows_;i++) {
867    if (fabs(array[i])>largest) {
868      largest=array[i];
869      iLargest=i;
870    }
871  }
872  printf("largest %g at row %d\n",array[iLargest],iLargest);
873#endif
874  abcFactorization_->updateFullColumn(*thisVector);
875#if 0
876  largest=0.0;
877  iLargest=-1;
878  for (int i=0;i<numberRows_;i++) {
879    if (fabs(array[i])>largest) {
880      largest=array[i];
881      iLargest=i;
882    }
883  }
884  printf("after largest %g at row %d\n",array[iLargest],iLargest);
885#endif
886#if 0
887  if (xxxxxx==0)
888  for (int i=0;i<numberRows_;i++)
889    printf("xx %d %.19g\n",i,array[i]);
890  if (numberIterations_>300)
891  exit(0);
892#endif
[2385]893  int numberRefinements = 0;
[1878]894  for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
895    int numberIn = thisVector->getNumElements();
[2385]896    const int *COIN_RESTRICT indexIn = thisVector->getIndices();
897    const double *COIN_RESTRICT arrayIn = thisVector->denseVector();
898    CoinAbcScatterToList(arrayIn, abcSolution_, indexIn, abcPivotVariable_, numberIn);
[1878]899    // check Ax == b  (for all)
900    abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, solutionBasic_);
901#if 0
902    if (xxxxxx==0)
903      for (int i=0;i<numberTotal_;i++)
904        printf("sol %d %.19g\n",i,abcSolution_[i]);
905    if (xxxxxx==0)
906      for (int i=0;i<numberRows_;i++)
907        printf("basic %d %.19g\n",i,solutionBasic_[i]);
908#endif
909    double multiplier = 131072.0;
[2385]910    largestPrimalError_ = CoinAbcMaximumAbsElementAndScale(solutionBasic_, multiplier, numberRows_);
911    double maxValue = 0.0;
912    for (int i = 0; i < numberRows_; i++) {
913      double value = fabs(solutionBasic_[i]);
914      if (value > maxValue) {
[1878]915#if 0
916        if (xxxxxx==0)
917          printf("larger %.19gg at row %d\n",value,i);
918        maxValue=value;
919#endif
920      }
921    }
922    if (largestPrimalError_ >= lastError) {
923      // restore
[2385]924      CoinIndexedVector *temp = thisVector;
[1878]925      thisVector = lastVector;
926      lastVector = temp;
927      //goodSolution = false;
928      break;
929    }
930    if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
931      // try and make better
932      numberRefinements++;
933      // save this
[2385]934      CoinIndexedVector *temp = thisVector;
[1878]935      thisVector = lastVector;
936      lastVector = temp;
[2385]937      int *COIN_RESTRICT indexOut = thisVector->getIndices();
[1878]938      int number = 0;
939      array = thisVector->denseVector();
940      thisVector->clear();
941      for (int iRow = 0; iRow < numberRows_; iRow++) {
[2385]942        double value = solutionBasic_[iRow];
943        if (value) {
944          array[iRow] = value;
945          indexOut[number++] = iRow;
946          solutionBasic_[iRow] = 0.0;
947        }
[1878]948      }
949      thisVector->setNumElements(number);
950      lastError = largestPrimalError_;
951      abcFactorization_->updateFullColumn(*thisVector);
[2385]952      double *previous = lastVector->denseVector();
[1878]953      number = 0;
[2385]954      multiplier = 1.0 / multiplier;
[1878]955      for (int iRow = 0; iRow < numberRows_; iRow++) {
[2385]956        double value = previous[iRow] + multiplier * array[iRow];
957        if (value) {
958          array[iRow] = value;
959          indexOut[number++] = iRow;
960        } else {
961          array[iRow] = 0.0;
962        }
[1878]963      }
964      thisVector->setNumElements(number);
965    } else {
966      break;
967    }
968  }
[2385]969
[1878]970  // solution as accurate as we are going to get
971  //if (!goodSolution) {
[2385]972  CoinAbcMemcpy(solutionBasic_, thisVector->denseVector(), numberRows_);
973  CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
[1878]974  arrayVector->clear();
975  previousVector->clear();
976  return numberRefinements;
977}
978// Moves basic stuff to basic area
[2385]979void AbcSimplex::moveToBasic(int which)
[1878]980{
[2385]981  if ((which & 1) != 0)
982    CoinAbcGatherFrom(abcSolution_, solutionBasic_, abcPivotVariable_, numberRows_);
983  if ((which & 2) != 0)
984    CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_);
985  if ((which & 4) != 0)
986    CoinAbcGatherFrom(abcLower_, lowerBasic_, abcPivotVariable_, numberRows_);
987  if ((which & 8) != 0)
988    CoinAbcGatherFrom(abcUpper_, upperBasic_, abcPivotVariable_, numberRows_);
[1878]989}
990// now dual side
[2385]991int AbcSimplex::computeDuals(double *givenDjs, CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
[1878]992{
[2385]993  double *COIN_RESTRICT array = arrayVector->denseVector();
994  int *COIN_RESTRICT index = arrayVector->getIndices();
[1878]995  int number = 0;
996  if (!givenDjs) {
997    for (int iRow = 0; iRow < numberRows_; iRow++) {
998      double value = costBasic_[iRow];
999      if (value) {
[2385]1000        array[iRow] = value;
1001        index[number++] = iRow;
[1878]1002      }
1003    }
1004  } else {
1005    // dual values pass - djs may not be zero
1006    for (int iRow = 0; iRow < numberRows_; iRow++) {
1007      int iPivot = abcPivotVariable_[iRow];
1008      // make sure zero if done
1009      if (!pivoted(iPivot))
[2385]1010        givenDjs[iPivot] = 0.0;
[1878]1011      double value = abcCost_[iPivot] - givenDjs[iPivot];
1012      if (value) {
[2385]1013        array[iRow] = value;
1014        index[number++] = iRow;
[1878]1015      }
1016    }
1017  }
1018  arrayVector->setNumElements(number);
1019  // Btran basic costs and get as accurate as possible
1020  double lastError = COIN_DBL_MAX;
[2385]1021  CoinIndexedVector *thisVector = arrayVector;
1022  CoinIndexedVector *lastVector = previousVector;
[1878]1023  abcFactorization_->updateFullColumnTranspose(*thisVector);
[2385]1024
1025  int numberRefinements = 0;
1026  for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
[1878]1027    // check basic reduced costs zero
1028    // put reduced cost of basic into djBasic_
[2385]1029    CoinAbcMemcpy(djBasic_, costBasic_, numberRows_);
1030    abcMatrix_->transposeTimesBasic(-1.0, thisVector->denseVector(), djBasic_);
1031    largestDualError_ = CoinAbcMaximumAbsElement(djBasic_, numberRows_);
[1878]1032    if (largestDualError_ >= lastError) {
1033      // restore
[2385]1034      CoinIndexedVector *temp = thisVector;
[1878]1035      thisVector = lastVector;
1036      lastVector = temp;
1037      break;
1038    }
1039    if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
[2385]1040      && !givenDjs) {
[1878]1041      numberRefinements++;
1042      // try and make better
1043      // save this
[2385]1044      CoinIndexedVector *temp = thisVector;
[1878]1045      thisVector = lastVector;
1046      lastVector = temp;
1047      array = thisVector->denseVector();
1048      double multiplier = 131072.0;
1049      //array=djBasic_*multiplier
[2385]1050      CoinAbcMultiplyAdd(djBasic_, numberRows_, multiplier, array, 0.0);
[1878]1051      lastError = largestDualError_;
[2385]1052      abcFactorization_->updateFullColumnTranspose(*thisVector);
[1878]1053#if ABC_DEBUG
1054      thisVector->checkClean();
1055#endif
1056      multiplier = 1.0 / multiplier;
[2385]1057      double *COIN_RESTRICT previous = lastVector->denseVector();
[1878]1058      // array = array*multiplier+previous
[2385]1059      CoinAbcMultiplyAdd(previous, numberRows_, 1.0, array, multiplier);
[1878]1060    } else {
1061      break;
1062    }
1063  }
1064  // now look at dual solution
[2385]1065  CoinAbcMemcpy(dual_, thisVector->denseVector(), numberRows_);
1066  CoinAbcMemset0(thisVector->denseVector(), numberRows_);
[1878]1067  thisVector->setNumElements(0);
1068  if (numberRefinements) {
[2385]1069    CoinAbcMemset0(lastVector->denseVector(), numberRows_);
[1878]1070    lastVector->setNumElements(0);
1071  }
[2385]1072  CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
1073  abcMatrix_->transposeTimesNonBasic(-1.0, dual_, abcDj_);
[1878]1074  // If necessary - override results
1075  if (givenDjs) {
1076    // restore accurate duals
1077    CoinMemcpyN(abcDj_, numberTotal_, givenDjs);
1078  }
1079  //arrayVector->clear();
1080  //previousVector->clear();
1081#if ABC_DEBUG
1082  arrayVector->checkClean();
1083  previousVector->checkClean();
1084#endif
1085  return numberRefinements;
1086}
1087
1088/* Factorizes using current basis.
1089   solveType - 1 iterating, 0 initial, -1 external
1090   - 2 then iterating but can throw out of basis
1091   If 10 added then in primal values pass
1092   Return codes are as from AbcSimplexFactorization unless initial factorization
1093   when total number of singularities is returned.
1094   Special case is numberRows_+1 -> all slack basis.
1095*/
[2385]1096int AbcSimplex::internalFactorize(int solveType)
[1878]1097{
[2385]1098  assert(status_);
1099
[1878]1100  bool valuesPass = false;
1101  if (solveType >= 10) {
1102    valuesPass = true;
1103    solveType -= 10;
1104  }
1105#if 0
1106  // Make sure everything is clean
1107  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1108    if(getInternalStatus(iSequence) == isFixed) {
1109      // double check fixed
1110      assert (abcUpper_[iSequence] == abcLower_[iSequence]);
1111      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<100.0*primalTolerance_);
1112    } else if (getInternalStatus(iSequence) == isFree) {
1113      assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
1114    } else if (getInternalStatus(iSequence) == atLowerBound) {
1115      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
1116    } else if (getInternalStatus(iSequence) == atUpperBound) {
1117      assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
1118    } else if (getInternalStatus(iSequence) == superBasic) {
1119      assert (!valuesPass);
1120    }
1121  }
1122#endif
1123#if 0 //ndef NDEBUG
1124  // Make sure everything is clean
1125  double sumOutside=0.0;
1126  int numberOutside=0;
1127  //double sumOutsideLarge=0.0;
1128  int numberOutsideLarge=0;
1129  double sumInside=0.0;
1130  int numberInside=0;
1131  //double sumInsideLarge=0.0;
1132  int numberInsideLarge=0;
1133  char rowcol[] = {'R', 'C'};
1134  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1135    if(getInternalStatus(iSequence) == isFixed) {
1136      // double check fixed
1137      assert (abcUpper_[iSequence] == abcLower_[iSequence]);
1138      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<primalTolerance_);
1139    } else if (getInternalStatus(iSequence) == isFree) {
1140      assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
1141    } else if (getInternalStatus(iSequence) == atLowerBound) {
1142      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
1143      if (abcSolution_[iSequence]<abcLower_[iSequence]) {
1144        numberOutside++;
[2385]1145#if ABC_NORMAL_DEBUG > 1
[1878]1146        if (handler_->logLevel()==63)
1147          printf("%c%d below by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
1148                 abcLower_[iSequence]-abcSolution_[iSequence]);
1149#endif
1150        sumOutside-=abcSolution_[iSequence]-abcLower_[iSequence];
1151        if (abcSolution_[iSequence]<abcLower_[iSequence]-primalTolerance_)
1152          numberOutsideLarge++;
1153      } else if (abcSolution_[iSequence]>abcLower_[iSequence]) {
1154        numberInside++;
1155        sumInside+=abcSolution_[iSequence]-abcLower_[iSequence];
1156        if (abcSolution_[iSequence]>abcLower_[iSequence]+primalTolerance_)
1157          numberInsideLarge++;
1158      }
1159    } else if (getInternalStatus(iSequence) == atUpperBound) {
1160      assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
1161      if (abcSolution_[iSequence]>abcUpper_[iSequence]) {
1162        numberOutside++;
[2385]1163#if ABC_NORMAL_DEBUG > 0
[1878]1164        if (handler_->logLevel()==63)
1165          printf("%c%d above by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
1166                 -(abcUpper_[iSequence]-abcSolution_[iSequence]));
1167#endif
1168        sumOutside+=abcSolution_[iSequence]-abcUpper_[iSequence];
1169        if (abcSolution_[iSequence]>abcUpper_[iSequence]+primalTolerance_)
1170          numberOutsideLarge++;
1171      } else if (abcSolution_[iSequence]<abcUpper_[iSequence]) {
1172        numberInside++;
1173        sumInside-=abcSolution_[iSequence]-abcUpper_[iSequence];
1174        if (abcSolution_[iSequence]<abcUpper_[iSequence]-primalTolerance_)
1175          numberInsideLarge++;
1176      }
1177    } else if (getInternalStatus(iSequence) == superBasic) {
1178      //assert (!valuesPass);
1179    }
1180  }
[2385]1181#if ABC_NORMAL_DEBUG > 0
[1878]1182  if (numberInside+numberOutside)
1183    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
1184           numberOutside,sumOutside,numberOutsideLarge,
1185           numberInside,sumInside,numberInsideLarge);
1186#endif
1187#endif
[2024]1188  // *** replace below by cleanStatus
[1878]1189  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[2385]1190    AbcSimplex::Status status = getInternalStatus(iSequence);
1191    if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1192      setInternalStatus(iSequence, isFixed);
[1878]1193  }
[2385]1194  if (numberIterations_ == baseIteration_ && !valuesPass) {
[1878]1195    double largeValue = this->largeValue();
[2385]1196    double *COIN_RESTRICT solution = abcSolution_;
[1878]1197    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[2385]1198      AbcSimplex::Status status = getInternalStatus(iSequence);
1199      if (status == superBasic) {
1200        double lower = abcLower_[iSequence];
1201        double upper = abcUpper_[iSequence];
1202        double value = solution[iSequence];
1203        AbcSimplex::Status thisStatus = isFree;
1204        if (lower > -largeValue || upper < largeValue) {
1205          if (lower != upper) {
1206            if (fabs(value - lower) < fabs(value - upper)) {
1207              thisStatus = AbcSimplex::atLowerBound;
1208              solution[iSequence] = lower;
1209            } else {
1210              thisStatus = AbcSimplex::atUpperBound;
1211              solution[iSequence] = upper;
1212            }
1213          } else {
1214            thisStatus = AbcSimplex::isFixed;
1215            solution[iSequence] = upper;
1216          }
1217          setInternalStatus(iSequence, thisStatus);
1218        }
[1878]1219      }
1220    }
1221  }
1222  int status = abcFactorization_->factorize(this, solveType, valuesPass);
1223  if (status) {
1224    handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1225      << status
1226      << CoinMessageEol;
1227    return -1;
1228  } else if (!solveType) {
1229    // Initial basis - return number of singularities
1230    int numberSlacks = 0;
1231    for (int iRow = 0; iRow < numberRows_; iRow++) {
1232      if (getInternalStatus(iRow) == basic)
[2385]1233        numberSlacks++;
[1878]1234    }
1235    status = CoinMax(numberSlacks - numberRows_, 0);
[1990]1236    if (status)
[2385]1237      printf("%d singularities\n", status);
[1878]1238    // special case if all slack
1239    if (numberSlacks == numberRows_) {
1240      status = numberRows_ + 1;
1241    }
1242  }
[2078]1243#ifdef ABC_USE_COIN_FACTORIZATION
1244  // sparse methods
1245  abcFactorization_->goSparse();
1246#endif
[1878]1247  return status;
1248}
[2024]1249// Make sure no superbasic etc
[2385]1250void AbcSimplex::cleanStatus(bool valuesPass)
[2024]1251{
1252  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[2385]1253    AbcSimplex::Status status = getInternalStatus(iSequence);
1254    if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1255      setInternalStatus(iSequence, isFixed);
[2024]1256  }
[2385]1257  if (numberIterations_ == baseIteration_ && !valuesPass) {
[2024]1258    double largeValue = this->largeValue();
[2385]1259    double *COIN_RESTRICT solution = abcSolution_;
[2024]1260    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[2385]1261      AbcSimplex::Status status = getInternalStatus(iSequence);
1262      if (status == superBasic) {
1263        double lower = abcLower_[iSequence];
1264        double upper = abcUpper_[iSequence];
1265        double value = solution[iSequence];
1266        AbcSimplex::Status thisStatus = isFree;
1267        if (lower > -largeValue || upper < largeValue) {
1268          if (lower != upper) {
1269            if (fabs(value - lower) < fabs(value - upper)) {
1270              thisStatus = AbcSimplex::atLowerBound;
1271              solution[iSequence] = lower;
1272            } else {
1273              thisStatus = AbcSimplex::atUpperBound;
1274              solution[iSequence] = upper;
1275            }
1276          } else {
1277            thisStatus = AbcSimplex::isFixed;
1278            solution[iSequence] = upper;
1279          }
1280          setInternalStatus(iSequence, thisStatus);
1281        }
[2024]1282      }
1283    }
1284  }
1285}
[1878]1286// Sets objectiveValue_ from rawObjectiveValue_
[2385]1287void AbcSimplex::setClpSimplexObjectiveValue()
[1878]1288{
[2385]1289  objectiveValue_ = rawObjectiveValue_ / (objectiveScale_ * rhsScale_);
[1878]1290  objectiveValue_ += objectiveOffset_;
1291}
1292/*
1293  This does basis housekeeping and does values for in/out variables.
1294  Can also decide to re-factorize
1295*/
[2385]1296int AbcSimplex::housekeeping()
[1878]1297{
1298#define DELAYED_UPDATE
1299#ifdef DELAYED_UPDATE
[2385]1300  if (algorithm_ < 0) {
[1878]1301    // modify dualout
1302    dualOut_ /= alpha_;
1303    dualOut_ *= -directionOut_;
1304    //double oldValue = valueIn_;
1305    if (directionIn_ == -1) {
1306      // as if from upper bound
1307      valueIn_ = upperIn_ + dualOut_;
1308#if 0 //def ABC_DEBUG
1309      printf("In from upper dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1310             dualOut_,movement_,movementOld,valueIn_,upperIn_+movement_);
1311#endif
1312    } else {
1313      // as if from lower bound
1314      valueIn_ = lowerIn_ + dualOut_;
1315#if 0 //def ABC_DEBUG
1316      printf("In from lower dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1317             dualOut_,movement_,movementOld,valueIn_,lowerIn_+movement_);
1318#endif
1319    }
1320    // outgoing
1321    if (directionOut_ > 0) {
1322      valueOut_ = lowerOut_;
1323    } else {
1324      valueOut_ = upperOut_;
1325    }
[2385]1326#if ABC_NORMAL_DEBUG > 3
[1878]1327    if (rawObjectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1328      printf("obj backwards %g %g\n", rawObjectiveValue_, oldobj);
1329#endif
1330  }
1331#endif
1332#if 0 //ndef NDEBUG
1333  {
1334    int numberFlagged=0;
1335    for (int iRow = 0; iRow < numberRows_; iRow++) {
1336      int iPivot = abcPivotVariable_[iRow];
1337      if (flagged(iPivot))
1338        numberFlagged++;
1339    }
1340    assert (numberFlagged==numberFlagged_);
1341  }
1342#endif
1343  // save value of incoming and outgoing
1344  numberIterations_++;
1345  changeMade_++; // something has happened
1346  // incoming variable
1347  if (handler_->logLevel() > 7) {
1348    //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
1349    handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
1350      << directionOut_
1351      << directionIn_ << theta_
1352      << dualOut_ << dualIn_ << alpha_
1353      << CoinMessageEol;
1354    if (getInternalStatus(sequenceIn_) == isFree) {
1355      handler_->message(CLP_SIMPLEX_FREEIN, messages_)
[2385]1356        << sequenceIn_
1357        << CoinMessageEol;
[1878]1358    }
1359  }
1360  // change of incoming
[2385]1361  char rowcol[] = { 'R', 'C' };
[1878]1362  if (abcUpper_[sequenceIn_] > 1.0e20 && abcLower_[sequenceIn_] < -1.0e20)
1363    progressFlag_ |= 2; // making real progress
1364#ifndef DELAYED_UPDATE
1365  abcSolution_[sequenceIn_] = valueIn_;
1366#endif
1367  if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] < 1.0e-12)
1368    progressFlag_ |= 1; // making real progress
1369  if (sequenceIn_ != sequenceOut_) {
1370    if (alphaAccuracy_ > 0.0) {
1371      double value = fabs(alpha_);
1372      if (value > 1.0)
[2385]1373        alphaAccuracy_ *= value;
[1878]1374      else
[2385]1375        alphaAccuracy_ /= value;
[1878]1376    }
1377    setInternalStatus(sequenceIn_, basic);
1378    if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] > 0) {
1379      // As Nonlinear costs may have moved bounds (to more feasible)
1380      // Redo using value
1381      if (fabs(valueOut_ - abcLower_[sequenceOut_]) < fabs(valueOut_ - abcUpper_[sequenceOut_])) {
[2385]1382        // going to lower
1383        setInternalStatus(sequenceOut_, atLowerBound);
[1878]1384      } else {
[2385]1385        // going to upper
1386        setInternalStatus(sequenceOut_, atUpperBound);
[1878]1387      }
1388    } else {
1389      // fixed
1390      setInternalStatus(sequenceOut_, isFixed);
1391    }
1392#ifndef DELAYED_UPDATE
1393    abcSolution_[sequenceOut_] = valueOut_;
1394#endif
[2385]1395#if PARTITION_ROW_COPY == 1
[1878]1396    // move in row copy
[2385]1397    abcMatrix_->inOutUseful(sequenceIn_, sequenceOut_);
[1878]1398#endif
1399  } else {
1400    //if (objective_->type()<2)
1401    //assert (fabs(theta_)>1.0e-13);
1402    // flip from bound to bound
1403    // As Nonlinear costs may have moved bounds (to more feasible)
1404    // Redo using value
1405    if (fabs(valueIn_ - abcLower_[sequenceIn_]) < fabs(valueIn_ - abcUpper_[sequenceIn_])) {
1406      // as if from upper bound
1407      setInternalStatus(sequenceIn_, atLowerBound);
1408    } else {
1409      // as if from lower bound
1410      setInternalStatus(sequenceIn_, atUpperBound);
1411    }
1412  }
1413  setClpSimplexObjectiveValue();
1414  if (handler_->logLevel() > 7) {
1415    //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
1416    handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
1417      << numberIterations_ << objectiveValue()
1418      << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
1419      << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
1420    handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
1421    handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
1422    handler_->message() << CoinMessageEol;
1423  }
1424#if 0
1425  if (numberIterations_ > 10000)
1426    printf(" it %d %g %c%d %c%d\n"
1427           , numberIterations_, objectiveValue()
1428           , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
1429           , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
1430#endif
1431  if (hitMaximumIterations())
1432    return 2;
1433  // check for small cycles
1434  int in = sequenceIn_;
1435  int out = sequenceOut_;
1436  int cycle = abcProgress_.cycle(in, out,
[2385]1437    directionIn_, directionOut_);
1438  if (cycle > 0) {
1439#if ABC_NORMAL_DEBUG > 0
[1878]1440    if (handler_->logLevel() >= 63)
1441      printf("Cycle of %d\n", cycle);
1442#endif
1443    // reset
1444    abcProgress_.startCheck();
1445    double random = randomNumberGenerator_.randomDouble();
[2385]1446    int extra = static_cast< int >(9.999 * random);
1447    int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
[1878]1448    if (abcFactorization_->pivots() > cycle) {
1449      forceFactorization_ = CoinMax(1, cycle - off[extra]);
1450    } else {
1451      // need to reject something
1452      int iSequence;
[2385]1453      if (algorithm_ < 0) {
1454        iSequence = sequenceIn_;
[1878]1455      } else {
[2385]1456        /* should be better if don't reject incoming
[1878]1457         as it is in basis */
[2385]1458        iSequence = sequenceOut_;
1459        // so can't happen immediately again
1460        sequenceOut_ = -1;
[1878]1461      }
1462      char x = isColumn(iSequence) ? 'C' : 'R';
1463      if (handler_->logLevel() >= 63)
[2385]1464        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1465          << x << sequenceWithin(iSequence)
1466          << CoinMessageEol;
[1878]1467      setFlagged(iSequence);
1468      //printf("flagging %d\n",iSequence);
1469    }
1470    return 1;
1471  }
[2385]1472  int invertNow = 0;
[1878]1473  // only time to re-factorize if one before real time
1474  // this is so user won't be surprised that maximumPivots has exact meaning
1475  int numberPivots = abcFactorization_->pivots();
[2385]1476  if (algorithm_ < 0)
[1878]1477    numberPivots++; // allow for update not done
1478  int maximumPivots = abcFactorization_->maximumPivots();
1479  int numberDense = 0; //abcFactorization_->numberDense();
[2385]1480  bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations());
1481  if (numberPivots == maximumPivots || maximumPivots < 2) {
[1878]1482    // If dense then increase
[2074]1483    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots && false) {
[1878]1484      abcFactorization_->maximumPivots(numberDense);
1485    }
1486    //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots);
[2385]1487#if CLP_FACTORIZATION_NEW_TIMING > 1
[2078]1488    abcFactorization_->statsRefactor('M');
1489#endif
[1878]1490    return 1;
1491  } else if ((abcFactorization_->timeToRefactorize() && !dontInvert)
[2385]1492    || invertNow) {
[1878]1493    //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots);
[2385]1494#if CLP_FACTORIZATION_NEW_TIMING > 1
[2078]1495    abcFactorization_->statsRefactor('T');
1496#endif
[1878]1497    return 1;
1498  } else if (forceFactorization_ > 0 &&
1499#ifndef DELAYED_UPDATE
[2385]1500    abcFactorization_->pivots()
[1878]1501#else
[2385]1502    abcFactorization_->pivots() + 1
1503#endif
1504      >= forceFactorization_) {
[1878]1505    // relax
1506    forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
1507    if (forceFactorization_ > abcFactorization_->maximumPivots())
1508      forceFactorization_ = -1; //off
[2385]1509      //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots);
1510#if CLP_FACTORIZATION_NEW_TIMING > 1
[2078]1511    abcFactorization_->statsRefactor('F');
1512#endif
[1878]1513    return 1;
1514  } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
1515    // A bit worried problem may be cycling - lets factorize at random intervals for a short period
1516    int numberTooManyIterations = numberIterations_ - 10 * (numberRows_ + (numberColumns_ >> 2));
1517    double random = randomNumberGenerator_.randomDouble();
[2385]1518    int window = numberTooManyIterations % 5000;
1519    if (window < 2 * maximumPivots)
1520      random = 0.2 * random + 0.8; // randomly re-factorize but not too soon
[1878]1521    else
[2385]1522      random = 1.0; // switch off if not in window of opportunity
[1878]1523    int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
1524    if (abcFactorization_->pivots() >= random * maxNumber) {
1525      //printf("ZZ cycling a %d\n",numberPivots);
1526      return 1;
[2385]1527    } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1000010 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
[1878]1528      //printf("ZZ cycling b %d\n",numberPivots);
1529      return 1;
1530    } else {
1531      // carry on iterating
1532      return 0;
1533    }
1534  } else {
1535    // carry on iterating
1536    return 0;
1537  }
1538}
[2385]1539// Swaps primal stuff
1540void AbcSimplex::swapPrimalStuff()
[1878]1541{
[2385]1542  if (sequenceOut_ < 0)
[1878]1543    return;
[2385]1544  assert(sequenceIn_ >= 0);
[1878]1545  abcSolution_[sequenceOut_] = valueOut_;
1546  abcSolution_[sequenceIn_] = valueIn_;
[2385]1547  solutionBasic_[pivotRow_] = valueIn_;
1548  if (algorithm_ < 0) {
[1878]1549    // and set bounds correctly
[2385]1550    static_cast< AbcSimplexDual * >(this)->originalBound(sequenceIn_);
1551    static_cast< AbcSimplexDual * >(this)->changeBound(sequenceOut_);
[1878]1552  } else {
[2385]1553#if ABC_NORMAL_DEBUG > 2
1554    if (handler_->logLevel() == 63)
[1878]1555      printf("Code swapPrimalStuff for primal\n");
1556#endif
1557  }
[2385]1558  if (pivotRow_ >= 0) { // may be flip in primal
1559    lowerBasic_[pivotRow_] = abcLower_[sequenceIn_];
1560    upperBasic_[pivotRow_] = abcUpper_[sequenceIn_];
1561    costBasic_[pivotRow_] = abcCost_[sequenceIn_];
[1878]1562    abcPivotVariable_[pivotRow_] = sequenceIn_;
1563  }
1564}
1565// Swaps dual stuff
[2385]1566void AbcSimplex::swapDualStuff(int lastSequenceOut, int lastDirectionOut)
[1878]1567{
1568  // outgoing
1569  // set dj to zero unless values pass
[2385]1570  if (lastSequenceOut >= 0) {
1571    if ((stateOfProblem_ & VALUES_PASS) == 0) {
[1878]1572      if (lastDirectionOut > 0) {
[2385]1573        abcDj_[lastSequenceOut] = theta_;
[1878]1574      } else {
[2385]1575        abcDj_[lastSequenceOut] = -theta_;
[1878]1576      }
1577    } else {
1578      if (lastDirectionOut > 0) {
[2385]1579        abcDj_[lastSequenceOut] -= theta_;
[1878]1580      } else {
[2385]1581        abcDj_[lastSequenceOut] += theta_;
[1878]1582      }
1583    }
1584  }
[2385]1585  if (sequenceIn_ >= 0) {
1586    abcDj_[sequenceIn_] = 0.0;
[1878]1587    //costBasic_[pivotRow_]=abcCost_[sequenceIn_];
1588  }
1589}
[2385]1590static void solveMany(int number, ClpSimplex **simplex)
[1878]1591{
[2385]1592  for (int i = 0; i < number - 1; i++)
[1878]1593    cilk_spawn simplex[i]->dual(0);
[2385]1594  simplex[number - 1]->dual(0);
[1878]1595  cilk_sync;
1596}
[2385]1597void AbcSimplex::crash(int type)
[1878]1598{
1599  int i;
[2385]1600  for (i = 0; i < numberRows_; i++) {
1601    if (getInternalStatus(i) != basic)
[1878]1602      break;
1603  }
[2385]1604  if (i < numberRows_)
[1878]1605    return;
1606  // all slack
[2385]1607  if (type == 3) {
[1878]1608    // decomposition crash
1609    if (!abcMatrix_->gotRowCopy())
1610      abcMatrix_->createRowCopy();
1611    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
[2385]1612    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
1613    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
[1878]1614    //const double * elementByRow = abcMatrix_->rowElements();
[2385]1615    const CoinBigIndex *rowStart = abcMatrix_->rowStart();
1616    const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
1617    const int *column = abcMatrix_->rowColumns();
1618    int *blockStart = usefulArray_[0].getIndices();
1619    int *columnBlock = blockStart + numberRows_;
1620    int *nextColumn = usefulArray_[1].getIndices();
1621    int *blockCount = nextColumn + numberColumns_;
1622    int direction[2] = { -1, 1 };
1623    int bestBreak = -1;
1624    double bestValue = 0.0;
1625    int iPass = 0;
1626    int halfway = (numberRows_ + 1) / 2;
1627    int firstMaster = -1;
1628    int lastMaster = -2;
1629    int numberBlocks = 0;
1630    int largestRows = 0;
1631    int largestColumns = 0;
1632    int numberEmpty = 0;
1633    int numberMaster = 0;
1634    int numberEmptyColumns = 0;
1635    int numberMasterColumns = 0;
1636    while (iPass < 2) {
1637      int increment = direction[iPass];
1638      int start = increment > 0 ? 0 : numberRows_ - 1;
1639      int stop = increment > 0 ? numberRows_ : -1;
1640      numberBlocks = 0;
1641      int thisBestBreak = -1;
1642      double thisBestValue = COIN_DBL_MAX;
1643      int numberRowsDone = 0;
1644      int numberMarkedColumns = 0;
1645      int maximumBlockSize = 0;
1646      for (int i = 0; i < numberRows_; i++) {
1647        blockStart[i] = -1;
1648        blockCount[i] = 0;
[1878]1649      }
[2385]1650      for (int i = 0; i < numberColumns_; i++) {
1651        columnBlock[i] = -1;
1652        nextColumn[i] = -1;
[1878]1653      }
[2385]1654      for (int iRow = start; iRow != stop; iRow += increment) {
1655        int iBlock = -1;
1656        for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1657          int iColumn = column[j] - maximumAbcNumberRows_;
1658          int whichColumnBlock = columnBlock[iColumn];
1659          if (whichColumnBlock >= 0) {
1660            // column marked
1661            if (iBlock < 0) {
1662              // put row in that block
1663              iBlock = whichColumnBlock;
1664            } else if (iBlock != whichColumnBlock) {
1665              // merge
1666              blockCount[iBlock] += blockCount[whichColumnBlock];
1667              blockCount[whichColumnBlock] = 0;
1668              int jColumn = blockStart[whichColumnBlock];
[1878]1669#ifndef NDEBUG
[2385]1670              int iiColumn = iColumn;
[1878]1671#endif
[2385]1672              while (jColumn >= 0) {
1673                assert(columnBlock[jColumn] == whichColumnBlock);
1674                columnBlock[jColumn] = iBlock;
[1878]1675#ifndef NDEBUG
[2385]1676                if (jColumn == iiColumn)
1677                  iiColumn = -1;
[1878]1678#endif
[2385]1679                iColumn = jColumn;
1680                jColumn = nextColumn[jColumn];
1681              }
[1878]1682#ifndef NDEBUG
[2385]1683              assert(iiColumn == -1);
[1878]1684#endif
[2385]1685              nextColumn[iColumn] = blockStart[iBlock];
1686              blockStart[iBlock] = blockStart[whichColumnBlock];
1687              blockStart[whichColumnBlock] = -1;
1688            }
1689          }
1690        }
1691        int n = numberMarkedColumns;
1692        if (iBlock < 0) {
1693          //new block
1694          if (rowEnd[iRow] > rowStart[iRow]) {
1695            numberBlocks++;
1696            iBlock = numberBlocks;
1697            int jColumn = column[rowStart[iRow]] - maximumAbcNumberRows_;
1698            columnBlock[jColumn] = iBlock;
1699            blockStart[iBlock] = jColumn;
1700            numberMarkedColumns++;
1701            for (CoinBigIndex j = rowStart[iRow] + 1; j < rowEnd[iRow]; j++) {
1702              int iColumn = column[j] - maximumAbcNumberRows_;
1703              columnBlock[iColumn] = iBlock;
1704              numberMarkedColumns++;
1705              nextColumn[jColumn] = iColumn;
1706              jColumn = iColumn;
1707            }
1708            blockCount[iBlock] = numberMarkedColumns - n;
1709          } else {
1710            // empty
1711          }
1712        } else {
1713          // put in existing block
1714          int jColumn = blockStart[iBlock];
1715          for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1716            int iColumn = column[j] - maximumAbcNumberRows_;
1717            assert(columnBlock[iColumn] < 0 || columnBlock[iColumn] == iBlock);
1718            if (columnBlock[iColumn] < 0) {
1719              columnBlock[iColumn] = iBlock;
1720              numberMarkedColumns++;
1721              nextColumn[iColumn] = jColumn;
1722              jColumn = iColumn;
1723            }
1724          }
1725          blockStart[iBlock] = jColumn;
1726          blockCount[iBlock] += numberMarkedColumns - n;
1727        }
1728        if (iBlock >= 0)
1729          maximumBlockSize = CoinMax(maximumBlockSize, blockCount[iBlock]);
1730        numberRowsDone++;
1731        if (thisBestValue * numberRowsDone > maximumBlockSize && numberRowsDone > halfway) {
1732          thisBestBreak = iRow;
1733          thisBestValue = static_cast< double >(maximumBlockSize) / static_cast< double >(numberRowsDone);
1734        }
[1878]1735      }
[2385]1736      if (thisBestBreak == stop)
1737        thisBestValue = COIN_DBL_MAX;
[1878]1738      iPass++;
[2385]1739      if (iPass == 1) {
1740        bestBreak = thisBestBreak;
1741        bestValue = thisBestValue;
[1878]1742      } else {
[2385]1743        if (bestValue < thisBestValue) {
1744          firstMaster = 0;
1745          lastMaster = bestBreak;
1746        } else {
1747          firstMaster = thisBestBreak; // ? +1
1748          lastMaster = numberRows_;
1749        }
[1878]1750      }
1751    }
[2385]1752    bool useful = false;
1753    if (firstMaster < lastMaster) {
1754      for (int i = 0; i < numberRows_; i++)
1755        blockStart[i] = -1;
1756      for (int i = firstMaster; i < lastMaster; i++)
1757        blockStart[i] = -2;
1758      for (int i = 0; i < numberColumns_; i++)
1759        columnBlock[i] = -1;
1760      int firstRow = 0;
1761      numberBlocks = -1;
[1878]1762      while (true) {
[2385]1763        for (; firstRow < numberRows_; firstRow++) {
1764          if (blockStart[firstRow] == -1)
1765            break;
1766        }
1767        if (firstRow == numberRows_)
1768          break;
1769        int nRows = 0;
1770        numberBlocks++;
1771        int numberStack = 1;
1772        blockCount[0] = firstRow;
1773        while (numberStack) {
1774          int iRow = blockCount[--numberStack];
1775          for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1776            int iColumn = column[j] - maximumAbcNumberRows_;
1777            int iBlock = columnBlock[iColumn];
1778            if (iBlock < 0) {
1779              columnBlock[iColumn] = numberBlocks;
1780              for (CoinBigIndex k = columnStart[iColumn];
1781                   k < columnStart[iColumn + 1]; k++) {
1782                int jRow = row[k];
1783                int rowBlock = blockStart[jRow];
1784                if (rowBlock == -1) {
1785                  nRows++;
1786                  blockStart[jRow] = numberBlocks;
1787                  blockCount[numberStack++] = jRow;
1788                }
1789              }
1790            }
1791          }
1792        }
1793        if (!nRows) {
1794          // empty!!
1795          numberBlocks--;
1796        }
1797        firstRow++;
[1878]1798      }
[2385]1799      // adjust
[1878]1800      numberBlocks++;
[2385]1801      for (int i = 0; i < numberBlocks; i++) {
1802        blockCount[i] = 0;
1803        nextColumn[i] = 0;
[1878]1804      }
1805      for (int iRow = 0; iRow < numberRows_; iRow++) {
[2385]1806        int iBlock = blockStart[iRow];
1807        if (iBlock >= 0) {
1808          blockCount[iBlock]++;
1809        } else {
1810          if (iBlock == -2)
1811            numberMaster++;
1812          else
1813            numberEmpty++;
1814        }
[1878]1815      }
1816      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
[2385]1817        int iBlock = columnBlock[iColumn];
1818        if (iBlock >= 0) {
1819          nextColumn[iBlock]++;
1820        } else {
1821          if (columnStart[iColumn + 1] > columnStart[iColumn])
1822            numberMasterColumns++;
1823          else
1824            numberEmptyColumns++;
1825        }
[1878]1826      }
[2385]1827      for (int i = 0; i < numberBlocks; i++) {
1828        if (blockCount[i] + nextColumn[i] > largestRows + largestColumns) {
1829          largestRows = blockCount[i];
1830          largestColumns = nextColumn[i];
1831        }
[1878]1832      }
[2385]1833      useful = true;
1834      if (numberMaster > halfway || largestRows * 3 > numberRows_)
1835        useful = false;
[1878]1836    }
1837    if (useful) {
[2385]1838#if ABC_NORMAL_DEBUG > 0
1839      if (logLevel() >= 2)
1840        printf("%d master rows %d <= < %d\n", lastMaster - firstMaster,
1841          firstMaster, lastMaster);
[1878]1842      printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n",
[2385]1843        useful ? "**Useful" : "NoGood",
1844        numberBlocks, largestRows, largestColumns, numberMaster, numberEmpty, numberRows_,
1845        numberMasterColumns, numberEmptyColumns, numberColumns_);
1846      if (logLevel() >= 2) {
1847        for (int i = 0; i < numberBlocks; i++)
1848          printf("Block %d has %d rows and %d columns\n",
1849            i, blockCount[i], nextColumn[i]);
[1878]1850      }
1851#endif
1852#define NUMBER_DW_BLOCKS 20
[2385]1853      int minSize1 = (numberRows_ - numberMaster + NUMBER_DW_BLOCKS - 1) / NUMBER_DW_BLOCKS;
1854      int minSize2 = (numberRows_ - numberMaster + 2 * NUMBER_DW_BLOCKS - 1) / (2 * NUMBER_DW_BLOCKS);
1855      int *backRow = usefulArray_[1].getIndices();
[1878]1856      // first sort
[2385]1857      for (int i = 0; i < numberBlocks; i++) {
1858        backRow[i] = -(3 * blockCount[i] + 0 * nextColumn[i]);
1859        nextColumn[i] = i;
[1878]1860      }
[2385]1861      CoinSort_2(backRow, backRow + numberBlocks, nextColumn);
[1878]1862      // keep if >minSize2 or sum >minSize1;
[2385]1863      int n = 0;
1864      for (int i = 0; i < numberBlocks; i++) {
1865        int originalBlock = nextColumn[i];
1866        if (blockCount[originalBlock] < minSize2)
1867          break;
1868        n++;
[1878]1869      }
[2385]1870      int size = minSize1;
1871      for (int i = n; i < numberBlocks; i++) {
1872        int originalBlock = nextColumn[i];
1873        size -= blockCount[originalBlock];
1874        if (size > 0 && i < numberBlocks - 1) {
1875          blockCount[originalBlock] = -1;
1876        } else {
1877          size = minSize1;
1878          n++;
1879        }
[1878]1880      }
[2385]1881      int n2 = numberBlocks;
1882      numberBlocks = n;
1883      for (int i = n2 - 1; i >= 0; i--) {
1884        int originalBlock = nextColumn[i];
1885        if (blockCount[originalBlock] > 0)
1886          n--;
1887        blockCount[originalBlock] = n;
[1878]1888      }
[2385]1889      assert(!n);
[1878]1890      for (int iRow = 0; iRow < numberRows_; iRow++) {
[2385]1891        int iBlock = blockStart[iRow];
1892        if (iBlock >= 0)
1893          blockStart[iRow] = blockCount[iBlock];
[1878]1894      }
1895      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
[2385]1896        int iBlock = columnBlock[iColumn];
1897        if (iBlock >= 0)
1898          columnBlock[iColumn] = blockCount[iBlock];
[1878]1899      }
1900      // stick to Clp for now
[2385]1901      ClpSimplex **simplex = new ClpSimplex *[numberBlocks];
1902      for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1903        int nRow = 0;
1904        int nColumn = 0;
1905        for (int iRow = 0; iRow < numberRows_; iRow++) {
1906          if (blockStart[iRow] == iBlock)
1907            blockCount[nRow++] = iRow;
1908        }
1909        for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1910          if (columnBlock[iColumn] == iBlock)
1911            nextColumn[nColumn++] = iColumn;
1912        }
1913        simplex[iBlock] = new ClpSimplex(this, nRow, blockCount, nColumn, nextColumn);
1914        simplex[iBlock]->setSpecialOptions(simplex[iBlock]->specialOptions() & (~65536));
1915        if (logLevel() < 2)
1916          simplex[iBlock]->setLogLevel(0);
[1878]1917      }
[2385]1918      solveMany(numberBlocks, simplex);
1919      int numberBasic = numberMaster;
1920      int numberStructurals = 0;
1921      for (int i = 0; i < numberMaster; i++)
1922        abcPivotVariable_[i] = i + firstMaster;
1923      for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1924        int nRow = 0;
1925        int nColumn = 0;
1926        // from Clp enum to Abc enum (and bound flip)
1927        unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
1928        unsigned char *COIN_RESTRICT getStatus = simplex[iBlock]->statusArray() + simplex[iBlock]->numberColumns();
1929        double *COIN_RESTRICT solutionSaved = solutionSaved_;
1930        double *COIN_RESTRICT lowerSaved = lowerSaved_;
1931        double *COIN_RESTRICT upperSaved = upperSaved_;
1932        for (int iRow = 0; iRow < numberRows_; iRow++) {
1933          if (blockStart[iRow] == iBlock) {
1934            unsigned char status = getStatus[nRow++] & 7;
1935            AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbcSlack[status]);
1936            if (status != ClpSimplex::basic) {
1937              double lowerValue = lowerSaved[iRow];
1938              double upperValue = upperSaved[iRow];
1939              if (lowerValue == -COIN_DBL_MAX) {
1940                if (upperValue == COIN_DBL_MAX) {
1941                  // free
1942                  abcStatus = isFree;
1943                } else {
1944                  abcStatus = atUpperBound;
1945                }
1946              } else if (upperValue == COIN_DBL_MAX) {
1947                abcStatus = atLowerBound;
1948              } else if (lowerValue == upperValue) {
1949                abcStatus = isFixed;
1950              }
1951              switch (abcStatus) {
1952              case isFixed:
1953              case atLowerBound:
1954                solutionSaved[iRow] = lowerValue;
1955                break;
1956              case atUpperBound:
1957                solutionSaved[iRow] = upperValue;
1958                break;
1959              default:
1960                break;
1961              }
1962            } else {
1963              // basic
1964              abcPivotVariable_[numberBasic++] = iRow;
1965            }
1966            internalStatus_[iRow] = abcStatus;
1967          }
1968        }
1969        // from Clp enum to Abc enum
1970        unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
1971        unsigned char *COIN_RESTRICT putStatus = internalStatus_ + maximumAbcNumberRows_;
1972        getStatus = simplex[iBlock]->statusArray();
1973        solutionSaved += maximumAbcNumberRows_;
1974        lowerSaved += maximumAbcNumberRows_;
1975        upperSaved += maximumAbcNumberRows_;
1976        int numberSaved = numberBasic;
1977        for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1978          if (columnBlock[iColumn] == iBlock) {
1979            unsigned char status = getStatus[nColumn++] & 7;
1980            AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbc[status]);
1981            if (status != ClpSimplex::basic) {
1982              double lowerValue = lowerSaved[iColumn];
1983              double upperValue = upperSaved[iColumn];
1984              if (lowerValue == -COIN_DBL_MAX) {
1985                if (upperValue == COIN_DBL_MAX) {
1986                  // free
1987                  abcStatus = isFree;
1988                } else {
1989                  abcStatus = atUpperBound;
1990                }
1991              } else if (upperValue == COIN_DBL_MAX) {
1992                abcStatus = atLowerBound;
1993              } else if (lowerValue == upperValue) {
1994                abcStatus = isFixed;
1995              } else if (abcStatus == isFree) {
1996                abcStatus = superBasic;
1997              }
1998              switch (abcStatus) {
1999              case isFixed:
2000              case atLowerBound:
2001                solutionSaved[iColumn] = lowerValue;
2002                break;
2003              case atUpperBound:
2004                solutionSaved[iColumn] = upperValue;
2005                break;
2006              default:
2007                break;
2008              }
2009            } else {
2010              // basic
2011              if (numberBasic < numberRows_)
2012                abcPivotVariable_[numberBasic++] = iColumn + maximumAbcNumberRows_;
2013              else
2014                abcStatus = superBasic;
2015            }
2016            putStatus[iColumn] = abcStatus;
2017          }
2018        }
2019        numberStructurals += numberBasic - numberSaved;
2020        delete simplex[iBlock];
[1878]2021      }
[2385]2022#if ABC_NORMAL_DEBUG > 0
2023      printf("%d structurals put into basis\n", numberStructurals);
[1878]2024#endif
[2385]2025      if (numberBasic < numberRows_) {
2026        for (int iRow = 0; iRow < numberRows_; iRow++) {
2027          AbcSimplex::Status status = getInternalStatus(iRow);
2028          if (status != AbcSimplex::basic) {
2029            abcPivotVariable_[numberBasic++] = iRow;
2030            setInternalStatus(iRow, basic);
2031            if (numberBasic == numberRows_)
2032              break;
2033          }
2034        }
[1878]2035      }
[2385]2036      assert(numberBasic == numberRows_);
[1878]2037#if 0
2038      int n2=0;
2039      for (int i=0;i<numberTotal_;i++) {
2040        if (getInternalStatus(i)==basic)
2041          n2++;
2042      }
2043      assert (n2==numberRows_);
2044      std::sort(abcPivotVariable_,abcPivotVariable_+numberRows_);
2045      n2=-1;
2046      for (int i=0;i<numberRows_;i++) {
2047        assert (abcPivotVariable_[i]>n2);
2048        n2=abcPivotVariable_[i];
2049      }
2050#endif
[2385]2051      delete[] simplex;
[1878]2052      return;
2053    } else {
2054      // try another
[2385]2055      type = 2;
[1878]2056    }
2057  }
[2385]2058  if (type == 1) {
2059    const double *linearObjective = abcCost_ + maximumAbcNumberRows_;
2060    double gamma = 0.0;
2061    int numberTotal = numberRows_ + numberColumns_;
2062    double *q = new double[numberTotal];
2063    double *v = q + numberColumns_;
2064    int *which = new int[numberTotal + 3 * numberRows_];
2065    int *ii = which + numberColumns_;
2066    int *r = ii + numberRows_;
2067    int *pivoted = r + numberRows_;
[1878]2068    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
[2385]2069      gamma = CoinMax(gamma, linearObjective[iColumn]);
[1878]2070    }
2071    for (int iRow = 0; iRow < numberRows_; iRow++) {
2072      double lowerBound = abcLower_[iRow];
2073      double upperBound = abcUpper_[iRow];
[2385]2074      pivoted[iRow] = -1;
2075      ii[iRow] = 0;
2076      r[iRow] = 0;
2077      v[iRow] = COIN_DBL_MAX;
2078      if (lowerBound == upperBound)
2079        continue;
2080      if (lowerBound <= 0.0 && upperBound >= 0.0) {
2081        pivoted[iRow] = iRow;
2082        ii[iRow] = 1;
2083        r[iRow] = 1;
[1878]2084      }
2085    }
[2385]2086    int nPossible = 0;
2087    int lastPossible = 0;
[1878]2088    double cMaxDiv;
2089    if (gamma)
[2385]2090      cMaxDiv = 1.0 / (1000.0 * gamma);
[1878]2091    else
[2385]2092      cMaxDiv = 1.0;
2093    const double *columnLower = abcLower_ + maximumAbcNumberRows_;
2094    const double *columnUpper = abcUpper_ + maximumAbcNumberRows_;
2095    for (int iPass = 0; iPass < 3; iPass++) {
[1878]2096      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
[2385]2097        double lowerBound = columnLower[iColumn];
2098        double upperBound = columnUpper[iColumn];
2099        if (lowerBound == upperBound)
2100          continue;
2101        double qValue;
2102        if (lowerBound > -1.0e20) {
2103          if (upperBound < 1.0e20) {
2104            // both
2105            qValue = lowerBound - upperBound;
2106            if (iPass != 2)
2107              qValue = COIN_DBL_MAX;
2108          } else {
2109            // just lower
2110            qValue = lowerBound;
2111            if (iPass != 1)
2112              qValue = COIN_DBL_MAX;
2113          }
2114        } else {
2115          if (upperBound < 1.0e20) {
2116            // just upper
2117            qValue = -upperBound;
2118            if (iPass != 1)
2119              qValue = COIN_DBL_MAX;
2120          } else {
2121            // free
2122            qValue = 0.0;
2123            if (iPass != 0)
2124              qValue = COIN_DBL_MAX;
2125          }
2126        }
2127        if (qValue != COIN_DBL_MAX) {
2128          which[nPossible] = iColumn;
2129          q[nPossible++] = qValue + linearObjective[iColumn] * cMaxDiv;
2130          ;
2131        }
[1878]2132      }
[2385]2133      CoinSort_2(q + lastPossible, q + nPossible, which + lastPossible);
2134      lastPossible = nPossible;
[1878]2135    }
[2385]2136    const double *element = abcMatrix_->getPackedMatrix()->getElements();
2137    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
[1878]2138    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
[2385]2139    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2140    int nPut = 0;
2141    for (int i = 0; i < nPossible; i++) {
2142      int iColumn = which[i];
2143      double maxAlpha = 0.0;
2144      int kRow = -1;
2145      double alternativeAlpha = 0.0;
2146      int jRow = -1;
2147      bool canTake = true;
2148      for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2149        int iRow = row[j];
2150        double alpha = fabs(element[j]);
2151        if (alpha > 0.01 * v[iRow]) {
2152          canTake = false;
2153        } else if (!ii[iRow] && alpha > alternativeAlpha) {
2154          alternativeAlpha = alpha;
2155          jRow = iRow;
2156        }
2157        if (!r[iRow] && alpha > maxAlpha) {
2158          maxAlpha = alpha;
2159          kRow = iRow;
2160        }
[1878]2161      }
2162      // only really works if scaled
[2385]2163      if (maxAlpha > 0.99) {
2164        pivoted[kRow] = iColumn + maximumAbcNumberRows_;
2165        v[kRow] = maxAlpha;
2166        ii[kRow] = 1;
2167        for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2168          int iRow = row[j];
2169          r[iRow]++;
2170        }
2171        nPut++;
2172      } else if (canTake && jRow >= 0) {
2173        pivoted[jRow] = iColumn + maximumAbcNumberRows_;
2174        v[jRow] = maxAlpha;
2175        ii[jRow] = 1;
2176        for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2177          int iRow = row[j];
2178          r[iRow]++;
2179        }
2180        nPut++;
[1878]2181      }
2182    }
[2385]2183    for (int iRow = 0; iRow < numberRows_; iRow++) {
2184      int iSequence = pivoted[iRow];
2185      if (iSequence >= 0 && iSequence < numberColumns_) {
2186        abcPivotVariable_[iRow] = iSequence;
2187        if (fabs(abcLower_[iRow]) < fabs(abcUpper_[iRow])) {
2188          setInternalStatus(iRow, atLowerBound);
2189          abcSolution_[iRow] = abcLower_[iRow];
2190          solutionSaved_[iRow] = abcLower_[iRow];
2191        } else {
2192          setInternalStatus(iRow, atUpperBound);
2193          abcSolution_[iRow] = abcUpper_[iRow];
2194          solutionSaved_[iRow] = abcUpper_[iRow];
2195        }
2196        setInternalStatus(iSequence, basic);
[1878]2197      }
2198    }
[2385]2199#if ABC_NORMAL_DEBUG > 0
2200    printf("%d put into basis\n", nPut);
[1878]2201#endif
[2385]2202    delete[] q;
2203    delete[] which;
2204  } else if (type == 2) {
[1878]2205    //return;
[2385]2206    int numberBad = 0;
2207    CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
[1878]2208    // Work on savedSolution and current
[2385]2209    int iVector = getAvailableArray();
[1878]2210#define ALLOW_BAD_DJS
2211#ifdef ALLOW_BAD_DJS
[2385]2212    double *modDj = usefulArray_[iVector].denseVector();
[1878]2213#endif
[2385]2214    for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
2215      double dj = abcDj_[iSequence];
2216      modDj[iSequence] = 0.0;
2217      if (getInternalStatus(iSequence) == atLowerBound) {
2218        if (dj < -dualTolerance_) {
2219          double costThru = -dj * (abcUpper_[iSequence] - abcLower_[iSequence]);
2220          if (costThru < dualBound_) {
2221            // flip
2222            setInternalStatus(iSequence, atUpperBound);
2223            solutionSaved_[iSequence] = abcUpper_[iSequence];
2224            abcSolution_[iSequence] = abcUpper_[iSequence];
2225          } else {
2226            numberBad++;
[1878]2227#ifdef ALLOW_BAD_DJS
[2385]2228            modDj[iSequence] = dj;
2229            dj = 0.0;
[1878]2230#else
[2385]2231            break;
[1878]2232#endif
[2385]2233          }
2234        } else {
2235          dj = CoinMax(dj, 0.0);
2236        }
2237      } else if (getInternalStatus(iSequence) == atLowerBound) {
2238        if (dj > dualTolerance_) {
2239          double costThru = dj * (abcUpper_[iSequence] - abcLower_[iSequence]);
2240          if (costThru < dualBound_) {
2241            assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10);
2242            // flip
2243            setInternalStatus(iSequence, atLowerBound);
2244            solutionSaved_[iSequence] = abcLower_[iSequence];
2245            abcSolution_[iSequence] = abcLower_[iSequence];
2246          } else {
2247            numberBad++;
[1878]2248#ifdef ALLOW_BAD_DJS
[2385]2249            modDj[iSequence] = dj;
2250            dj = 0.0;
[1878]2251#else
[2385]2252            break;
[1878]2253#endif
[2385]2254          }
2255        } else {
2256          dj = CoinMin(dj, 0.0);
2257        }
[1878]2258      } else {
[2385]2259        if (fabs(dj) < dualTolerance_) {
2260          dj = 0.0;
2261        } else {
2262          numberBad++;
[1878]2263#ifdef ALLOW_BAD_DJS
[2385]2264          modDj[iSequence] = dj;
2265          dj = 0.0;
[1878]2266#else
[2385]2267          break;
[1878]2268#endif
[2385]2269        }
[1878]2270      }
[2385]2271      abcDj_[iSequence] = dj;
[1878]2272    }
2273#ifndef ALLOW_BAD_DJS
2274    if (numberBad) {
2275      //CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_);
[2385]2276      return;
[1878]2277    }
2278#endif
2279    if (!abcMatrix_->gotRowCopy())
2280      abcMatrix_->createRowCopy();
2281    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
[2385]2282    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_;
2283    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2284    const double *elementByRow = abcMatrix_->rowElements();
2285    const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2286    const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2287    const int *column = abcMatrix_->rowColumns();
2288    CoinAbcMemset0(solutionBasic_, numberRows_);
2289    CoinAbcMemcpy(lowerBasic_, abcLower_, numberRows_);
2290    CoinAbcMemcpy(upperBasic_, abcUpper_, numberRows_);
2291    abcMatrix_->timesIncludingSlacks(-1.0, solutionSaved_, solutionBasic_);
2292    const double multiplier[] = { 1.0, -1.0 };
2293    int nBasic = 0;
2294    int *index = usefulArray_[iVector].getIndices();
2295    int iVector2 = getAvailableArray();
2296    int *index2 = usefulArray_[iVector2].getIndices();
2297    double *sort = usefulArray_[iVector2].denseVector();
2298    int average = CoinMax(5, rowEnd[numberRows_ - 1] / (8 * numberRows_));
2299    int nPossible = 0;
2300    if (numberRows_ > 10000) {
[1878]2301      // allow more
[2385]2302      for (int iRow = 0; iRow < numberRows_; iRow++) {
2303        double solutionValue = solutionBasic_[iRow];
2304        double infeasibility = 0.0;
2305        double lowerValue = lowerBasic_[iRow];
2306        double upperValue = upperBasic_[iRow];
2307        if (solutionValue < lowerValue - primalTolerance_) {
2308          infeasibility = -(lowerValue - solutionValue);
2309        } else if (solutionValue > upperValue + primalTolerance_) {
2310          infeasibility = upperValue - solutionValue;
2311        }
2312        int length = rowEnd[iRow] - rowStart[iRow];
2313        if (infeasibility)
2314          index2[nPossible++] = length;
[1878]2315      }
[2385]2316      std::sort(index2, index2 + nPossible);
[1878]2317      // see how much we need to get numberRows/10 or nPossible/3
[2385]2318      average = CoinMax(average, index2[CoinMin(numberRows_ / 10, nPossible / 3)]);
2319      nPossible = 0;
[1878]2320    }
[2385]2321    for (int iRow = 0; iRow < numberRows_; iRow++) {
2322      double solutionValue = solutionBasic_[iRow];
2323      double infeasibility = 0.0;
2324      double lowerValue = lowerBasic_[iRow];
2325      double upperValue = upperBasic_[iRow];
2326      if (solutionValue < lowerValue - primalTolerance_) {
2327        infeasibility = -(lowerValue - solutionValue);
2328      } else if (solutionValue > upperValue + primalTolerance_) {
2329        infeasibility = upperValue - solutionValue;
[1878]2330      }
[2385]2331      int length = rowEnd[iRow] - rowStart[iRow];
2332      if (infeasibility && length < average) {
2333        index2[nPossible] = iRow;
2334        sort[nPossible++] = 1.0e5 * infeasibility - iRow;
2335        //sort[nPossible++]=1.0e9*length-iRow;//infeasibility;
[1878]2336      }
2337    }
[2385]2338    CoinSort_2(sort, sort + nPossible, index2);
2339    for (int iWhich = 0; iWhich < nPossible; iWhich++) {
[1878]2340      int iRow = index2[iWhich];
[2385]2341      sort[iWhich] = 0.0;
[1878]2342      if (abcDj_[iRow])
[2385]2343        continue; // marked as invalid
2344      double solutionValue = solutionBasic_[iRow];
2345      double infeasibility = 0.0;
2346      double lowerValue = lowerBasic_[iRow];
2347      double upperValue = upperBasic_[iRow];
2348      if (solutionValue < lowerValue - primalTolerance_) {
2349        infeasibility = lowerValue - solutionValue;
2350      } else if (solutionValue > upperValue + primalTolerance_) {
2351        infeasibility = upperValue - solutionValue;
[1878]2352      }
[2385]2353      assert(infeasibility);
2354      double direction = infeasibility > 0 ? 1.0 : -1.0;
[1878]2355      infeasibility *= direction;
[2385]2356      int whichColumn = -1;
2357      double upperTheta = 1.0e30;
2358      int whichColumn2 = -1;
2359      double upperTheta2 = 1.0e30;
2360      double costThru = 0.0;
2361      int nThru = 0;
2362      for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2363        int iSequence = column[j];
2364        assert(iSequence >= maximumAbcNumberRows_);
2365        double dj = abcDj_[iSequence];
2366        double tableauValue = -elementByRow[j] * direction;
2367        unsigned char iStatus = internalStatus_[iSequence] & 7;
2368        if ((iStatus & 4) == 0) {
2369          double mult = multiplier[iStatus];
2370          double alpha = tableauValue * mult;
2371          double oldValue = dj * mult;
2372          assert(oldValue > -1.0e-2);
2373          double value = oldValue - upperTheta * alpha;
2374          if (value < 0.0) {
2375            upperTheta2 = upperTheta;
2376            whichColumn2 = whichColumn;
2377            costThru = alpha * (abcUpper_[iSequence] - abcLower_[iSequence]);
2378            nThru = 0;
2379            upperTheta = oldValue / alpha;
2380            whichColumn = iSequence;
2381          } else if (oldValue - upperTheta2 * alpha < 0.0) {
2382            costThru += alpha * (abcUpper_[iSequence] - abcLower_[iSequence]);
2383            index[nThru++] = iSequence;
2384          }
2385        } else if (iStatus < 6) {
2386          upperTheta = -1.0;
2387          upperTheta2 = elementByRow[j];
2388          whichColumn = iSequence;
2389          break;
2390        }
[1878]2391      }
[2385]2392      if (whichColumn < 0)
2393        continue;
2394      if (upperTheta != -1.0) {
2395        assert(upperTheta >= 0.0);
2396        if (costThru < infeasibility && whichColumn2 >= 0) {
2397          index[nThru++] = whichColumn;
2398          for (int i = 0; i < nThru; i++) {
2399            int iSequence = index[i];
2400            assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10);
2401            // flip and use previous
2402            if (getInternalStatus(iSequence) == atLowerBound) {
2403              setInternalStatus(iSequence, atUpperBound);
2404              solutionSaved_[iSequence] = abcUpper_[iSequence];
2405              abcSolution_[iSequence] = abcUpper_[iSequence];
2406            } else {
2407              setInternalStatus(iSequence, atLowerBound);
2408              solutionSaved_[iSequence] = abcLower_[iSequence];
2409              abcSolution_[iSequence] = abcLower_[iSequence];
2410            }
2411          }
2412          whichColumn = whichColumn2;
2413          upperTheta = upperTheta2;
2414        }
[1878]2415      } else {
[2385]2416        // free coming in
2417        upperTheta = (abcDj_[whichColumn] * direction) / upperTheta2;
[1878]2418      }
2419      // update djs
2420      upperTheta *= -direction;
[2385]2421      for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2422        int iSequence = column[j];
2423        double dj = abcDj_[iSequence];
2424        double tableauValue = elementByRow[j];
2425        dj -= upperTheta * tableauValue;
2426        unsigned char iStatus = internalStatus_[iSequence] & 7;
2427        if ((iStatus & 4) == 0) {
2428          if (!iStatus) {
2429            assert(dj > -1.0e-2);
2430            dj = CoinMax(dj, 0.0);
[1878]2431#ifdef ALLOW_BAD_DJS
[2385]2432            if (numberBad && modDj[iSequence]) {
2433              double bad = modDj[iSequence];
2434              if (dj + bad >= 0.0) {
2435                numberBad--;
2436                modDj[iSequence] = 0.0;
2437                dj += bad;
2438              } else {
2439                modDj[iSequence] += dj;
2440                dj = 0.0;
2441              }
2442            }
[1878]2443#endif
[2385]2444          } else {
2445            assert(dj < 1.0e-2);
2446            dj = CoinMin(dj, 0.0);
[1878]2447#ifdef ALLOW_BAD_DJS
[2385]2448            if (numberBad && modDj[iSequence]) {
2449              double bad = modDj[iSequence];
2450              if (dj + bad <= 0.0) {
2451                numberBad--;
2452                modDj[iSequence] = 0.0;
2453                dj += bad;
2454              } else {
2455                modDj[iSequence] += dj;
2456                dj = 0.0;
2457              }
2458            }
[1878]2459#endif
[2385]2460          }
2461        } else if (iStatus < 6) {
2462          assert(fabs(dj) < 1.0e-4);
2463          dj = 0.0;
2464        }
2465        abcDj_[iSequence] = dj;
[1878]2466      }
2467      // do basis
[2385]2468      if (direction > 0.0) {
2469        if (upperBasic_[iRow] > lowerBasic_[iRow])
2470          setInternalStatus(iRow, atLowerBound);
2471        else
2472          setInternalStatus(iRow, isFixed);
2473        solutionSaved_[iRow] = abcLower_[iRow];
2474        abcSolution_[iRow] = abcLower_[iRow];
[1878]2475      } else {
[2385]2476        if (upperBasic_[iRow] > lowerBasic_[iRow])
2477          setInternalStatus(iRow, atUpperBound);
2478        else
2479          setInternalStatus(iRow, isFixed);
2480        solutionSaved_[iRow] = abcUpper_[iRow];
2481        abcSolution_[iRow] = abcUpper_[iRow];
[1878]2482      }
[2385]2483      setInternalStatus(whichColumn, basic);
2484      abcPivotVariable_[iRow] = whichColumn;
[1878]2485      nBasic++;
2486      // mark rows
[2385]2487      for (CoinBigIndex j = columnStart[whichColumn]; j < columnStart[whichColumn + 1]; j++) {
2488        int jRow = row[j];
2489        abcDj_[jRow] = 1.0;
[1878]2490      }
2491    }
2492#ifdef ALLOW_BAD_DJS
[2385]2493    CoinAbcMemset0(modDj + maximumAbcNumberRows_, numberColumns_);
[1878]2494#endif
2495    setAvailableArray(iVector);
2496    setAvailableArray(iVector2);
[2385]2497#if ABC_NORMAL_DEBUG > 0
2498    printf("dual crash put %d in basis\n", nBasic);
[1878]2499#endif
2500  } else {
[2385]2501    assert((stateOfProblem_ & VALUES_PASS2) != 0);
[1878]2502    // The idea is to put as many likely variables into basis as possible
[2385]2503    int n = 0;
2504    int iVector = getAvailableArray();
2505    int *index = usefulArray_[iVector].getIndices();
2506    double *array = usefulArray_[iVector].denseVector();
2507    int iVector2 = getAvailableArray();
2508    int *index2 = usefulArray_[iVector].getIndices();
2509    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[1878]2510      double dj = djSaved_[iSequence];
2511      double value = solution_[iSequence];
2512      double lower = abcLower_[iSequence];
2513      double upper = abcUpper_[iSequence];
[2385]2514      double gapUp = CoinMin(1.0e3, upper - value);
2515      assert(gapUp >= -1.0e-3);
2516      gapUp = CoinMax(gapUp, 0.0);
2517      double gapDown = CoinMin(1.0e3, value - lower);
2518      assert(gapDown >= -1.0e-3);
2519      gapDown = CoinMax(gapDown, 0.0);
2520      double measure = (CoinMin(gapUp, gapDown) + 1.0e-6) / (fabs(dj) + 1.0e-6);
2521      if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) {
2522        // set to ub
2523        setInternalStatus(iSequence, atUpperBound);
2524        solutionSaved_[iSequence] = abcUpper_[iSequence];
2525        abcSolution_[iSequence] = abcUpper_[iSequence];
2526      } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) {
2527        // set to lb
2528        setInternalStatus(iSequence, atLowerBound);
2529        solutionSaved_[iSequence] = abcLower_[iSequence];
2530        abcSolution_[iSequence] = abcLower_[iSequence];
2531      } else if (upper > lower) {
2532        // set to nearest
2533        if (gapUp < gapDown) {
2534          // set to ub
2535          setInternalStatus(iSequence, atUpperBound);
2536          solutionSaved_[iSequence] = abcUpper_[iSequence];
2537          abcSolution_[iSequence] = abcUpper_[iSequence];
2538        } else {
2539          // set to lb
2540          setInternalStatus(iSequence, atLowerBound);
2541          solutionSaved_[iSequence] = abcLower_[iSequence];
2542          abcSolution_[iSequence] = abcLower_[iSequence];
2543        }
2544        array[n] = -measure;
2545        index[n++] = iSequence;
[1878]2546      }
2547    }
2548    // set slacks basic
[2385]2549    memset(internalStatus_, 6, numberRows_);
2550    CoinSort_2(solution_, solution_ + n, index);
2551    CoinAbcMemset0(array, n);
2552    for (int i = 0; i < numberRows_; i++)
2553      index2[i] = 0;
[1878]2554    // first put all possible slacks in
[2385]2555    int n2 = 0;
2556    for (int i = 0; i < n; i++) {
2557      int iSequence = index[i];
2558      if (iSequence < numberRows_) {
2559        index2[iSequence] = numberRows_;
[1878]2560      } else {
[2385]2561        index[n2++] = iSequence;
[1878]2562      }
2563    }
[2385]2564    n = n2;
2565    int numberIn = 0;
2566    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_;
[1878]2567    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
[2385]2568    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
[1878]2569    if (!abcMatrix_->gotRowCopy())
2570      abcMatrix_->createRowCopy();
2571    //const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2572    //const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
[2385]2573    for (int i = 0; i < n; i++) {
2574      int iSequence = index[i];
2575      int bestRow = -1;
2576      int bestCount = numberRows_ + 1;
2577      for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) {
2578        int iRow = row[j];
2579        if (!abcMatrix_->gotRowCopy())
2580          abcMatrix_->createRowCopy();
2581        const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2582        const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2583        if (!index2[iRow]) {
2584          int length = rowEnd[iRow] - rowStart[iRow];
2585          if (length < bestCount) {
2586            bestCount = length;
2587            bestRow = iRow;
2588          }
2589        }
[1878]2590      }
[2385]2591      if (bestRow >= 0) {
2592        numberIn++;
2593        for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) {
2594          int iRow = row[j];
2595          index2[iRow]++;
2596        }
2597        setInternalStatus(iSequence, basic);
2598        abcPivotVariable_[bestRow] = iSequence;
2599        double dj = djSaved_[bestRow];
2600        double value = solution_[bestRow];
2601        double lower = abcLower_[bestRow];
2602        double upper = abcUpper_[bestRow];
2603        double gapUp = CoinMax(CoinMin(1.0e3, upper - value), 0.0);
2604        double gapDown = CoinMax(CoinMin(1.0e3, value - lower), 0.0);
2605        //double measure = (CoinMin(gapUp,gapDown)+1.0e-6)/(fabs(dj)+1.0e-6);
2606        if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) {
2607          // set to ub
2608          setInternalStatus(bestRow, atUpperBound);
2609          solutionSaved_[bestRow] = abcUpper_[bestRow];
2610          abcSolution_[bestRow] = abcUpper_[bestRow];
2611        } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) {
2612          // set to lb
2613          setInternalStatus(bestRow, atLowerBound);
2614          solutionSaved_[bestRow] = abcLower_[bestRow];
2615          abcSolution_[bestRow] = abcLower_[bestRow];
2616        } else if (upper > lower) {
2617          // set to nearest
2618          if (gapUp < gapDown) {
2619            // set to ub
2620            setInternalStatus(bestRow, atUpperBound);
2621            solutionSaved_[bestRow] = abcUpper_[bestRow];
2622            abcSolution_[bestRow] = abcUpper_[bestRow];
2623          } else {
2624            // set to lb
2625            setInternalStatus(bestRow, atLowerBound);
2626            solutionSaved_[bestRow] = abcLower_[bestRow];
2627            abcSolution_[bestRow] = abcLower_[bestRow];
2628          }
2629        }
[1878]2630      }
2631    }
[2385]2632#if ABC_NORMAL_DEBUG > 0
2633    printf("Idiot crash put %d in basis\n", numberIn);
[1878]2634#endif
2635    setAvailableArray(iVector);
2636    setAvailableArray(iVector2);
[2385]2637    delete[] solution_;
2638    solution_ = NULL;
[1878]2639  }
2640}
2641/* Puts more stuff in basis
2642   1 bit set - do even if basis exists
2643   2 bit set - don't bother staying triangular
2644*/
[2385]2645void AbcSimplex::putStuffInBasis(int type)
[1878]2646{
2647  int i;
[2385]2648  for (i = 0; i < numberRows_; i++) {
2649    if (getInternalStatus(i) != basic)
[1878]2650      break;
2651  }
[2385]2652  if (i < numberRows_ && (type & 1) == 0)
[1878]2653    return;
[2385]2654  int iVector = getAvailableArray();
[1878]2655  // Column part is mustColumnIn
[2385]2656  int *COIN_RESTRICT mustRowOut = usefulArray_[iVector].getIndices();
[1878]2657  if (!abcMatrix_->gotRowCopy())
2658    abcMatrix_->createRowCopy();
[2385]2659  const double *elementByRow = abcMatrix_->rowElements();
2660  const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2661  const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2662  const int *column = abcMatrix_->rowColumns();
2663  for (int i = 0; i < numberTotal_; i++)
2664    mustRowOut[i] = -1;
2665  int nPossible = 0;
[1878]2666  // find equality rows with effective nonzero rhs
[2385]2667  for (int iRow = 0; iRow < numberRows_; iRow++) {
2668    if (abcUpper_[iRow] > abcLower_[iRow] || getInternalStatus(iRow) != basic) {
2669      mustRowOut[iRow] = -2;
[1878]2670      continue;
2671    }
[2385]2672    int chooseColumn[2] = { -1, -1 };
2673    for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2674      int iColumn = column[j];
2675      if (elementByRow[j] > 0.0) {
2676        if (chooseColumn[0] == -1)
2677          chooseColumn[0] = iColumn;
2678        else
2679          chooseColumn[0] = -2;
[1878]2680      } else {
[2385]2681        if (chooseColumn[1] == -1)
2682          chooseColumn[1] = iColumn;
2683        else
2684          chooseColumn[1] = -2;
[1878]2685      }
2686    }
[2385]2687    for (int iTry = 0; iTry < 2; iTry++) {
2688      int jColumn = chooseColumn[iTry];
2689      if (jColumn >= 0 && getInternalStatus(jColumn) != basic) {
2690        // see if has to be basic
2691        double lowerValue = -abcUpper_[iRow]; // check swap
2692        double upperValue = -abcLower_[iRow];
2693        int lowerInf = 0;
2694        int upperInf = 0;
2695        double alpha = 0.0;
2696        for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2697          int iColumn = column[j];
2698          if (iColumn != jColumn) {
2699            if (abcLower_[iColumn] > -1.0e20)
2700              lowerValue -= abcLower_[iColumn] * elementByRow[j];
2701            else
2702              lowerInf++;
2703            if (abcUpper_[iColumn] < 1.0e20)
2704              upperValue -= abcUpper_[iColumn] * elementByRow[j];
2705            else
2706              upperInf++;
2707          } else {
2708            alpha = elementByRow[j];
2709          }
2710        }
2711        // find values column must lie between (signs again)
2712        if (upperInf)
2713          upperValue = COIN_DBL_MAX;
2714        else
2715          upperValue /= alpha;
2716        if (lowerInf)
2717          lowerValue = -COIN_DBL_MAX;
2718        else
2719          lowerValue /= alpha;
2720        if (iTry) {
2721          // swap
2722          double temp = lowerValue;
2723          lowerValue = upperValue;
2724          upperValue = temp;
2725        }
2726        if (lowerValue > abcLower_[jColumn] + 10.0 * primalTolerance_ && upperValue < abcUpper_[jColumn] - 10.0 * primalTolerance_) {
2727          nPossible++;
2728          if (mustRowOut[jColumn] >= 0) {
2729            // choose one ???
2730            //printf("Column %d already selected on row %d now on %d\n",
2731            //     jColumn,mustRowOut[jColumn],iRow);
2732            continue;
2733          }
2734          mustRowOut[jColumn] = iRow;
2735          mustRowOut[iRow] = jColumn;
2736        }
[1878]2737      }
2738    }
2739  }
2740  if (nPossible) {
[2385]2741#if ABC_NORMAL_DEBUG > 0
2742    printf("%d possible candidates\n", nPossible);
[1878]2743#endif
[2385]2744    if ((type & 2) == 0) {
[1878]2745      // triangular
[2385]2746      int iVector2 = getAvailableArray();
2747      int *COIN_RESTRICT counts = usefulArray_[iVector2].getIndices();
2748      CoinAbcMemset0(counts, numberRows_);
2749      for (int iRow = 0; iRow < numberRows_; iRow++) {
2750        int n = 0;
2751        for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2752          int iColumn = column[j];
2753          if (getInternalStatus(iColumn) == basic)
2754            n++;
2755        }
2756        counts[iRow] = n;
[1878]2757      }
[2385]2758      const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts()
2759        - maximumAbcNumberRows_;
2760      const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2761      for (int iRow = 0; iRow < numberRows_; iRow++) {
2762        if (!counts[iRow]) {
2763          int iColumn = mustRowOut[iRow];
2764          if (iColumn >= 0) {
2765            setInternalStatus(iRow, isFixed);
2766            solutionSaved_[iRow] = abcLower_[iRow];
2767            setInternalStatus(iColumn, basic);
2768            for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2769              int jRow = row[j];
2770              counts[jRow]++;
2771            }
2772          }
2773        }
[1878]2774      }
2775      setAvailableArray(iVector2);
2776    } else {
2777      // all
[2385]2778      for (int iRow = 0; iRow < numberRows_; iRow++) {
2779        int iColumn = mustRowOut[iRow];
2780        if (iColumn >= 0) {
2781          setInternalStatus(iRow, isFixed);
2782          solutionSaved_[iRow] = abcLower_[iRow];
2783          setInternalStatus(iColumn, basic);
2784        }
[1878]2785      }
2786    }
2787    // redo pivot array
[2385]2788    int numberBasic = 0;
2789    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2790      if (getInternalStatus(iSequence) == basic)
2791        abcPivotVariable_[numberBasic++] = iSequence;
[1878]2792    }
[2385]2793    assert(numberBasic == numberRows_);
[1878]2794  }
2795  setAvailableArray(iVector);
2796}
2797// Computes nonbasic cost and total cost
[2385]2798void AbcSimplex::computeObjective()
[1878]2799{
2800  sumNonBasicCosts_ = 0.0;
2801  rawObjectiveValue_ = 0.0;
2802  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
[2385]2803    double value = abcSolution_[iSequence] * abcCost_[iSequence];
[1878]2804    rawObjectiveValue_ += value;
[2385]2805    if (getInternalStatus(iSequence) != basic)
[1878]2806      sumNonBasicCosts_ += value;
2807  }
2808  setClpSimplexObjectiveValue();
2809}
2810// This sets largest infeasibility and most infeasible
[2385]2811void AbcSimplex::checkPrimalSolution(bool justBasic)
[1878]2812{
[2385]2813  rawObjectiveValue_ = CoinAbcInnerProduct(costBasic_, numberRows_, solutionBasic_);
2814  //rawObjectiveValue_ += sumNonBasicCosts_;
[1878]2815  setClpSimplexObjectiveValue();
2816  // now look at primal solution
2817  sumPrimalInfeasibilities_ = 0.0;
2818  numberPrimalInfeasibilities_ = 0;
2819  double primalTolerance = primalTolerance_;
2820  double relaxedTolerance = primalTolerance_;
2821  // we can't really trust infeasibilities if there is primal error
[2385]2822  double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 5.0 * primalTolerance_));
[1878]2823  // allow tolerance at least slightly bigger than standard
[2385]2824  relaxedTolerance = relaxedTolerance + error;
[1878]2825  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
[2385]2826  const double *COIN_RESTRICT lowerBasic = lowerBasic_;
2827  const double *COIN_RESTRICT upperBasic = upperBasic_;
2828  const double *COIN_RESTRICT solutionBasic = solutionBasic_;
[1878]2829  if (justBasic) {
2830    for (int iRow = 0; iRow < numberRows_; iRow++) {
2831      double infeasibility = 0.0;
2832      if (solutionBasic[iRow] > upperBasic[iRow]) {
[2385]2833        infeasibility = solutionBasic[iRow] - upperBasic[iRow];
[1878]2834      } else if (solutionBasic[iRow] < lowerBasic[iRow]) {
[2385]2835        infeasibility = lowerBasic[iRow] - solutionBasic[iRow];
[1878]2836      }
2837      if (infeasibility > primalTolerance) {
[2385]2838        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2839        if (infeasibility > relaxedTolerance)
2840          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2841        numberPrimalInfeasibilities_++;
[1878]2842      }
2843    }
2844  } else {
2845    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2846      double infeasibility = 0.0;
2847      if (abcSolution_[iSequence] > abcUpper_[iSequence]) {
[2385]2848        infeasibility = abcSolution_[iSequence] - abcUpper_[iSequence];
[1878]2849      } else if (abcSolution_[iSequence] < abcLower_[iSequence]) {
[2385]2850        infeasibility = abcLower_[iSequence] - abcSolution_[iSequence];
[1878]2851      }
2852      if (infeasibility > primalTolerance) {
[2385]2853        //assert (getInternalStatus(iSequence)==basic);
2854        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2855        if (infeasibility > relaxedTolerance)
2856          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2857        numberPrimalInfeasibilities_++;
[1878]2858      }
2859    }
2860  }
2861}
[2385]2862void AbcSimplex::checkDualSolution()
[1878]2863{
[2385]2864
[1878]2865  sumDualInfeasibilities_ = 0.0;
2866  numberDualInfeasibilities_ = 0;
2867  int firstFreePrimal = -1;
2868  int firstFreeDual = -1;
2869  int numberSuperBasicWithDj = 0;
2870  bestPossibleImprovement_ = 0.0;
2871  // we can't really trust infeasibilities if there is dual error
[2385]2872  double error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_));
[1878]2873  // allow tolerance at least slightly bigger than standard
[2385]2874  double relaxedTolerance = currentDualTolerance_ + error;
[1878]2875  // allow bigger tolerance for possible improvement
2876  double possTolerance = 5.0 * relaxedTolerance;
2877  sumOfRelaxedDualInfeasibilities_ = 0.0;
[2385]2878  int numberNeedToMove = 0;
[1878]2879  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2880    if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
2881      // not basic
2882      double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
2883      double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
2884      double value = abcDj_[iSequence];
2885      if (distanceUp > primalTolerance_) {
[2385]2886        // Check if "free"
2887        if (distanceDown <= primalTolerance_) {
2888          // should not be negative
2889          if (value < 0.0) {
2890            value = -value;
2891            if (value > currentDualTolerance_) {
2892              sumDualInfeasibilities_ += value - currentDualTolerance_;
2893              if (value > possTolerance)
2894                bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2895              if (value > relaxedTolerance)
2896                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2897              numberDualInfeasibilities_++;
2898            }
2899          }
2900        } else {
2901          // free
2902          value = fabs(value);
2903          if (value > 1.0 * relaxedTolerance) {
2904            numberSuperBasicWithDj++;
2905            if (firstFreeDual < 0)
2906              firstFreeDual = iSequence;
2907          }
2908          if (value > currentDualTolerance_ || getInternalStatus(iSequence) != AbcSimplex::isFree) {
2909            numberNeedToMove++;
2910            if (firstFreePrimal < 0)
2911              firstFreePrimal = iSequence;
2912            sumDualInfeasibilities_ += value - currentDualTolerance_;
2913            if (value > possTolerance)
2914              bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2915            if (value > relaxedTolerance)
2916              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2917            numberDualInfeasibilities_++;
2918          }
2919        }
[1878]2920      } else if (distanceDown > primalTolerance_) {
[2385]2921        // should not be positive
2922        if (value > 0.0) {
2923          if (value > currentDualTolerance_) {
2924            sumDualInfeasibilities_ += value - currentDualTolerance_;
2925            if (value > possTolerance)
2926              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2927            if (value > relaxedTolerance)
2928              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2929            numberDualInfeasibilities_++;
2930          }
2931        }
[1878]2932      }
2933    }
2934  }
[2385]2935  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj;
[1878]2936  if (algorithm_ < 0 && firstFreeDual >= 0) {
2937    // dual
2938    firstFree_ = firstFreeDual;
[2385]2939  } else if (numberSuperBasicWithDj || numberNeedToMove) {
[1878]2940    //(abcProgress_.lastIterationNumber(0) <= 0)) {
2941    firstFree_ = firstFreePrimal;
2942    if (!sumDualInfeasibilities_)
[2385]2943      sumDualInfeasibilities_ = 1.0e-8;
[1878]2944  }
2945}
2946/* This sets sum and number of infeasibilities (Dual and Primal) */
[2385]2947void AbcSimplex::checkBothSolutions()
[1878]2948{
[2385]2949  // old way
[1878]2950  checkPrimalSolution(true);
2951  checkDualSolution();
2952}
2953/*
2954  Unpacks one column of the matrix into indexed array
2955*/
[2385]2956void AbcSimplex::unpack(CoinIndexedVector &rowArray, int sequence) const
2957{
2958  abcMatrix_->unpack(rowArray, sequence);
[1878]2959}
2960// Sets row pivot choice algorithm in dual
[2385]2961void AbcSimplex::setDualRowPivotAlgorithm(AbcDualRowPivot &choice)
[1878]2962{
2963  delete abcDualRowPivot_;
2964  abcDualRowPivot_ = choice.clone(true);
2965  abcDualRowPivot_->setModel(this);
2966}
2967// Sets column pivot choice algorithm in primal
[2385]2968void AbcSimplex::setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot &choice)
[1878]2969{
2970  delete abcPrimalColumnPivot_;
2971  abcPrimalColumnPivot_ = choice.clone(true);
2972  abcPrimalColumnPivot_->setModel(this);
2973}
[2385]2974void AbcSimplex::setFactorization(AbcSimplexFactorization &factorization)
[1878]2975{
2976  if (abcFactorization_)
2977    abcFactorization_->setFactorization(factorization);
2978  else
2979    abcFactorization_ = new AbcSimplexFactorization(factorization,
[2385]2980      numberRows_);
[1878]2981}
2982
2983// Swaps factorization
2984AbcSimplexFactorization *
[2385]2985AbcSimplex::swapFactorization(AbcSimplexFactorization *factorization)
[1878]2986{
[2385]2987  AbcSimplexFactorization *swap = abcFactorization_;
[1878]2988  abcFactorization_ = factorization;
2989  return swap;
2990}
2991
2992/* Tightens primal bounds to make dual faster.  Unless
2993   fixed, bounds are slightly looser than they could be.
2994   This is to make dual go faster and is probably not needed
2995   with a presolve.  Returns non-zero if problem infeasible
2996*/
[2385]2997int AbcSimplex::tightenPrimalBounds()
[1878]2998{
[2385]2999  int tightenType = 1;
3000  if (maximumIterations() >= 1000000 && maximumIterations() < 1000010)
3001    tightenType = maximumIterations() - 1000000;
[1878]3002  if (!tightenType)
3003    return 0;
3004  if (integerType_) {
[2385]3005#if ABC_NORMAL_DEBUG > 0
[1878]3006    printf("Redo tighten to relax by 1 for integers (and don't be shocked by infeasibility)\n");
3007#endif
3008    return 0;
3009  }
[2385]3010  // This needs to work on scaled matrix - then replace
[1878]3011  // Get a row copy in standard format
3012  CoinPackedMatrix copy;
3013  copy.setExtraGap(0.0);
3014  copy.setExtraMajor(0.0);
3015  copy.reverseOrderedCopyOf(*abcMatrix_->matrix());
3016  // get matrix data pointers
[2385]3017  const int *COIN_RESTRICT column = copy.getIndices();
3018  const CoinBigIndex *COIN_RESTRICT rowStart = copy.getVectorStarts();
3019  const int *COIN_RESTRICT rowLength = copy.getVectorLengths();
3020  const double *COIN_RESTRICT element = copy.getElements();
[1878]3021  int numberChanged = 1, iPass = 0;
3022  double large = largeValue(); // treat bounds > this as infinite
3023#ifndef NDEBUG
3024  double large2 = 1.0e10 * large;
3025#endif
3026  int numberInfeasible = 0;
3027  int totalTightened = 0;
[2385]3028
[1878]3029  double tolerance = primalTolerance();
[2385]3030
[1878]3031  // A copy of bounds is up at top
3032  translate(0); // move down
[2385]3033
[1878]3034#define MAXPASS 10
[2385]3035
[1878]3036  // Loop round seeing if we can tighten bounds
3037  // Would be faster to have a stack of possible rows
3038  // and we put altered rows back on stack
3039  int numberCheck = -1;
3040  // temp swap signs so can use existing code
[2385]3041  double *COIN_RESTRICT rowLower = abcLower_;
3042  double *COIN_RESTRICT rowUpper = abcUpper_;
3043  for (int iRow = 0; iRow < numberRows_; iRow++) {
3044    double value = -rowLower[iRow];
3045    rowLower[iRow] = -rowUpper[iRow];
3046    rowUpper[iRow] = value;
[1878]3047  }
3048  // Keep lower and upper for rows
3049  int whichArray[5];
[2385]3050  for (int i = 0; i < 5; i++)
3051    whichArray[i] = getAvailableArray();
3052  double *COIN_RESTRICT minRhs = usefulArray_[whichArray[0]].denseVector();
3053  double *COIN_RESTRICT maxRhs = usefulArray_[whichArray[1]].denseVector();
3054  int *COIN_RESTRICT changedList = usefulArray_[whichArray[0]].getIndices();
3055  int *COIN_RESTRICT changedListNext = usefulArray_[whichArray[1]].getIndices();
3056  char *COIN_RESTRICT changed = reinterpret_cast< char * >(usefulArray_[whichArray[2]].getIndices());
[1878]3057  // -1 no infinite, -2 more than one infinite >=0 column
[2385]3058  int *COIN_RESTRICT whichInfiniteDown = usefulArray_[whichArray[3]].getIndices();
3059  int *COIN_RESTRICT whichInfiniteUp = usefulArray_[whichArray[3]].getIndices();
3060  int numberToDoNext = 0;
3061  for (int iRow = 0; iRow < numberRows_; iRow++) {
[1878]3062    if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
[2385]3063      changedListNext[numberToDoNext++] = iRow;
3064      whichInfiniteDown[iRow] = -1;
3065      whichInfiniteUp[iRow] = -1;
[1878]3066    } else {
[2385]3067      minRhs[iRow] = -COIN_DBL_MAX;
3068      maxRhs[iRow] = COIN_DBL_MAX;
3069      whichInfiniteDown[iRow] = -2;
3070      whichInfiniteUp[iRow] = -2;
[1878]3071    }
3072  }
[2385]3073  const int *COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
3074  const CoinBigIndex *COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
3075  const double *COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
3076
3077  double *COIN_RESTRICT columnLower = abcLower_ + maximumAbcNumberRows_;
3078  double *COIN_RESTRICT columnUpper = abcUpper_ + maximumAbcNumberRows_;
3079  while (numberChanged > numberCheck) {
3080    int numberToDo = numberToDoNext;
3081    numberToDoNext = 0;
3082    int *COIN_RESTRICT temp = changedList;
3083    changedList = changedListNext;
3084    changedListNext = temp;
3085    CoinAbcMemset0(changed, numberRows_);
3086
[1878]3087    numberChanged = 0; // Bounds tightened this pass
[2385]3088
3089    if (iPass == MAXPASS)
3090      break;
[1878]3091    iPass++;
[2385]3092    for (int k = 0; k < numberToDo; k++) {
3093      int iRow = changedList[k];
[1878]3094      // possible row
3095      int infiniteUpper = 0;
3096      int infiniteLower = 0;
[2385]3097      int firstInfiniteUpper = -1;
3098      int firstInfiniteLower = -1;
[1878]3099      double maximumUp = 0.0;
3100      double maximumDown = 0.0;
3101      double newBound;
3102      CoinBigIndex rStart = rowStart[iRow];
3103      CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3104      CoinBigIndex j;
3105      // Compute possible lower and upper ranges
[2385]3106
[1878]3107      for (j = rStart; j < rEnd; ++j) {
[2385]3108        double value = element[j];
3109        int iColumn = column[j];
3110        if (value > 0.0) {
3111          if (columnUpper[iColumn] >= large) {
3112            firstInfiniteUpper = iColumn;
3113            ++infiniteUpper;
3114          } else {
3115            maximumUp += columnUpper[iColumn] * value;
3116          }
3117          if (columnLower[iColumn] <= -large) {
3118            firstInfiniteLower = iColumn;
3119            ++infiniteLower;
3120          } else {
3121            maximumDown += columnLower[iColumn] * value;
3122          }
3123        } else if (value < 0.0) {
3124          if (columnUpper[iColumn] >= large) {
3125            firstInfiniteLower = iColumn;
3126            ++infiniteLower;
3127          } else {
3128            maximumDown += columnUpper[iColumn] * value;
3129          }
3130          if (columnLower[iColumn] <= -large) {
3131            firstInfiniteUpper = iColumn;
3132            ++infiniteUpper;
3133          } else {
3134            maximumUp += columnLower[iColumn] * value;
3135          }
3136        }
[1878]3137      }
3138      // Build in a margin of error
3139      maximumUp += 1.0e-8 * fabs(maximumUp);
3140      maximumDown -= 1.0e-8 * fabs(maximumDown);
3141      double maxUp = maximumUp + infiniteUpper * 1.0e31;
3142      double maxDown = maximumDown - infiniteLower * 1.0e31;
[2385]3143      minRhs[iRow] = infiniteLower ? -COIN_DBL_MAX : maximumDown;
3144      maxRhs[iRow] = infiniteUpper ? COIN_DBL_MAX : maximumUp;
3145      if (infiniteLower == 0)
3146        whichInfiniteDown[iRow] = -1;
3147      else if (infiniteLower == 1)
3148        whichInfiniteDown[iRow] = firstInfiniteLower;
[1878]3149      else
[2385]3150        whichInfiniteDown[iRow] = -2;
3151      if (infiniteUpper == 0)
3152        whichInfiniteUp[iRow] = -1;
3153      else if (infiniteUpper == 1)
3154        whichInfiniteUp[iRow] = firstInfiniteUpper;
[1878]3155      else
[2385]3156        whichInfiniteUp[iRow] = -2;
3157      if (maxUp <= rowUpper[iRow] + tolerance && maxDown >= rowLower[iRow] - tolerance) {
3158
3159        // Row is redundant - make totally free
3160        // NO - can't do this for postsolve
3161        // rowLower[iRow]=-COIN_DBL_MAX;
3162        // rowUpper[iRow]=COIN_DBL_MAX;
3163        //printf("Redundant row in presolveX %d\n",iRow);
3164
[1878]3165      } else {
[2385]3166        if (maxUp < rowLower[iRow] - 100.0 * tolerance || maxDown > rowUpper[iRow] + 100.0 * tolerance) {
3167          // problem is infeasible - exit at once
3168          numberInfeasible++;
3169          break;
3170        }
3171        double lower = rowLower[iRow];
3172        double upper = rowUpper[iRow];
3173        for (j = rStart; j < rEnd; ++j) {
3174          double value = element[j];
3175          int iColumn = column[j];
3176          double nowLower = columnLower[iColumn];
3177          double nowUpper = columnUpper[iColumn];
3178          if (value > 0.0) {
3179            // positive value
3180            if (lower > -large) {
3181              if (!infiniteUpper) {
3182                assert(nowUpper < large2);
3183                newBound = nowUpper + (lower - maximumUp) / value;
3184                // relax if original was large
3185                if (fabs(maximumUp) > 1.0e8)
3186                  newBound -= 1.0e-12 * fabs(maximumUp);
3187              } else if (infiniteUpper == 1 && nowUpper > large) {
3188                newBound = (lower - maximumUp) / value;
3189                // relax if original was large
3190                if (fabs(maximumUp) > 1.0e8)
3191                  newBound -= 1.0e-12 * fabs(maximumUp);
3192              } else {
3193                newBound = -COIN_DBL_MAX;
3194              }
3195              if (newBound > nowLower + 1.0e-12 && newBound > -large) {
3196                // Tighten the lower bound
3197                numberChanged++;
3198                // check infeasible (relaxed)
3199                if (nowUpper < newBound) {
3200                  if (nowUpper - newBound < -100.0 * tolerance)
3201                    numberInfeasible++;
3202                  else
3203                    newBound = nowUpper;
3204                }
3205                columnLower[iColumn] = newBound;
3206                for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3207                  int jRow = row[j1];
3208                  if (!changed[jRow]) {
3209                    changed[jRow] = 1;
3210                    changedListNext[numberToDoNext++] = jRow;
3211                  }
3212                }
3213                // adjust
3214                double now;
3215                if (nowLower < -large) {
3216                  now = 0.0;
3217                  infiniteLower--;
3218                } else {
3219                  now = nowLower;
3220                }
3221                maximumDown += (newBound - now) * value;
3222                nowLower = newBound;
3223              }
3224            }
3225            if (upper < large) {
3226              if (!infiniteLower) {
3227                assert(nowLower > -large2);
3228                newBound = nowLower + (upper - maximumDown) / value;
3229                // relax if original was large
3230                if (fabs(maximumDown) > 1.0e8)
3231                  newBound += 1.0e-12 * fabs(maximumDown);
3232              } else if (infiniteLower == 1 && nowLower < -large) {
3233                newBound = (upper - maximumDown) / value;
3234                // relax if original was large
3235                if (fabs(maximumDown) > 1.0e8)
3236                  newBound += 1.0e-12 * fabs(maximumDown);
3237              } else {
3238                newBound = COIN_DBL_MAX;
3239              }
3240              if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3241                // Tighten the upper bound
3242                numberChanged++;
3243                // check infeasible (relaxed)
3244                if (nowLower > newBound) {
3245                  if (newBound - nowLower < -100.0 * tolerance)
3246                    numberInfeasible++;
3247                  else
3248                    newBound = nowLower;
3249                }
3250                columnUpper[iColumn] = newBound;
3251                for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3252                  int jRow = row[j1];
3253                  if (!changed[jRow]) {
3254                    changed[jRow] = 1;
3255                    changedListNext[numberToDoNext++] = jRow;
3256                  }
3257                }
3258                // adjust
3259                double now;
3260                if (nowUpper > large) {
3261                  now = 0.0;
3262                  infiniteUpper--;
3263                } else {
3264                  now = nowUpper;
3265                }
3266                maximumUp += (newBound - now) * value;
3267                nowUpper = newBound;
3268              }
3269            }
3270          } else {
3271            // negative value
3272            if (lower > -large) {
3273              if (!infiniteUpper) {
3274                assert(nowLower < large2);
3275                newBound = nowLower + (lower - maximumUp) / value;
3276                // relax if original was large
3277                if (fabs(maximumUp) > 1.0e8)
3278                  newBound += 1.0e-12 * fabs(maximumUp);
3279              } else if (infiniteUpper == 1 && nowLower < -large) {
3280                newBound = (lower - maximumUp) / value;
3281                // relax if original was large
3282                if (fabs(maximumUp) > 1.0e8)
3283                  newBound += 1.0e-12 * fabs(maximumUp);
3284              } else {
3285                newBound = COIN_DBL_MAX;
3286              }
3287              if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3288                // Tighten the upper bound
3289                numberChanged++;
3290                // check infeasible (relaxed)
3291                if (nowLower > newBound) {
3292                  if (newBound - nowLower < -100.0 * tolerance)
3293                    numberInfeasible++;
3294                  else
3295                    newBound = nowLower;
3296                }
3297                columnUpper[iColumn] = newBound;
3298                for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3299                  int jRow = row[j1];
3300                  if (!changed[jRow]) {
3301                    changed[jRow] = 1;
3302                    changedListNext[numberToDoNext++] = jRow;
3303                  }
3304                }
3305                // adjust
3306                double now;
3307                if (nowUpper > large) {
3308                  now = 0.0;
3309                  infiniteLower--;
3310                } else {
3311                  now = nowUpper;
3312                }
3313                maximumDown += (newBound