source: stable/1.17/Clp/src/AbcSimplex.cpp @ 2431

Last change on this file since 2431 was 2431, checked in by stefan, 3 months ago

sync with trunk

  • Property svn:keywords set to Id
File size: 226.0 KB
Line 
1/* $Id: AbcSimplex.cpp 2431 2019-03-15 15:56:51Z stefan $ */
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
12#if SLIM_CLP == 2
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//#############################################################################
41AbcSimplex::AbcSimplex(bool emptyMessages)
42  :
43
44  ClpSimplex(emptyMessages)
45{
46  gutsOfDelete(0);
47  gutsOfInitialize(0, 0, true);
48  assert(maximumAbcNumberRows_ >= 0);
49  //printf("XX %x AbcSimplex constructor\n",this);
50}
51
52//-----------------------------------------------------------------------------
53
54AbcSimplex::~AbcSimplex()
55{
56  //printf("XX %x AbcSimplex destructor\n",this);
57  gutsOfDelete(1);
58}
59// Copy constructor.
60AbcSimplex::AbcSimplex(const AbcSimplex &rhs)
61  : ClpSimplex(rhs)
62{
63  //printf("XX %x AbcSimplex constructor from %x\n",this,&rhs);
64  gutsOfDelete(0);
65  gutsOfInitialize(numberRows_, numberColumns_, false);
66  gutsOfCopy(rhs);
67  assert(maximumAbcNumberRows_ >= 0);
68}
69#include "ClpDualRowSteepest.hpp"
70#include "ClpPrimalColumnSteepest.hpp"
71#include "ClpFactorization.hpp"
72// Copy constructor from model
73AbcSimplex::AbcSimplex(const ClpSimplex &rhs)
74  : ClpSimplex(rhs)
75{
76  //printf("XX %x AbcSimplex constructor from ClpSimplex\n",this);
77  gutsOfDelete(0);
78  gutsOfInitialize(numberRows_, numberColumns_, true);
79  gutsOfResize(numberRows_, numberColumns_);
80  // Set up row/column selection and factorization type
81  ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(rhs.dualRowPivot());
82  if (pivot) {
83    AbcDualRowSteepest steep(pivot->mode());
84    setDualRowPivotAlgorithm(steep);
85  } else {
86    AbcDualRowDantzig dantzig;
87    setDualRowPivotAlgorithm(dantzig);
88  }
89  ClpPrimalColumnSteepest *pivotColumn = dynamic_cast< ClpPrimalColumnSteepest * >(rhs.primalColumnPivot());
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();
99  factorization()->synchronize(rhs.factorization(), this);
100  //factorization_->setGoDenseThreshold(rhs.factorization()->goDenseThreshold());
101  //factorization_->setGoSmallThreshold(rhs.factorization()->goSmallThreshold());
102  translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
103  assert(maximumAbcNumberRows_ >= 0);
104}
105// Assignment operator. This copies the data
106AbcSimplex &
107AbcSimplex::operator=(const AbcSimplex &rhs)
108{
109  if (this != &rhs) {
110    gutsOfDelete(1);
111    ClpSimplex::operator=(rhs);
112    gutsOfCopy(rhs);
113    assert(maximumAbcNumberRows_ >= 0);
114  }
115  return *this;
116}
117// fills in perturbationSaved_ from start with 0.5+random
118void AbcSimplex::fillPerturbation(int start, int number)
119{
120  double *array = perturbationSaved_ + start;
121  for (int i = 0; i < number; i++)
122    array[i] = 0.5 + 0.5 * randomNumberGenerator_.randomDouble();
123}
124// Sets up all extra pointers
125void AbcSimplex::setupPointers(int maxRows, int maxColumns)
126{
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;
146}
147// Copies all saved versions to working versions and may do something for perturbation
148void AbcSimplex::copyFromSaved(int which)
149{
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_);
163  }
164}
165void AbcSimplex::gutsOfCopy(const AbcSimplex &rhs)
166{
167#ifdef ABC_INHERIT
168  abcSimplex_ = this;
169#endif
170  assert(numberRows_ == rhs.numberRows_);
171  assert(numberColumns_ == rhs.numberColumns_);
172  numberTotal_ = rhs.numberTotal_;
173  maximumNumberTotal_ = rhs.maximumNumberTotal_;
174  // special options here to be safe
175  specialOptions_ = rhs.specialOptions_;
176  //assert (maximumInternalRows_ >= numberRows_);
177  //assert (maximumInternalColumns_ >= numberColumns_);
178  // Copy all scalars
179  CoinAbcMemcpy(reinterpret_cast< int * >(&sumNonBasicCosts_),
180    reinterpret_cast< const int * >(&rhs.sumNonBasicCosts_),
181    (&swappedAlgorithm_ - reinterpret_cast< int * >(&sumNonBasicCosts_)) + 1);
182  // could add 2 so can go off end
183  int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_;
184  internalStatus_ = ClpCopyOfArray(rhs.internalStatus_, sizeArray + maximumNumberTotal_);
185  abcLower_ = ClpCopyOfArray(rhs.abcLower_, sizeArray);
186  abcUpper_ = ClpCopyOfArray(rhs.abcUpper_, sizeArray);
187  abcCost_ = ClpCopyOfArray(rhs.abcCost_, sizeArray + maximumNumberTotal_);
188  abcDj_ = ClpCopyOfArray(rhs.abcDj_, sizeArray);
189
190  abcSolution_ = ClpCopyOfArray(rhs.abcSolution_, sizeArray + maximumNumberTotal_);
191  abcPerturbation_ = ClpCopyOfArray(rhs.abcPerturbation_, sizeArray);
192  abcPivotVariable_ = ClpCopyOfArray(rhs.abcPivotVariable_, maximumAbcNumberRows_);
193  //fromExternal_ = ClpCopyOfArray(rhs.fromExternal_,sizeArray);
194  //toExternal_ = ClpCopyOfArray(rhs.toExternal_,sizeArray);
195  scaleFromExternal_ = ClpCopyOfArray(rhs.scaleFromExternal_, sizeArray);
196  offset_ = ClpCopyOfArray(rhs.offset_, sizeArray);
197  setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
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
237void AbcSimplex::gutsOfDelete(int type)
238{
239  if (type) {
240    delete[] abcLower_;
241    delete[] abcUpper_;
242    delete[] abcCost_;
243    delete[] abcDj_;
244    delete[] abcSolution_;
245    //delete [] fromExternal_ ;
246    //delete [] toExternal_ ;
247    delete[] scaleFromExternal_;
248    //delete [] scaleToExternal_ ;
249    delete[] offset_;
250    delete[] internalStatus_;
251    delete[] abcPerturbation_;
252    delete abcMatrix_;
253    delete abcFactorization_;
254#ifdef EARLY_FACTORIZE
255    delete abcEarlyFactorization_;
256#endif
257    delete[] abcPivotVariable_;
258    delete abcDualRowPivot_;
259    delete abcPrimalColumnPivot_;
260    delete abcBaseModel_;
261    delete clpModel_;
262    delete abcNonLinearCost_;
263  }
264  CoinAbcMemset0(reinterpret_cast< char * >(&scaleToExternal_),
265    reinterpret_cast< char * >(&usefulArray_[0]) - reinterpret_cast< char * >(&scaleToExternal_));
266}
267template < class T >
268T *newArray(T * /*nullArray*/, int size)
269{
270  if (size) {
271    T *arrayNew = new T[size];
272#ifndef NDEBUG
273    memset(arrayNew, 'A', (size) * sizeof(T));
274#endif
275    return arrayNew;
276  } else {
277    return NULL;
278  }
279}
280template < class T >
281T *resizeArray(T *array, int oldSize1, int oldSize2, int newSize1, int newSize2, int extra)
282{
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];
289#ifndef NDEBUG
290      memset(arrayNew, 'A', (2 * newTotal + newSize1 + extra) * sizeof(T));
291#endif
292      CoinAbcMemcpy(arrayNew, array, oldSize1);
293      CoinAbcMemcpy(arrayNew + newSize1, array + oldSize1, oldSize2);
294      // and second half
295      CoinAbcMemcpy(arrayNew + newTotal, array, oldSize1 + oldTotal);
296      CoinAbcMemcpy(arrayNew + newSize1 + newTotal, array + oldSize1 + oldTotal, oldSize2);
297      delete[] array;
298    } else {
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];
306    }
307    return arrayNew;
308  } else {
309    return array;
310  }
311}
312// Initializes arrays
313void AbcSimplex::gutsOfInitialize(int numberRows, int numberColumns, bool doMore)
314{
315#ifdef ABC_INHERIT
316  abcSimplex_ = this;
317#endif
318  // Zero all
319  CoinAbcMemset0(&sumNonBasicCosts_, (reinterpret_cast< double * >(&usefulArray_[0]) - &sumNonBasicCosts_));
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;
340  numberRefinements_ = 1;
341  changeMade_ = 1;
342  forceFactorization_ = -1;
343  if (perturbation_ < 50 || (perturbation_ > 60 && perturbation_ < 100))
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();
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_;
362  if (doMore) {
363    // say Steepest pricing
364    abcDualRowPivot_ = new AbcDualRowSteepest();
365    abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
366    internalStatus_ = newArray((unsigned char *)NULL,
367      sizeArray + maximumNumberTotal_);
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_);
379    // Fill perturbation array
380    setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
381    fillPerturbation(0, maximumNumberTotal_);
382  }
383  // get an empty factorization so we can set tolerances etc
384  getEmptyFactorization();
385  for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
386    usefulArray_[i].reserve(CoinMax(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200), 2 * numberRows_));
387  //savedStatus_ = internalStatus_+maximumNumberTotal_;
388  //startPermanentArrays();
389}
390// resizes arrays
391void AbcSimplex::gutsOfResize(int numberRows, int numberColumns)
392{
393  if (!abcSolution_) {
394    numberRows_ = 0;
395    numberColumns_ = 0;
396    numberTotal_ = 0;
397    maximumAbcNumberRows_ = 0;
398    maximumAbcNumberColumns_ = 0;
399    maximumNumberTotal_ = 0;
400  }
401  if (numberRows == numberRows_ && numberColumns == numberColumns_)
402    return;
403  // can do on state bit
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_;
411  //fromExternal_ = resizeArray(fromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1);
412  //toExternal_ = resizeArray(toExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
413  scaleFromExternal_ = resizeArray(scaleFromExternal_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
414  //scaleToExternal_ = resizeArray(scaleToExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
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);
426  // Fill gaps in perturbation array
427  fillPerturbation(maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_);
428  fillPerturbation(newSize1 + maximumAbcNumberColumns_, newSize2 - maximumAbcNumberColumns_);
429  // Clean gap
430  //CoinIotaN(fromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
431  //CoinIotaN(toExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
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_;
445  abcPivotVariable_ = new int[maximumAbcNumberRows_];
446  for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
447    usefulArray_[i].reserve(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200));
448}
449void AbcSimplex::translate(int type)
450{
451  // clear bottom bits
452  stateOfProblem_ &= ~(VALUES_PASS - 1);
453  if ((type & DO_SCALE_AND_MATRIX) != 0) {
454    //stateOfProblem_ |= DO_SCALE_AND_MATRIX;
455    stateOfProblem_ |= DO_SCALE_AND_MATRIX;
456    delete abcMatrix_;
457    abcMatrix_ = new AbcMatrix(*matrix());
458    abcMatrix_->setModel(this);
459    abcMatrix_->scale(scalingFlag_ ? 0 : -1);
460  }
461  if ((type & DO_STATUS) != 0 && (type & DO_BASIS_AND_ORDER) == 0) {
462    // from Clp enum to Abc enum (and bound flip)
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_;
466    int i;
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];
473      else
474        break;
475    }
476    if (!i) {
477      // from Clp enum to Abc enum
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;
489      }
490      if (i)
491        type |= DO_BASIS_AND_ORDER;
492    } else {
493      type |= DO_BASIS_AND_ORDER;
494    }
495    stateOfProblem_ |= DO_STATUS;
496  }
497  if ((type & DO_BASIS_AND_ORDER) != 0) {
498    permuteIn();
499    permuteBasis();
500    stateOfProblem_ |= DO_BASIS_AND_ORDER;
501    type &= ~DO_SOLUTION;
502  }
503  if ((type & DO_SOLUTION) != 0) {
504    permuteBasis();
505    stateOfProblem_ |= DO_SOLUTION;
506  }
507  if ((type & DO_JUST_BOUNDS) != 0) {
508    stateOfProblem_ |= DO_JUST_BOUNDS;
509  }
510  if (!type) {
511    // just move stuff down
512    CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_);
513    CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_);
514    CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
515  }
516}
517#ifdef ABC_SPRINT
518// Overwrite to create sub problem (just internal arrays) - save full stuff
519AbcSimplex *
520AbcSimplex::createSubProblem(int numberColumns, const int *whichColumn)
521{
522  int numberFullColumns = numberColumns_;
523  numberColumns_ = numberColumns;
524  AbcSimplex *subProblem = this;
525  AbcSimplex *fullProblem = reinterpret_cast< AbcSimplex * >(new char[sizeof(AbcSimplex)]);
526  memset(fullProblem, 0, sizeof(AbcSimplex));
527  fullProblem->maximumAbcNumberRows_ = maximumAbcNumberRows_;
528  fullProblem->numberColumns_ = numberFullColumns;
529  fullProblem->numberTotal_ = numberTotal_;
530  fullProblem->maximumNumberTotal_ = maximumNumberTotal_;
531  fullProblem->numberTotalWithoutFixed_ = numberTotalWithoutFixed_;
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_;
544  fullProblem->abcNonLinearCost_ = abcNonLinearCost_;
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_;
550  subProblem->internalStatus_ = newArray((unsigned char *)NULL,
551    sizeArray + subProblem->maximumNumberTotal_);
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_);
563  subProblem->setupPointers(maximumAbcNumberRows_, numberColumns);
564  // could use arrays - but for now be safe
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;
573    int iPivot = fullProblem->abcPivotVariable_[i];
574    if (iPivot < numberRows_) {
575      subProblem->abcPivotVariable_[i] = iPivot;
576    } else {
577      assert(backward[iPivot - numberRows_] >= 0);
578      subProblem->abcPivotVariable_[i] = backward[iPivot - numberRows_];
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  }
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];
619  }
620  subProblem->abcNonLinearCost_ = new AbcNonLinearCost(subProblem);
621  subProblem->abcNonLinearCost_->checkInfeasibilities(0.0);
622  subProblem->abcMatrix_ = new AbcMatrix(*fullProblem->abcMatrix_, numberRows_, whichRow,
623    numberColumns, whichColumn);
624  subProblem->abcMatrix_->setModel(subProblem);
625  subProblem->abcMatrix_->rebalance();
626  subProblem->abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
627  subProblem->abcPrimalColumnPivot_->saveWeights(subProblem, 2);
628  delete[] backward;
629  // swap
630  return fullProblem;
631}
632// Restore stuff from sub problem (and delete sub problem)
633void AbcSimplex::restoreFromSubProblem(AbcSimplex *fullProblem, const int *whichColumn)
634{
635  AbcSimplex *subProblem = this;
636  for (int i = 0; i < numberRows_; i++) {
637    int iPivot = subProblem->abcPivotVariable_[i];
638    if (iPivot < numberRows_) {
639      fullProblem->abcPivotVariable_[i] = iPivot;
640    } else {
641      fullProblem->abcPivotVariable_[i] = whichColumn[iPivot - numberRows_] + numberRows_;
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_;
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];
682  }
683  delete[] subProblem->internalStatus_;
684  delete[] subProblem->abcPerturbation_;
685  delete subProblem->abcMatrix_;
686  delete[] subProblem->abcLower_;
687  delete[] subProblem->abcUpper_;
688  delete[] subProblem->abcCost_;
689  delete[] subProblem->abcSolution_;
690  delete[] subProblem->abcDj_;
691  delete subProblem->abcPrimalColumnPivot_;
692  delete[] subProblem->scaleFromExternal_;
693  delete[] subProblem->offset_;
694  delete[] subProblem->abcPivotVariable_;
695  delete[] subProblem->reversePivotVariable_;
696  delete subProblem->abcNonLinearCost_;
697  numberColumns_ = fullProblem->numberColumns_;
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);
714  // ? redo nonlinearcost
715  abcNonLinearCost_ = fullProblem->abcNonLinearCost_;
716  abcNonLinearCost_->refresh();
717  delete[] reinterpret_cast< char * >(fullProblem);
718}
719#endif
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*/
725void AbcSimplex::setupDualValuesPass(const double *fakeDuals,
726  const double *fakePrimals,
727  int type)
728{
729  // allslack basis
730  memset(internalStatus_, 6, numberRows_);
731  // temp
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;
738      }
739    }
740    if (allEqual) {
741      // just modify costs
742      transposeTimes(-1.0, fakeDuals, objective());
743      return;
744    }
745  }
746  // compute unscaled djs
747  CoinAbcMemcpy(djSaved_ + maximumAbcNumberRows_, objective(), numberColumns_);
748  matrix_->transposeTimes(-1.0, fakeDuals, djSaved_ + maximumAbcNumberRows_);
749  // save fake solution
750  assert(solution_);
751  //solution_ = new double[numberTotal_];
752  CoinAbcMemset0(solution_, numberRows_);
753  CoinAbcMemcpy(solution_ + maximumAbcNumberRows_, fakePrimals, numberColumns_);
754  // adjust
755  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++)
756    solution_[iSequence] -= offset_[iSequence];
757  matrix_->times(-1.0, solution_ + maximumAbcNumberRows_, solution_);
758  double direction = optimizationDirection_;
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];
763    solution_[iRow] *= rowScale[iRow];
764  }
765  const double *COIN_RESTRICT columnScale = scaleToExternal_;
766  for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++)
767    djSaved_[iColumn] *= direction * columnScale[iColumn];
768  // Set solution values
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        }
811      } else {
812        solution[iColumn] = lower[iColumn];
813        setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound);
814        setStatus(iColumn, ClpSimplex::atLowerBound);
815      }
816    } else if (thisUpper < 1.0e30) {
817      solution[iColumn] = upper[iColumn];
818      setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound);
819      setStatus(iColumn, ClpSimplex::atUpperBound);
820    } else {
821      // free
822      solution[iColumn] = thisValue * inverseColumnScale[iColumn];
823      setInternalStatus(iColumn + maximumAbcNumberRows_, isFree);
824      setStatus(iColumn, ClpSimplex::isFree);
825    }
826  }
827  switch (type) {
828  case 1:
829    stateOfProblem_ |= VALUES_PASS;
830    break;
831  case 2:
832    stateOfProblem_ |= VALUES_PASS2;
833    delete[] solution_;
834    solution_ = NULL;
835    break;
836  case 3:
837    stateOfProblem_ |= VALUES_PASS2;
838    break;
839  }
840}
841//#############################################################################
842int AbcSimplex::computePrimals(CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
843{
844
845  arrayVector->clear();
846  previousVector->clear();
847  // accumulate non basic stuff
848  double *COIN_RESTRICT array = arrayVector->denseVector();
849  CoinAbcScatterZeroTo(abcSolution_, abcPivotVariable_, numberRows_);
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;
860  CoinIndexedVector *thisVector = arrayVector;
861  CoinIndexedVector *lastVector = previousVector;
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
893  int numberRefinements = 0;
894  for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
895    int numberIn = thisVector->getNumElements();
896    const int *COIN_RESTRICT indexIn = thisVector->getIndices();
897    const double *COIN_RESTRICT arrayIn = thisVector->denseVector();
898    CoinAbcScatterToList(arrayIn, abcSolution_, indexIn, abcPivotVariable_, numberIn);
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;
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) {
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
924      CoinIndexedVector *temp = thisVector;
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
934      CoinIndexedVector *temp = thisVector;
935      thisVector = lastVector;
936      lastVector = temp;
937      int *COIN_RESTRICT indexOut = thisVector->getIndices();
938      int number = 0;
939      array = thisVector->denseVector();
940      thisVector->clear();
941      for (int iRow = 0; iRow < numberRows_; iRow++) {
942        double value = solutionBasic_[iRow];
943        if (value) {
944          array[iRow] = value;
945          indexOut[number++] = iRow;
946          solutionBasic_[iRow] = 0.0;
947        }
948      }
949      thisVector->setNumElements(number);
950      lastError = largestPrimalError_;
951      abcFactorization_->updateFullColumn(*thisVector);
952      double *previous = lastVector->denseVector();
953      number = 0;
954      multiplier = 1.0 / multiplier;
955      for (int iRow = 0; iRow < numberRows_; iRow++) {
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        }
963      }
964      thisVector->setNumElements(number);
965    } else {
966      break;
967    }
968  }
969
970  // solution as accurate as we are going to get
971  //if (!goodSolution) {
972  CoinAbcMemcpy(solutionBasic_, thisVector->denseVector(), numberRows_);
973  CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
974  arrayVector->clear();
975  previousVector->clear();
976  return numberRefinements;
977}
978// Moves basic stuff to basic area
979void AbcSimplex::moveToBasic(int which)
980{
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_);
989}
990// now dual side
991int AbcSimplex::computeDuals(double *givenDjs, CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
992{
993  double *COIN_RESTRICT array = arrayVector->denseVector();
994  int *COIN_RESTRICT index = arrayVector->getIndices();
995  int number = 0;
996  if (!givenDjs) {
997    for (int iRow = 0; iRow < numberRows_; iRow++) {
998      double value = costBasic_[iRow];
999      if (value) {
1000        array[iRow] = value;
1001        index[number++] = iRow;
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))
1010        givenDjs[iPivot] = 0.0;
1011      double value = abcCost_[iPivot] - givenDjs[iPivot];
1012      if (value) {
1013        array[iRow] = value;
1014        index[number++] = iRow;
1015      }
1016    }
1017  }
1018  arrayVector->setNumElements(number);
1019  // Btran basic costs and get as accurate as possible
1020  double lastError = COIN_DBL_MAX;
1021  CoinIndexedVector *thisVector = arrayVector;
1022  CoinIndexedVector *lastVector = previousVector;
1023  abcFactorization_->updateFullColumnTranspose(*thisVector);
1024
1025  int numberRefinements = 0;
1026  for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1027    // check basic reduced costs zero
1028    // put reduced cost of basic into djBasic_
1029    CoinAbcMemcpy(djBasic_, costBasic_, numberRows_);
1030    abcMatrix_->transposeTimesBasic(-1.0, thisVector->denseVector(), djBasic_);
1031    largestDualError_ = CoinAbcMaximumAbsElement(djBasic_, numberRows_);
1032    if (largestDualError_ >= lastError) {
1033      // restore
1034      CoinIndexedVector *temp = thisVector;
1035      thisVector = lastVector;
1036      lastVector = temp;
1037      break;
1038    }
1039    if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
1040      && !givenDjs) {
1041      numberRefinements++;
1042      // try and make better
1043      // save this
1044      CoinIndexedVector *temp = thisVector;
1045      thisVector = lastVector;
1046      lastVector = temp;
1047      array = thisVector->denseVector();
1048      double multiplier = 131072.0;
1049      //array=djBasic_*multiplier
1050      CoinAbcMultiplyAdd(djBasic_, numberRows_, multiplier, array, 0.0);
1051      lastError = largestDualError_;
1052      abcFactorization_->updateFullColumnTranspose(*thisVector);
1053#if ABC_DEBUG
1054      thisVector->checkClean();
1055#endif
1056      multiplier = 1.0 / multiplier;
1057      double *COIN_RESTRICT previous = lastVector->denseVector();
1058      // array = array*multiplier+previous
1059      CoinAbcMultiplyAdd(previous, numberRows_, 1.0, array, multiplier);
1060    } else {
1061      break;
1062    }
1063  }
1064  // now look at dual solution
1065  CoinAbcMemcpy(dual_, thisVector->denseVector(), numberRows_);
1066  CoinAbcMemset0(thisVector->denseVector(), numberRows_);
1067  thisVector->setNumElements(0);
1068  if (numberRefinements) {
1069    CoinAbcMemset0(lastVector->denseVector(), numberRows_);
1070    lastVector->setNumElements(0);
1071  }
1072  CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
1073  abcMatrix_->transposeTimesNonBasic(-1.0, dual_, abcDj_);
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*/
1096int AbcSimplex::internalFactorize(int solveType)
1097{
1098  assert(status_);
1099
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++;
1145#if ABC_NORMAL_DEBUG > 1
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++;
1163#if ABC_NORMAL_DEBUG > 0
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  }
1181#if ABC_NORMAL_DEBUG > 0
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
1188  // *** replace below by cleanStatus
1189  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1190    AbcSimplex::Status status = getInternalStatus(iSequence);
1191    if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1192      setInternalStatus(iSequence, isFixed);
1193  }
1194  if (numberIterations_ == baseIteration_ && !valuesPass) {
1195    double largeValue = this->largeValue();
1196    double *COIN_RESTRICT solution = abcSolution_;
1197    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
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        }
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)
1233        numberSlacks++;
1234    }
1235    status = CoinMax(numberSlacks - numberRows_, 0);
1236    if (status)
1237      printf("%d singularities\n", status);
1238    // special case if all slack
1239    if (numberSlacks == numberRows_) {
1240      status = numberRows_ + 1;
1241    }
1242  }
1243#ifdef ABC_USE_COIN_FACTORIZATION
1244  // sparse methods
1245  abcFactorization_->goSparse();
1246#endif
1247  return status;
1248}
1249// Make sure no superbasic etc
1250void AbcSimplex::cleanStatus(bool valuesPass)
1251{
1252  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1253    AbcSimplex::Status status = getInternalStatus(iSequence);
1254    if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1255      setInternalStatus(iSequence, isFixed);
1256  }
1257  if (numberIterations_ == baseIteration_ && !valuesPass) {
1258    double largeValue = this->largeValue();
1259    double *COIN_RESTRICT solution = abcSolution_;
1260    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
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        }
1282      }
1283    }
1284  }
1285}
1286// Sets objectiveValue_ from rawObjectiveValue_
1287void AbcSimplex::setClpSimplexObjectiveValue()
1288{
1289  objectiveValue_ = rawObjectiveValue_ / (objectiveScale_ * rhsScale_);
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*/
1296int AbcSimplex::housekeeping()
1297{
1298#define DELAYED_UPDATE
1299#ifdef DELAYED_UPDATE
1300  if (algorithm_ < 0) {
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    }
1326#if ABC_NORMAL_DEBUG > 3
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_)
1356        << sequenceIn_
1357        << CoinMessageEol;
1358    }
1359  }
1360  // change of incoming
1361  char rowcol[] = { 'R', 'C' };
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)
1373        alphaAccuracy_ *= value;
1374      else
1375        alphaAccuracy_ /= value;
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_])) {
1382        // going to lower
1383        setInternalStatus(sequenceOut_, atLowerBound);
1384      } else {
1385        // going to upper
1386        setInternalStatus(sequenceOut_, atUpperBound);
1387      }
1388    } else {
1389      // fixed
1390      setInternalStatus(sequenceOut_, isFixed);
1391    }
1392#ifndef DELAYED_UPDATE
1393    abcSolution_[sequenceOut_] = valueOut_;
1394#endif
1395#if PARTITION_ROW_COPY == 1
1396    // move in row copy
1397    abcMatrix_->inOutUseful(sequenceIn_, sequenceOut_);
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,
1437    directionIn_, directionOut_);
1438  if (cycle > 0) {
1439#if ABC_NORMAL_DEBUG > 0
1440    if (handler_->logLevel() >= 63)
1441      printf("Cycle of %d\n", cycle);
1442#endif
1443    // reset
1444    abcProgress_.startCheck();
1445    double random = randomNumberGenerator_.randomDouble();
1446    int extra = static_cast< int >(9.999 * random);
1447    int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
1448    if (abcFactorization_->pivots() > cycle) {
1449      forceFactorization_ = CoinMax(1, cycle - off[extra]);
1450    } else {
1451      // need to reject something
1452      int iSequence;
1453      if (algorithm_ < 0) {
1454        iSequence = sequenceIn_;
1455      } else {
1456        /* should be better if don't reject incoming
1457         as it is in basis */
1458        iSequence = sequenceOut_;
1459        // so can't happen immediately again
1460        sequenceOut_ = -1;
1461      }
1462      char x = isColumn(iSequence) ? 'C' : 'R';
1463      if (handler_->logLevel() >= 63)
1464        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1465          << x << sequenceWithin(iSequence)
1466          << CoinMessageEol;
1467      setFlagged(iSequence);
1468      //printf("flagging %d\n",iSequence);
1469    }
1470    return 1;
1471  }
1472  int invertNow = 0;
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();
1476  if (algorithm_ < 0)
1477    numberPivots++; // allow for update not done
1478  int maximumPivots = abcFactorization_->maximumPivots();
1479  int numberDense = 0; //abcFactorization_->numberDense();
1480  bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations());
1481  if (numberPivots == maximumPivots || maximumPivots < 2) {
1482    // If dense then increase
1483    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots && false) {
1484      abcFactorization_->maximumPivots(numberDense);
1485    }
1486    //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots);
1487#if CLP_FACTORIZATION_NEW_TIMING > 1
1488    abcFactorization_->statsRefactor('M');
1489#endif
1490    return 1;
1491  } else if ((abcFactorization_->timeToRefactorize() && !dontInvert)
1492    || invertNow) {
1493    //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots);
1494#if CLP_FACTORIZATION_NEW_TIMING > 1
1495    abcFactorization_->statsRefactor('T');
1496#endif
1497    return 1;
1498  } else if (forceFactorization_ > 0 &&
1499#ifndef DELAYED_UPDATE
1500    abcFactorization_->pivots()
1501#else
1502    abcFactorization_->pivots() + 1
1503#endif
1504      >= forceFactorization_) {
1505    // relax
1506    forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
1507    if (forceFactorization_ > abcFactorization_->maximumPivots())
1508      forceFactorization_ = -1; //off
1509      //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots);
1510#if CLP_FACTORIZATION_NEW_TIMING > 1
1511    abcFactorization_->statsRefactor('F');
1512#endif
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();
1518    int window = numberTooManyIterations % 5000;
1519    if (window < 2 * maximumPivots)
1520      random = 0.2 * random + 0.8; // randomly re-factorize but not too soon
1521    else
1522      random = 1.0; // switch off if not in window of opportunity
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;
1527    } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1000010 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
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}
1539// Swaps primal stuff
1540void AbcSimplex::swapPrimalStuff()
1541{
1542  if (sequenceOut_ < 0)
1543    return;
1544  assert(sequenceIn_ >= 0);
1545  abcSolution_[sequenceOut_] = valueOut_;
1546  abcSolution_[sequenceIn_] = valueIn_;
1547  solutionBasic_[pivotRow_] = valueIn_;
1548  if (algorithm_ < 0) {
1549    // and set bounds correctly
1550    static_cast< AbcSimplexDual * >(this)->originalBound(sequenceIn_);
1551    static_cast< AbcSimplexDual * >(this)->changeBound(sequenceOut_);
1552  } else {
1553#if ABC_NORMAL_DEBUG > 2
1554    if (handler_->logLevel() == 63)
1555      printf("Code swapPrimalStuff for primal\n");
1556#endif
1557  }
1558  if (pivotRow_ >= 0) { // may be flip in primal
1559    lowerBasic_[pivotRow_] = abcLower_[sequenceIn_];
1560    upperBasic_[pivotRow_] = abcUpper_[sequenceIn_];
1561    costBasic_[pivotRow_] = abcCost_[sequenceIn_];
1562    abcPivotVariable_[pivotRow_] = sequenceIn_;
1563  }
1564}
1565// Swaps dual stuff
1566void AbcSimplex::swapDualStuff(int lastSequenceOut, int lastDirectionOut)
1567{
1568  // outgoing
1569  // set dj to zero unless values pass
1570  if (lastSequenceOut >= 0) {
1571    if ((stateOfProblem_ & VALUES_PASS) == 0) {
1572      if (lastDirectionOut > 0) {
1573        abcDj_[lastSequenceOut] = theta_;
1574      } else {
1575        abcDj_[lastSequenceOut] = -theta_;
1576      }
1577    } else {
1578      if (lastDirectionOut > 0) {
1579        abcDj_[lastSequenceOut] -= theta_;
1580      } else {
1581        abcDj_[lastSequenceOut] += theta_;
1582      }
1583    }
1584  }
1585  if (sequenceIn_ >= 0) {
1586    abcDj_[sequenceIn_] = 0.0;
1587    //costBasic_[pivotRow_]=abcCost_[sequenceIn_];
1588  }
1589}
1590static void solveMany(int number, ClpSimplex **simplex)
1591{
1592  for (int i = 0; i < number - 1; i++)
1593    cilk_spawn simplex[i]->dual(0);
1594  simplex[number - 1]->dual(0);
1595  cilk_sync;
1596}
1597void AbcSimplex::crash(int type)
1598{
1599  int i;
1600  for (i = 0; i < numberRows_; i++) {
1601    if (getInternalStatus(i) != basic)
1602      break;
1603  }
1604  if (i < numberRows_)
1605    return;
1606  // all slack
1607  if (type == 3) {
1608    // decomposition crash
1609    if (!abcMatrix_->gotRowCopy())
1610      abcMatrix_->createRowCopy();
1611    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
1612    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
1613    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
1614    //const double * elementByRow = abcMatrix_->rowElements();
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;
1649      }
1650      for (int i = 0; i < numberColumns_; i++) {
1651        columnBlock[i] = -1;
1652        nextColumn[i] = -1;
1653      }
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];
1669#ifndef NDEBUG
1670              int iiColumn = iColumn;
1671#endif
1672              while (jColumn >= 0) {
1673                assert(columnBlock[jColumn] == whichColumnBlock);
1674                columnBlock[jColumn] = iBlock;
1675#ifndef NDEBUG
1676                if (jColumn == iiColumn)
1677                  iiColumn = -1;
1678#endif
1679                iColumn = jColumn;
1680                jColumn = nextColumn[jColumn];
1681              }
1682#ifndef NDEBUG
1683              assert(iiColumn == -1);
1684#endif
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        }
1735      }
1736      if (thisBestBreak == stop)
1737        thisBestValue = COIN_DBL_MAX;
1738      iPass++;
1739      if (iPass == 1) {
1740        bestBreak = thisBestBreak;
1741        bestValue = thisBestValue;
1742      } else {
1743        if (bestValue < thisBestValue) {
1744          firstMaster = 0;
1745          lastMaster = bestBreak;
1746        } else {
1747          firstMaster = thisBestBreak; // ? +1
1748          lastMaster = numberRows_;
1749        }
1750      }
1751    }
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;
1762      while (true) {
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++;
1798      }
1799      // adjust
1800      numberBlocks++;
1801      for (int i = 0; i < numberBlocks; i++) {
1802        blockCount[i] = 0;
1803        nextColumn[i] = 0;
1804      }
1805      for (int iRow = 0; iRow < numberRows_; iRow++) {
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        }
1815      }
1816      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
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        }
1826      }
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        }
1832      }
1833      useful = true;
1834      if (numberMaster > halfway || largestRows * 3 > numberRows_)
1835        useful = false;
1836    }
1837    if (useful) {
1838#if ABC_NORMAL_DEBUG > 0
1839      if (logLevel() >= 2)
1840        printf("%d master rows %d <= < %d\n", lastMaster - firstMaster,
1841          firstMaster, lastMaster);
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",
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]);
1850      }
1851#endif
1852#define NUMBER_DW_BLOCKS 20
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();
1856      // first sort
1857      for (int i = 0; i < numberBlocks; i++) {
1858        backRow[i] = -(3 * blockCount[i] + 0 * nextColumn[i]);
1859        nextColumn[i] = i;
1860      }
1861      CoinSort_2(backRow, backRow + numberBlocks, nextColumn);
1862      // keep if >minSize2 or sum >minSize1;
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++;
1869      }
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        }
1880      }
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;
1888      }
1889      assert(!n);
1890      for (int iRow = 0; iRow < numberRows_; iRow++) {
1891        int iBlock = blockStart[iRow];
1892        if (iBlock >= 0)
1893          blockStart[iRow] = blockCount[iBlock];
1894      }
1895      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1896        int iBlock = columnBlock[iColumn];
1897        if (iBlock >= 0)
1898          columnBlock[iColumn] = blockCount[iBlock];
1899      }
1900      // stick to Clp for now
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);
1917      }
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];
2021      }
2022#if ABC_NORMAL_DEBUG > 0
2023      printf("%d structurals put into basis\n", numberStructurals);
2024#endif
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        }
2035      }
2036      assert(numberBasic == numberRows_);
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
2051      delete[] simplex;
2052      return;
2053    } else {
2054      // try another
2055      type = 2;
2056    }
2057  }
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_;
2068    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2069      gamma = CoinMax(gamma, linearObjective[iColumn]);
2070    }
2071    for (int iRow = 0; iRow < numberRows_; iRow++) {
2072      double lowerBound = abcLower_[iRow];
2073      double upperBound = abcUpper_[iRow];
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;
2084      }
2085    }
2086    int nPossible = 0;
2087    int lastPossible = 0;
2088    double cMaxDiv;
2089    if (gamma)
2090      cMaxDiv = 1.0 / (1000.0 * gamma);
2091    else
2092      cMaxDiv = 1.0;
2093    const double *columnLower = abcLower_ + maximumAbcNumberRows_;
2094    const double *columnUpper = abcUpper_ + maximumAbcNumberRows_;
2095    for (int iPass = 0; iPass < 3; iPass++) {
2096      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
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        }
2132      }
2133      CoinSort_2(q + lastPossible, q + nPossible, which + lastPossible);
2134      lastPossible = nPossible;
2135    }
2136    const double *element = abcMatrix_->getPackedMatrix()->getElements();
2137    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
2138    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
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        }
2161      }
2162      // only really works if scaled
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++;
2181      }
2182    }
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);
2197      }
2198    }
2199#if ABC_NORMAL_DEBUG > 0
2200    printf("%d put into basis\n", nPut);
2201#endif
2202    delete[] q;
2203    delete[] which;
2204  } else if (type == 2) {
2205    //return;
2206    int numberBad = 0;
2207    CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
2208    // Work on savedSolution and current
2209    int iVector = getAvailableArray();
2210#define ALLOW_BAD_DJS
2211#ifdef ALLOW_BAD_DJS
2212    double *modDj = usefulArray_[iVector].denseVector();
2213#endif
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++;
2227#ifdef ALLOW_BAD_DJS
2228            modDj[iSequence] = dj;
2229            dj = 0.0;
2230#else
2231            break;
2232#endif
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++;
2248#ifdef ALLOW_BAD_DJS
2249            modDj[iSequence] = dj;
2250            dj = 0.0;
2251#else
2252            break;
2253#endif
2254          }
2255        } else {
2256          dj = CoinMin(dj, 0.0);
2257        }
2258      } else {
2259        if (fabs(dj) < dualTolerance_) {
2260          dj = 0.0;
2261        } else {
2262          numberBad++;
2263#ifdef ALLOW_BAD_DJS
2264          modDj[iSequence] = dj;
2265          dj = 0.0;
2266#else
2267          break;
2268#endif
2269        }
2270      }
2271      abcDj_[iSequence] = dj;
2272    }
2273#ifndef ALLOW_BAD_DJS
2274    if (numberBad) {
2275      //CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_);
2276      return;
2277    }
2278#endif
2279    if (!abcMatrix_->gotRowCopy())
2280      abcMatrix_->createRowCopy();
2281    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
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) {
2301      // allow more
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;
2315      }
2316      std::sort(index2, index2 + nPossible);
2317      // see how much we need to get numberRows/10 or nPossible/3
2318      average = CoinMax(average, index2[CoinMin(numberRows_ / 10, nPossible / 3)]);
2319      nPossible = 0;
2320    }
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;
2330      }
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;
2336      }
2337    }
2338    CoinSort_2(sort, sort + nPossible, index2);
2339    for (int iWhich = 0; iWhich < nPossible; iWhich++) {
2340      int iRow = index2[iWhich];
2341      sort[iWhich] = 0.0;
2342      if (abcDj_[iRow])
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;
2352      }
2353      assert(infeasibility);
2354      double direction = infeasibility > 0 ? 1.0 : -1.0;
2355      infeasibility *= direction;
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        }
2391      }
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        }
2415      } else {
2416        // free coming in
2417        upperTheta = (abcDj_[whichColumn] * direction) / upperTheta2;
2418      }
2419      // update djs
2420      upperTheta *= -direction;
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);
2431#ifdef ALLOW_BAD_DJS
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            }
2443#endif
2444          } else {
2445            assert(dj < 1.0e-2);
2446            dj = CoinMin(dj, 0.0);
2447#ifdef ALLOW_BAD_DJS
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            }
2459#endif
2460          }
2461        } else if (iStatus < 6) {
2462          assert(fabs(dj) < 1.0e-4);
2463          dj = 0.0;
2464        }
2465        abcDj_[iSequence] = dj;
2466      }
2467      // do basis
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];
2475      } else {
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];
2482      }
2483      setInternalStatus(whichColumn, basic);
2484      abcPivotVariable_[iRow] = whichColumn;
2485      nBasic++;
2486      // mark rows
2487      for (CoinBigIndex j = columnStart[whichColumn]; j < columnStart[whichColumn + 1]; j++) {
2488        int jRow = row[j];
2489        abcDj_[jRow] = 1.0;
2490      }
2491    }
2492#ifdef ALLOW_BAD_DJS
2493    CoinAbcMemset0(modDj + maximumAbcNumberRows_, numberColumns_);
2494#endif
2495    setAvailableArray(iVector);
2496    setAvailableArray(iVector2);
2497#if ABC_NORMAL_DEBUG > 0
2498    printf("dual crash put %d in basis\n", nBasic);
2499#endif
2500  } else {
2501    assert((stateOfProblem_ & VALUES_PASS2) != 0);
2502    // The idea is to put as many likely variables into basis as possible
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++) {
2510      double dj = djSaved_[iSequence];
2511      double value = solution_[iSequence];
2512      double lower = abcLower_[iSequence];
2513      double upper = abcUpper_[iSequence];
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;
2546      }
2547    }
2548    // set slacks basic
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;
2554    // first put all possible slacks in
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_;
2560      } else {
2561        index[n2++] = iSequence;
2562      }
2563    }
2564    n = n2;
2565    int numberIn = 0;
2566    const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_;
2567    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
2568    const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2569    if (!abcMatrix_->gotRowCopy())
2570      abcMatrix_->createRowCopy();
2571    //const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2572    //const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
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        }
2590      }
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        }
2630      }
2631    }
2632#if ABC_NORMAL_DEBUG > 0
2633    printf("Idiot crash put %d in basis\n", numberIn);
2634#endif
2635    setAvailableArray(iVector);
2636    setAvailableArray(iVector2);
2637    delete[] solution_;
2638    solution_ = NULL;
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*/
2645void AbcSimplex::putStuffInBasis(int type)
2646{
2647  int i;
2648  for (i = 0; i < numberRows_; i++) {
2649    if (getInternalStatus(i) != basic)
2650      break;
2651  }
2652  if (i < numberRows_ && (type & 1) == 0)
2653    return;
2654  int iVector = getAvailableArray();
2655  // Column part is mustColumnIn
2656  int *COIN_RESTRICT mustRowOut = usefulArray_[iVector].getIndices();
2657  if (!abcMatrix_->gotRowCopy())
2658    abcMatrix_->createRowCopy();
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;
2666  // find equality rows with effective nonzero rhs
2667  for (int iRow = 0; iRow < numberRows_; iRow++) {
2668    if (abcUpper_[iRow] > abcLower_[iRow] || getInternalStatus(iRow) != basic) {
2669      mustRowOut[iRow] = -2;
2670      continue;
2671    }
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;
2680      } else {
2681        if (chooseColumn[1] == -1)
2682          chooseColumn[1] = iColumn;
2683        else
2684          chooseColumn[1] = -2;
2685      }
2686    }
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        }
2737      }
2738    }
2739  }
2740  if (nPossible) {
2741#if ABC_NORMAL_DEBUG > 0
2742    printf("%d possible candidates\n", nPossible);
2743#endif
2744    if ((type & 2) == 0) {
2745      // triangular
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;
2757      }
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        }
2774      }
2775      setAvailableArray(iVector2);
2776    } else {
2777      // all
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        }
2785      }
2786    }
2787    // redo pivot array
2788    int numberBasic = 0;
2789    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2790      if (getInternalStatus(iSequence) == basic)
2791        abcPivotVariable_[numberBasic++] = iSequence;
2792    }
2793    assert(numberBasic == numberRows_);
2794  }
2795  setAvailableArray(iVector);
2796}
2797// Computes nonbasic cost and total cost
2798void AbcSimplex::computeObjective()
2799{
2800  sumNonBasicCosts_ = 0.0;
2801  rawObjectiveValue_ = 0.0;
2802  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2803    double value = abcSolution_[iSequence] * abcCost_[iSequence];
2804    rawObjectiveValue_ += value;
2805    if (getInternalStatus(iSequence) != basic)
2806      sumNonBasicCosts_ += value;
2807  }
2808  setClpSimplexObjectiveValue();
2809}
2810// This sets largest infeasibility and most infeasible
2811void AbcSimplex::checkPrimalSolution(bool justBasic)
2812{
2813  rawObjectiveValue_ = CoinAbcInnerProduct(costBasic_, numberRows_, solutionBasic_);
2814  //rawObjectiveValue_ += sumNonBasicCosts_;
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
2822  double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 5.0 * primalTolerance_));
2823  // allow tolerance at least slightly bigger than standard
2824  relaxedTolerance = relaxedTolerance + error;
2825  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2826  const double *COIN_RESTRICT lowerBasic = lowerBasic_;
2827  const double *COIN_RESTRICT upperBasic = upperBasic_;
2828  const double *COIN_RESTRICT solutionBasic = solutionBasic_;
2829  if (justBasic) {
2830    for (int iRow = 0; iRow < numberRows_; iRow++) {
2831      double infeasibility = 0.0;
2832      if (solutionBasic[iRow] > upperBasic[iRow]) {
2833        infeasibility = solutionBasic[iRow] - upperBasic[iRow];
2834      } else if (solutionBasic[iRow] < lowerBasic[iRow]) {
2835        infeasibility = lowerBasic[iRow] - solutionBasic[iRow];
2836      }
2837      if (infeasibility > primalTolerance) {
2838        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2839        if (infeasibility > relaxedTolerance)
2840          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2841        numberPrimalInfeasibilities_++;
2842      }
2843    }
2844  } else {
2845    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2846      double infeasibility = 0.0;
2847      if (abcSolution_[iSequence] > abcUpper_[iSequence]) {
2848        infeasibility = abcSolution_[iSequence] - abcUpper_[iSequence];
2849      } else if (abcSolution_[iSequence] < abcLower_[iSequence]) {
2850        infeasibility = abcLower_[iSequence] - abcSolution_[iSequence];
2851      }
2852      if (infeasibility > primalTolerance) {
2853        //assert (getInternalStatus(iSequence)==basic);
2854        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2855        if (infeasibility > relaxedTolerance)
2856          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2857        numberPrimalInfeasibilities_++;
2858      }
2859    }
2860  }
2861}
2862void AbcSimplex::checkDualSolution()
2863{
2864
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
2872  double error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_));
2873  // allow tolerance at least slightly bigger than standard
2874  double relaxedTolerance = currentDualTolerance_ + error;
2875  // allow bigger tolerance for possible improvement
2876  double possTolerance = 5.0 * relaxedTolerance;
2877  sumOfRelaxedDualInfeasibilities_ = 0.0;
2878  int numberNeedToMove = 0;
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_) {
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        }
2920      } else if (distanceDown > primalTolerance_) {
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        }
2932      }
2933    }
2934  }
2935  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj;
2936  if (algorithm_ < 0 && firstFreeDual >= 0) {
2937    // dual
2938    firstFree_ = firstFreeDual;
2939  } else if (numberSuperBasicWithDj || numberNeedToMove) {
2940    //(abcProgress_.lastIterationNumber(0) <= 0)) {
2941    firstFree_ = firstFreePrimal;
2942    if (!sumDualInfeasibilities_)
2943      sumDualInfeasibilities_ = 1.0e-8;
2944  }
2945}
2946/* This sets sum and number of infeasibilities (Dual and Primal) */
2947void AbcSimplex::checkBothSolutions()
2948{
2949  // old way
2950  checkPrimalSolution(true);
2951  checkDualSolution();
2952}
2953/*
2954  Unpacks one column of the matrix into indexed array
2955*/
2956void AbcSimplex::unpack(CoinIndexedVector &rowArray, int sequence) const
2957{
2958  abcMatrix_->unpack(rowArray, sequence);
2959}
2960// Sets row pivot choice algorithm in dual
2961void AbcSimplex::setDualRowPivotAlgorithm(AbcDualRowPivot &choice)
2962{
2963  delete abcDualRowPivot_;
2964  abcDualRowPivot_ = choice.clone(true);
2965  abcDualRowPivot_->setModel(this);
2966}
2967// Sets column pivot choice algorithm in primal
2968void AbcSimplex::setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot &choice)
2969{
2970  delete abcPrimalColumnPivot_;
2971  abcPrimalColumnPivot_ = choice.clone(true);
2972  abcPrimalColumnPivot_->setModel(this);
2973}
2974void AbcSimplex::setFactorization(AbcSimplexFactorization &factorization)
2975{
2976  if (abcFactorization_)
2977    abcFactorization_->setFactorization(factorization);
2978  else
2979    abcFactorization_ = new AbcSimplexFactorization(factorization,
2980      numberRows_);
2981}
2982
2983// Swaps factorization
2984AbcSimplexFactorization *
2985AbcSimplex::swapFactorization(AbcSimplexFactorization *factorization)
2986{
2987  AbcSimplexFactorization *swap = abcFactorization_;
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*/
2997int AbcSimplex::tightenPrimalBounds()
2998{
2999  int tightenType = 1;
3000  if (maximumIterations() >= 1000000 && maximumIterations() < 1000010)
3001    tightenType = maximumIterations() - 1000000;
3002  if (!tightenType)
3003    return 0;
3004  if (integerType_) {
3005#if ABC_NORMAL_DEBUG > 0
3006    printf("Redo tighten to relax by 1 for integers (and don't be shocked by infeasibility)\n");
3007#endif
3008    return 0;
3009  }
3010  // This needs to work on scaled matrix - then replace
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
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();
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;
3028
3029  double tolerance = primalTolerance();
3030
3031  // A copy of bounds is up at top
3032  translate(0); // move down
3033
3034#define MAXPASS 10
3035
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
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;
3047  }
3048  // Keep lower and upper for rows
3049  int whichArray[5];
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());
3057  // -1 no infinite, -2 more than one infinite >=0 column
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++) {
3062    if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3063      changedListNext[numberToDoNext++] = iRow;
3064      whichInfiniteDown[iRow] = -1;
3065      whichInfiniteUp[iRow] = -1;
3066    } else {
3067      minRhs[iRow] = -COIN_DBL_MAX;
3068      maxRhs[iRow] = COIN_DBL_MAX;
3069      whichInfiniteDown[iRow] = -2;
3070      whichInfiniteUp[iRow] = -2;
3071    }
3072  }
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
3087    numberChanged = 0; // Bounds tightened this pass
3088
3089    if (iPass == MAXPASS)
3090      break;
3091    iPass++;
3092    for (int k = 0; k < numberToDo; k++) {
3093      int iRow = changedList[k];
3094      // possible row
3095      int infiniteUpper = 0;
3096      int infiniteLower = 0;
3097      int firstInfiniteUpper = -1;
3098      int firstInfiniteLower = -1;
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
3106
3107      for (j = rStart; j < rEnd; ++j) {
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        }
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;
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;
3149      else
3150        whichInfiniteDown[iRow] = -2;
3151      if (infiniteUpper == 0)
3152        whichInfiniteUp[iRow] = -1;
3153      else if (infiniteUpper == 1)
3154        whichInfiniteUp[iRow] = firstInfiniteUpper;
3155      else
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
3165      } else {
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 - now) * value;
3314                nowUpper = newBound;
3315              }
3316            }
3317            if (upper < large) {
3318              if (!infiniteLower) {
3319                assert(nowUpper < large2);
3320                newBound = nowUpper + (upper - maximumDown) / value;
3321                // relax if original was large
3322                if (fabs(maximumDown) > 1.0e8)
3323                  newBound -= 1.0e-12 * fabs(maximumDown);
3324              } else if (infiniteLower == 1 && nowUpper > large) {
3325                newBound = (upper - maximumDown) / value;
3326                // relax if original was large
3327                if (fabs(maximumDown) > 1.0e8)
3328                  newBound -= 1.0e-12 * fabs(maximumDown);
3329              } else {
3330                newBound = -COIN_DBL_MAX;
3331              }
3332              if (newBound > nowLower + 1.0e-12 && newBound > -large) {
3333                // Tighten the lower bound
3334                numberChanged++;
3335                // check infeasible (relaxed)
3336                if (nowUpper < newBound) {
3337                  if (nowUpper - newBound < -100.0 * tolerance)
3338                    numberInfeasible++;
3339                  else
3340                    newBound = nowUpper;
3341                }
3342                columnLower[iColumn] = newBound;
3343                for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3344                  int jRow = row[j1];
3345                  if (!changed[jRow]) {
3346                    changed[jRow] = 1;
3347                    changedListNext[numberToDoNext++] = jRow;
3348                  }
3349                }
3350                // adjust
3351                double now;
3352                if (nowLower < -large) {
3353                  now = 0.0;
3354                  infiniteUpper--;
3355                } else {
3356                  now = nowLower;
3357                }
3358                maximumUp += (newBound - now) * value;
3359                nowLower = newBound;
3360              }
3361            }
3362          }
3363        }
3364      }
3365    }
3366#if 0
3367    // see if we can tighten doubletons with infinite bounds
3368    if (iPass==1) {
3369      const double * COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
3370      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
3371        CoinBigIndex start = columnStart[iColumn];
3372        if (start+2==columnStart[iColumn+1]&&
3373          (columnLower[iColumn]<-1.0e30||columnUpper[iColumn])) {
3374          int iRow0=row[start];
3375          int iRow1=row[start+1];
3376          printf("col %d row0 %d el0 %g whichDown0 %d whichUp0 %d row1 %d el1 %g whichDown1 %d whichUp1 %d\n",iColumn,
3377                 iRow0,elementByColumn[start],whichInfiniteDown[iRow0],whichInfiniteUp[iRow0],
3378                 iRow1,elementByColumn[start+1],whichInfiniteDown[iRow1],whichInfiniteUp[iRow1]);
3379        }
3380      }
3381    }
3382#endif
3383    totalTightened += numberChanged;
3384    if (iPass == 1)
3385      numberCheck = numberChanged >> 4;
3386    if (numberInfeasible)
3387      break;
3388  }
3389  const double *COIN_RESTRICT saveLower = abcLower_ + maximumNumberTotal_ + maximumAbcNumberRows_;
3390  const double *COIN_RESTRICT saveUpper = abcUpper_ + maximumNumberTotal_ + maximumAbcNumberRows_;
3391  if (!numberInfeasible) {
3392    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_)
3393      << totalTightened
3394      << CoinMessageEol;
3395    int numberLowerChanged = 0;
3396    int numberUpperChanged = 0;
3397    if (!integerType_) {
3398      // Set bounds slightly loose
3399      // tighten rows as well
3400      //#define PRINT_TIGHTEN_PROGRESS 0
3401      for (int iRow = 0; iRow < numberRows_; iRow++) {
3402        assert(maxRhs[iRow] >= rowLower[iRow]);
3403        assert(minRhs[iRow] <= rowUpper[iRow]);
3404        rowLower[iRow] = CoinMax(rowLower[iRow], minRhs[iRow]);
3405        rowUpper[iRow] = CoinMin(rowUpper[iRow], maxRhs[iRow]);
3406#if PRINT_TIGHTEN_PROGRESS > 3
3407        if (handler_->logLevel() > 5)
3408          printf("Row %d min %g max %g lower %g upper %g\n",
3409            iRow, minRhs[iRow], maxRhs[iRow], rowLower[iRow], rowUpper[iRow]);
3410#endif
3411      }
3412#if 0
3413    double useTolerance = 1.0e-4;
3414    double multiplier = 100.0;
3415    // To do this we need to compute min and max JUST for no costs?
3416    const double * COIN_RESTRICT cost = abcCost_+maximumAbcNumberRows_;
3417    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3418      if (saveUpper[iColumn] > saveLower[iColumn] + useTolerance) {
3419        double awayFromLower=0.0;
3420        double awayFromUpper=0.0;
3421        //double gap=columnUpper[iColumn]-columnLower[iColumn];
3422        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3423          int iRow=row[j];
3424          double value=elementByColumn[j];
3425#if PRINT_TIGHTEN_PROGRESS > 3
3426          if (handler_->logLevel()>5)
3427            printf("row %d element %g\n",iRow,value);
3428#endif
3429          if (value>0.0) {
3430            if (minRhs[iRow]+value*awayFromLower<rowLower[iRow]) {
3431              if (minRhs[iRow]>-1.0e10&&awayFromLower<1.0e10)
3432                awayFromLower = CoinMax(awayFromLower,(rowLower[iRow]-minRhs[iRow])/value);
3433              else
3434                awayFromLower=COIN_DBL_MAX;
3435            }
3436            if (maxRhs[iRow]-value*awayFromUpper>rowUpper[iRow]) {
3437              if (maxRhs[iRow]<1.0e10&&awayFromUpper<1.0e10)
3438                awayFromUpper = CoinMax(awayFromUpper,(maxRhs[iRow]-rowUpper[iRow])/value);
3439              else
3440                awayFromUpper=COIN_DBL_MAX;
3441            }
3442          } else {
3443            if (maxRhs[iRow]+value*awayFromLower>rowUpper[iRow]) {
3444              if (maxRhs[iRow]<1.0e10&&awayFromLower<1.0e10)
3445                awayFromLower = CoinMax(awayFromLower,(rowUpper[iRow]-maxRhs[iRow])/value);
3446              else
3447                awayFromLower=COIN_DBL_MAX;
3448            }
3449            if (minRhs[iRow]-value*awayFromUpper<rowLower[iRow]) {
3450              if (minRhs[iRow]>-1.0e10&&awayFromUpper<1.0e10)
3451                awayFromUpper = CoinMin(awayFromUpper,(minRhs[iRow]-rowLower[iRow])/value);
3452              else
3453                awayFromUpper=COIN_DBL_MAX;
3454            }
3455          }
3456        }
3457        // Might have to go as high as
3458        double upper = CoinMin(columnLower[iColumn]+awayFromLower,columnUpper[iColumn]);
3459        // and as low as
3460        double lower = CoinMax(columnUpper[iColumn]-awayFromUpper,columnLower[iColumn]);
3461        // but be sensible
3462        double gap=0.999*(CoinMin(columnUpper[iColumn]-columnLower[iColumn],1.0e10));
3463        if (awayFromLower>gap||awayFromUpper>gap) {
3464          if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3465#if PRINT_TIGHTEN_PROGRESS > 1
3466            printf("YYchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3467                   columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3468#endif
3469            upper=columnUpper[iColumn];
3470          }
3471          if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3472#if PRINT_TIGHTEN_PROGRESS > 1
3473            printf("YYchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3474                   columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3475#endif
3476            lower=columnLower[iColumn];
3477          }
3478        }
3479        if (lower>upper) {
3480#if PRINT_TIGHTEN_PROGRESS
3481          printf("XXchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3482                 columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3483          printf("XXchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3484                 columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3485#endif
3486          lower=columnLower[iColumn];
3487          upper=columnUpper[iColumn];
3488        }
3489        //upper=CoinMax(upper,0.0);
3490        //upper=CoinMax(upper,0.0);
3491        if (cost[iColumn]>=0.0&&awayFromLower<1.0e10&&columnLower[iColumn]>-1.0e10) {
3492          // doesn't want to be too positive
3493          if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3494#if PRINT_TIGHTEN_PROGRESS > 2
3495            if (handler_->logLevel()>1)
3496              printf("change upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3497                     columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3498#endif
3499            double difference=columnUpper[iColumn]-upper;
3500            for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3501              int iRow=row[j];
3502              double value=elementByColumn[j];
3503              if (value>0.0) {
3504                assert (difference*value>=0.0);
3505                maxRhs[iRow]-=difference*value;
3506              } else {
3507                assert (difference*value<=0.0);
3508                minRhs[iRow]-=difference*value;
3509              }
3510            }
3511            columnUpper[iColumn]=upper;
3512            numberUpperChanged++;
3513          }
3514        }
3515        if (cost[iColumn]<=0.0&&awayFromUpper<1.0e10&&columnUpper[iColumn]<1.0e10) {
3516          // doesn't want to be too negative
3517          if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3518#if PRINT_TIGHTEN_PROGRESS > 2
3519            if (handler_->logLevel()>1)
3520              printf("change lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3521                     columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3522#endif
3523            double difference=columnLower[iColumn]-lower;
3524            for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3525              int iRow=row[j];
3526              double value=elementByColumn[j];
3527              if (value>0.0) {
3528                assert (difference*value<=0.0);
3529                minRhs[iRow]-=difference*value;
3530              } else {
3531                assert (difference*value>=0.0);
3532                maxRhs[iRow]-=difference*value;
3533              }
3534            }
3535            columnLower[iColumn]=lower;
3536            numberLowerChanged++;
3537          }
3538        }
3539        // Make large bounds stay infinite
3540        if (saveUpper[iColumn] > 1.0e30 && columnUpper[iColumn] > 1.0e10) {
3541          columnUpper[iColumn] = COIN_DBL_MAX;
3542        }
3543        if (saveLower[iColumn] < -1.0e30 && columnLower[iColumn] < -1.0e10) {
3544          columnLower[iColumn] = -COIN_DBL_MAX;
3545        }
3546        if (columnUpper[iColumn] - columnLower[iColumn] < useTolerance + 1.0e-8) {
3547          // relax enough so will have correct dj
3548          columnLower[iColumn] = CoinMax(saveLower[iColumn],
3549                                          columnLower[iColumn] - multiplier * useTolerance);
3550          columnUpper[iColumn] = CoinMin(saveUpper[iColumn],
3551                                          columnUpper[iColumn] + multiplier * useTolerance);
3552        } else {
3553          if (columnUpper[iColumn] < saveUpper[iColumn]) {
3554            // relax a bit
3555            columnUpper[iColumn] = CoinMin(columnUpper[iColumn] + multiplier * useTolerance,
3556                                            saveUpper[iColumn]);
3557          }
3558          if (columnLower[iColumn] > saveLower[iColumn]) {
3559            // relax a bit
3560            columnLower[iColumn] = CoinMax(columnLower[iColumn] - multiplier * useTolerance,
3561                                            saveLower[iColumn]);
3562          }
3563        }
3564        if (0) {
3565          // debug
3566          // tighten rows as well
3567          for (int iRow = 0; iRow < numberRows_; iRow++) {
3568            if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3569              // possible row
3570              int infiniteUpper = 0;
3571              int infiniteLower = 0;
3572              double maximumUp = 0.0;
3573              double maximumDown = 0.0;
3574              CoinBigIndex rStart = rowStart[iRow];
3575              CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3576              CoinBigIndex j;
3577              // Compute possible lower and upper ranges
3578              for (j = rStart; j < rEnd; ++j) {
3579                double value = element[j];
3580                int iColumn = column[j];
3581                if (value > 0.0) {
3582                  if (columnUpper[iColumn] >= large) {
3583                    ++infiniteUpper;
3584                  } else {
3585                    maximumUp += columnUpper[iColumn] * value;
3586                  }
3587                  if (columnLower[iColumn] <= -large) {
3588                    ++infiniteLower;
3589                  } else {
3590                    maximumDown += columnLower[iColumn] * value;
3591                  }
3592                } else if (value < 0.0) {
3593                  if (columnUpper[iColumn] >= large) {
3594                    ++infiniteLower;
3595                  } else {
3596                    maximumDown += columnUpper[iColumn] * value;
3597                  }
3598                  if (columnLower[iColumn] <= -large) {
3599                    ++infiniteUpper;
3600                  } else {
3601                    maximumUp += columnLower[iColumn] * value;
3602                  }
3603                }
3604              }
3605              // Build in a margin of error
3606              double maxUp = maximumUp+1.0e-8 * fabs(maximumUp);
3607              double maxDown = maximumDown-1.0e-8 * fabs(maximumDown);
3608              double rLower=rowLower[iRow];
3609              double rUpper=rowUpper[iRow];
3610              if (!infiniteUpper&&maxUp < rUpper - tolerance) {
3611                if (rUpper!=COIN_DBL_MAX||maximumUp<1.0e5)
3612                  rUpper=maximumUp;
3613              }
3614              if (!infiniteLower&&maxDown > rLower + tolerance) {
3615                if (rLower!=-COIN_DBL_MAX||maximumDown>-1.0e5)
3616                  rLower=maximumDown;
3617              }
3618              if (infiniteUpper)
3619                maxUp=COIN_DBL_MAX;
3620              if (fabs(maxUp-maxRhs[iRow])>1.0e-3*(1.0+fabs(maxUp)))
3621                printf("bad Maxup row %d maxUp %g maxRhs %g\n",iRow,maxUp,maxRhs[iRow]);
3622              maxRhs[iRow]=maxUp;
3623              if (infiniteLower)
3624                maxDown=-COIN_DBL_MAX;
3625              if (fabs(maxDown-minRhs[iRow])>1.0e-3*(1.0+fabs(maxDown)))
3626                printf("bad Maxdown row %d maxDown %g minRhs %g\n",iRow,maxDown,minRhs[iRow]);
3627              minRhs[iRow]=maxDown;
3628              assert(rLower<=rUpper);
3629            }
3630          }
3631            // end debug
3632        }
3633      }
3634    }
3635    if (tightenType>1) {
3636      // relax
3637      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3638        // relax enough so will have correct dj
3639        double lower=saveLower[iColumn];
3640        double upper=saveUpper[iColumn];
3641        columnLower[iColumn] = CoinMax(lower,columnLower[iColumn] - multiplier);
3642        columnUpper[iColumn] = CoinMin(upper,columnUpper[iColumn] + multiplier);
3643      }
3644    }
3645#endif
3646    } else {
3647#if ABC_NORMAL_DEBUG > 0
3648      printf("Redo tighten to relax by 1 for integers\n");
3649#endif
3650    }
3651#if ABC_NORMAL_DEBUG > 0
3652    printf("Tighten - phase1 %d bounds changed, phase2 %d lower, %d upper\n",
3653      totalTightened, numberLowerChanged, numberUpperChanged);
3654#endif
3655    if (integerType_) {
3656      const double *columnScale = scaleToExternal_ + maximumAbcNumberRows_;
3657      const double *inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
3658      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3659        if (integerType_[iColumn]) {
3660          double value;
3661          double lower = columnLower[iColumn] * inverseColumnScale[iColumn];
3662          double upper = columnUpper[iColumn] * inverseColumnScale[iColumn];
3663          value = floor(lower + 0.5);
3664          if (fabs(value - lower) > primalTolerance_)
3665            value = ceil(lower);
3666          columnLower_[iColumn] = value;
3667          columnLower[iColumn] = value * columnScale[iColumn];
3668          value = floor(upper + 0.5);
3669          if (fabs(value - upper) > primalTolerance_)
3670            value = floor(upper);
3671          columnUpper_[iColumn] = value;
3672          columnUpper[iColumn] = value * columnScale[iColumn];
3673          if (columnLower_[iColumn] > columnUpper_[iColumn])
3674            numberInfeasible++;
3675        }
3676      }
3677      if (numberInfeasible) {
3678        handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3679          << numberInfeasible
3680          << CoinMessageEol;
3681        // restore column bounds
3682        CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3683        CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3684      }
3685    }
3686  } else {
3687    handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3688      << numberInfeasible
3689      << CoinMessageEol;
3690    // restore column bounds
3691    CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3692    CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3693  }
3694  if (!numberInfeasible) {
3695    // tighten rows as well
3696    for (int iRow = 0; iRow < numberRows_; iRow++) {
3697      if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3698        // possible row
3699        int infiniteUpper = 0;
3700        int infiniteLower = 0;
3701        double maximumUp = 0.0;
3702        double maximumDown = 0.0;
3703        CoinBigIndex rStart = rowStart[iRow];
3704        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3705        CoinBigIndex j;
3706        // Compute possible lower and upper ranges
3707        for (j = rStart; j < rEnd; ++j) {
3708          double value = element[j];
3709          int iColumn = column[j];
3710          if (value > 0.0) {
3711            if (columnUpper[iColumn] >= large) {
3712              ++infiniteUpper;
3713            } else {
3714              maximumUp += columnUpper[iColumn] * value;
3715            }
3716            if (columnLower[iColumn] <= -large) {
3717              ++infiniteLower;
3718            } else {
3719              maximumDown += columnLower[iColumn] * value;
3720            }
3721          } else if (value < 0.0) {
3722            if (columnUpper[iColumn] >= large) {
3723              ++infiniteLower;
3724            } else {
3725              maximumDown += columnUpper[iColumn] * value;
3726            }
3727            if (columnLower[iColumn] <= -large) {
3728              ++infiniteUpper;
3729            } else {
3730              maximumUp += columnLower[iColumn] * value;
3731            }
3732          }
3733        }
3734        // Build in a margin of error
3735        double maxUp = maximumUp + 1.0e-8 * fabs(maximumUp);
3736        double maxDown = maximumDown - 1.0e-8 * fabs(maximumDown);
3737        if (!infiniteUpper && maxUp < rowUpper[iRow] - tolerance) {
3738          if (rowUpper[iRow] != COIN_DBL_MAX || maximumUp < 1.0e5)
3739            rowUpper[iRow] = maximumUp;
3740        }
3741        if (!infiniteLower && maxDown > rowLower[iRow] + tolerance) {
3742          if (rowLower[iRow] != -COIN_DBL_MAX || maximumDown > -1.0e5)
3743            rowLower[iRow] = maximumDown;
3744        }
3745        assert(rowLower[iRow] <= rowUpper[iRow]);
3746        if (rowUpper[iRow] - rowLower[iRow] < 1.0e-12) {
3747          if (fabs(rowLower[iRow]) > fabs(rowUpper[iRow]))
3748            rowLower[iRow] = rowUpper[iRow];
3749          else
3750            rowUpper[iRow] = rowLower[iRow];
3751        }
3752      }
3753    }
3754    // flip back
3755    for (int iRow = 0; iRow < numberRows_; iRow++) {
3756      double value = -rowLower[iRow];
3757      rowLower[iRow] = -rowUpper[iRow];
3758      rowUpper[iRow] = value;
3759    }
3760    // move back - relaxed unless integer
3761    if (!integerType_) {
3762      double *COIN_RESTRICT lowerSaved = abcLower_ + maximumNumberTotal_;
3763      double *COIN_RESTRICT upperSaved = abcUpper_ + maximumNumberTotal_;
3764      double tolerance = 2.0 * primalTolerance_;
3765      for (int i = 0; i < numberTotal_; i++) {
3766        double lower = abcLower_[i];
3767        double upper = abcUpper_[i];
3768        if (lower != upper) {
3769          lower = CoinMax(lower - tolerance, lowerSaved[i]);
3770          upper = CoinMin(upper + tolerance, upperSaved[i]);
3771        }
3772        abcLower_[i] = lower;
3773        lowerSaved[i] = lower;
3774        abcUpper_[i] = upper;
3775        upperSaved[i] = upper;
3776      }
3777    } else {
3778      CoinAbcMemcpy(abcLower_ + maximumNumberTotal_, abcLower_, numberTotal_);
3779      CoinAbcMemcpy(abcUpper_ + maximumNumberTotal_, abcUpper_, numberTotal_);
3780    }
3781#if 0
3782    for (int i=0;i<numberTotal_;i++) {
3783      if (abcSolution_[i+maximumNumberTotal_]>abcUpper_[i+maximumNumberTotal_]+1.0e-5)
3784        printf ("above %d %g %g\n",i,
3785                abcSolution_[i+maximumNumberTotal_],abcUpper_[i+maximumNumberTotal_]);
3786      if (abcSolution_[i+maximumNumberTotal_]<abcLower_[i+maximumNumberTotal_]-1.0e-5)
3787        printf ("below %d %g %g\n",i,
3788                abcSolution_[i+maximumNumberTotal_],abcLower_[i+maximumNumberTotal_]);
3789    }
3790#endif
3791    permuteBasis();
3792    // move solution
3793    CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_);
3794  }
3795  CoinAbcMemset0(minRhs, numberRows_);
3796  CoinAbcMemset0(maxRhs, numberRows_);
3797  for (int i = 0; i < 5; i++)
3798    setAvailableArray(whichArray[i]);
3799  return (numberInfeasible);
3800}
3801// dual
3802#include "AbcSimplexDual.hpp"
3803#ifdef ABC_INHERIT
3804// Returns 0 if dual can be skipped
3805int ClpSimplex::doAbcDual()
3806{
3807  if ((abcState_ & CLP_ABC_WANTED) == 0) {
3808    return 1; // not wanted
3809  } else {
3810    assert(abcSimplex_);
3811    int crashState = (abcState_ >> 8) & 7;
3812    if (crashState && crashState < 4) {
3813      switch (crashState) {
3814      case 1:
3815        crash(1000.0, 1);
3816        break;
3817      case 2:
3818        crash(abcSimplex_->dualBound(), 0);
3819        break;
3820      case 3:
3821        crash(0.0, 3);
3822        break;
3823      }
3824    }
3825    // this and abcSimplex_ are approximately same pointer
3826    if ((abcState_ & CLP_ABC_FULL_DONE) == 0) {
3827      abcSimplex_->gutsOfResize(numberRows_, numberColumns_);
3828      abcSimplex_->translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
3829      //abcSimplex_->permuteIn();
3830      int maximumPivotsAbc = CoinMin(abcSimplex_->factorization()->maximumPivots(), numberColumns_ + 200);
3831      abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
3832      if (numberColumns_)
3833        abcSimplex_->tightenPrimalBounds();
3834    }
3835    if ((abcSimplex_->stateOfProblem() & VALUES_PASS2) != 0) {
3836      if (solution_)
3837        crashState = 6;
3838    }
3839    if (crashState && crashState > 3) {
3840      abcSimplex_->crash(crashState - 2);
3841    }
3842    if (crashState && crashState < 4) {
3843      abcSimplex_->putStuffInBasis(1 + 2);
3844    }
3845    int returnCode = abcSimplex_->doAbcDual();
3846    //set to force crossover test returnCode=1;
3847    abcSimplex_->permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
3848    if (returnCode) {
3849      // fix as this model is all messed up e.g. scaling
3850      scalingFlag_ = abs(scalingFlag_);
3851      //delete [] rowScale_;
3852      //delete [] columnScale_;
3853      rowScale_ = NULL;
3854      columnScale_ = NULL;
3855#if 0
3856      delete [] savedRowScale_;
3857      delete [] savedColumnScale_;
3858      savedRowScale_ = NULL;
3859      savedColumnScale_ = NULL;
3860      inverseRowScale_ = NULL;
3861      inverseColumnScale_ = NULL;
3862#endif
3863      if (perturbation_ == 50)
3864        perturbation_ = 51;
3865        //perturbation_=100;
3866#if ABC_NORMAL_DEBUG > 0
3867      printf("Bad exit with status of %d\n", problemStatus_);
3868#endif
3869      //problemStatus_=10;
3870      int status = problemStatus_;
3871      //problemStatus_=-1;
3872      if (status != 10) {
3873        dual(); // Do ClpSimplexDual
3874      } else {
3875        moreSpecialOptions_ |= 256;
3876        primal(1);
3877      }
3878    } else {
3879      // double check objective
3880      //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
3881      perturbation_ = 100;
3882      moreSpecialOptions_ |= 256;
3883      //primal(1);
3884    }
3885    return returnCode;
3886  }
3887}
3888#endif
3889#include "AbcSimplexPrimal.hpp"
3890// Do dual (return 1 if cleanup needed)
3891int AbcSimplex::doAbcDual()
3892{
3893  if (objective_) {
3894    objective_->setActivated(0);
3895  } else {
3896    // create dummy stuff
3897    assert(!numberColumns_);
3898    if (!numberRows_)
3899      problemStatus_ = 0; // say optimal
3900    return 0;
3901  }
3902  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
3903      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
3904      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
3905      additional data and have no destructor or (non-default) constructor.
3906     
3907      This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
3908      in primal and dual.
3909     
3910      As far as I can see this is perfectly safe.
3911  */
3912  algorithm_ = -1;
3913  baseIteration_ = numberIterations_;
3914  if (!abcMatrix_->gotRowCopy())
3915    abcMatrix_->createRowCopy();
3916#ifdef EARLY_FACTORIZE
3917  if (maximumIterations() > 1000000 && maximumIterations() < 1000999) {
3918    numberEarly_ = maximumIterations() - 1000000;
3919#if ABC_NORMAL_DEBUG > 0
3920    printf("Setting numberEarly_ to %d\n", numberEarly_);
3921#endif
3922    numberEarly_ += numberEarly_ << 16;
3923  }
3924#endif
3925  minimumThetaMovement_ = 1.0e-10;
3926  if (fabs(infeasibilityCost_ - 1.0e10) < 1.1e9
3927    && fabs(infeasibilityCost_ - 1.0e10) > 1.0e5 && false) {
3928    minimumThetaMovement_ = 1.0e-3 / fabs(infeasibilityCost_ - 1.0e10);
3929    printf("trying minimum theta movement of %g\n", minimumThetaMovement_);
3930    infeasibilityCost_ = 1.0e10;
3931  }
3932  int savePerturbation = perturbation_;
3933  static_cast< AbcSimplexDual * >(this)->dual();
3934  minimumThetaMovement_ = 1.0e-10;
3935  if ((specialOptions_ & 2048) != 0 && problemStatus_ == 10 && !numberPrimalInfeasibilities_
3936    && sumDualInfeasibilities_ < 1000.0 * dualTolerance_)
3937    problemStatus_ = 0; // ignore
3938  onStopped(); // set secondary status if stopped
3939  if (problemStatus_ == 1 && numberFlagged_) {
3940    problemStatus_ = 10;
3941    static_cast< AbcSimplexPrimal * >(this)->unflag();
3942  }
3943  if (perturbation_ < 101) {
3944    //printf("Perturbed djs?\n");
3945    double *save = abcCost_;
3946    abcCost_ = costSaved_;
3947    computeInternalObjectiveValue();
3948    abcCost_ = save;
3949  }
3950  int totalNumberIterations = numberIterations_;
3951  if (problemStatus_ == 10 && (moreSpecialOptions_ & 32768) != 0 && sumDualInfeasibilities_ < 0.1) {
3952    problemStatus_ = 0;
3953  }
3954#ifndef TRY_ABC_GUS
3955  if (problemStatus_ == 10) {
3956    int savePerturbation = perturbation_;
3957    perturbation_ = 100;
3958    copyFromSaved(2);
3959    minimumThetaMovement_ = 1.0e-12;
3960    baseIteration_ = numberIterations_;
3961    static_cast< AbcSimplexPrimal * >(this)->primal(0);
3962    totalNumberIterations += numberIterations_ - baseIteration_;
3963    perturbation_ = savePerturbation;
3964  }
3965#else
3966  int iPass = 0;
3967  while (problemStatus_ == 10 && minimumThetaMovement_ > 0.99999e-12) {
3968    iPass++;
3969    if (initialSumInfeasibilities_ >= 0.0) {
3970      if (savePerturbation >= 100) {
3971        perturbation_ = 100;
3972      } else {
3973        if (savePerturbation == 50)
3974          perturbation_ = CoinMax(51, HEAVY_PERTURBATION - 4 - iPass); //smaller
3975        else
3976          perturbation_ = CoinMax(51, savePerturbation - 1 - iPass); //smaller
3977      }
3978    } else {
3979      perturbation_ = savePerturbation;
3980      abcFactorization_->setPivots(100000); // force factorization
3981      initialSumInfeasibilities_ = 1.23456789e-5;
3982      // clean pivots
3983      int numberBasic = 0;
3984      int *pivot = pivotVariable();
3985      for (int i = 0; i < numberTotal_; i++) {
3986        if (getInternalStatus(i) == basic)
3987          pivot[numberBasic++] = i;
3988      }
3989      assert(numberBasic == numberRows_);
3990    }
3991    copyFromSaved(14);
3992    // Say second call
3993    if (iPass > 3)
3994      moreSpecialOptions_ |= 256;
3995    if (iPass > 2)
3996      perturbation_ = 101;
3997    baseIteration_ = numberIterations_;
3998    static_cast< AbcSimplexPrimal * >(this)->primal(0);
3999    totalNumberIterations += numberIterations_ - baseIteration_;
4000    if (problemStatus_ == 10) {
4001      // reduce
4002      minimumThetaMovement_ *= 1.0e-1;
4003      copyFromSaved(14);
4004      baseIteration_ = numberIterations_;
4005      static_cast< AbcSimplexDual * >(this)->dual();
4006      totalNumberIterations += numberIterations_ - baseIteration_;
4007    }
4008    perturbation_ = savePerturbation;
4009  }
4010#endif
4011  numberIterations_ = totalNumberIterations;
4012  return (problemStatus_ < 0 || problemStatus_ == 4 || problemStatus_ == 10) ? 1 : 0;
4013}
4014int AbcSimplex::dual()
4015{
4016  if (!abcMatrix_->gotRowCopy())
4017    abcMatrix_->createRowCopy();
4018  assert(!numberFlagged_);
4019  baseIteration_ = numberIterations_;
4020  //double savedPivotTolerance = abcFactorization_->pivotTolerance();
4021  if (objective_) {
4022    objective_->setActivated(0);
4023  } else {
4024    // create dummy stuff
4025    assert(!numberColumns_);
4026    if (!numberRows_)
4027      problemStatus_ = 0; // say optimal
4028    return 0;
4029  }
4030  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
4031      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
4032      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
4033      additional data and have no destructor or (non-default) constructor.
4034     
4035      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4036      in primal and dual.
4037     
4038      As far as I can see this is perfectly safe.
4039  */
4040  minimumThetaMovement_ = 1.0e-12;
4041  int returnCode = static_cast< AbcSimplexDual * >(this)->dual();
4042  if ((specialOptions_ & 2048) != 0 && problemStatus_ == 10 && !numberPrimalInfeasibilities_
4043    && sumDualInfeasibilities_ < 1000.0 * dualTolerance_ && perturbation_ >= 100)
4044    problemStatus_ = 0; // ignore
4045  if (problemStatus_ == 10) {
4046#if 1 //ABC_NORMAL_DEBUG>0
4047    printf("return and use ClpSimplexPrimal\n");
4048    abort();
4049#else
4050    //printf("Cleaning up with primal\n");
4051    //lastAlgorithm=1;
4052    int savePerturbation = perturbation_;
4053    int saveLog = handler_->logLevel();
4054    //handler_->setLogLevel(63);
4055    perturbation_ = 101;
4056    bool denseFactorization = initialDenseFactorization();
4057    // It will be safe to allow dense
4058    setInitialDenseFactorization(true);
4059    // Allow for catastrophe
4060    int saveMax = intParam_[ClpMaxNumIteration];
4061    if (numberIterations_) {
4062      // normal
4063      if (intParam_[ClpMaxNumIteration] > 100000 + numberIterations_)
4064        intParam_[ClpMaxNumIteration]
4065          = numberIterations_ + 1000 + 2 * numberRows_ + numberColumns_;
4066    } else {
4067      // Not normal allow more
4068      baseIteration_ += 2 * (numberRows_ + numberColumns_);
4069    }
4070    // check which algorithms allowed
4071    int dummy;
4072    if (problemStatus_ == 10 && saveObjective == objective_)
4073      startFinishOptions |= 2;
4074    baseIteration_ = numberIterations_;
4075    // Say second call
4076    moreSpecialOptions_ |= 256;
4077    minimumThetaMovement_ = 1.0e-12;
4078    returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(1, startFinishOptions);
4079    // Say not second call
4080    moreSpecialOptions_ &= ~256;
4081    baseIteration_ = 0;
4082    if (saveObjective != objective_) {
4083      // We changed objective to see if infeasible
4084      delete objective_;
4085      objective_ = saveObjective;
4086      if (!problemStatus_) {
4087        // carry on
4088        returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(1, startFinishOptions);
4089      }
4090    }
4091    if (problemStatus_ == 3 && numberIterations_ < saveMax) {
4092      // flatten solution and try again
4093      for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4094        if (getInternalStatus(iSequence) != basic) {
4095          setInternalStatus(iSequence, superBasic);
4096          // but put to bound if close
4097          if (fabs(rowActivity_[iSequence] - rowLower_[iSequence])
4098            <= primalTolerance_) {
4099            abcSolution_[iSequence] = abcLower_[iSequence];
4100            setInternalStatus(iSequence, atLowerBound);
4101          } else if (fabs(rowActivity_[iSequence] - rowUpper_[iSequence])
4102            <= primalTolerance_) {
4103            abcSolution_[iSequence] = abcUpper_[iSequence];
4104            setInternalStatus(iSequence, atUpperBound);
4105          }
4106        }
4107      }
4108      problemStatus_ = -1;
4109      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 2 * numberRows_ + numberColumns_, saveMax);
4110      perturbation_ = savePerturbation;
4111      baseIteration_ = numberIterations_;
4112      // Say second call
4113      moreSpecialOptions_ |= 256;
4114      returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(0, startFinishOptions);
4115      // Say not second call
4116      moreSpecialOptions_ &= ~256;
4117      baseIteration_ = 0;
4118      computeObjectiveValue();
4119      // can't rely on djs either
4120      memset(reducedCost_, 0, numberColumns_ * sizeof(double));
4121    }
4122    intParam_[ClpMaxNumIteration] = saveMax;
4123
4124    setInitialDenseFactorization(denseFactorization);
4125    perturbation_ = savePerturbation;
4126    if (problemStatus_ == 10) {
4127      if (!numberPrimalInfeasibilities_)
4128        problemStatus_ = 0;
4129      else
4130        problemStatus_ = 4;
4131    }
4132    handler_->setLogLevel(saveLog);
4133#endif
4134  }
4135  onStopped(); // set secondary status if stopped
4136  return returnCode;
4137}
4138#ifdef ABC_INHERIT
4139// primal
4140// Returns 0 if primal can be skipped
4141int ClpSimplex::doAbcPrimal(int ifValuesPass)
4142{
4143  if ((abcState_ & CLP_ABC_WANTED) == 0) {
4144    return 1; // not wanted
4145  } else {
4146    assert(abcSimplex_);
4147    if (ifValuesPass)
4148      abcSimplex_->setStateOfProblem(abcSimplex_->stateOfProblem() | VALUES_PASS);
4149    int crashState = (abcState_ >> 8) & 7;
4150    if (crashState && crashState < 4) {
4151      switch (crashState) {
4152      case 1:
4153        crash(1000.0, 1);
4154        break;
4155      case 2:
4156        crash(abcSimplex_->dualBound(), 0);
4157        break;
4158      case 3:
4159        crash(0.0, 3);
4160        break;
4161      }
4162    }
4163    // this and abcSimplex_ are approximately same pointer
4164    if ((abcState_ & CLP_ABC_FULL_DONE) == 0) {
4165      abcSimplex_->gutsOfResize(numberRows_, numberColumns_);
4166      abcSimplex_->translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
4167      //abcSimplex_->permuteIn();
4168      int maximumPivotsAbc = CoinMin(abcSimplex_->factorization()->maximumPivots(), numberColumns_ + 200);
4169      abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
4170      abcSimplex_->copyFromSaved(1);
4171    }
4172    if ((abcSimplex_->stateOfProblem() & VALUES_PASS2) != 0) {
4173      if (solution_)
4174        crashState = 6;
4175    }
4176    if (crashState && crashState > 3) {
4177      abcSimplex_->crash(crashState - 2);
4178    }
4179    if (crashState && crashState < 4) {
4180      abcSimplex_->putStuffInBasis(1 + 2);
4181    }
4182    int returnCode = abcSimplex_->doAbcPrimal(ifValuesPass);
4183    //set to force crossover test returnCode=1;
4184    abcSimplex_->permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
4185    if (returnCode && problemStatus_ < 3) {
4186      // fix as this model is all messed up e.g. scaling
4187      scalingFlag_ = abs(scalingFlag_);
4188      //delete [] rowScale_;
4189      //delete [] columnScale_;
4190      rowScale_ = NULL;
4191      columnScale_ = NULL;
4192#if 0
4193     delete [] savedRowScale_;
4194     delete [] savedColumnScale_;
4195     savedRowScale_ = NULL;
4196     savedColumnScale_ = NULL;
4197     inverseRowScale_ = NULL;
4198     inverseColumnScale_ = NULL;
4199#endif
4200      if (perturbation_ == 50)
4201        perturbation_ = 51;
4202        //perturbation_=100;
4203#if ABC_NORMAL_DEBUG > 0
4204      if (problemStatus_)
4205        printf("Bad exit with status of %d\n", problemStatus_);
4206#endif
4207      //problemStatus_=10;
4208      int status = problemStatus_;
4209      // should not be using
4210      for (int i = 0; i < 6; i++) {
4211        if (rowArray_[i])
4212          rowArray_[i]->clear();
4213        if (columnArray_[i])
4214          columnArray_[i]->clear();
4215      }
4216      //problemStatus_=-1;
4217      if (status < 3 && status != 0) {
4218        if (status != 10) {
4219          primal(); // Do ClpSimplexPrimal
4220        } else {
4221          moreSpecialOptions_ |= 256;
4222          dual();
4223        }
4224      }
4225    } else {
4226      // double check objective
4227      //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
4228    }
4229    numberIterations_ = abcSimplex_->numberIterations();
4230    return returnCode;
4231  }
4232}
4233#endif
4234// Do primal (return 1 if cleanup needed)
4235int AbcSimplex::doAbcPrimal(int ifValuesPass)
4236{
4237  if (objective_) {
4238    objective_->setActivated(0);
4239  } else {
4240    // create dummy stuff
4241    assert(!numberColumns_);
4242    if (!numberRows_)
4243      problemStatus_ = 0; // say optimal
4244    return 0;
4245  }
4246  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
4247      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
4248      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
4249      additional data and have no destructor or (non-default) constructor.
4250     
4251      This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
4252      in primal and dual.
4253     
4254      As far as I can see this is perfectly safe.
4255  */
4256  algorithm_ = -1;
4257  int savePerturbation = perturbation_;
4258  baseIteration_ = numberIterations_;
4259  if (!abcMatrix_->gotRowCopy())
4260    abcMatrix_->createRowCopy();
4261#ifdef EARLY_FACTORIZE
4262  if (maximumIterations() > 1000000 && maximumIterations() < 1000999) {
4263    numberEarly_ = maximumIterations() - 1000000;
4264#if ABC_NORMAL_DEBUG > 0
4265    printf("Setting numberEarly_ to %d\n", numberEarly_);
4266#endif
4267    numberEarly_ += numberEarly_ << 16;
4268  }
4269#endif
4270  // add values pass options
4271  minimumThetaMovement_ = 1.0e-12;
4272  if ((specialOptions_ & 8192) != 0)
4273    minimumThetaMovement_ = 1.0e-15;
4274  int returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(ifValuesPass);
4275  stateOfProblem_ &= ~VALUES_PASS;
4276#ifndef TRY_ABC_GUS
4277  if (returnCode == 10) {
4278    copyFromSaved(14);
4279    int savePerturbation = perturbation_;
4280    perturbation_ = 51;
4281    minimumThetaMovement_ = 1.0e-12;
4282    returnCode = static_cast< AbcSimplexDual * >(this)->dual();
4283    perturbation_ = savePerturbation;
4284  }
4285#else
4286  minimumThetaMovement_ = 1.0e-12;
4287  int iPass = 0;
4288  while (problemStatus_ == 10 && minimumThetaMovement_ > 1.0e-15) {
4289    iPass++;
4290    if (minimumThetaMovement_ == 1.0e-12)
4291      perturbation_ = CoinMin(savePerturbation, 55);
4292    else
4293      perturbation_ = 100;
4294    copyFromSaved(14);
4295    // reduce ?
4296    // Say second call
4297    // Say second call
4298    if (iPass > 3)
4299      moreSpecialOptions_ |= 256;
4300    if (iPass > 2)
4301      perturbation_ = 101;
4302    baseIteration_ = numberIterations_;
4303    // make sure no superbasic
4304    cleanStatus();
4305    if ((specialOptions_ & 8192) == 0)
4306      static_cast< AbcSimplexDual * >(this)->dual();
4307    baseIteration_ = numberIterations_;
4308    if (problemStatus_ == 10) {
4309      if (initialSumInfeasibilities_ >= 0.0) {
4310        if (savePerturbation >= 100) {
4311          perturbation_ = 100;
4312        } else {
4313          if (savePerturbation == 50)
4314            perturbation_ = CoinMax(51, HEAVY_PERTURBATION - 4 - iPass); //smaller
4315          else
4316            perturbation_ = CoinMax(51, savePerturbation - 1 - iPass); //smaller
4317        }
4318      } else {
4319        specialOptions_ |= 8192; // stop going to dual
4320        perturbation_ = savePerturbation;
4321        abcFactorization_->setPivots(100000); // force factorization
4322        initialSumInfeasibilities_ = 1.23456789e-5;
4323        // clean pivots
4324        int numberBasic = 0;
4325        int *pivot = pivotVariable();
4326        for (int i = 0; i < numberTotal_; i++) {
4327          if (getInternalStatus(i) == basic)
4328            pivot[numberBasic++] = i;
4329        }
4330        assert(numberBasic == numberRows_);
4331      }
4332      if (iPass > 2)
4333        perturbation_ = 100;
4334      copyFromSaved(14);
4335      baseIteration_ = numberIterations_;
4336      static_cast< AbcSimplexPrimal * >(this)->primal(0);
4337    }
4338    minimumThetaMovement_ *= 0.5;
4339    perturbation_ = savePerturbation;
4340  }
4341#endif
4342  return returnCode;
4343}
4344/* Sets up all slack basis and resets solution to
4345   as it was after initial load or readMps */
4346void AbcSimplex::allSlackBasis()
4347{
4348  // set column status to one nearest zero
4349  CoinFillN(internalStatus_, numberRows_, static_cast< unsigned char >(basic));
4350  const double *COIN_RESTRICT lower = abcLower_;
4351  const double *COIN_RESTRICT upper = abcUpper_;
4352  double *COIN_RESTRICT solution = abcSolution_;
4353  // But set value to zero if lb <0.0 and ub>0.0
4354  for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
4355    if (lower[i] >= 0.0) {
4356      solution[i] = lower[i];
4357      setInternalStatus(i, atLowerBound);
4358    } else if (upper[i] <= 0.0) {
4359      solution[i] = upper[i];
4360      setInternalStatus(i, atUpperBound);
4361    } else if (lower[i] < -1.0e20 && upper[i] > 1.0e20) {
4362      // free
4363      solution[i] = 0.0;
4364      setInternalStatus(i, isFree);
4365    } else if (fabs(lower[i]) < fabs(upper[i])) {
4366      solution[i] = 0.0;
4367      setInternalStatus(i, atLowerBound);
4368    } else {
4369      solution[i] = 0.0;
4370      setInternalStatus(i, atUpperBound);
4371    }
4372  }
4373}
4374#if 0 //ndef SLIM_NOIO
4375// Read an mps file from the given filename
4376int
4377AbcSimplex::readMps(const char *filename,
4378                    bool keepNames,
4379                    bool ignoreErrors)
4380{
4381  int status = ClpSimplex::readMps(filename, keepNames, ignoreErrors);
4382  ClpSimplex * thisModel=this;
4383  gutsOfResize(thisModel->numberRows(),thisModel->numberColumns());
4384  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4385  return status;
4386}
4387// Read GMPL files from the given filenames
4388int
4389AbcSimplex::readGMPL(const char *filename, const char * dataName,
4390                     bool keepNames)
4391{
4392  int status = ClpSimplex::readGMPL(filename, dataName, keepNames);
4393  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4394  return status;
4395}
4396// Read file in LP format from file with name filename.
4397int
4398AbcSimplex::readLp(const char *filename, const double epsilon )
4399{
4400  FILE *fp = fopen(filename, "r");
4401 
4402  if(!fp) {
4403    printf("### ERROR: AbcSimplex::readLp():  Unable to open file %s for reading\n",
4404           filename);
4405    return(1);
4406  }
4407  CoinLpIO m;
4408  m.readLp(fp, epsilon);
4409  fclose(fp);
4410 
4411  // set problem name
4412  setStrParam(ClpProbName, m.getProblemName());
4413  // no errors
4414  loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
4415              m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
4416 
4417  if (m.integerColumns()) {
4418    integerType_ = new char[numberColumns_];
4419    CoinMemcpyN(m.integerColumns(), numberColumns_, integerType_);
4420  } else {
4421    integerType_ = NULL;
4422  }
4423  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4424  unsigned int maxLength = 0;
4425  int iRow;
4426  rowNames_ = std::vector<std::string> ();
4427  columnNames_ = std::vector<std::string> ();
4428  rowNames_.reserve(numberRows_);
4429  for (iRow = 0; iRow < numberRows_; iRow++) {
4430    const char * name = m.rowName(iRow);
4431    if (name) {
4432      maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4433      rowNames_.push_back(name);
4434    } else {
4435      rowNames_.push_back("");
4436    }
4437  }
4438 
4439  int iColumn;
4440  columnNames_.reserve(numberColumns_);
4441  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4442    const char * name = m.columnName(iColumn);
4443    if (name) {
4444      maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4445      columnNames_.push_back(name);
4446    } else {
4447      columnNames_.push_back("");
4448    }
4449  }
4450  lengthNames_ = static_cast<int> (maxLength);
4451 
4452  return 0;
4453}
4454#endif
4455// Factorization frequency
4456int AbcSimplex::factorizationFrequency() const
4457{
4458  if (abcFactorization_)
4459    return abcFactorization_->maximumPivots();
4460  else
4461    return -1;
4462}
4463void AbcSimplex::setFactorizationFrequency(int value)
4464{
4465  if (abcFactorization_)
4466    abcFactorization_->maximumPivots(value);
4467}
4468
4469/* Get a clean factorization - i.e. throw out singularities
4470   may do more later */
4471int AbcSimplex::cleanFactorization(int ifValuesPass)
4472{
4473  int status = internalFactorize(ifValuesPass ? 10 : 0);
4474  if (status < 0) {
4475    return 1; // some error
4476  } else {
4477    firstFree_ = 0;
4478    return 0;
4479  }
4480}
4481// Save data
4482ClpDataSave
4483AbcSimplex::saveData()
4484{
4485  ClpDataSave saved;
4486  saved.dualBound_ = dualBound_;
4487  saved.infeasibilityCost_ = infeasibilityCost_;
4488  saved.pivotTolerance_ = abcFactorization_->pivotTolerance();
4489  saved.zeroFactorizationTolerance_ = abcFactorization_->zeroTolerance();
4490  saved.zeroSimplexTolerance_ = zeroTolerance_;
4491  saved.perturbation_ = perturbation_;
4492  saved.forceFactorization_ = forceFactorization_;
4493  saved.acceptablePivot_ = acceptablePivot_;
4494  saved.objectiveScale_ = objectiveScale_;
4495  // Progress indicator
4496  abcProgress_.fillFromModel(this);
4497  return saved;
4498}
4499// Restore data
4500void AbcSimplex::restoreData(ClpDataSave saved)
4501{
4502  //abcFactorization_->sparseThreshold(saved.sparseThreshold_);
4503  abcFactorization_->pivotTolerance(saved.pivotTolerance_);
4504  abcFactorization_->zeroTolerance(saved.zeroFactorizationTolerance_);
4505  zeroTolerance_ = saved.zeroSimplexTolerance_;
4506  perturbation_ = saved.perturbation_;
4507  infeasibilityCost_ = saved.infeasibilityCost_;
4508  dualBound_ = saved.dualBound_;
4509  forceFactorization_ = saved.forceFactorization_;
4510  objectiveScale_ = saved.objectiveScale_;
4511  acceptablePivot_ = saved.acceptablePivot_;
4512}
4513// To flag a variable
4514void AbcSimplex::setFlagged(int sequence)
4515{
4516  assert(sequence >= 0 && sequence < numberTotal_);
4517  if (sequence >= 0) {
4518    internalStatus_[sequence] |= 64;
4519    // use for real disaster lastFlaggedIteration_ = numberIterations_;
4520    numberFlagged_++;
4521  }
4522}
4523#ifndef NDEBUG
4524// For errors to make sure print to screen
4525// only called in debug mode
4526static void indexError(int index,
4527  std::string methodName)
4528{
4529  std::cerr << "Illegal index " << index << " in AbcSimplex::" << methodName << std::endl;
4530  throw CoinError("Illegal index", methodName, "AbcSimplex");
4531}
4532#endif
4533// After modifying first copy refreshes second copy and marks as updated
4534void AbcSimplex::refreshCosts()
4535{
4536  whatsChanged_ &= ~OBJECTIVE_SAME;
4537  CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
4538}
4539void AbcSimplex::refreshLower(unsigned int type)
4540{
4541  if (type)
4542    whatsChanged_ &= type;
4543  CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_);
4544}
4545void AbcSimplex::refreshUpper(unsigned int type)
4546{
4547  if (type)
4548    whatsChanged_ &= type;
4549  CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_);
4550}
4551/* Set an objective function coefficient */
4552void AbcSimplex::setObjectiveCoefficient(int elementIndex, double elementValue)
4553{
4554#ifndef NDEBUG
4555  if (elementIndex < 0 || elementIndex >= numberColumns_) {
4556    indexError(elementIndex, "setObjectiveCoefficient");
4557  }
4558#endif
4559  if (objective()[elementIndex] != elementValue) {
4560    objective()[elementIndex] = elementValue;
4561    // work arrays exist - update as well
4562    whatsChanged_ &= ~OBJECTIVE_SAME;
4563    double direction = optimizationDirection_ * objectiveScale_;
4564    costSaved_[elementIndex + maximumAbcNumberRows_] = direction * elementValue
4565      * columnUseScale_[elementIndex + maximumAbcNumberRows_];
4566  }
4567}
4568/* Set a single row lower bound<br>
4569   Use -DBL_MAX for -infinity. */
4570void AbcSimplex::setRowLower(int elementIndex, double elementValue)
4571{
4572#ifndef NDEBUG
4573  int n = numberRows_;
4574  if (elementIndex < 0 || elementIndex >= n) {
4575    indexError(elementIndex, "setRowLower");
4576  }
4577#endif
4578  if (elementValue < -1.0e27)
4579    elementValue = -COIN_DBL_MAX;
4580  if (rowLower_[elementIndex] != elementValue) {
4581    rowLower_[elementIndex] = elementValue;
4582    // work arrays exist - update as well
4583    whatsChanged_ &= ~ROW_LOWER_SAME;
4584    if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4585      lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4586    } else {
4587      lowerSaved_[elementIndex] = elementValue * rhsScale_
4588        * rowUseScale_[elementIndex];
4589    }
4590  }
4591}
4592
4593/* Set a single row upper bound<br>
4594   Use DBL_MAX for infinity. */
4595void AbcSimplex::setRowUpper(int elementIndex, double elementValue)
4596{
4597#ifndef NDEBUG
4598  int n = numberRows_;
4599  if (elementIndex < 0 || elementIndex >= n) {
4600    indexError(elementIndex, "setRowUpper");
4601  }
4602#endif
4603  if (elementValue > 1.0e27)
4604    elementValue = COIN_DBL_MAX;
4605  if (rowUpper_[elementIndex] != elementValue) {
4606    rowUpper_[elementIndex] = elementValue;
4607    // work arrays exist - update as well
4608    whatsChanged_ &= ~ROW_UPPER_SAME;
4609    if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4610      upperSaved_[elementIndex] = COIN_DBL_MAX;
4611    } else {
4612      upperSaved_[elementIndex] = elementValue * rhsScale_
4613        * rowUseScale_[elementIndex];
4614    }
4615  }
4616}
4617
4618/* Set a single row lower and upper bound */
4619void AbcSimplex::setRowBounds(int elementIndex,
4620  double lowerValue, double upperValue)
4621{
4622#ifndef NDEBUG
4623  int n = numberRows_;
4624  if (elementIndex < 0 || elementIndex >= n) {
4625    indexError(elementIndex, "setRowBounds");
4626  }
4627#endif
4628  if (lowerValue < -1.0e27)
4629    lowerValue = -COIN_DBL_MAX;
4630  if (upperValue > 1.0e27)
4631    upperValue = COIN_DBL_MAX;
4632  //CoinAssert (upperValue>=lowerValue);
4633  if (rowLower_[elementIndex] != lowerValue) {
4634    rowLower_[elementIndex] = lowerValue;
4635    // work arrays exist - update as well
4636    whatsChanged_ &= ~ROW_LOWER_SAME;
4637    if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4638      lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4639    } else {
4640      lowerSaved_[elementIndex] = lowerValue * rhsScale_
4641        * rowUseScale_[elementIndex];
4642    }
4643  }
4644  if (rowUpper_[elementIndex] != upperValue) {
4645    rowUpper_[elementIndex] = upperValue;
4646    // work arrays exist - update as well
4647    whatsChanged_ &= ~ROW_UPPER_SAME;
4648    if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4649      upperSaved_[elementIndex] = COIN_DBL_MAX;
4650    } else {
4651      upperSaved_[elementIndex] = upperValue * rhsScale_
4652        * rowUseScale_[elementIndex];
4653    }
4654  }
4655}
4656void AbcSimplex::setRowSetBounds(const int *indexFirst,
4657  const int *indexLast,
4658  const double *boundList)
4659{
4660#ifndef NDEBUG
4661  int n = numberRows_;
4662#endif
4663  int numberChanged = 0;
4664  const int *saveFirst = indexFirst;
4665  while (indexFirst != indexLast) {
4666    const int iRow = *indexFirst++;
4667#ifndef NDEBUG
4668    if (iRow < 0 || iRow >= n) {
4669      indexError(iRow, "setRowSetBounds");
4670    }
4671#endif
4672    double lowerValue = *boundList++;
4673    double upperValue = *boundList++;
4674    if (lowerValue < -1.0e27)
4675      lowerValue = -COIN_DBL_MAX;
4676    if (upperValue > 1.0e27)
4677      upperValue = COIN_DBL_MAX;
4678    //CoinAssert (upperValue>=lowerValue);
4679    if (rowLower_[iRow] != lowerValue) {
4680      rowLower_[iRow] = lowerValue;
4681      whatsChanged_ &= ~ROW_LOWER_SAME;
4682      numberChanged++;
4683    }
4684    if (rowUpper_[iRow] != upperValue) {
4685      rowUpper_[iRow] = upperValue;
4686      whatsChanged_ &= ~ROW_UPPER_SAME;
4687      numberChanged++;
4688    }
4689  }
4690  if (numberChanged) {
4691    indexFirst = saveFirst;
4692    while (indexFirst != indexLast) {
4693      const int iRow = *indexFirst++;
4694      if (rowLower_[iRow] == -COIN_DBL_MAX) {
4695        lowerSaved_[iRow] = -COIN_DBL_MAX;
4696      } else {
4697        lowerSaved_[iRow] = rowLower_[iRow] * rhsScale_
4698          * rowUseScale_[iRow];
4699      }
4700      if (rowUpper_[iRow] == COIN_DBL_MAX) {
4701        upperSaved_[iRow] = COIN_DBL_MAX;
4702      } else {
4703        upperSaved_[iRow] = rowUpper_[iRow] * rhsScale_
4704          * rowUseScale_[iRow];
4705      }
4706    }
4707  }
4708}
4709//-----------------------------------------------------------------------------
4710/* Set a single column lower bound<br>
4711   Use -DBL_MAX for -infinity. */
4712void AbcSimplex::setColumnLower(int elementIndex, double elementValue)
4713{
4714#ifndef NDEBUG
4715  int n = numberColumns_;
4716  if (elementIndex < 0 || elementIndex >= n) {
4717    indexError(elementIndex, "setColumnLower");
4718  }
4719#endif
4720  if (elementValue < -1.0e27)
4721    elementValue = -COIN_DBL_MAX;
4722  if (columnLower_[elementIndex] != elementValue) {
4723    columnLower_[elementIndex] = elementValue;
4724    // work arrays exist - update as well
4725    whatsChanged_ &= ~COLUMN_LOWER_SAME;
4726    double value;
4727    if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4728      value = -COIN_DBL_MAX;
4729    } else {
4730      value = elementValue * rhsScale_
4731        * inverseColumnUseScale_[elementIndex];
4732    }
4733    lowerSaved_[elementIndex + maximumAbcNumberRows_] = value;
4734  }
4735}
4736
4737/* Set a single column upper bound<br>
4738   Use DBL_MAX for infinity. */
4739void AbcSimplex::setColumnUpper(int elementIndex, double elementValue)
4740{
4741#ifndef NDEBUG
4742  int n = numberColumns_;
4743  if (elementIndex < 0 || elementIndex >= n) {
4744    indexError(elementIndex, "setColumnUpper");
4745  }
4746#endif
4747  if (elementValue > 1.0e27)
4748    elementValue = COIN_DBL_MAX;
4749  if (columnUpper_[elementIndex] != elementValue) {
4750    columnUpper_[elementIndex] = elementValue;
4751    // work arrays exist - update as well
4752    whatsChanged_ &= ~COLUMN_UPPER_SAME;
4753    double value;
4754    if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4755      value = COIN_DBL_MAX;
4756    } else {
4757      value = elementValue * rhsScale_
4758        * inverseColumnUseScale_[elementIndex];
4759    }
4760    upperSaved_[elementIndex + maximumAbcNumberRows_] = value;
4761  }
4762}
4763
4764/* Set a single column lower and upper bound */
4765void AbcSimplex::setColumnBounds(int elementIndex,
4766  double lowerValue, double upperValue)
4767{
4768#ifndef NDEBUG
4769  int n = numberColumns_;
4770  if (elementIndex < 0 || elementIndex >= n) {
4771    indexError(elementIndex, "setColumnBounds");
4772  }
4773#endif
4774  if (lowerValue < -1.0e27)
4775    lowerValue = -COIN_DBL_MAX;
4776  if (columnLower_[elementIndex] != lowerValue) {
4777    columnLower_[elementIndex] = lowerValue;
4778    // work arrays exist - update as well
4779    whatsChanged_ &= ~COLUMN_LOWER_SAME;
4780    if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4781      lowerSaved_[elementIndex + maximumAbcNumberRows_] = -COIN_DBL_MAX;
4782    } else {
4783      lowerSaved_[elementIndex + maximumAbcNumberRows_] = lowerValue * rhsScale_
4784        * inverseColumnUseScale_[elementIndex];
4785    }
4786  }
4787  if (upperValue > 1.0e27)
4788    upperValue = COIN_DBL_MAX;
4789  if (columnUpper_[elementIndex] != upperValue) {
4790    columnUpper_[elementIndex] = upperValue;
4791    // work arrays exist - update as well
4792    whatsChanged_ &= ~COLUMN_UPPER_SAME;
4793    if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4794      upperSaved_[elementIndex + maximumAbcNumberRows_] = COIN_DBL_MAX;
4795    } else {
4796      upperSaved_[elementIndex + maximumAbcNumberRows_] = upperValue * rhsScale_
4797        * inverseColumnUseScale_[elementIndex];
4798    }
4799  }
4800}
4801void AbcSimplex::setColumnSetBounds(const int *indexFirst,
4802  const int *indexLast,
4803  const double *boundList)
4804{
4805#ifndef NDEBUG
4806  int n = numberColumns_;
4807#endif
4808  int numberChanged = 0;
4809  const int *saveFirst = indexFirst;
4810  while (indexFirst != indexLast) {
4811    const int iColumn = *indexFirst++;
4812#ifndef NDEBUG
4813    if (iColumn < 0 || iColumn >= n) {
4814      indexError(iColumn, "setColumnSetBounds");
4815    }
4816#endif
4817    double lowerValue = *boundList++;
4818    double upperValue = *boundList++;
4819    if (lowerValue < -1.0e27)
4820      lowerValue = -COIN_DBL_MAX;
4821    if (upperValue > 1.0e27)
4822      upperValue = COIN_DBL_MAX;
4823    //CoinAssert (upperValue>=lowerValue);
4824    if (columnLower_[iColumn] != lowerValue) {
4825      columnLower_[iColumn] = lowerValue;
4826      whatsChanged_ &= ~COLUMN_LOWER_SAME;
4827      numberChanged++;
4828    }
4829    if (columnUpper_[iColumn] != upperValue) {
4830      columnUpper_[iColumn] = upperValue;
4831      whatsChanged_ &= ~COLUMN_UPPER_SAME;
4832      numberChanged++;
4833    }
4834  }
4835  if (numberChanged) {
4836    indexFirst = saveFirst;
4837    while (indexFirst != indexLast) {
4838      const int iColumn = *indexFirst++;
4839      if (columnLower_[iColumn] == -COIN_DBL_MAX) {
4840        lowerSaved_[iColumn + maximumAbcNumberRows_] = -COIN_DBL_MAX;
4841      } else {
4842        lowerSaved_[iColumn + maximumAbcNumberRows_] = columnLower_[iColumn] * rhsScale_
4843          * inverseColumnUseScale_[iColumn];
4844      }
4845      if (columnUpper_[iColumn] == COIN_DBL_MAX) {
4846        upperSaved_[iColumn + maximumAbcNumberRows_] = COIN_DBL_MAX;
4847      } else {
4848        upperSaved_[iColumn + maximumAbcNumberRows_] = columnUpper_[iColumn] * rhsScale_
4849          * inverseColumnUseScale_[iColumn];
4850      }
4851    }
4852  }
4853}
4854#ifndef SLIM_CLP
4855#include "CoinWarmStartBasis.hpp"
4856
4857// Returns a basis (to be deleted by user)
4858CoinWarmStartBasis *
4859AbcSimplex::getBasis() const
4860{
4861  CoinWarmStartBasis *basis = new CoinWarmStartBasis();
4862  basis->setSize(numberColumns_, numberRows_);
4863
4864  unsigned char lookup[8] = { 2, 3, 255, 255, 0, 0, 1, 3 };
4865  for (int iRow = 0; iRow < numberRows_; iRow++) {
4866    int iStatus = getInternalStatus(iRow);
4867    iStatus = lookup[iStatus];
4868    basis->setArtifStatus(iRow, static_cast< CoinWarmStartBasis::Status >(iStatus));
4869  }
4870  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4871    int iStatus = getInternalStatus(iColumn + maximumAbcNumberRows_);
4872    iStatus = lookup[iStatus];
4873    basis->setStructStatus(iColumn, static_cast< CoinWarmStartBasis::Status >(iStatus));
4874  }
4875  return basis;
4876}
4877#endif
4878// Compute objective value from solution
4879void AbcSimplex::computeObjectiveValue(bool useInternalArrays)
4880{
4881  const double *obj = objective();
4882  if (!useInternalArrays) {
4883    objectiveValue_ = 0.0;
4884    for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
4885      double value = columnActivity_[iSequence];
4886      objectiveValue_ += value * obj[iSequence];
4887    }
4888    // But remember direction as we are using external objective
4889    objectiveValue_ *= optimizationDirection_;
4890  } else {
4891    rawObjectiveValue_ = 0.0;
4892    for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4893      double value = abcSolution_[iSequence];
4894      rawObjectiveValue_ += value * abcCost_[iSequence];
4895    }
4896    setClpSimplexObjectiveValue();
4897  }
4898}
4899// Compute minimization objective value from internal solution
4900double
4901AbcSimplex::computeInternalObjectiveValue()
4902{
4903  rawObjectiveValue_ = 0.0;
4904  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4905    double value = abcSolution_[iSequence];
4906    rawObjectiveValue_ += value * abcCost_[iSequence];
4907  }
4908  setClpSimplexObjectiveValue();
4909  return objectiveValue_ - dblParam_[ClpObjOffset];
4910}
4911// If user left factorization frequency then compute
4912void AbcSimplex::defaultFactorizationFrequency()
4913{
4914  if (factorizationFrequency() == 200) {
4915    // User did not touch preset
4916    const int cutoff1 = 10000;
4917    const int cutoff2 = 100000;
4918    const int base = 75;
4919    const int freq0 = 50;
4920    const int freq1 = 150;
4921    const int maximum = 10000;
4922    int frequency;
4923    if (numberRows_ < cutoff1)
4924      frequency = base + numberRows_ / freq0;
4925    else
4926      frequency = base + cutoff1 / freq0 + (numberRows_ - cutoff1) / freq1;
4927    setFactorizationFrequency(CoinMin(maximum, frequency));
4928  }
4929}
4930// Gets clean and emptyish factorization
4931AbcSimplexFactorization *
4932AbcSimplex::getEmptyFactorization()
4933{
4934  if ((specialOptions_ & 65536) == 0) {
4935    assert(!abcFactorization_);
4936    abcFactorization_ = new AbcSimplexFactorization();
4937  } else if (!abcFactorization_) {
4938    abcFactorization_ = new AbcSimplexFactorization();
4939  }
4940  return abcFactorization_;
4941}
4942// Resizes rim part of model
4943void AbcSimplex::resize(int newNumberRows, int newNumberColumns)
4944{
4945  assert(maximumAbcNumberRows_ >= 0);
4946  // save
4947  int numberRows = numberRows_;
4948  int numberColumns = numberColumns_;
4949  ClpSimplex::resize(newNumberRows, newNumberColumns);
4950  // back out
4951  numberRows_ = numberRows;
4952  numberColumns_ = numberColumns;
4953  gutsOfResize(newNumberRows, newNumberColumns);
4954}
4955// Return true if the objective limit test can be relied upon
4956bool AbcSimplex::isObjectiveLimitTestValid() const
4957{
4958  if (problemStatus_ == 0) {
4959    return true;
4960  } else if (problemStatus_ == 1) {
4961    // ok if dual
4962    return (algorithm_ < 0);
4963  } else if (problemStatus_ == 2) {
4964    // ok if primal
4965    return (algorithm_ > 0);
4966  } else {
4967    return false;
4968  }
4969}
4970/*
4971  Permutes in from ClpModel data - assumes scale factors done
4972  and AbcMatrix exists but is in original order (including slacks)
4973  For now just add basicArray at end
4974  ==
4975  But could partition into
4976  normal (i.e. reasonable lower/upper)
4977  abnormal - free, odd bounds
4978  fixed
4979  ==
4980  sets a valid pivotVariable
4981  Slacks always shifted by offset
4982  Fixed variables always shifted by offset
4983  Recode to allow row objective so can use pi from idiot etc
4984*/
4985void AbcSimplex::permuteIn()
4986{
4987  // for now only if maximumAbcNumberRows_==numberRows_
4988  //assert(maximumAbcNumberRows_==numberRows_);
4989  numberTotal_ = maximumAbcNumberRows_ + numberColumns_;
4990  double direction = optimizationDirection_;
4991  // all this could use cilk
4992  // move primal stuff to end
4993  const double *COIN_RESTRICT rowScale = scaleFromExternal_;
4994  const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
4995  const double *COIN_RESTRICT columnScale = scaleToExternal_ + maximumAbcNumberRows_;
4996  const double *COIN_RESTRICT inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
4997  double *COIN_RESTRICT offsetRhs = offsetRhs_;
4998  double *COIN_RESTRICT saveLower = lowerSaved_ + maximumAbcNumberRows_;
4999  double *COIN_RESTRICT saveUpper = upperSaved_ + maximumAbcNumberRows_;
5000  double *COIN_RESTRICT saveCost = costSaved_ + maximumAbcNumberRows_;
5001  double *COIN_RESTRICT saveSolution = solutionSaved_ + maximumAbcNumberRows_;
5002  CoinAbcMemset0(offsetRhs, maximumAbcNumberRows_);
5003  const double *COIN_RESTRICT objective = this->objective();
5004  objectiveOffset_ = 0.0;
5005  double *COIN_RESTRICT offset = offset_ + maximumAbcNumberRows_;
5006  const int *COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
5007  const CoinBigIndex *COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
5008  const double *COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
5009  largestGap_ = 1.0e-12;
5010  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5011    double scale = inverseColumnScale[iColumn];
5012    double lowerValue = columnLower_[iColumn];
5013    double upperValue = columnUpper_[iColumn];
5014    double thisOffset = 0.0;
5015#if 1 //ndef CLP_USER_DRIVEN
5016    double lowerValue2 = fabs(lowerValue);
5017    double upperValue2 = fabs(upperValue);
5018    if (CoinMin(lowerValue2, upperValue2) < 1000.0) {
5019      // move to zero
5020      if (lowerValue2 > upperValue2)
5021        thisOffset = upperValue;
5022      else
5023        thisOffset = lowerValue;
5024    }
5025#endif
5026    offset[iColumn] = thisOffset;
5027    if (thisOffset) {
5028      objectiveOffset_ += thisOffset * objective[iColumn] * optimizationDirection_;
5029      double scaledOffset = thisOffset * scale;
5030      for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
5031        int iRow = row[j];
5032        offsetRhs[iRow] += scaledOffset * elementByColumn[j];
5033      }
5034    }
5035    lowerValue -= thisOffset;
5036    if (lowerValue > -1.0e30)
5037      lowerValue *= scale;
5038    saveLower[iColumn] = lowerValue;
5039    upperValue -= thisOffset;
5040    if (upperValue < 1.0e30)
5041      upperValue *= scale;
5042    saveUpper[iColumn] = upperValue;
5043    largestGap_ = CoinMax(largestGap_, upperValue - lowerValue);
5044    saveSolution[iColumn] = scale * (columnActivity_[iColumn] - thisOffset);
5045    saveCost[iColumn] = objective[iColumn] * direction * columnScale[iColumn];
5046  }
5047  CoinAbcMemset0(offset_, maximumAbcNumberRows_);
5048  saveLower -= maximumAbcNumberRows_;
5049  saveUpper -= maximumAbcNumberRows_;
5050  saveCost -= maximumAbcNumberRows_;
5051  saveSolution -= maximumAbcNumberRows_;
5052  for (int iRow = 0; iRow < numberRows_; iRow++) {
5053    // note switch of lower to upper
5054    double upperValue = -rowLower_[iRow];
5055    double lowerValue = -rowUpper_[iRow];
5056    double thisOffset = offsetRhs_[iRow];
5057    double scale = rowScale[iRow];
5058    if (lowerValue > -1.0e30)
5059      lowerValue = lowerValue * scale + thisOffset;
5060    saveLower[iRow] = lowerValue;
5061    if (upperValue < 1.0e30)
5062      upperValue = upperValue * scale + thisOffset;
5063    saveUpper[iRow] = upperValue;
5064    largestGap_ = CoinMax(largestGap_, upperValue - lowerValue);
5065    saveCost[iRow] = 0.0;
5066    dual_[iRow] *= direction * inverseRowScale[iRow];
5067    saveSolution[iRow] = 0.0; // not necessary
5068  }
5069  dualBound_ = CoinMin(dualBound_, largestGap_);
5070  // Compute rhsScale_ and objectiveScale_
5071  double minValue = COIN_DBL_MAX;
5072  double maxValue = 0.0;
5073  CoinAbcMinMaxAbsNormalValues(costSaved_ + maximumAbcNumberRows_, numberTotal_ - maximumAbcNumberRows_, minValue, maxValue);
5074  // scale to 1000.0 ?
5075  if (minValue && false) {
5076    objectiveScale_ = 1000.0 / sqrt(minValue * maxValue);
5077    objectiveScale_ = CoinMin(1.0, 1000.0 / maxValue);
5078#ifndef NDEBUG
5079    double smallestNormal = COIN_DBL_MAX;
5080    double smallestAny = COIN_DBL_MAX;
5081    double largestAny = 0.0;
5082    for (int i = 0; i < numberTotal_; i++) {
5083      double value = fabs(costSaved_[i]);
5084      if (value) {
5085        if (value > 1.0e-8)
5086          smallestNormal = CoinMin(smallestNormal, value);
5087        smallestAny = CoinMin(smallestAny, value);
5088        largestAny = CoinMax(largestAny, value);
5089      }
5090    }
5091    printf("objectiveScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
5092      objectiveScale_, minValue, smallestNormal, smallestAny, maxValue, largestAny);
5093#endif
5094  } else {
5095    //objectiveScale_=1.0;
5096  }
5097  CoinAbcScale(costSaved_, objectiveScale_, numberTotal_);
5098  minValue = COIN_DBL_MAX;
5099  maxValue = 0.0;
5100  CoinAbcMinMaxAbsNormalValues(lowerSaved_, numberTotal_, minValue, maxValue);
5101  CoinAbcMinMaxAbsNormalValues(upperSaved_, numberTotal_, minValue, maxValue);
5102  // scale to 100.0 ?
5103  if (minValue && false) {
5104    rhsScale_ = 100.0 / sqrt(minValue * maxValue);
5105#ifndef NDEBUG
5106    double smallestNormal = COIN_DBL_MAX;
5107    double smallestAny = COIN_DBL_MAX;
5108    double largestAny = 0.0;
5109    for (int i = 0; i < numberTotal_; i++) {
5110      double value = fabs(lowerSaved_[i]);
5111      if (value && value != COIN_DBL_MAX) {
5112        if (value > 1.0e-8)
5113          smallestNormal = CoinMin(smallestNormal, value);
5114        smallestAny = CoinMin(smallestAny, value);
5115        largestAny = CoinMax(largestAny, value);
5116      }
5117    }
5118    for (int i = 0; i < numberTotal_; i++) {
5119      double value = fabs(upperSaved_[i]);
5120      if (value && value != COIN_DBL_MAX) {
5121        if (value > 1.0e-8)
5122          smallestNormal = CoinMin(smallestNormal, value);
5123        smallestAny = CoinMin(smallestAny, value);
5124        largestAny = CoinMax(largestAny, value);
5125      }
5126    }
5127    printf("rhsScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
5128      rhsScale_, minValue, smallestNormal, smallestAny, maxValue, largestAny);
5129#endif
5130  } else {
5131    rhsScale_ = 1.0;
5132  }
5133  CoinAbcScaleNormalValues(lowerSaved_, rhsScale_, 1.0e-13, numberTotal_);
5134  CoinAbcScaleNormalValues(upperSaved_, rhsScale_, 1.0e-13, numberTotal_);
5135  // copy
5136  CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_);
5137  CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_);
5138  CoinAbcMemcpy(abcSolution_, abcSolution_ + maximumNumberTotal_, numberTotal_);
5139  CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
5140}
5141void AbcSimplex::permuteBasis()
5142{
5143  assert(abcPivotVariable_ || (!numberRows_ && !numberColumns_));
5144  int numberBasic = 0;
5145  // from Clp enum to Abc enum (and bound flip)
5146  unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
5147  unsigned char *COIN_RESTRICT getStatus = status_ + numberColumns_;
5148  double *COIN_RESTRICT solutionSaved = solutionSaved_;
5149  double *COIN_RESTRICT lowerSaved = lowerSaved_;
5150  double *COIN_RESTRICT upperSaved = upperSaved_;
5151  bool ordinaryVariables = true;
5152  bool valuesPass = (stateOfProblem_ & VALUES_PASS) != 0;
5153  if (valuesPass) {
5154    // get solution
5155    CoinAbcMemset0(abcSolution_, numberRows_);
5156    abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, abcSolution_);
5157    //double * temp = new double[numberRows_];
5158    //memset(temp,0,numberRows_*sizeof(double));
5159    //matrix_->times(1.0,columnActivity_,temp);
5160    CoinAbcMemcpy(solutionSaved_, abcSolution_, numberRows_);
5161    int n = 0;
5162    for (int i = 0; i < numberTotal_; i++) {
5163      if (solutionSaved_[i] > upperSaved_[i] + 1.0e-5)
5164        n++;
5165      else if (solutionSaved_[i] < lowerSaved_[i] - 1.0e-5)
5166        n++;
5167    }
5168#if ABC_NORMAL_DEBUG > 0
5169    if (n)
5170      printf("%d infeasibilities\n", n);
5171#endif
5172  }
5173  // dual at present does not like superbasic
5174  for (int iRow = 0; iRow < numberRows_; iRow++) {
5175    unsigned char status = getStatus[iRow] & 7;
5176    AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbcSlack[status]);
5177    if (status != ClpSimplex::basic) {
5178      double lowerValue = lowerSaved[iRow];
5179      double upperValue = upperSaved[iRow];
5180      if (lowerValue == -COIN_DBL_MAX) {
5181        if (upperValue == COIN_DBL_MAX) {
5182          // free
5183          abcStatus = isFree;
5184          ordinaryVariables = false;
5185        } else {
5186          abcStatus = atUpperBound;
5187        }
5188      } else if (upperValue == COIN_DBL_MAX) {
5189        abcStatus = atLowerBound;
5190      } else if (lowerValue == upperValue) {
5191        abcStatus = isFixed;
5192      } else if (abcStatus == isFixed) {
5193        double value = solutionSaved[iRow];
5194        if (value - lowerValue < upperValue - value)
5195          abcStatus = atLowerBound;
5196        else
5197          abcStatus = atUpperBound;
5198      }
5199      if (valuesPass) {
5200        double value = solutionSaved[iRow];
5201        if (value > lowerValue + primalTolerance_ && value < upperValue - primalTolerance_ && (abcStatus == atLowerBound || abcStatus == atUpperBound))
5202          abcStatus = superBasic;
5203      }
5204      switch (abcStatus) {
5205      case isFixed:
5206      case atLowerBound:
5207        solutionSaved[iRow] = lowerValue;
5208        break;
5209      case atUpperBound:
5210        solutionSaved[iRow] = upperValue;
5211        break;
5212      default:
5213        ordinaryVariables = false;
5214        break;
5215      }
5216    } else {
5217      // basic
5218      abcPivotVariable_[numberBasic++] = iRow;
5219    }
5220    internalStatus_[iRow] = abcStatus;
5221  }
5222  // from Clp enum to Abc enum
5223  unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
5224  unsigned char *COIN_RESTRICT putStatus = internalStatus_ + maximumAbcNumberRows_;
5225  getStatus = status_;
5226  solutionSaved += maximumAbcNumberRows_;
5227  lowerSaved += maximumAbcNumberRows_;
5228  upperSaved += maximumAbcNumberRows_;
5229  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5230    unsigned char status = getStatus[iColumn] & 7;
5231    AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbc[status]);
5232    if (status != ClpSimplex::basic) {
5233      double lowerValue = lowerSaved[iColumn];
5234      double upperValue = upperSaved[iColumn];
5235      if (lowerValue == -COIN_DBL_MAX) {
5236        if (upperValue == COIN_DBL_MAX) {
5237          // free
5238          abcStatus = isFree;
5239          ordinaryVariables = false;
5240        } else {
5241          abcStatus = atUpperBound;
5242        }
5243      } else if (upperValue == COIN_DBL_MAX) {
5244        abcStatus = atLowerBound;
5245      } else if (lowerValue == upperValue) {
5246        abcStatus = isFixed;
5247      } else if (abcStatus == isFixed) {
5248        double value = solutionSaved[iColumn];
5249        if (value - lowerValue < upperValue - value)
5250          abcStatus = atLowerBound;
5251        else
5252          abcStatus = atUpperBound;
5253      } else if (abcStatus == isFree) {
5254        abcStatus = superBasic;
5255      }
5256      if (valuesPass && (abcStatus == atLowerBound || abcStatus == atUpperBound)) {
5257        double value = solutionSaved[iColumn];
5258        if (value > lowerValue + primalTolerance_) {
5259          if (value < upperValue - primalTolerance_) {
5260            abcStatus = superBasic;
5261          } else {
5262            abcStatus = atUpperBound;
5263          }
5264        } else {
5265          abcStatus = atLowerBound;
5266        }
5267      }
5268      switch (abcStatus) {
5269      case isFixed:
5270      case atLowerBound:
5271        solutionSaved[iColumn] = lowerValue;
5272        break;
5273      case atUpperBound:
5274        solutionSaved[iColumn] = upperValue;
5275        break;
5276      default:
5277        ordinaryVariables = false;
5278        break;
5279      }
5280    } else {
5281      // basic
5282      if (numberBasic < numberRows_)
5283        abcPivotVariable_[numberBasic++] = iColumn + maximumAbcNumberRows_;
5284      else
5285        abcStatus = superBasic;
5286    }
5287    putStatus[iColumn] = abcStatus;
5288  }
5289  ordinaryVariables_ = ordinaryVariables ? 1 : 0;
5290  if (numberBasic < numberRows_) {
5291    for (int iRow = 0; iRow < numberRows_; iRow++) {
5292      AbcSimplex::Status status = getInternalStatus(iRow);
5293      if (status != AbcSimplex::basic) {
5294        abcPivotVariable_[numberBasic++] = iRow;
5295        setInternalStatus(iRow, basic);
5296        if (numberBasic == numberRows_)
5297          break;
5298      }
5299    }
5300  }
5301  // copy
5302  CoinAbcMemcpy(internalStatus_ + maximumNumberTotal_, internalStatus_, numberTotal_);
5303}
5304// Permutes out - bit settings same as stateOfProblem
5305void AbcSimplex::permuteOut(int whatsWanted)
5306{
5307  assert(whatsWanted);
5308  if ((whatsWanted & ALL_STATUS_OK) != 0 && (stateOfProblem_ & ALL_STATUS_OK) == 0) {
5309    stateOfProblem_ |= ALL_STATUS_OK;
5310    // from Abc enum to Clp enum (and bound flip)
5311    unsigned char lookupToClpSlack[8] = { 2, 3, 255, 255, 0, 0, 1, 5 };
5312    unsigned char *COIN_RESTRICT putStatus = status_ + numberColumns_;
5313    const unsigned char *COIN_RESTRICT getStatus = internalStatus_;
5314    for (int iRow = 0; iRow < numberRows_; iRow++) {
5315      putStatus[iRow] = lookupToClpSlack[getStatus[iRow] & 7];
5316    }
5317    // from Abc enum to Clp enum
5318    unsigned char lookupToClp[8] = { 3, 2, 255, 255, 0, 0, 1, 5 };
5319    putStatus = status_;
5320    getStatus = internalStatus_ + maximumAbcNumberRows_;
5321    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5322      putStatus[iColumn] = lookupToClp[getStatus[iColumn] & 7];
5323    }
5324  }
5325  const double *COIN_RESTRICT rowScale = scaleFromExternal_;
5326  const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
5327  const double *COIN_RESTRICT columnScale = scaleToExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5328  const double *COIN_RESTRICT inverseColumnScale = scaleFromExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5329  double scaleC = 1.0 / objectiveScale_;
5330  double scaleR = 1.0 / rhsScale_;
5331  int numberPrimalScaled = 0;
5332  int numberPrimalUnscaled = 0;
5333  int numberDualScaled = 0;
5334  int numberDualUnscaled = 0;
5335  bool filledInSolution = false;
5336  if ((whatsWanted & ROW_PRIMAL_OK) != 0 && (stateOfProblem_ & ROW_PRIMAL_OK) == 0) {
5337    stateOfProblem_ |= ROW_PRIMAL_OK;
5338    if (!filledInSolution) {
5339      filledInSolution = true;
5340      CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
5341      CoinAbcScatterZeroTo(abcDj_, abcPivotVariable_, numberRows_);
5342    }
5343    // Collect infeasibilities
5344    double *COIN_RESTRICT putSolution = rowActivity_;
5345    const double *COIN_RESTRICT rowLower = rowLower_;
5346    const double *COIN_RESTRICT rowUpper = rowUpper_;
5347    const double *COIN_RESTRICT abcLower = abcLower_;
5348    const double *COIN_RESTRICT abcUpper = abcUpper_;
5349    const double *COIN_RESTRICT abcSolution = abcSolution_;
5350    const double *COIN_RESTRICT offsetRhs = offsetRhs_;
5351    for (int i = 0; i < numberRows_; i++) {
5352      double scaleFactor = inverseRowScale[i];
5353      double valueScaled = abcSolution[i];
5354      double lowerScaled = abcLower[i];
5355      double upperScaled = abcUpper[i];
5356      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5357        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
5358          numberPrimalScaled++;
5359        else
5360          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5361      }
5362      double value = (offsetRhs[i] - valueScaled * scaleR) * scaleFactor;
5363      putSolution[i] = value;
5364      if (value < rowLower[i] - primalTolerance_)
5365        numberPrimalUnscaled++;
5366      else if (value > rowUpper[i] + primalTolerance_)
5367        numberPrimalUnscaled++;
5368    }
5369  }
5370  if ((whatsWanted & COLUMN_PRIMAL_OK) != 0 && (stateOfProblem_ & COLUMN_PRIMAL_OK) == 0) {
5371    stateOfProblem_ |= COLUMN_PRIMAL_OK;
5372    // Collect infeasibilities
5373    if (!filledInSolution) {
5374      filledInSolution = true;
5375      CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
5376      CoinAbcScatterZeroTo(abcDj_, abcPivotVariable_, numberRows_);
5377    }
5378    // Collect infeasibilities
5379    double *COIN_RESTRICT putSolution = columnActivity_ - maximumAbcNumberRows_;
5380    const double *COIN_RESTRICT columnLower = columnLower_ - maximumAbcNumberRows_;
5381    const double *COIN_RESTRICT columnUpper = columnUpper_ - maximumAbcNumberRows_;
5382    for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5383      double scaleFactor = columnScale[i];
5384      double valueScaled = abcSolution_[i];
5385      double lowerScaled = abcLower_[i];
5386      double upperScaled = abcUpper_[i];
5387      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5388        if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
5389          numberPrimalScaled++;
5390        else
5391          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5392      }
5393      double value = (valueScaled * scaleR) * scaleFactor + offset_[i];
5394      putSolution[i] = value;
5395      if (value < columnLower[i] - primalTolerance_)
5396        numberPrimalUnscaled++;
5397      else if (value > columnUpper[i] + primalTolerance_)
5398        numberPrimalUnscaled++;
5399    }
5400  }
5401  if ((whatsWanted & COLUMN_DUAL_OK) != 0 && (stateOfProblem_ & COLUMN_DUAL_OK) == 0) {
5402    // Get fixed djs
5403    CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
5404    abcMatrix_->transposeTimesAll(dual_, abcDj_);
5405    stateOfProblem_ |= COLUMN_DUAL_OK;
5406    // Collect infeasibilities
5407    double *COIN_RESTRICT putDj = reducedCost_ - maximumAbcNumberRows_;
5408    for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5409      double scaleFactor = inverseColumnScale[i];
5410      double valueDual = abcDj_[i];
5411      double value = abcSolution_[i];
5412      bool aboveLower = value > abcLower_[i] + primalTolerance_;
5413      bool belowUpper = value < abcUpper_[i] - primalTolerance_;
5414      if (aboveLower && valueDual > dualTolerance_)
5415        numberDualScaled++;
5416      if (belowUpper && valueDual < -dualTolerance_)
5417        numberDualScaled++;
5418      valueDual *= scaleFactor * scaleC;
5419      putDj[i] = valueDual;
5420      if (aboveLower && valueDual > dualTolerance_)
5421        numberDualUnscaled++;
5422      if (belowUpper && valueDual < -dualTolerance_)
5423        numberDualUnscaled++;
5424    }
5425  }
5426  if ((whatsWanted & ROW_DUAL_OK) != 0 && (stateOfProblem_ & ROW_DUAL_OK) == 0) {
5427    stateOfProblem_ |= ROW_DUAL_OK;
5428    // Collect infeasibilities
5429    for (int i = 0; i < numberRows_; i++) {
5430      double scaleFactor = rowScale[i];
5431      double valueDual = abcDj_[i]; // ? +- and direction
5432      double value = abcSolution_[i];
5433      bool aboveLower = value > abcLower_[i] + primalTolerance_;
5434      bool belowUpper = value < abcUpper_[i] - primalTolerance_;
5435      if (aboveLower && valueDual > dualTolerance_)
5436        numberDualScaled++;
5437      if (belowUpper && valueDual < -dualTolerance_)
5438        numberDualScaled++;
5439      valueDual *= scaleFactor * scaleC;
5440      dual_[i] = -(valueDual - abcCost_[i]); // sign
5441      if (aboveLower && valueDual > dualTolerance_)
5442        numberDualUnscaled++;
5443      if (belowUpper && valueDual < -dualTolerance_)
5444        numberDualUnscaled++;
5445    }
5446  }
5447  // do after djs
5448  if (!problemStatus_ && (!secondaryStatus_ || secondaryStatus_ == 2 || secondaryStatus_ == 3)) {
5449    // See if we need to set secondary status
5450    if (numberPrimalUnscaled) {
5451      if (numberDualUnscaled || secondaryStatus_ == 3)
5452        secondaryStatus_ = 4;
5453      else
5454        secondaryStatus_ = 2;
5455    } else if (numberDualUnscaled) {
5456      if (secondaryStatus_ == 0)
5457        secondaryStatus_ = 3;
5458      else
5459        secondaryStatus_ = 4;
5460    }
5461  }
5462  if (scalingFlag_) {
5463    if (problemStatus_ == 2) {
5464      for (int i = 0; i < numberColumns_; i++) {
5465        ray_[i] *= columnScale[i];
5466      }
5467    } else if (problemStatus_ == 1 && ray_) {
5468      for (int i = 0; i < numberRows_; i++) {
5469        ray_[i] *= rowScale[i];
5470      }
5471    }
5472  }
5473}
5474