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

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

fix a few problems with Aboca

File size: 195.6 KB
Line 
1/* $Id: AbcSimplex.cpp 1448 2011-06-19 15:34:41Z 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//#undef NDEBUG
7#if CLP_HAS_ABC
8
9#include "ClpConfig.h"
10
11#include "CoinPragma.hpp"
12#include <math.h>
13//#define ABC_DEBUG 2
14
15#if SLIM_CLP==2
16#define SLIM_NOIO
17#endif
18#include "CoinHelperFunctions.hpp"
19#include "CoinFloatEqual.hpp"
20#include "ClpSimplex.hpp"
21#include "AbcSimplex.hpp"
22#include "AbcSimplexDual.hpp"
23#include "AbcSimplexFactorization.hpp"
24#include "AbcNonLinearCost.hpp"
25#include "CoinAbcCommon.hpp"
26#include "AbcMatrix.hpp"
27#include "CoinIndexedVector.hpp"
28#include "AbcDualRowDantzig.hpp"
29#include "AbcDualRowSteepest.hpp"
30#include "AbcPrimalColumnDantzig.hpp"
31#include "AbcPrimalColumnSteepest.hpp"
32#include "ClpMessage.hpp"
33#include "ClpEventHandler.hpp"
34#include "ClpLinearObjective.hpp"
35#include "CoinAbcHelperFunctions.hpp"
36#include "CoinModel.hpp"
37#include "CoinLpIO.hpp"
38#include <cfloat>
39
40#include <string>
41#include <stdio.h>
42#include <iostream>
43//#############################################################################
44AbcSimplex::AbcSimplex (bool emptyMessages) :
45 
46  ClpSimplex(emptyMessages)
47{
48  gutsOfDelete(0);
49  gutsOfInitialize(0,0,true);
50  assert (maximumAbcNumberRows_>=0);
51  //printf("XX %x AbcSimplex constructor\n",this);
52}
53
54//-----------------------------------------------------------------------------
55
56AbcSimplex::~AbcSimplex ()
57{
58  //printf("XX %x AbcSimplex destructor\n",this);
59  gutsOfDelete(1);
60}
61// Copy constructor.
62AbcSimplex::AbcSimplex(const AbcSimplex &rhs) :
63  ClpSimplex(rhs)
64{
65  //printf("XX %x AbcSimplex constructor from %x\n",this,&rhs);
66  gutsOfDelete(0);
67  gutsOfInitialize(numberRows_,numberColumns_,false);
68  gutsOfCopy(rhs);
69  assert (maximumAbcNumberRows_>=0);
70}
71#include "ClpDualRowSteepest.hpp"
72#include "ClpPrimalColumnSteepest.hpp"
73#include "ClpFactorization.hpp"
74// Copy constructor from model
75AbcSimplex::AbcSimplex(const ClpSimplex &rhs) :
76  ClpSimplex(rhs)
77{
78  //printf("XX %x AbcSimplex constructor from ClpSimplex\n",this);
79  gutsOfDelete(0);
80  gutsOfInitialize(numberRows_,numberColumns_,true);
81  gutsOfResize(numberRows_,numberColumns_);
82  // Set up row/column selection and factorization type
83  ClpDualRowSteepest * pivot = dynamic_cast<ClpDualRowSteepest *>(rhs.dualRowPivot());
84  if (pivot) {
85    AbcDualRowSteepest steep(pivot->mode());
86    setDualRowPivotAlgorithm(steep);
87  } else {
88    AbcDualRowDantzig dantzig;
89    setDualRowPivotAlgorithm(dantzig);
90  }
91  ClpPrimalColumnSteepest * pivotColumn = dynamic_cast<ClpPrimalColumnSteepest *>(rhs.primalColumnPivot());
92  if (pivotColumn) {
93    AbcPrimalColumnSteepest steep(pivotColumn->mode());
94    setPrimalColumnPivotAlgorithm(steep);
95  } else {
96    AbcPrimalColumnDantzig dantzig;
97    setPrimalColumnPivotAlgorithm(dantzig);
98  }
99  //if (rhs.factorization()->)
100  //factorization_->forceOtherFactorization();
101  factorization()->synchronize(rhs.factorization(),this);
102  //factorization_->setGoDenseThreshold(rhs.factorization()->goDenseThreshold());
103  //factorization_->setGoSmallThreshold(rhs.factorization()->goSmallThreshold());
104  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
105  assert (maximumAbcNumberRows_>=0);
106}
107// Assignment operator. This copies the data
108AbcSimplex &
109AbcSimplex::operator=(const AbcSimplex & rhs)
110{
111  if (this != &rhs) {
112    gutsOfDelete(1);
113    ClpSimplex::operator=(rhs);
114    gutsOfCopy(rhs);
115    assert (maximumAbcNumberRows_>=0);
116  }
117  return *this;
118}
119// fills in perturbationSaved_ from start with 0.5+random
120void 
121AbcSimplex::fillPerturbation(int start, int number)
122{
123  double * array = perturbationSaved_+start;
124  for (int i=0;i<number;i++) 
125    array[i] = 0.5+0.5*randomNumberGenerator_.randomDouble();
126}
127// Sets up all extra pointers
128void 
129AbcSimplex::setupPointers(int maxRows,int maxColumns)
130{
131  int numberTotal=maxRows+maxColumns;
132  scaleToExternal_ = scaleFromExternal_+numberTotal;
133  tempArray_=offset_+numberTotal;
134  lowerSaved_=abcLower_+numberTotal;
135  upperSaved_=abcUpper_+numberTotal;
136  costSaved_=abcCost_+numberTotal;
137  solutionSaved_=abcSolution_+numberTotal;
138  djSaved_=abcDj_+numberTotal;
139  internalStatusSaved_= internalStatus_+numberTotal;
140  perturbationSaved_=abcPerturbation_+numberTotal;
141  offsetRhs_=tempArray_+numberTotal;
142  lowerBasic_=lowerSaved_+numberTotal;
143  upperBasic_=upperSaved_+numberTotal;
144  costBasic_=costSaved_+numberTotal;
145  solutionBasic_=solutionSaved_+numberTotal;
146  djBasic_=djSaved_+numberTotal;
147  perturbationBasic_=perturbationSaved_+numberTotal;
148  columnUseScale_ = scaleToExternal_+maxRows;
149  inverseColumnUseScale_ = scaleFromExternal_+maxRows;
150}
151// Copies all saved versions to working versions and may do something for perturbation
152void 
153AbcSimplex::copyFromSaved(int which)
154{
155  if ((which&1)!=0)
156    CoinAbcMemcpy(abcSolution_,solutionSaved_,maximumNumberTotal_);
157  if ((which&2)!=0)
158    CoinAbcMemcpy(abcCost_,costSaved_,maximumNumberTotal_);
159  if ((which&4)!=0)
160    CoinAbcMemcpy(abcLower_,lowerSaved_,maximumNumberTotal_);
161  if ((which&8)!=0)
162    CoinAbcMemcpy(abcUpper_,upperSaved_,maximumNumberTotal_);
163  if ((which&16)!=0)
164    CoinAbcMemcpy(abcDj_,djSaved_,maximumNumberTotal_);
165  if ((which&32)!=0) {
166    CoinAbcMemcpy(abcLower_,abcPerturbation_,numberTotal_);
167    CoinAbcMemcpy(abcUpper_,abcPerturbation_+numberTotal_,numberTotal_);
168  }
169}
170void
171AbcSimplex::gutsOfCopy(const AbcSimplex & rhs)
172{
173  abcSimplex_=this;
174  assert (numberRows_ == rhs.numberRows_);
175  assert (numberColumns_ == rhs.numberColumns_);
176  numberTotal_ = rhs.numberTotal_;
177  maximumNumberTotal_ = rhs.maximumNumberTotal_;
178  // special options here to be safe
179  specialOptions_=rhs.specialOptions_;
180  //assert (maximumInternalRows_ >= numberRows_);
181  //assert (maximumInternalColumns_ >= numberColumns_);
182  // Copy all scalars
183  CoinAbcMemcpy(reinterpret_cast<int *>(&sumNonBasicCosts_),
184            reinterpret_cast<const int *>(&rhs.sumNonBasicCosts_),
185            (&swappedAlgorithm_-reinterpret_cast<int *>(&sumNonBasicCosts_))+1);
186  // could add 2 so can go off end
187  int sizeArray=2*maximumNumberTotal_+maximumAbcNumberRows_;
188  internalStatus_ = ClpCopyOfArray(rhs.internalStatus_,sizeArray+maximumNumberTotal_);
189  abcLower_ = ClpCopyOfArray(rhs.abcLower_, sizeArray);
190  abcUpper_ = ClpCopyOfArray(rhs.abcUpper_, sizeArray);
191  abcCost_ = ClpCopyOfArray(rhs.abcCost_, sizeArray);
192  abcDj_ = ClpCopyOfArray(rhs.abcDj_, sizeArray);
193  abcSolution_ = ClpCopyOfArray(rhs.abcSolution_, sizeArray+maximumNumberTotal_);
194  abcPerturbation_ = ClpCopyOfArray(rhs.abcPerturbation_,sizeArray);
195  abcPivotVariable_ = ClpCopyOfArray(rhs.abcPivotVariable_,maximumAbcNumberRows_);
196  //fromExternal_ = ClpCopyOfArray(rhs.fromExternal_,sizeArray);
197  //toExternal_ = ClpCopyOfArray(rhs.toExternal_,sizeArray);
198  scaleFromExternal_ = ClpCopyOfArray(rhs.scaleFromExternal_,sizeArray);
199  offset_ = ClpCopyOfArray(rhs.offset_,sizeArray);
200  setupPointers(maximumAbcNumberRows_,maximumAbcNumberColumns_);
201  if (rhs.abcMatrix_)
202    abcMatrix_ = new AbcMatrix(*rhs.abcMatrix_);
203  else
204    abcMatrix_ = NULL;
205  for (int i = 0; i < ABC_NUMBER_USEFUL; i++) {
206    usefulArray_[i] = rhs.usefulArray_[i];
207  }
208  if (rhs.abcFactorization_) {
209    setFactorization(*rhs.abcFactorization_);
210  } else {
211    delete abcFactorization_;
212    abcFactorization_ = NULL;
213  }
214#ifdef EARLY_FACTORIZE
215  delete abcEarlyFactorization_;
216  if (rhs.abcEarlyFactorization_) {
217    abcEarlyFactorization_ = new AbcSimplexFactorization(*rhs.abcEarlyFactorization_);
218  } else {
219    abcEarlyFactorization_ = NULL;
220  }
221#endif
222  abcDualRowPivot_ = rhs.abcDualRowPivot_->clone(true);
223  abcDualRowPivot_->setModel(this);
224  abcPrimalColumnPivot_ = rhs.abcPrimalColumnPivot_->clone(true);
225  abcPrimalColumnPivot_->setModel(this);
226  if (rhs.abcBaseModel_) {
227    abcBaseModel_ = new AbcSimplex(*rhs.abcBaseModel_);
228  } else {
229    abcBaseModel_ = NULL;
230  }
231  if (rhs.clpModel_) {
232    clpModel_ = new ClpSimplex(*rhs.clpModel_);
233  } else {
234    clpModel_ = NULL;
235  }
236  abcProgress_ = rhs.abcProgress_;
237  solveType_ = rhs.solveType_;
238}
239// type == 0 nullify, 1 also delete
240void
241AbcSimplex::gutsOfDelete(int type)
242{
243  if (type) {
244    delete [] abcLower_;
245    delete [] abcUpper_;
246    delete [] abcCost_;
247    delete [] abcDj_;
248    delete [] abcSolution_;
249    //delete [] fromExternal_ ;
250    //delete [] toExternal_ ;
251    delete [] scaleFromExternal_ ;
252    //delete [] scaleToExternal_ ;
253    delete [] offset_ ;
254    delete [] internalStatus_ ;
255    delete [] abcPerturbation_ ;
256    delete abcMatrix_ ;
257    delete abcFactorization_;
258#ifdef EARLY_FACTORIZE
259    delete abcEarlyFactorization_;
260#endif
261    delete [] abcPivotVariable_;
262    delete abcDualRowPivot_;
263    delete abcPrimalColumnPivot_;
264    delete abcBaseModel_;
265    delete clpModel_;
266    delete abcNonLinearCost_;
267  }
268  CoinAbcMemset0(reinterpret_cast<char *>(&scaleToExternal_),
269             reinterpret_cast<char *>(&usefulArray_[0])-reinterpret_cast<char *>(&scaleToExternal_));
270}
271template <class T> T *
272newArray(T * /*nullArray*/, int size)
273{
274  if (size) {
275    T * arrayNew = new T[size];
276#ifndef NDEBUG
277    memset(arrayNew,'A',(size)*sizeof(T));
278#endif
279    return arrayNew;
280  } else {
281    return NULL;
282  }
283}
284template <class T> T *
285resizeArray( T * array, int oldSize1, int oldSize2, int newSize1, int newSize2,int extra)
286{
287  if ((array||!oldSize1)&&(newSize1!=oldSize1||newSize2!=oldSize2)) {
288    int newTotal=newSize1+newSize2;
289    int oldTotal=oldSize1+oldSize2;
290    T * arrayNew;
291    if (newSize1>oldSize1||newSize2>oldSize2) {
292      arrayNew = new T[2*newTotal+newSize1+extra];
293#ifndef NDEBUG
294      memset(arrayNew,'A',(2*newTotal+newSize1+extra)*sizeof(T));
295#endif
296      CoinAbcMemcpy(arrayNew,array,oldSize1);
297      CoinAbcMemcpy(arrayNew+newSize1,array+oldSize1,oldSize2);
298      // and second half
299      CoinAbcMemcpy(arrayNew+newTotal,array,oldSize1+oldTotal);
300      CoinAbcMemcpy(arrayNew+newSize1+newTotal,array+oldSize1+oldTotal,oldSize2);
301      delete [] array;
302    } else {
303      arrayNew=array;
304      for (int i=0;i<newSize2;i++)
305        array[newSize1+i]=array[oldSize1+i];
306      for (int i=0;i<newSize1;i++)
307        array[newTotal+i]=array[oldTotal+i];
308      for (int i=0;i<newSize2;i++)
309        array[newSize1+newTotal+i]=array[oldSize1+oldTotal+i];
310    }
311    return arrayNew;
312  } else {
313    return array;
314  }
315}
316// Initializes arrays
317void 
318AbcSimplex::gutsOfInitialize(int numberRows,int numberColumns,bool doMore)
319{
320  abcSimplex_=this;
321  // Zero all
322  CoinAbcMemset0(&sumNonBasicCosts_,(reinterpret_cast<double *>(&usefulArray_[0])-&sumNonBasicCosts_));
323  zeroTolerance_ = 1.0e-13;
324  bestObjectiveValue_ = -COIN_DBL_MAX;
325  primalToleranceToGetOptimal_ = -1.0;
326  primalTolerance_ = 1.0e-6;
327  //dualTolerance_ = 1.0e-6;
328  alphaAccuracy_ = -1.0;
329  upperIn_ = -COIN_DBL_MAX;
330  lowerOut_ = -1;
331  valueOut_ = -1;
332  upperOut_ = -1;
333  dualOut_ = -1;
334  acceptablePivot_ = 1.0e-8;
335  //dualBound_=1.0e9;
336  sequenceIn_ = -1;
337  directionIn_ = -1;
338  sequenceOut_ = -1;
339  directionOut_ = -1;
340  pivotRow_ = -1;
341  lastGoodIteration_ = -100;
342  numberPrimalInfeasibilities_ = 100;
343  numberRefinements_=1;
344  changeMade_ = 1;
345  forceFactorization_ = -1;
346  if (perturbation_<50||(perturbation_>60&&perturbation_<100))
347    perturbation_ = 50;
348  lastBadIteration_ = -999999;
349  lastFlaggedIteration_ = -999999; // doesn't seem to be used
350  firstFree_ = -1;
351  incomingInfeasibility_ = 1.0;
352  allowedInfeasibility_ = 10.0;
353  solveType_ = 1; // say simplex based life form
354  //specialOptions_|=65536;
355  //ClpModel::startPermanentArrays();
356  maximumInternalRows_ =0;
357  maximumInternalColumns_ =0;
358  numberRows_=numberRows;
359  numberColumns_=numberColumns;
360  numberTotal_=numberRows_+numberColumns_;
361  maximumAbcNumberRows_=numberRows;
362  maximumAbcNumberColumns_=numberColumns;
363  maximumNumberTotal_=numberTotal_;
364  int sizeArray=2*maximumNumberTotal_+maximumAbcNumberRows_;
365  if (doMore) {
366    // say Steepest pricing
367    abcDualRowPivot_ = new AbcDualRowSteepest();
368    abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
369    internalStatus_ = newArray(reinterpret_cast<unsigned char *>(NULL),sizeArray+maximumNumberTotal_);
370    abcLower_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
371    abcUpper_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
372    abcCost_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
373    abcDj_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
374    abcSolution_ = newArray(reinterpret_cast<double *>(NULL),sizeArray+maximumNumberTotal_);
375    //fromExternal_ = newArray(reinterpret_cast<int *>(NULL),sizeArray);
376    //toExternal_ = newArray(reinterpret_cast<int *>(NULL),sizeArray);
377    scaleFromExternal_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
378    offset_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
379    abcPerturbation_ = newArray(reinterpret_cast<double *>(NULL),sizeArray);
380    abcPivotVariable_ = newArray(reinterpret_cast<int *>(NULL),maximumAbcNumberRows_);
381    // Fill perturbation array
382    setupPointers(maximumAbcNumberRows_,maximumAbcNumberColumns_);
383    fillPerturbation(0,maximumNumberTotal_);
384  }
385  // get an empty factorization so we can set tolerances etc
386  getEmptyFactorization();
387  for (int i=0;i<ABC_NUMBER_USEFUL;i++) 
388    usefulArray_[i].reserve(CoinMax(CoinMax(numberTotal_,maximumAbcNumberRows_+200),2*numberRows_));
389  //savedStatus_ = internalStatus_+maximumNumberTotal_;
390  //startPermanentArrays();
391}
392// resizes arrays
393void 
394AbcSimplex::gutsOfResize(int numberRows,int numberColumns)
395{
396  if (!abcSolution_) {
397    numberRows_=0;
398    numberColumns_=0;
399    numberTotal_=0;
400    maximumAbcNumberRows_=0;
401    maximumAbcNumberColumns_=0;
402    maximumNumberTotal_=0;
403  }
404  if (numberRows==numberRows_&&numberColumns==numberColumns_)
405    return;
406  // can do on state bit
407  int newSize1=CoinMax(numberRows,maximumAbcNumberRows_);
408  if ((stateOfProblem_&ADD_A_BIT)!=0&&numberRows>maximumAbcNumberRows_)
409    newSize1=CoinMax(numberRows,maximumAbcNumberRows_+CoinMin(100,numberRows_/10));
410  int newSize2=CoinMax(numberColumns,maximumAbcNumberColumns_);
411  numberRows_=numberRows;
412  numberColumns_=numberColumns;
413  numberTotal_=numberRows_+numberColumns_;
414  //fromExternal_ = resizeArray(fromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1);
415  //toExternal_ = resizeArray(toExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
416  scaleFromExternal_ = resizeArray(scaleFromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
417  //scaleToExternal_ = resizeArray(scaleToExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
418  internalStatus_ = resizeArray(internalStatus_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,numberTotal_);
419  abcLower_ = resizeArray(abcLower_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
420  abcUpper_ = resizeArray(abcUpper_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
421  abcCost_ = resizeArray(abcCost_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
422  abcDj_ = resizeArray(abcDj_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
423  abcSolution_ = resizeArray(abcSolution_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,numberTotal_);
424  abcPerturbation_ = resizeArray(abcPerturbation_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
425  offset_ = resizeArray(offset_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,0);
426  setupPointers(newSize1,newSize2);
427  // Fill gaps in perturbation array
428  fillPerturbation(maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_);
429  fillPerturbation(newSize1+maximumAbcNumberColumns_,newSize2-maximumAbcNumberColumns_);
430  // Clean gap
431  //CoinIotaN(fromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
432  //CoinIotaN(toExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
433  CoinFillN(scaleFromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,1.0);
434  CoinFillN(scaleToExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,1.0);
435  CoinFillN(internalStatusSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,static_cast<unsigned char>(1));
436  CoinFillN(lowerSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,-COIN_DBL_MAX);
437  CoinFillN(upperSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,COIN_DBL_MAX);
438  CoinFillN(costSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,0.0);
439  CoinFillN(djSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,0.0);
440  CoinFillN(solutionSaved_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,0.0);
441  CoinFillN(offset_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,0.0);
442  maximumAbcNumberRows_=newSize1;
443  maximumAbcNumberColumns_=newSize2;
444  maximumNumberTotal_=newSize1+newSize2;
445  delete [] abcPivotVariable_;
446  abcPivotVariable_ = new int[maximumAbcNumberRows_];
447  for (int i=0;i<ABC_NUMBER_USEFUL;i++) 
448    usefulArray_[i].reserve(CoinMax(numberTotal_,maximumAbcNumberRows_+200));
449}
450void
451AbcSimplex::translate(int type)
452{
453  // clear bottom bits
454  stateOfProblem_ &= ~(VALUES_PASS-1);
455  if ((type&DO_SCALE_AND_MATRIX)!=0) {
456    //stateOfProblem_ |= DO_SCALE_AND_MATRIX;
457    stateOfProblem_ |= DO_SCALE_AND_MATRIX;
458    delete abcMatrix_;
459    abcMatrix_=new AbcMatrix(*matrix());
460    abcMatrix_->setModel(this);
461    abcMatrix_->scale(scalingFlag_ ? 0 : -1);
462  }
463  if ((type&DO_STATUS)!=0&&(type&DO_BASIS_AND_ORDER)==0) {
464    // from Clp enum to Abc enum (and bound flip)
465    unsigned char lookupToAbcSlack[6]={4,6,0/*1*/,1/*0*/,5,7};
466    unsigned char *  COIN_RESTRICT statusAbc=internalStatus_;
467    const unsigned char *  COIN_RESTRICT statusClp=status_;
468    int i;
469    for (i=numberRows_-1;i>=0;i--) {
470      unsigned char status=statusClp[i]&7;
471      bool basicClp=status==1;
472      bool basicAbc=(statusAbc[i]&7)==4;
473      if (basicClp==basicAbc)
474        statusAbc[i]=lookupToAbcSlack[status];
475      else
476        break;
477    }
478    if (!i) {
479      // from Clp enum to Abc enum
480      unsigned char lookupToAbc[6]={4,6,1,0,5,7};
481      statusAbc+=maximumAbcNumberRows_;
482      statusClp+=maximumAbcNumberRows_;
483      for (i=numberColumns_-1;i>=0;i--) {
484        unsigned char status=statusClp[i]&7;
485        bool basicClp=status==1;
486        bool basicAbc=(statusAbc[i]&7)==4;
487        if (basicClp==basicAbc)
488          statusAbc[i]=lookupToAbc[status];
489        else
490          break;
491      }
492      if (i)
493        type |=DO_BASIS_AND_ORDER;
494    } else {
495      type |=DO_BASIS_AND_ORDER;
496    }
497    stateOfProblem_ |= DO_STATUS;
498  }
499  if ((type&DO_BASIS_AND_ORDER)!=0) {
500    permuteIn();
501    permuteBasis();
502    stateOfProblem_ |= DO_BASIS_AND_ORDER;
503    type &= ~DO_SOLUTION;
504  }
505  if ((type&DO_SOLUTION)!=0) {
506    permuteBasis();
507    stateOfProblem_ |= DO_SOLUTION;
508  }
509  if ((type&DO_JUST_BOUNDS)!=0) {
510    stateOfProblem_ |= DO_JUST_BOUNDS;
511  }
512  if (!type) {
513    // just move stuff down
514    CoinAbcMemcpy(abcLower_,abcLower_+maximumNumberTotal_,numberTotal_);
515    CoinAbcMemcpy(abcUpper_,abcUpper_+maximumNumberTotal_,numberTotal_);
516    CoinAbcMemcpy(abcCost_,abcCost_+maximumNumberTotal_,numberTotal_);
517  }
518}
519/* Sets dual values pass djs using unscaled duals
520   type 1 - values pass
521   type 2 - just use as infeasibility weights
522   type 3 - as 2 but crash
523*/
524void 
525AbcSimplex::setupDualValuesPass(const double * fakeDuals,
526                                const double * fakePrimals,
527                                int type)
528{
529  // allslack basis
530  memset(internalStatus_,6,numberRows_);
531  // temp
532  if (type==1) {
533    bool allEqual=true;
534    for (int i=0;i<numberRows_;i++) {
535      if (rowUpper_[i]>rowLower_[i]) {
536        allEqual=false;
537        break;
538      }
539    }
540    if (allEqual) {
541      // just modify costs
542      transposeTimes(-1.0,fakeDuals,objective());
543      return;
544    }
545  }
546  // compute unscaled djs
547  CoinAbcMemcpy(djSaved_+maximumAbcNumberRows_,objective(),numberColumns_);
548  matrix_->transposeTimes(-1.0,fakeDuals,djSaved_+maximumAbcNumberRows_);
549  // save fake solution
550  assert (solution_);
551  //solution_ = new double[numberTotal_];
552  CoinAbcMemset0(solution_,numberRows_);
553  CoinAbcMemcpy(solution_+maximumAbcNumberRows_,fakePrimals,numberColumns_);
554  // adjust
555  for (int iSequence=maximumAbcNumberRows_;iSequence<numberTotal_;iSequence++)
556    solution_[iSequence]-=offset_[iSequence];
557  matrix_->times(-1.0,solution_+maximumAbcNumberRows_,solution_);
558  double direction = optimizationDirection_;
559  const double *  COIN_RESTRICT rowScale=scaleFromExternal_;
560  const double *  COIN_RESTRICT inverseRowScale=scaleToExternal_;
561  for (int iRow=0;iRow<numberRows_;iRow++) {
562    djSaved_[iRow]=direction*fakeDuals[iRow]*inverseRowScale[iRow];
563    solution_[iRow] *= rowScale[iRow];
564  }
565  const double *  COIN_RESTRICT columnScale=scaleToExternal_;
566  for (int iColumn=maximumAbcNumberRows_;iColumn<numberTotal_;iColumn++)
567    djSaved_[iColumn]*=direction*columnScale[iColumn];
568  // Set solution values
569  const double * lower = abcLower_+maximumAbcNumberRows_;
570  const double * upper = abcUpper_+maximumAbcNumberRows_;
571  double * solution = abcSolution_+maximumAbcNumberRows_;
572  const double * djs = djSaved_+maximumAbcNumberRows_;
573  const double * inverseColumnScale=scaleFromExternal_+maximumAbcNumberRows_;
574  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
575    double thisValue=fakePrimals[iColumn];
576    double thisLower=columnLower_[iColumn];
577    double thisUpper=columnUpper_[iColumn];
578    double thisDj=djs[iColumn];
579    solution_[iColumn+maximumAbcNumberRows_]=solution_[iColumn+maximumAbcNumberRows_]*
580      inverseColumnScale[iColumn];
581    if (thisLower>-1.0e30) {
582      if (thisUpper<1.0e30) {
583        double gapUp=thisUpper-thisValue;
584        double gapDown=thisValue-thisLower;
585        bool goUp;
586        if (gapUp>gapDown&&thisDj>-dualTolerance_) {
587          goUp=false;
588        } else if (gapUp<gapDown&&thisDj<dualTolerance_) {
589          goUp=true;
590        } else {
591          // wild guess
592          double badUp;
593          double badDown;
594          if (gapUp>gapDown) {
595            badUp=gapUp*dualTolerance_;
596            badDown=-gapDown*thisDj;
597          } else {
598            badUp=gapUp*thisDj;
599            badDown=gapDown*dualTolerance_;
600          }
601          goUp=badDown>badUp;
602        }
603        if (goUp) {
604          solution[iColumn]=upper[iColumn];
605          setInternalStatus(iColumn+maximumAbcNumberRows_,atUpperBound);
606          setStatus(iColumn,ClpSimplex::atUpperBound);
607        } else {
608          solution[iColumn]=lower[iColumn];
609          setInternalStatus(iColumn+maximumAbcNumberRows_,atLowerBound);
610          setStatus(iColumn,ClpSimplex::atLowerBound);
611        }
612      } else {
613        solution[iColumn]=lower[iColumn];
614        setInternalStatus(iColumn+maximumAbcNumberRows_,atLowerBound);
615        setStatus(iColumn,ClpSimplex::atLowerBound);
616      }
617    } else if (thisUpper<1.0e30) {
618      solution[iColumn]=upper[iColumn];
619      setInternalStatus(iColumn+maximumAbcNumberRows_,atUpperBound);
620      setStatus(iColumn,ClpSimplex::atUpperBound);
621    } else {
622      // free
623      solution[iColumn]=thisValue*inverseColumnScale[iColumn];
624      setInternalStatus(iColumn+maximumAbcNumberRows_,isFree);
625      setStatus(iColumn,ClpSimplex::isFree);
626    }
627  }
628  switch (type) {
629  case 1:
630    stateOfProblem_ |= VALUES_PASS;
631    break;
632  case 2:
633    stateOfProblem_ |= VALUES_PASS2;
634    delete [] solution_;
635    solution_=NULL;
636    break;
637  case 3:
638    stateOfProblem_ |= VALUES_PASS2;
639    break;
640  }
641}
642//#############################################################################
643int
644AbcSimplex::computePrimals(CoinIndexedVector * arrayVector, CoinIndexedVector * previousVector)
645{
646 
647  arrayVector->clear();
648  previousVector->clear();
649  // accumulate non basic stuff
650 
651  double *  COIN_RESTRICT array = arrayVector->denseVector();
652  CoinAbcScatterZeroTo(abcSolution_,abcPivotVariable_,numberRows_);
653  abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, array);
654#if 0
655  static int xxxxxx=0;
656  if (xxxxxx==0)
657  for (int i=0;i<numberRows_;i++)
658    printf("%d %.18g\n",i,array[i]);
659#endif
660  //arrayVector->scan(0,numberRows_,zeroTolerance_);
661  // Ftran adjusted RHS and iterate to improve accuracy
662  double lastError = COIN_DBL_MAX;
663  CoinIndexedVector * thisVector = arrayVector;
664  CoinIndexedVector * lastVector = previousVector;
665  //if (arrayVector->getNumElements())
666#if 0
667  double largest=0.0;
668  int iLargest=-1;
669  for (int i=0;i<numberRows_;i++) {
670    if (fabs(array[i])>largest) {
671      largest=array[i];
672      iLargest=i;
673    }
674  }
675  printf("largest %g at row %d\n",array[iLargest],iLargest);
676#endif
677  abcFactorization_->updateFullColumn(*thisVector);
678#if 0
679  largest=0.0;
680  iLargest=-1;
681  for (int i=0;i<numberRows_;i++) {
682    if (fabs(array[i])>largest) {
683      largest=array[i];
684      iLargest=i;
685    }
686  }
687  printf("after largest %g at row %d\n",array[iLargest],iLargest);
688#endif
689#if 0
690  if (xxxxxx==0)
691  for (int i=0;i<numberRows_;i++)
692    printf("xx %d %.19g\n",i,array[i]);
693  if (numberIterations_>300)
694  exit(0);
695#endif
696  int numberRefinements=0;
697  for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
698    int numberIn = thisVector->getNumElements();
699    const int *  COIN_RESTRICT indexIn = thisVector->getIndices();
700    const double *  COIN_RESTRICT arrayIn = thisVector->denseVector();
701    CoinAbcScatterToList(arrayIn,abcSolution_,indexIn,abcPivotVariable_,numberIn);
702    // check Ax == b  (for all)
703    abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, solutionBasic_);
704#if 0
705    if (xxxxxx==0)
706      for (int i=0;i<numberTotal_;i++)
707        printf("sol %d %.19g\n",i,abcSolution_[i]);
708    if (xxxxxx==0)
709      for (int i=0;i<numberRows_;i++)
710        printf("basic %d %.19g\n",i,solutionBasic_[i]);
711#endif
712    double multiplier = 131072.0;
713    largestPrimalError_ = CoinAbcMaximumAbsElementAndScale(solutionBasic_,multiplier,numberRows_);
714    double maxValue=0.0;
715    for (int i=0;i<numberRows_;i++) {
716      double value=fabs(solutionBasic_[i]);
717      if (value>maxValue) {
718#if 0
719        if (xxxxxx==0)
720          printf("larger %.19gg at row %d\n",value,i);
721        maxValue=value;
722#endif
723      }
724    }
725    if (largestPrimalError_ >= lastError) {
726      // restore
727      CoinIndexedVector * temp = thisVector;
728      thisVector = lastVector;
729      lastVector = temp;
730      //goodSolution = false;
731      break;
732    }
733    if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
734      // try and make better
735      numberRefinements++;
736      // save this
737      CoinIndexedVector * temp = thisVector;
738      thisVector = lastVector;
739      lastVector = temp;
740      int *  COIN_RESTRICT indexOut = thisVector->getIndices();
741      int number = 0;
742      array = thisVector->denseVector();
743      thisVector->clear();
744      for (int iRow = 0; iRow < numberRows_; iRow++) {
745        double value = solutionBasic_[iRow];
746        if (value) {
747          array[iRow] = value;
748          indexOut[number++] = iRow;
749          solutionBasic_[iRow] = 0.0;
750        }
751      }
752      thisVector->setNumElements(number);
753      lastError = largestPrimalError_;
754      abcFactorization_->updateFullColumn(*thisVector);
755      double * previous = lastVector->denseVector();
756      number = 0;
757      multiplier=1.0/multiplier;
758      for (int iRow = 0; iRow < numberRows_; iRow++) {
759        double value = previous[iRow] + multiplier * array[iRow];
760        if (value) {
761          array[iRow] = value;
762          indexOut[number++] = iRow;
763        } else {
764          array[iRow] = 0.0;
765        }
766      }
767      thisVector->setNumElements(number);
768    } else {
769      break;
770    }
771  }
772 
773  // solution as accurate as we are going to get
774  //if (!goodSolution) {
775  CoinAbcMemcpy(solutionBasic_,thisVector->denseVector(),numberRows_);
776  CoinAbcScatterTo(solutionBasic_,abcSolution_,abcPivotVariable_,numberRows_);
777  arrayVector->clear();
778  previousVector->clear();
779  return numberRefinements;
780}
781// Moves basic stuff to basic area
782void 
783AbcSimplex::moveToBasic(int which)
784{
785  if ((which&1)!=0)
786    CoinAbcGatherFrom(abcSolution_,solutionBasic_,abcPivotVariable_,numberRows_);
787  if ((which&2)!=0)
788    CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);
789  if ((which&4)!=0)
790    CoinAbcGatherFrom(abcLower_,lowerBasic_,abcPivotVariable_,numberRows_);
791  if ((which&8)!=0)
792    CoinAbcGatherFrom(abcUpper_,upperBasic_,abcPivotVariable_,numberRows_);
793}
794// now dual side
795int
796AbcSimplex::computeDuals(double * givenDjs, CoinIndexedVector * arrayVector, CoinIndexedVector * previousVector)
797{
798  double *  COIN_RESTRICT array = arrayVector->denseVector();
799  int *  COIN_RESTRICT index = arrayVector->getIndices();
800  int number = 0;
801  if (!givenDjs) {
802    for (int iRow = 0; iRow < numberRows_; iRow++) {
803      double value = costBasic_[iRow];
804      if (value) {
805        array[iRow] = value;
806        index[number++] = iRow;
807      }
808    }
809  } else {
810    // dual values pass - djs may not be zero
811    for (int iRow = 0; iRow < numberRows_; iRow++) {
812      int iPivot = abcPivotVariable_[iRow];
813      // make sure zero if done
814      if (!pivoted(iPivot))
815        givenDjs[iPivot] = 0.0;
816      double value = abcCost_[iPivot] - givenDjs[iPivot];
817      if (value) {
818        array[iRow] = value;
819        index[number++] = iRow;
820      }
821    }
822  }
823  arrayVector->setNumElements(number);
824  // Btran basic costs and get as accurate as possible
825  double lastError = COIN_DBL_MAX;
826  CoinIndexedVector * thisVector = arrayVector;
827  CoinIndexedVector * lastVector = previousVector;
828  abcFactorization_->updateFullColumnTranspose(*thisVector);
829 
830  int numberRefinements=0;
831  for (int iRefine = 0; iRefine < numberRefinements_+1; iRefine++) {
832    // check basic reduced costs zero
833    // put reduced cost of basic into djBasic_
834    CoinAbcMemcpy(djBasic_,costBasic_,numberRows_);
835    abcMatrix_->transposeTimesBasic(-1.0,thisVector->denseVector(),djBasic_);
836    largestDualError_ = CoinAbcMaximumAbsElement(djBasic_,numberRows_);
837    if (largestDualError_ >= lastError) {
838      // restore
839      CoinIndexedVector * temp = thisVector;
840      thisVector = lastVector;
841      lastVector = temp;
842      break;
843    }
844    if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
845        && !givenDjs) {
846      numberRefinements++;
847      // try and make better
848      // save this
849      CoinIndexedVector * temp = thisVector;
850      thisVector = lastVector;
851      lastVector = temp;
852      array = thisVector->denseVector();
853      double multiplier = 131072.0;
854      //array=djBasic_*multiplier
855      CoinAbcMultiplyAdd(djBasic_,numberRows_,multiplier,array,0.0);
856      lastError = largestDualError_;
857      abcFactorization_->updateFullColumnTranspose( *thisVector);
858#if ABC_DEBUG
859      thisVector->checkClean();
860#endif
861      multiplier = 1.0 / multiplier;
862      double *  COIN_RESTRICT previous = lastVector->denseVector();
863      // array = array*multiplier+previous
864      CoinAbcMultiplyAdd(previous,numberRows_,1.0,array,multiplier);
865    } else {
866      break;
867    }
868  }
869  // now look at dual solution
870  CoinAbcMemcpy(dual_,thisVector->denseVector(),numberRows_);
871  CoinAbcMemset0(thisVector->denseVector(),numberRows_);
872  thisVector->setNumElements(0);
873  if (numberRefinements) {
874    CoinAbcMemset0(lastVector->denseVector(),numberRows_);
875    lastVector->setNumElements(0);
876  }
877  CoinAbcMemcpy(abcDj_,abcCost_,numberTotal_);
878  abcMatrix_->transposeTimesNonBasic(-1.0, dual_,abcDj_);
879  // If necessary - override results
880  if (givenDjs) {
881    // restore accurate duals
882    CoinMemcpyN(abcDj_, numberTotal_, givenDjs);
883  }
884  //arrayVector->clear();
885  //previousVector->clear();
886#if ABC_DEBUG
887  arrayVector->checkClean();
888  previousVector->checkClean();
889#endif
890  return numberRefinements;
891}
892
893/* Factorizes using current basis.
894   solveType - 1 iterating, 0 initial, -1 external
895   - 2 then iterating but can throw out of basis
896   If 10 added then in primal values pass
897   Return codes are as from AbcSimplexFactorization unless initial factorization
898   when total number of singularities is returned.
899   Special case is numberRows_+1 -> all slack basis.
900*/
901int AbcSimplex::internalFactorize ( int solveType)
902{
903  assert (status_);
904 
905  bool valuesPass = false;
906  if (solveType >= 10) {
907    valuesPass = true;
908    solveType -= 10;
909  }
910#if 0
911  // Make sure everything is clean
912  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
913    if(getInternalStatus(iSequence) == isFixed) {
914      // double check fixed
915      assert (abcUpper_[iSequence] == abcLower_[iSequence]);
916      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<100.0*primalTolerance_);
917    } else if (getInternalStatus(iSequence) == isFree) {
918      assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
919    } else if (getInternalStatus(iSequence) == atLowerBound) {
920      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
921    } else if (getInternalStatus(iSequence) == atUpperBound) {
922      assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
923    } else if (getInternalStatus(iSequence) == superBasic) {
924      assert (!valuesPass);
925    }
926  }
927#endif
928#if 0 //ndef NDEBUG
929  // Make sure everything is clean
930  double sumOutside=0.0;
931  int numberOutside=0;
932  //double sumOutsideLarge=0.0;
933  int numberOutsideLarge=0;
934  double sumInside=0.0;
935  int numberInside=0;
936  //double sumInsideLarge=0.0;
937  int numberInsideLarge=0;
938  char rowcol[] = {'R', 'C'};
939  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
940    if(getInternalStatus(iSequence) == isFixed) {
941      // double check fixed
942      assert (abcUpper_[iSequence] == abcLower_[iSequence]);
943      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<primalTolerance_);
944    } else if (getInternalStatus(iSequence) == isFree) {
945      assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
946    } else if (getInternalStatus(iSequence) == atLowerBound) {
947      assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
948      if (abcSolution_[iSequence]<abcLower_[iSequence]) {
949        numberOutside++;
950#if ABC_NORMAL_DEBUG>1
951        if (handler_->logLevel()==63)
952          printf("%c%d below by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
953                 abcLower_[iSequence]-abcSolution_[iSequence]);
954#endif
955        sumOutside-=abcSolution_[iSequence]-abcLower_[iSequence];
956        if (abcSolution_[iSequence]<abcLower_[iSequence]-primalTolerance_)
957          numberOutsideLarge++;
958      } else if (abcSolution_[iSequence]>abcLower_[iSequence]) {
959        numberInside++;
960        sumInside+=abcSolution_[iSequence]-abcLower_[iSequence];
961        if (abcSolution_[iSequence]>abcLower_[iSequence]+primalTolerance_)
962          numberInsideLarge++;
963      }
964    } else if (getInternalStatus(iSequence) == atUpperBound) {
965      assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
966      if (abcSolution_[iSequence]>abcUpper_[iSequence]) {
967        numberOutside++;
968#if ABC_NORMAL_DEBUG>0
969        if (handler_->logLevel()==63)
970          printf("%c%d above by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
971                 -(abcUpper_[iSequence]-abcSolution_[iSequence]));
972#endif
973        sumOutside+=abcSolution_[iSequence]-abcUpper_[iSequence];
974        if (abcSolution_[iSequence]>abcUpper_[iSequence]+primalTolerance_)
975          numberOutsideLarge++;
976      } else if (abcSolution_[iSequence]<abcUpper_[iSequence]) {
977        numberInside++;
978        sumInside-=abcSolution_[iSequence]-abcUpper_[iSequence];
979        if (abcSolution_[iSequence]<abcUpper_[iSequence]-primalTolerance_)
980          numberInsideLarge++;
981      }
982    } else if (getInternalStatus(iSequence) == superBasic) {
983      //assert (!valuesPass);
984    }
985  }
986#if ABC_NORMAL_DEBUG>0
987  if (numberInside+numberOutside)
988    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
989           numberOutside,sumOutside,numberOutsideLarge,
990           numberInside,sumInside,numberInsideLarge);
991#endif
992#endif
993  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
994    AbcSimplex::Status status=getInternalStatus(iSequence);
995    if (status!= basic&&status!=isFixed&&abcUpper_[iSequence] == abcLower_[iSequence])
996      setInternalStatus(iSequence,isFixed);
997  }
998  if (numberIterations_==baseIteration_&&!valuesPass) {
999    double largeValue = this->largeValue();
1000    double * COIN_RESTRICT solution = abcSolution_;
1001    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1002      AbcSimplex::Status status=getInternalStatus(iSequence);
1003      if (status== superBasic) {
1004        double lower = abcLower_[iSequence];
1005        double upper = abcUpper_[iSequence];
1006        double value = solution[iSequence];
1007        AbcSimplex::Status thisStatus=isFree;
1008        if (lower > -largeValue || upper < largeValue) {
1009          if (lower!=upper) {
1010            if (fabs(value - lower) < fabs(value - upper)) {
1011              thisStatus=AbcSimplex::atLowerBound;
1012              solution[iSequence] = lower;
1013            } else {
1014              thisStatus= AbcSimplex::atUpperBound;
1015              solution[iSequence] = upper;
1016            }
1017          } else {
1018            thisStatus= AbcSimplex::isFixed;
1019            solution[iSequence] = upper;
1020          }
1021          setInternalStatus(iSequence,thisStatus);
1022        }
1023      }
1024    }
1025  }
1026  int status = abcFactorization_->factorize(this, solveType, valuesPass);
1027  if (status) {
1028    handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1029      << status
1030      << CoinMessageEol;
1031    return -1;
1032  } else if (!solveType) {
1033    // Initial basis - return number of singularities
1034    int numberSlacks = 0;
1035    for (int iRow = 0; iRow < numberRows_; iRow++) {
1036      if (getInternalStatus(iRow) == basic)
1037        numberSlacks++;
1038    }
1039    status = CoinMax(numberSlacks - numberRows_, 0);
1040    // special case if all slack
1041    if (numberSlacks == numberRows_) {
1042      status = numberRows_ + 1;
1043    }
1044  }
1045  return status;
1046}
1047// Sets objectiveValue_ from rawObjectiveValue_
1048void 
1049AbcSimplex::setClpSimplexObjectiveValue()
1050{
1051  objectiveValue_ = rawObjectiveValue_/(objectiveScale_ * rhsScale_);
1052  objectiveValue_ += objectiveOffset_;
1053}
1054/*
1055  This does basis housekeeping and does values for in/out variables.
1056  Can also decide to re-factorize
1057*/
1058int
1059AbcSimplex::housekeeping()
1060{
1061#define DELAYED_UPDATE
1062#ifdef DELAYED_UPDATE
1063  if (algorithm_<0) {
1064    // modify dualout
1065    dualOut_ /= alpha_;
1066    dualOut_ *= -directionOut_;
1067    //double oldValue = valueIn_;
1068    if (directionIn_ == -1) {
1069      // as if from upper bound
1070      valueIn_ = upperIn_ + dualOut_;
1071#if 0 //def ABC_DEBUG
1072      printf("In from upper dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1073             dualOut_,movement_,movementOld,valueIn_,upperIn_+movement_);
1074#endif
1075    } else {
1076      // as if from lower bound
1077      valueIn_ = lowerIn_ + dualOut_;
1078#if 0 //def ABC_DEBUG
1079      printf("In from lower dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1080             dualOut_,movement_,movementOld,valueIn_,lowerIn_+movement_);
1081#endif
1082    }
1083    // outgoing
1084    if (directionOut_ > 0) {
1085      valueOut_ = lowerOut_;
1086    } else {
1087      valueOut_ = upperOut_;
1088    }
1089#if ABC_NORMAL_DEBUG>3
1090    if (rawObjectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1091      printf("obj backwards %g %g\n", rawObjectiveValue_, oldobj);
1092#endif
1093  }
1094#endif
1095#if 0 //ndef NDEBUG
1096  {
1097    int numberFlagged=0;
1098    for (int iRow = 0; iRow < numberRows_; iRow++) {
1099      int iPivot = abcPivotVariable_[iRow];
1100      if (flagged(iPivot))
1101        numberFlagged++;
1102    }
1103    assert (numberFlagged==numberFlagged_);
1104  }
1105#endif
1106  // save value of incoming and outgoing
1107  numberIterations_++;
1108  changeMade_++; // something has happened
1109  // incoming variable
1110  if (handler_->logLevel() > 7) {
1111    //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
1112    handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
1113      << directionOut_
1114      << directionIn_ << theta_
1115      << dualOut_ << dualIn_ << alpha_
1116      << CoinMessageEol;
1117    if (getInternalStatus(sequenceIn_) == isFree) {
1118      handler_->message(CLP_SIMPLEX_FREEIN, messages_)
1119        << sequenceIn_
1120        << CoinMessageEol;
1121    }
1122  }
1123  // change of incoming
1124  char rowcol[] = {'R', 'C'};
1125  if (abcUpper_[sequenceIn_] > 1.0e20 && abcLower_[sequenceIn_] < -1.0e20)
1126    progressFlag_ |= 2; // making real progress
1127#ifndef DELAYED_UPDATE
1128  abcSolution_[sequenceIn_] = valueIn_;
1129#endif
1130  if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] < 1.0e-12)
1131    progressFlag_ |= 1; // making real progress
1132  if (sequenceIn_ != sequenceOut_) {
1133    if (alphaAccuracy_ > 0.0) {
1134      double value = fabs(alpha_);
1135      if (value > 1.0)
1136        alphaAccuracy_ *= value;
1137      else
1138        alphaAccuracy_ /= value;
1139    }
1140    setInternalStatus(sequenceIn_, basic);
1141    if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] > 0) {
1142      // As Nonlinear costs may have moved bounds (to more feasible)
1143      // Redo using value
1144      if (fabs(valueOut_ - abcLower_[sequenceOut_]) < fabs(valueOut_ - abcUpper_[sequenceOut_])) {
1145        // going to lower
1146        setInternalStatus(sequenceOut_, atLowerBound);
1147      } else {
1148        // going to upper
1149        setInternalStatus(sequenceOut_, atUpperBound);
1150      }
1151    } else {
1152      // fixed
1153      setInternalStatus(sequenceOut_, isFixed);
1154    }
1155#ifndef DELAYED_UPDATE
1156    abcSolution_[sequenceOut_] = valueOut_;
1157#endif
1158#if PARTITION_ROW_COPY==1
1159    // move in row copy
1160    abcMatrix_->inOutUseful(sequenceIn_,sequenceOut_);
1161#endif
1162  } else {
1163    //if (objective_->type()<2)
1164    //assert (fabs(theta_)>1.0e-13);
1165    // flip from bound to bound
1166    // As Nonlinear costs may have moved bounds (to more feasible)
1167    // Redo using value
1168    if (fabs(valueIn_ - abcLower_[sequenceIn_]) < fabs(valueIn_ - abcUpper_[sequenceIn_])) {
1169      // as if from upper bound
1170      setInternalStatus(sequenceIn_, atLowerBound);
1171    } else {
1172      // as if from lower bound
1173      setInternalStatus(sequenceIn_, atUpperBound);
1174    }
1175  }
1176  setClpSimplexObjectiveValue();
1177  if (handler_->logLevel() > 7) {
1178    //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
1179    handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
1180      << numberIterations_ << objectiveValue()
1181      << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
1182      << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
1183    handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
1184    handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
1185    handler_->message() << CoinMessageEol;
1186  }
1187#if 0
1188  if (numberIterations_ > 10000)
1189    printf(" it %d %g %c%d %c%d\n"
1190           , numberIterations_, objectiveValue()
1191           , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
1192           , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
1193#endif
1194  if (hitMaximumIterations())
1195    return 2;
1196  // check for small cycles
1197  int in = sequenceIn_;
1198  int out = sequenceOut_;
1199  int cycle = abcProgress_.cycle(in, out,
1200                              directionIn_, directionOut_);
1201  if (cycle > 0 ) {
1202#if ABC_NORMAL_DEBUG>0
1203    if (handler_->logLevel() >= 63)
1204      printf("Cycle of %d\n", cycle);
1205#endif
1206    // reset
1207    abcProgress_.startCheck();
1208    double random = randomNumberGenerator_.randomDouble();
1209    int extra = static_cast<int> (9.999 * random);
1210    int off[] = {1, 1, 1, 1, 2, 2, 2, 3, 3, 4};
1211    if (abcFactorization_->pivots() > cycle) {
1212      forceFactorization_ = CoinMax(1, cycle - off[extra]);
1213    } else {
1214      // need to reject something
1215      int iSequence;
1216      if (algorithm_<0) {
1217        iSequence = sequenceIn_;
1218      } else {
1219        /* should be better if don't reject incoming
1220         as it is in basis */
1221        iSequence = sequenceOut_;
1222        // so can't happen immediately again
1223        sequenceOut_=-1;
1224      }
1225      char x = isColumn(iSequence) ? 'C' : 'R';
1226      if (handler_->logLevel() >= 63)
1227        handler_->message(CLP_SIMPLEX_FLAG, messages_)
1228          << x << sequenceWithin(iSequence)
1229          << CoinMessageEol;
1230      setFlagged(iSequence);
1231      //printf("flagging %d\n",iSequence);
1232    }
1233    return 1;
1234  }
1235  int invertNow=0;
1236  // only time to re-factorize if one before real time
1237  // this is so user won't be surprised that maximumPivots has exact meaning
1238  int numberPivots = abcFactorization_->pivots();
1239  if (algorithm_<0)
1240    numberPivots++; // allow for update not done
1241  int maximumPivots = abcFactorization_->maximumPivots();
1242  int numberDense = 0; //abcFactorization_->numberDense();
1243  bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 >
1244                     2 * maximumIterations());
1245  if (numberPivots == maximumPivots ||
1246      maximumPivots < 2) {
1247    // If dense then increase
1248    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots) {
1249      abcFactorization_->maximumPivots(numberDense);
1250    }
1251    //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots);
1252    return 1;
1253  } else if ((abcFactorization_->timeToRefactorize() && !dontInvert)
1254             ||invertNow) {
1255    //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots);
1256    return 1;
1257  } else if (forceFactorization_ > 0 &&
1258#ifndef DELAYED_UPDATE
1259             abcFactorization_->pivots() 
1260#else
1261             abcFactorization_->pivots()+1
1262#endif
1263             >= forceFactorization_) {
1264    // relax
1265    forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
1266    if (forceFactorization_ > abcFactorization_->maximumPivots())
1267      forceFactorization_ = -1; //off
1268    //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots);
1269    return 1;
1270  } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
1271    // A bit worried problem may be cycling - lets factorize at random intervals for a short period
1272    int numberTooManyIterations = numberIterations_ - 10 * (numberRows_ + (numberColumns_ >> 2));
1273    double random = randomNumberGenerator_.randomDouble();
1274    int window = numberTooManyIterations%5000;
1275    if (window<2*maximumPivots)
1276      random = 0.2*random+0.8; // randomly re-factorize but not too soon
1277    else
1278      random=1.0; // switch off if not in window of opportunity
1279    int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
1280    if (abcFactorization_->pivots() >= random * maxNumber) {
1281      //printf("ZZ cycling a %d\n",numberPivots);
1282      return 1;
1283    } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) &&
1284               numberIterations_ < 1000010 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
1285      //printf("ZZ cycling b %d\n",numberPivots);
1286      return 1;
1287    } else {
1288      // carry on iterating
1289      return 0;
1290    }
1291  } else {
1292    // carry on iterating
1293    return 0;
1294  }
1295}
1296// Swaps primal stuff
1297void 
1298AbcSimplex::swapPrimalStuff()
1299{
1300  if (sequenceOut_<0)
1301    return;
1302  assert (sequenceIn_>=0);
1303  abcSolution_[sequenceOut_] = valueOut_;
1304  abcSolution_[sequenceIn_] = valueIn_;
1305  solutionBasic_[pivotRow_]=valueIn_;
1306  if (algorithm_<0) {
1307    // and set bounds correctly
1308    static_cast<AbcSimplexDual *> (this)->originalBound(sequenceIn_);
1309    static_cast<AbcSimplexDual *> (this)->changeBound(sequenceOut_);
1310  } else {
1311#if ABC_NORMAL_DEBUG>2
1312    if (handler_->logLevel()==63)
1313      printf("Code swapPrimalStuff for primal\n");
1314#endif
1315  }
1316  if (pivotRow_>=0) { // may be flip in primal
1317    lowerBasic_[pivotRow_]=abcLower_[sequenceIn_];
1318    upperBasic_[pivotRow_]=abcUpper_[sequenceIn_];
1319    costBasic_[pivotRow_]=abcCost_[sequenceIn_];
1320    abcPivotVariable_[pivotRow_] = sequenceIn_;
1321  }
1322}
1323// Swaps dual stuff
1324void 
1325AbcSimplex::swapDualStuff(int lastSequenceOut,int lastDirectionOut)
1326{
1327  // outgoing
1328  // set dj to zero unless values pass
1329  if (lastSequenceOut>=0) {
1330    if ((stateOfProblem_&VALUES_PASS)==0) {
1331      if (lastDirectionOut > 0) {
1332        abcDj_[lastSequenceOut] = theta_;
1333      } else {
1334        abcDj_[lastSequenceOut] = -theta_;
1335      }
1336    } else {
1337      if (lastDirectionOut > 0) {
1338        abcDj_[lastSequenceOut] -= theta_;
1339      } else {
1340        abcDj_[lastSequenceOut] += theta_;
1341      }
1342    }
1343  }
1344  if (sequenceIn_>=0) {
1345    abcDj_[sequenceIn_]=0.0;
1346    //costBasic_[pivotRow_]=abcCost_[sequenceIn_];
1347  }
1348}
1349static void solveMany(int number,ClpSimplex ** simplex)
1350{
1351  for (int i=0;i<number-1;i++)
1352    cilk_spawn simplex[i]->dual(0);
1353  simplex[number-1]->dual(0);
1354  cilk_sync;
1355}
1356void 
1357AbcSimplex::crash (int type)
1358{
1359  int i;
1360  for (i=0;i<numberRows_;i++) {
1361    if (getInternalStatus(i)!=basic)
1362      break;
1363  }
1364  if (i<numberRows_)
1365    return;
1366  // all slack
1367  if (type==3) {
1368    // decomposition crash
1369    if (!abcMatrix_->gotRowCopy())
1370      abcMatrix_->createRowCopy();
1371    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
1372    const CoinBigIndex * columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
1373    const int * row = abcMatrix_->getPackedMatrix()->getIndices();
1374    //const double * elementByRow = abcMatrix_->rowElements();
1375    const CoinBigIndex * rowStart = abcMatrix_->rowStart();
1376    const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
1377    const int * column = abcMatrix_->rowColumns();
1378    int * blockStart = usefulArray_[0].getIndices();
1379    int * columnBlock = blockStart+numberRows_;
1380    int * nextColumn = usefulArray_[1].getIndices();
1381    int * blockCount = nextColumn+numberColumns_;
1382    int direction[2]={-1,1};
1383    int bestBreak=-1;
1384    double bestValue=0.0;
1385    int iPass=0;
1386    int halfway=(numberRows_+1)/2;
1387    int firstMaster=-1;
1388    int lastMaster=-2;
1389    int numberBlocks=0;
1390    int largestRows=0;
1391    int largestColumns=0;
1392    int numberEmpty=0;
1393    int numberMaster=0;
1394    int numberEmptyColumns=0;
1395    int numberMasterColumns=0;
1396    while (iPass<2) {
1397      int increment=direction[iPass];
1398      int start= increment>0 ? 0 : numberRows_-1;
1399      int stop=increment>0 ? numberRows_ : -1;
1400      numberBlocks=0;
1401      int thisBestBreak=-1;
1402      double thisBestValue=COIN_DBL_MAX;
1403      int numberRowsDone=0;
1404      int numberMarkedColumns=0;
1405      int maximumBlockSize=0;
1406      for (int i=0;i<numberRows_;i++) { 
1407        blockStart[i]=-1;
1408        blockCount[i]=0;
1409      }
1410      for (int i=0;i<numberColumns_;i++) {
1411        columnBlock[i]=-1;
1412        nextColumn[i]=-1;
1413      }
1414      for (int iRow=start;iRow!=stop;iRow+=increment) {
1415        int iBlock = -1;
1416        for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
1417          int iColumn=column[j]-maximumAbcNumberRows_;
1418          int whichColumnBlock=columnBlock[iColumn];
1419          if (whichColumnBlock>=0) {
1420            // column marked
1421            if (iBlock<0) {
1422              // put row in that block
1423              iBlock=whichColumnBlock;
1424            } else if (iBlock!=whichColumnBlock) {
1425              // merge
1426              blockCount[iBlock]+=blockCount[whichColumnBlock];
1427              blockCount[whichColumnBlock]=0;
1428              int jColumn=blockStart[whichColumnBlock];
1429#ifndef NDEBUG
1430              int iiColumn=iColumn;
1431#endif
1432              while (jColumn>=0) {
1433                assert (columnBlock[jColumn]==whichColumnBlock);
1434                columnBlock[jColumn]=iBlock;
1435#ifndef NDEBUG
1436                if (jColumn==iiColumn)
1437                  iiColumn=-1;
1438#endif
1439                iColumn=jColumn;
1440                jColumn=nextColumn[jColumn];
1441              }
1442#ifndef NDEBUG
1443              assert (iiColumn==-1);
1444#endif
1445              nextColumn[iColumn]=blockStart[iBlock];
1446              blockStart[iBlock]=blockStart[whichColumnBlock];
1447              blockStart[whichColumnBlock]=-1;
1448            }
1449          }
1450        }
1451        int n=numberMarkedColumns;
1452        if (iBlock<0) {
1453          //new block
1454          if (rowEnd[iRow]>rowStart[iRow]) {
1455            numberBlocks++;
1456            iBlock=numberBlocks;
1457            int jColumn=column[rowStart[iRow]]-maximumAbcNumberRows_;
1458            columnBlock[jColumn]=iBlock;
1459            blockStart[iBlock]=jColumn;
1460            numberMarkedColumns++;
1461            for (CoinBigIndex j=rowStart[iRow]+1;j<rowEnd[iRow];j++) {
1462              int iColumn=column[j]-maximumAbcNumberRows_;
1463              columnBlock[iColumn]=iBlock;
1464              numberMarkedColumns++;
1465              nextColumn[jColumn]=iColumn;
1466              jColumn=iColumn;
1467            }
1468            blockCount[iBlock]=numberMarkedColumns-n;
1469          } else {
1470            // empty
1471          }
1472        } else {
1473          // put in existing block
1474          int jColumn=blockStart[iBlock];
1475          for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
1476            int iColumn=column[j]-maximumAbcNumberRows_;
1477            assert (columnBlock[iColumn]<0||columnBlock[iColumn]==iBlock);
1478            if (columnBlock[iColumn]<0) {
1479              columnBlock[iColumn]=iBlock;
1480              numberMarkedColumns++;
1481              nextColumn[iColumn]=jColumn;
1482              jColumn=iColumn;
1483            }
1484          }
1485          blockStart[iBlock]=jColumn;
1486          blockCount[iBlock]+=numberMarkedColumns-n;
1487        }
1488        if (iBlock>=0)
1489          maximumBlockSize=CoinMax(maximumBlockSize,blockCount[iBlock]);
1490        numberRowsDone++;
1491        if (thisBestValue*numberRowsDone > maximumBlockSize&&numberRowsDone>halfway) { 
1492          thisBestBreak=iRow;
1493          thisBestValue=static_cast<double>(maximumBlockSize)/
1494            static_cast<double>(numberRowsDone);
1495        }
1496      }
1497      if (thisBestBreak==stop)
1498        thisBestValue=COIN_DBL_MAX;
1499      iPass++;
1500      if (iPass==1) {
1501        bestBreak=thisBestBreak;
1502        bestValue=thisBestValue;
1503      } else {
1504        if (bestValue<thisBestValue) {
1505          firstMaster=0;
1506          lastMaster=bestBreak;
1507        } else {
1508          firstMaster=thisBestBreak; // ? +1
1509          lastMaster=numberRows_;
1510        }
1511      }
1512    }
1513    bool useful=false;
1514    if (firstMaster<lastMaster) {
1515      for (int i=0;i<numberRows_;i++) 
1516        blockStart[i]=-1;
1517      for (int i=firstMaster;i<lastMaster;i++)
1518        blockStart[i]=-2;
1519      for (int i=0;i<numberColumns_;i++)
1520        columnBlock[i]=-1;
1521      int firstRow=0;
1522      numberBlocks=-1;
1523      while (true) {
1524        for (;firstRow<numberRows_;firstRow++) {
1525          if (blockStart[firstRow]==-1)
1526            break;
1527        }
1528        if (firstRow==numberRows_)
1529          break;
1530        int nRows=0;
1531        numberBlocks++;
1532        int numberStack=1;
1533        blockCount[0] = firstRow;
1534        while (numberStack) {
1535          int iRow=blockCount[--numberStack];
1536          for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
1537            int iColumn=column[j]-maximumAbcNumberRows_;
1538            int iBlock=columnBlock[iColumn];
1539            if (iBlock<0) {
1540              columnBlock[iColumn]=numberBlocks;
1541              for (CoinBigIndex k=columnStart[iColumn];
1542                   k<columnStart[iColumn+1];k++) {
1543                int jRow=row[k];
1544                int rowBlock=blockStart[jRow];
1545                if (rowBlock==-1) {
1546                  nRows++;
1547                  blockStart[jRow]=numberBlocks;
1548                  blockCount[numberStack++]=jRow;
1549                }
1550              }
1551            }
1552          }
1553        }
1554        if (!nRows) {
1555          // empty!!
1556          numberBlocks--;
1557        }
1558        firstRow++;
1559      }
1560      // adjust
1561      numberBlocks++;
1562      for (int i=0;i<numberBlocks;i++) {
1563        blockCount[i]=0;
1564        nextColumn[i]=0;
1565      }
1566      for (int iRow = 0; iRow < numberRows_; iRow++) {
1567        int iBlock=blockStart[iRow];
1568        if (iBlock>=0) {
1569          blockCount[iBlock]++;
1570        } else {
1571          if (iBlock==-2)
1572            numberMaster++;
1573          else
1574            numberEmpty++;
1575        }
1576      }
1577      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1578        int iBlock=columnBlock[iColumn];
1579        if (iBlock>=0) {
1580          nextColumn[iBlock]++;
1581        } else {
1582          if (columnStart[iColumn+1]>columnStart[iColumn])
1583            numberMasterColumns++;
1584          else
1585            numberEmptyColumns++;
1586        }
1587      }
1588      for (int i=0;i<numberBlocks;i++) {
1589        if (blockCount[i]+nextColumn[i]>largestRows+largestColumns) {
1590          largestRows=blockCount[i];
1591          largestColumns=nextColumn[i];
1592        }
1593      }
1594      useful=true;
1595      if (numberMaster>halfway||largestRows*3>numberRows_)
1596        useful=false;
1597    }
1598    if (useful) {
1599#if ABC_NORMAL_DEBUG>0
1600      if (logLevel()>=2)
1601        printf("%d master rows %d <= < %d\n",lastMaster-firstMaster,
1602               firstMaster,lastMaster);
1603      printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n",
1604             useful ? "**Useful" : "NoGood",
1605             numberBlocks,largestRows,largestColumns,numberMaster,numberEmpty,numberRows_,
1606             numberMasterColumns,numberEmptyColumns,numberColumns_);
1607      if (logLevel()>=2) {
1608        for (int i=0;i<numberBlocks;i++) 
1609          printf("Block %d has %d rows and %d columns\n",
1610               i,blockCount[i],nextColumn[i]);
1611      }
1612#endif
1613#define NUMBER_DW_BLOCKS 20
1614      int minSize1=(numberRows_-numberMaster+NUMBER_DW_BLOCKS-1)/NUMBER_DW_BLOCKS;
1615      int minSize2=(numberRows_-numberMaster+2*NUMBER_DW_BLOCKS-1)/(2*NUMBER_DW_BLOCKS);
1616      int * backRow = usefulArray_[1].getIndices();
1617      // first sort
1618      for (int i=0;i<numberBlocks;i++) {
1619        backRow[i]=-(3*blockCount[i]+0*nextColumn[i]);
1620        nextColumn[i]=i;
1621      }
1622      CoinSort_2(backRow,backRow+numberBlocks,nextColumn);
1623      // keep if >minSize2 or sum >minSize1;
1624      int n=0;
1625      for (int i=0;i<numberBlocks;i++) {
1626        int originalBlock=nextColumn[i];
1627        if (blockCount[originalBlock]<minSize2)
1628          break;
1629        n++;
1630      }
1631      int size=minSize1;
1632      for (int i=n;i<numberBlocks;i++) {
1633        int originalBlock=nextColumn[i];
1634        size-=blockCount[originalBlock];
1635        if (size>0&&i<numberBlocks-1) {
1636          blockCount[originalBlock]=-1;
1637        } else {
1638          size=minSize1;
1639          n++;
1640        }
1641      }
1642      int n2=numberBlocks;
1643      numberBlocks=n;
1644      for (int i=n2-1;i>=0;i--) {
1645        int originalBlock=nextColumn[i];
1646        if (blockCount[originalBlock]>0)
1647          n--;
1648        blockCount[originalBlock]=n;
1649      }
1650      assert (!n);
1651      for (int iRow = 0; iRow < numberRows_; iRow++) {
1652        int iBlock=blockStart[iRow];
1653        if (iBlock>=0) 
1654          blockStart[iRow]=blockCount[iBlock];
1655      }
1656      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1657        int iBlock=columnBlock[iColumn];
1658        if (iBlock>=0) 
1659          columnBlock[iColumn]=blockCount[iBlock];
1660      }
1661      // stick to Clp for now
1662      ClpSimplex ** simplex = new ClpSimplex * [numberBlocks]; 
1663      for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1664        int nRow=0;
1665        int nColumn=0;
1666        for (int iRow=0;iRow<numberRows_;iRow++) {
1667          if (blockStart[iRow]==iBlock) 
1668            blockCount[nRow++]=iRow;
1669        }
1670        for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1671          if (columnBlock[iColumn]==iBlock) 
1672            nextColumn[nColumn++]=iColumn;
1673        }
1674        simplex[iBlock]=new ClpSimplex(this,nRow,blockCount,nColumn,nextColumn);
1675        simplex[iBlock]->setSpecialOptions(simplex[iBlock]->specialOptions()&(~65536));
1676        if (logLevel()<2)
1677          simplex[iBlock]->setLogLevel(0);
1678      }
1679      solveMany(numberBlocks,simplex);
1680      int numberBasic=numberMaster;
1681      int numberStructurals=0;
1682      for (int i=0;i<numberMaster;i++)
1683        abcPivotVariable_[i]=i+firstMaster;
1684      for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
1685        int nRow=0;
1686        int nColumn=0;
1687        // from Clp enum to Abc enum (and bound flip)
1688        unsigned char lookupToAbcSlack[6]={4,6,0/*1*/,1/*0*/,5,7};
1689        unsigned char *  COIN_RESTRICT getStatus = simplex[iBlock]->statusArray()+
1690          simplex[iBlock]->numberColumns();
1691        double *  COIN_RESTRICT solutionSaved=solutionSaved_;
1692        double *  COIN_RESTRICT lowerSaved=lowerSaved_;
1693        double *  COIN_RESTRICT upperSaved=upperSaved_;
1694        for (int iRow=0;iRow<numberRows_;iRow++) {
1695          if (blockStart[iRow]==iBlock) {
1696            unsigned char status=getStatus[nRow++]&7;
1697            AbcSimplex::Status abcStatus=static_cast<AbcSimplex::Status>(lookupToAbcSlack[status]);
1698            if (status!=ClpSimplex::basic) {
1699              double lowerValue=lowerSaved[iRow];
1700              double upperValue=upperSaved[iRow];
1701              if (lowerValue==-COIN_DBL_MAX) {
1702                if(upperValue==COIN_DBL_MAX) {
1703                  // free
1704                  abcStatus=isFree;
1705                } else {
1706                  abcStatus=atUpperBound;
1707                }
1708              } else if (upperValue==COIN_DBL_MAX) {
1709                abcStatus=atLowerBound;
1710              } else if (lowerValue==upperValue) {
1711                abcStatus=isFixed;
1712              }
1713              switch (abcStatus) {
1714              case isFixed:
1715              case atLowerBound:
1716                solutionSaved[iRow]=lowerValue;
1717                break;
1718              case atUpperBound:
1719                solutionSaved[iRow]=upperValue;
1720                break;
1721              default:
1722                break;
1723              }
1724            } else {
1725              // basic
1726              abcPivotVariable_[numberBasic++]=iRow;
1727            }
1728            internalStatus_[iRow]=abcStatus;
1729          }
1730        }
1731        // from Clp enum to Abc enum
1732        unsigned char lookupToAbc[6]={4,6,1,0,5,7};
1733        unsigned char *  COIN_RESTRICT putStatus=internalStatus_+maximumAbcNumberRows_;
1734        getStatus = simplex[iBlock]->statusArray();
1735        solutionSaved+= maximumAbcNumberRows_;
1736        lowerSaved+= maximumAbcNumberRows_;
1737        upperSaved+= maximumAbcNumberRows_;
1738        int numberSaved=numberBasic;
1739        for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1740          if (columnBlock[iColumn]==iBlock) {
1741            unsigned char status=getStatus[nColumn++]&7;
1742            AbcSimplex::Status abcStatus=static_cast<AbcSimplex::Status>(lookupToAbc[status]);
1743            if (status!=ClpSimplex::basic) {
1744              double lowerValue=lowerSaved[iColumn];
1745              double upperValue=upperSaved[iColumn];
1746              if (lowerValue==-COIN_DBL_MAX) {
1747                if(upperValue==COIN_DBL_MAX) {
1748                  // free
1749                  abcStatus=isFree;
1750                } else {
1751                  abcStatus=atUpperBound;
1752                }
1753              } else if (upperValue==COIN_DBL_MAX) {
1754                abcStatus=atLowerBound;
1755              } else if (lowerValue==upperValue) {
1756                abcStatus=isFixed;
1757              } else if (abcStatus==isFree) {
1758                abcStatus=superBasic;
1759              }
1760              switch (abcStatus) {
1761              case isFixed:
1762              case atLowerBound:
1763                solutionSaved[iColumn]=lowerValue;
1764                break;
1765              case atUpperBound:
1766                solutionSaved[iColumn]=upperValue;
1767                break;
1768            default:
1769              break;
1770              }
1771            } else {
1772              // basic
1773              if (numberBasic<numberRows_) 
1774                abcPivotVariable_[numberBasic++]=iColumn+maximumAbcNumberRows_;
1775              else
1776                abcStatus=superBasic;
1777            }
1778            putStatus[iColumn]=abcStatus;
1779          }
1780        }
1781        numberStructurals+=numberBasic-numberSaved;
1782        delete simplex[iBlock];
1783      }
1784#if ABC_NORMAL_DEBUG>0
1785      printf("%d structurals put into basis\n",numberStructurals);
1786#endif
1787      if (numberBasic<numberRows_) {
1788        for (int iRow=0;iRow<numberRows_;iRow++) {
1789          AbcSimplex::Status status=getInternalStatus(iRow);
1790          if (status!=AbcSimplex::basic) {
1791            abcPivotVariable_[numberBasic++]=iRow;
1792            setInternalStatus(iRow,basic);
1793            if (numberBasic==numberRows_)
1794              break;
1795          }
1796        }
1797      }
1798      assert (numberBasic==numberRows_);
1799#if 0
1800      int n2=0;
1801      for (int i=0;i<numberTotal_;i++) {
1802        if (getInternalStatus(i)==basic)
1803          n2++;
1804      }
1805      assert (n2==numberRows_);
1806      std::sort(abcPivotVariable_,abcPivotVariable_+numberRows_);
1807      n2=-1;
1808      for (int i=0;i<numberRows_;i++) {
1809        assert (abcPivotVariable_[i]>n2);
1810        n2=abcPivotVariable_[i];
1811      }
1812#endif
1813      delete [] simplex;
1814      return;
1815    } else {
1816      // try another
1817      type=2;
1818    }
1819  }
1820  if (type==1) {
1821    const double * linearObjective = abcCost_+maximumAbcNumberRows_;
1822    double gamma=0.0;
1823    int numberTotal=numberRows_+numberColumns_;
1824    double * q = new double [numberTotal];
1825    double * v = q+numberColumns_;
1826    int * which = new int [numberTotal+3*numberRows_];
1827    int * ii = which+numberColumns_;
1828    int * r = ii+numberRows_;
1829    int * pivoted = r+numberRows_;
1830    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1831      gamma=CoinMax(gamma,linearObjective[iColumn]);
1832    }
1833    for (int iRow = 0; iRow < numberRows_; iRow++) {
1834      double lowerBound = abcLower_[iRow];
1835      double upperBound = abcUpper_[iRow];
1836      pivoted[iRow]=-1;
1837      ii[iRow]=0;
1838      r[iRow]=0;
1839      v[iRow]=COIN_DBL_MAX;
1840      if (lowerBound==upperBound)
1841        continue;
1842      if (lowerBound<=0.0&&upperBound>=0.0) {
1843        pivoted[iRow]=iRow;
1844        ii[iRow]=1;
1845        r[iRow]=1;
1846      }
1847    }
1848    int nPossible=0;
1849    int lastPossible=0;
1850    double cMaxDiv;
1851    if (gamma)
1852      cMaxDiv=1.0/(1000.0*gamma);
1853    else
1854      cMaxDiv=1.0;
1855    const double * columnLower = abcLower_+maximumAbcNumberRows_;
1856    const double * columnUpper = abcUpper_+maximumAbcNumberRows_;
1857    for (int iPass=0;iPass<3;iPass++) {
1858      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1859        double lowerBound = columnLower[iColumn];
1860        double upperBound = columnUpper[iColumn];
1861        if (lowerBound==upperBound)
1862          continue;
1863        double qValue; 
1864        if (lowerBound > -1.0e20 ) {
1865          if (upperBound < 1.0e20) {
1866            // both
1867            qValue=lowerBound-upperBound;
1868            if (iPass!=2)
1869              qValue=COIN_DBL_MAX;
1870          } else {
1871            // just lower
1872            qValue=lowerBound;
1873            if (iPass!=1)
1874              qValue=COIN_DBL_MAX;
1875          }
1876        } else {
1877          if (upperBound < 1.0e20) {
1878            // just upper
1879            qValue=-upperBound;
1880            if (iPass!=1)
1881              qValue=COIN_DBL_MAX;
1882          } else {
1883            // free
1884            qValue=0.0;
1885            if (iPass!=0)
1886              qValue=COIN_DBL_MAX;
1887          }
1888        }
1889        if (qValue!=COIN_DBL_MAX) {
1890          which[nPossible]=iColumn;
1891          q[nPossible++]=qValue+linearObjective[iColumn]*cMaxDiv;;
1892        }
1893      }
1894      CoinSort_2(q+lastPossible,q+nPossible,which+lastPossible);
1895      lastPossible=nPossible;
1896    }
1897    const double * element = abcMatrix_->getPackedMatrix()->getElements();
1898    const CoinBigIndex * columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
1899    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
1900    const int * row = abcMatrix_->getPackedMatrix()->getIndices();
1901    int nPut=0;
1902    for (int i=0;i<nPossible;i++) {
1903      int iColumn=which[i];
1904      double maxAlpha=0.0;
1905      int kRow=-1;
1906      double alternativeAlpha=0.0;
1907      int jRow=-1;
1908      bool canTake=true;
1909      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1910        int iRow=row[j];
1911        double alpha=fabs(element[j]);
1912        if (alpha>0.01*v[iRow]) {
1913          canTake=false;
1914        } else if (!ii[iRow]&&alpha>alternativeAlpha) {
1915          alternativeAlpha=alpha;
1916          jRow=iRow;
1917        }
1918        if (!r[iRow]&&alpha>maxAlpha) {
1919          maxAlpha=alpha;
1920          kRow=iRow;
1921        }
1922      }
1923      // only really works if scaled
1924      if (maxAlpha>0.99) {
1925        pivoted[kRow]=iColumn+maximumAbcNumberRows_;
1926        v[kRow]=maxAlpha;
1927        ii[kRow]=1;
1928        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1929          int iRow=row[j];
1930          r[iRow]++;
1931        }
1932        nPut++;
1933      } else if (canTake&&jRow>=0) {
1934        pivoted[jRow]=iColumn+maximumAbcNumberRows_;
1935        v[jRow]=maxAlpha;
1936        ii[jRow]=1;
1937        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
1938          int iRow=row[j];
1939          r[iRow]++;
1940        }
1941        nPut++;
1942      }
1943    }
1944    for (int iRow=0;iRow<numberRows_;iRow++) {
1945      int iSequence=pivoted[iRow];
1946      if (iSequence>=0&&iSequence<numberColumns_) {
1947        abcPivotVariable_[iRow]=iSequence;
1948        if (fabs(abcLower_[iRow])<fabs(abcUpper_[iRow])) {
1949          setInternalStatus(iRow,atLowerBound);
1950          abcSolution_[iRow]=abcLower_[iRow];
1951          solutionSaved_[iRow]=abcLower_[iRow];
1952        } else {
1953          setInternalStatus(iRow,atUpperBound);
1954          abcSolution_[iRow]=abcUpper_[iRow];
1955          solutionSaved_[iRow]=abcUpper_[iRow];
1956        }
1957        setInternalStatus(iSequence,basic);
1958      }
1959    }
1960#if ABC_NORMAL_DEBUG>0
1961    printf("%d put into basis\n",nPut);
1962#endif
1963    delete [] q;
1964    delete [] which;
1965  } else if (type==2) {
1966    //return;
1967    int numberBad=0;
1968    CoinAbcMemcpy(abcDj_,abcCost_,numberTotal_);
1969    // Work on savedSolution and current
1970    int iVector=getAvailableArray();
1971#define ALLOW_BAD_DJS
1972#ifdef ALLOW_BAD_DJS
1973    double * modDj = usefulArray_[iVector].denseVector();
1974#endif
1975    for (int iSequence=maximumAbcNumberRows_;iSequence<numberTotal_;iSequence++) {
1976      double dj=abcDj_[iSequence];
1977      modDj[iSequence]=0.0;
1978      if (getInternalStatus(iSequence)==atLowerBound) {
1979        if (dj<-dualTolerance_) {
1980          double costThru=-dj*(abcUpper_[iSequence]-abcLower_[iSequence]);
1981          if (costThru<dualBound_) {
1982            // flip
1983            setInternalStatus(iSequence,atUpperBound);
1984            solutionSaved_[iSequence]=abcUpper_[iSequence];
1985            abcSolution_[iSequence]=abcUpper_[iSequence];
1986          } else {
1987            numberBad++;
1988#ifdef ALLOW_BAD_DJS
1989            modDj[iSequence]=dj;
1990            dj=0.0;
1991#else
1992            break;
1993#endif
1994          }
1995        } else {
1996          dj=CoinMax(dj,0.0);
1997        }
1998      } else if (getInternalStatus(iSequence)==atLowerBound) {
1999        if (dj>dualTolerance_) {
2000          double costThru=dj*(abcUpper_[iSequence]-abcLower_[iSequence]);
2001          if (costThru<dualBound_) {
2002            assert (abcUpper_[iSequence]-abcLower_[iSequence]<1.0e10);
2003            // flip
2004            setInternalStatus(iSequence,atLowerBound);
2005            solutionSaved_[iSequence]=abcLower_[iSequence];
2006            abcSolution_[iSequence]=abcLower_[iSequence];
2007          } else {
2008            numberBad++;
2009#ifdef ALLOW_BAD_DJS
2010            modDj[iSequence]=dj;
2011            dj=0.0;
2012#else
2013            break;
2014#endif
2015          }
2016        } else {
2017          dj=CoinMin(dj,0.0);
2018        }
2019      } else {
2020        if (fabs(dj)<dualTolerance_) {
2021          dj=0.0;
2022        } else {
2023          numberBad++;
2024#ifdef ALLOW_BAD_DJS
2025          modDj[iSequence]=dj;
2026          dj=0.0;
2027#else
2028          break;
2029#endif
2030        }
2031      }
2032      abcDj_[iSequence]=dj;
2033    }
2034#ifndef ALLOW_BAD_DJS
2035    if (numberBad) {
2036      //CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_);
2037      return ;
2038    }
2039#endif
2040    if (!abcMatrix_->gotRowCopy())
2041      abcMatrix_->createRowCopy();
2042    //const double * element = abcMatrix_->getPackedMatrix()->getElements();
2043    const CoinBigIndex * columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts()-maximumAbcNumberRows_;
2044    const int * row = abcMatrix_->getPackedMatrix()->getIndices();
2045    const double * elementByRow = abcMatrix_->rowElements();
2046    const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2047    const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
2048    const int * column = abcMatrix_->rowColumns();
2049    CoinAbcMemset0(solutionBasic_,numberRows_);
2050    CoinAbcMemcpy(lowerBasic_,abcLower_,numberRows_);
2051    CoinAbcMemcpy(upperBasic_,abcUpper_,numberRows_);
2052    abcMatrix_->timesIncludingSlacks(-1.0,solutionSaved_,solutionBasic_);
2053    const double multiplier[] = { 1.0, -1.0};
2054    int nBasic=0;
2055    int * index=usefulArray_[iVector].getIndices();
2056    int iVector2=getAvailableArray();
2057    int * index2=usefulArray_[iVector2].getIndices();
2058    double * sort =usefulArray_[iVector2].denseVector();
2059    int average = CoinMax(5,rowEnd[numberRows_-1]/(8*numberRows_));
2060    int nPossible=0;
2061    if (numberRows_>10000) {
2062      // allow more
2063      for (int iRow=0;iRow<numberRows_;iRow++) {
2064        double solutionValue=solutionBasic_[iRow];
2065        double infeasibility=0.0;
2066        double lowerValue=lowerBasic_[iRow];
2067        double upperValue=upperBasic_[iRow];
2068        if (solutionValue<lowerValue-primalTolerance_) {
2069          infeasibility=-(lowerValue-solutionValue);
2070        } else if (solutionValue>upperValue+primalTolerance_) {
2071          infeasibility=upperValue-solutionValue;
2072        }
2073        int length=rowEnd[iRow]-rowStart[iRow];
2074        if (infeasibility) 
2075          index2[nPossible++]=length;
2076      }
2077      std::sort(index2,index2+nPossible);
2078      // see how much we need to get numberRows/10 or nPossible/3
2079      average=CoinMax(average,index2[CoinMin(numberRows_/10,nPossible/3)]);
2080      nPossible=0;
2081    }
2082    for (int iRow=0;iRow<numberRows_;iRow++) {
2083      double solutionValue=solutionBasic_[iRow];
2084      double infeasibility=0.0;
2085      double lowerValue=lowerBasic_[iRow];
2086      double upperValue=upperBasic_[iRow];
2087      if (solutionValue<lowerValue-primalTolerance_) {
2088        infeasibility=-(lowerValue-solutionValue);
2089      } else if (solutionValue>upperValue+primalTolerance_) {
2090        infeasibility=upperValue-solutionValue;
2091      }
2092      int length=rowEnd[iRow]-rowStart[iRow];
2093      if (infeasibility&&length<average) {
2094        index2[nPossible]=iRow;
2095        sort[nPossible++]=1.0e5*infeasibility-iRow;
2096        //sort[nPossible++]=1.0e9*length-iRow;//infeasibility;
2097      }
2098    }
2099    CoinSort_2(sort,sort+nPossible,index2);
2100    for (int iWhich=0;iWhich<nPossible;iWhich++) {
2101      int iRow = index2[iWhich];
2102      sort[iWhich]=0.0;
2103      if (abcDj_[iRow])
2104        continue; // marked as invalid
2105      double solutionValue=solutionBasic_[iRow];
2106      double infeasibility=0.0;
2107      double lowerValue=lowerBasic_[iRow];
2108      double upperValue=upperBasic_[iRow];
2109      if (solutionValue<lowerValue-primalTolerance_) {
2110        infeasibility=lowerValue-solutionValue;
2111      } else if (solutionValue>upperValue+primalTolerance_) {
2112        infeasibility=upperValue-solutionValue;
2113      }
2114      assert (infeasibility);
2115      double direction = infeasibility >0 ? 1.0 : -1.0;
2116      infeasibility *= direction;
2117      int whichColumn=-1;
2118      double upperTheta=1.0e30;
2119      int whichColumn2=-1;
2120      double upperTheta2=1.0e30;
2121      double costThru=0.0;
2122      int nThru=0;
2123      for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
2124        int iSequence=column[j];
2125        assert (iSequence>=maximumAbcNumberRows_);
2126        double dj = abcDj_[iSequence];
2127        double tableauValue=-elementByRow[j]*direction;
2128        unsigned char iStatus=internalStatus_[iSequence]&7;
2129        if ((iStatus&4)==0) {
2130          double mult = multiplier[iStatus];
2131          double alpha = tableauValue * mult;
2132          double oldValue = dj * mult;
2133          assert (oldValue>-1.0e-2);
2134          double value = oldValue - upperTheta * alpha;
2135          if (value < 0.0) {
2136            upperTheta2=upperTheta;
2137            whichColumn2=whichColumn;
2138            costThru=alpha*(abcUpper_[iSequence]-abcLower_[iSequence]);
2139            nThru=0;
2140            upperTheta = oldValue / alpha;
2141            whichColumn=iSequence;
2142          } else if (oldValue-upperTheta2*alpha<0.0) {
2143            costThru+=alpha*(abcUpper_[iSequence]-abcLower_[iSequence]);
2144            index[nThru++]=iSequence;
2145          }
2146        } else if (iStatus<6) {
2147          upperTheta=-1.0;
2148          upperTheta2=elementByRow[j];
2149          whichColumn=iSequence;
2150          break;
2151        } 
2152      }
2153      if (whichColumn<0)
2154        continue;
2155      if (upperTheta!=-1.0) {
2156        assert (upperTheta>=0.0);
2157        if (costThru<infeasibility&&whichColumn2>=0) {
2158          index[nThru++]=whichColumn;
2159          for (int i=0;i<nThru;i++) {
2160            int iSequence=index[i];
2161            assert (abcUpper_[iSequence]-abcLower_[iSequence]<1.0e10);
2162            // flip and use previous
2163            if (getInternalStatus(iSequence)==atLowerBound) {
2164              setInternalStatus(iSequence,atUpperBound);
2165              solutionSaved_[iSequence]=abcUpper_[iSequence];
2166              abcSolution_[iSequence]=abcUpper_[iSequence];
2167            } else {
2168              setInternalStatus(iSequence,atLowerBound);
2169              solutionSaved_[iSequence]=abcLower_[iSequence];
2170              abcSolution_[iSequence]=abcLower_[iSequence];
2171            }
2172          }
2173          whichColumn=whichColumn2;
2174          upperTheta=upperTheta2;
2175        }
2176      } else {
2177        // free coming in
2178        upperTheta=(abcDj_[whichColumn]*direction)/upperTheta2;
2179      }
2180      // update djs
2181      upperTheta *= -direction;
2182      for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
2183        int iSequence=column[j];
2184        double dj = abcDj_[iSequence];
2185        double tableauValue=elementByRow[j];
2186        dj -= upperTheta*tableauValue;
2187        unsigned char iStatus=internalStatus_[iSequence]&7;
2188        if ((iStatus&4)==0) {
2189          if (!iStatus) {
2190            assert (dj>-1.0e-2);
2191            dj = CoinMax(dj,0.0);
2192#ifdef ALLOW_BAD_DJS
2193            if (numberBad&&modDj[iSequence]) {
2194              double bad=modDj[iSequence];
2195              if (dj+bad>=0.0) {
2196                numberBad--;
2197                modDj[iSequence]=0.0;
2198                dj += bad;
2199              } else {
2200                modDj[iSequence]+=dj;
2201                dj =0.0;
2202              }
2203            }
2204#endif
2205          } else {
2206            assert (dj<1.0e-2);
2207            dj = CoinMin(dj,0.0);
2208#ifdef ALLOW_BAD_DJS
2209            if (numberBad&&modDj[iSequence]) {
2210              double bad=modDj[iSequence];
2211              if (dj+bad<=0.0) {
2212                numberBad--;
2213                modDj[iSequence]=0.0;
2214                dj += bad;
2215              } else {
2216                modDj[iSequence]+=dj;
2217                dj =0.0;
2218              }
2219            }
2220#endif
2221          }
2222        } else if (iStatus<6) {
2223          assert (fabs(dj)<1.0e-4);
2224          dj = 0.0;
2225        }
2226        abcDj_[iSequence]=dj;
2227      }
2228      // do basis
2229      if (direction>0.0) {
2230        if (upperBasic_[iRow]>lowerBasic_[iRow]) 
2231          setInternalStatus(iRow,atLowerBound);
2232        else
2233          setInternalStatus(iRow,isFixed);
2234        solutionSaved_[iRow]=abcLower_[iRow];
2235        abcSolution_[iRow]=abcLower_[iRow];
2236      } else {
2237        if (upperBasic_[iRow]>lowerBasic_[iRow]) 
2238          setInternalStatus(iRow,atUpperBound);
2239        else
2240          setInternalStatus(iRow,isFixed);
2241        solutionSaved_[iRow]=abcUpper_[iRow];
2242        abcSolution_[iRow]=abcUpper_[iRow];
2243      }
2244      setInternalStatus(whichColumn,basic);
2245      abcPivotVariable_[iRow]=whichColumn;
2246      nBasic++;
2247      // mark rows
2248      for (CoinBigIndex j=columnStart[whichColumn];j<columnStart[whichColumn+1];j++) {
2249        int jRow=row[j];
2250        abcDj_[jRow]=1.0;
2251      }
2252    }
2253#ifdef ALLOW_BAD_DJS
2254    CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_);
2255#endif
2256    setAvailableArray(iVector);
2257    setAvailableArray(iVector2);
2258#if ABC_NORMAL_DEBUG>0
2259    printf("dual crash put %d in basis\n",nBasic);
2260#endif
2261  } else {
2262    assert ((stateOfProblem_&VALUES_PASS2)!=0);
2263    // The idea is to put as many likely variables into basis as possible
2264    int n=0;
2265    int iVector=getAvailableArray();
2266    int * index=usefulArray_[iVector].getIndices();
2267    double * array = usefulArray_[iVector].denseVector();
2268    int iVector2=getAvailableArray();
2269    int * index2=usefulArray_[iVector].getIndices();
2270    for (int iSequence=0;iSequence<numberTotal_;iSequence++) {
2271      double dj = djSaved_[iSequence];
2272      double value = solution_[iSequence];
2273      double lower = abcLower_[iSequence];
2274      double upper = abcUpper_[iSequence];
2275      double gapUp=CoinMin(1.0e3,upper-value);
2276      assert (gapUp>=-1.0e-3);
2277      gapUp=CoinMax(gapUp,0.0);
2278      double gapDown=CoinMin(1.0e3,value-lower);
2279      assert (gapDown>=-1.0e-3);
2280      gapDown=CoinMax(gapDown,0.0);
2281      double measure = (CoinMin(gapUp,gapDown)+1.0e-6)/(fabs(dj)+1.0e-6);
2282      if (gapUp<primalTolerance_*10.0&&dj<dualTolerance_) {
2283        // set to ub
2284        setInternalStatus(iSequence,atUpperBound);
2285        solutionSaved_[iSequence]=abcUpper_[iSequence];
2286        abcSolution_[iSequence]=abcUpper_[iSequence];
2287      } else if (gapDown<primalTolerance_*10.0&&dj>-dualTolerance_) {
2288        // set to lb
2289        setInternalStatus(iSequence,atLowerBound);
2290        solutionSaved_[iSequence]=abcLower_[iSequence];
2291        abcSolution_[iSequence]=abcLower_[iSequence];
2292      } else if (upper>lower) {
2293        // set to nearest
2294        if (gapUp<gapDown) {
2295          // set to ub
2296          setInternalStatus(iSequence,atUpperBound);
2297          solutionSaved_[iSequence]=abcUpper_[iSequence];
2298          abcSolution_[iSequence]=abcUpper_[iSequence];
2299        } else {
2300          // set to lb
2301          setInternalStatus(iSequence,atLowerBound);
2302          solutionSaved_[iSequence]=abcLower_[iSequence];
2303          abcSolution_[iSequence]=abcLower_[iSequence];
2304        }
2305        array[n]=-measure;
2306        index[n++]=iSequence;
2307      }
2308    }
2309    // set slacks basic
2310    memset(internalStatus_,6,numberRows_);
2311    CoinSort_2(solution_,solution_+n,index);
2312    CoinAbcMemset0(array,n);
2313    for (int i=0;i<numberRows_;i++)
2314      index2[i]=0;
2315    // first put all possible slacks in
2316    int n2=0;
2317    for (int i=0;i<n;i++) {
2318      int iSequence=index[i];
2319      if (iSequence<numberRows_) {
2320        index2[iSequence]=numberRows_;
2321      } else {
2322        index[n2++]=iSequence;
2323      }
2324    }
2325    n=n2;
2326    int numberIn=0;
2327    const CoinBigIndex * columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts()-
2328      maximumAbcNumberRows_;
2329    //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
2330    const int * row = abcMatrix_->getPackedMatrix()->getIndices();
2331    if (!abcMatrix_->gotRowCopy())
2332      abcMatrix_->createRowCopy();
2333    //const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2334    //const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
2335    for (int i=0;i<n;i++) {
2336      int iSequence=index[i];
2337      int bestRow=-1;
2338      int bestCount=numberRows_+1;
2339      for (CoinBigIndex j=columnStart[iSequence];j<columnStart[iSequence+1];j++) {
2340        int iRow=row[j];
2341    if (!abcMatrix_->gotRowCopy())
2342      abcMatrix_->createRowCopy();
2343    const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2344    const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
2345        if (!index2[iRow]) {
2346          int length=rowEnd[iRow]-rowStart[iRow];
2347          if (length<bestCount) {
2348            bestCount=length;
2349            bestRow=iRow;
2350          }
2351        }
2352      }
2353      if (bestRow>=0) {
2354        numberIn++;
2355        for (CoinBigIndex j=columnStart[iSequence];j<columnStart[iSequence+1];j++) {
2356          int iRow=row[j];
2357          index2[iRow]++;
2358        }
2359        setInternalStatus(iSequence,basic);
2360        abcPivotVariable_[bestRow]=iSequence;
2361        double dj = djSaved_[bestRow];
2362        double value = solution_[bestRow];
2363        double lower = abcLower_[bestRow];
2364        double upper = abcUpper_[bestRow];
2365        double gapUp=CoinMax(CoinMin(1.0e3,upper-value),0.0);
2366        double gapDown=CoinMax(CoinMin(1.0e3,value-lower),0.0);
2367        //double measure = (CoinMin(gapUp,gapDown)+1.0e-6)/(fabs(dj)+1.0e-6);
2368        if (gapUp<primalTolerance_*10.0&&dj<dualTolerance_) {
2369          // set to ub
2370          setInternalStatus(bestRow,atUpperBound);
2371          solutionSaved_[bestRow]=abcUpper_[bestRow];
2372          abcSolution_[bestRow]=abcUpper_[bestRow];
2373        } else if (gapDown<primalTolerance_*10.0&&dj>-dualTolerance_) {
2374          // set to lb
2375          setInternalStatus(bestRow,atLowerBound);
2376          solutionSaved_[bestRow]=abcLower_[bestRow];
2377          abcSolution_[bestRow]=abcLower_[bestRow];
2378        } else if (upper>lower) {
2379          // set to nearest
2380          if (gapUp<gapDown) {
2381            // set to ub
2382            setInternalStatus(bestRow,atUpperBound);
2383            solutionSaved_[bestRow]=abcUpper_[bestRow];
2384            abcSolution_[bestRow]=abcUpper_[bestRow];
2385          } else {
2386            // set to lb
2387            setInternalStatus(bestRow,atLowerBound);
2388            solutionSaved_[bestRow]=abcLower_[bestRow];
2389            abcSolution_[bestRow]=abcLower_[bestRow];
2390          }
2391        }
2392      }
2393    }
2394#if ABC_NORMAL_DEBUG>0
2395    printf("Idiot crash put %d in basis\n",numberIn);
2396#endif
2397    setAvailableArray(iVector);
2398    setAvailableArray(iVector2);
2399    delete [] solution_;
2400    solution_=NULL;
2401  }
2402}
2403/* Puts more stuff in basis
2404   1 bit set - do even if basis exists
2405   2 bit set - don't bother staying triangular
2406*/
2407void 
2408AbcSimplex::putStuffInBasis(int type)
2409{
2410  int i;
2411  for (i=0;i<numberRows_;i++) {
2412    if (getInternalStatus(i)!=basic)
2413      break;
2414  }
2415  if (i<numberRows_&&(type&1)==0)
2416    return;
2417  int iVector=getAvailableArray();
2418  // Column part is mustColumnIn
2419  int * COIN_RESTRICT mustRowOut = usefulArray_[iVector].getIndices();
2420  if (!abcMatrix_->gotRowCopy())
2421    abcMatrix_->createRowCopy();
2422  const double * elementByRow = abcMatrix_->rowElements();
2423  const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2424  const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
2425  const int * column = abcMatrix_->rowColumns();
2426  for (int i=0;i<numberTotal_;i++)
2427    mustRowOut[i]=-1;
2428  int nPossible=0;
2429  // find equality rows with effective nonzero rhs
2430  for (int iRow=0;iRow<numberRows_;iRow++) {
2431    if (abcUpper_[iRow]>abcLower_[iRow]||getInternalStatus(iRow)!=basic) {
2432      mustRowOut[iRow]=-2;
2433      continue;
2434    }
2435    int chooseColumn[2]={-1,-1};
2436    for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
2437      int iColumn=column[j];
2438      if (elementByRow[j]>0.0) {
2439        if (chooseColumn[0]==-1)
2440          chooseColumn[0]=iColumn;
2441        else
2442          chooseColumn[0]=-2;
2443      } else {
2444        if (chooseColumn[1]==-1)
2445          chooseColumn[1]=iColumn;
2446        else
2447          chooseColumn[1]=-2;
2448      }
2449    }
2450    for (int iTry=0;iTry<2;iTry++) {
2451      int jColumn=chooseColumn[iTry];
2452      if (jColumn>=0&&getInternalStatus(jColumn)!=basic) {
2453        // see if has to be basic
2454        double lowerValue=-abcUpper_[iRow]; // check swap
2455        double upperValue=-abcLower_[iRow];
2456        int lowerInf=0;
2457        int upperInf=0;
2458        double alpha=0.0;
2459        for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
2460          int iColumn=column[j];
2461          if (iColumn!=jColumn) {
2462            if (abcLower_[iColumn]>-1.0e20)
2463              lowerValue -= abcLower_[iColumn]*elementByRow[j];
2464            else
2465              lowerInf ++;
2466            if (abcUpper_[iColumn]<1.0e20)
2467              upperValue -= abcUpper_[iColumn]*elementByRow[j];
2468            else
2469              upperInf ++;
2470          } else {
2471            alpha=elementByRow[j];
2472          }
2473        }
2474        // find values column must lie between (signs again)
2475        if (upperInf)
2476          upperValue=COIN_DBL_MAX;
2477        else
2478          upperValue /=alpha; 
2479        if (lowerInf)
2480          lowerValue=-COIN_DBL_MAX;
2481        else
2482          lowerValue /=alpha;
2483        if (iTry) {
2484          // swap
2485          double temp=lowerValue;
2486          lowerValue=upperValue;
2487          upperValue=temp;
2488        }
2489        if (lowerValue>abcLower_[jColumn]+10.0*primalTolerance_&&
2490            upperValue<abcUpper_[jColumn]-10.0*primalTolerance_) {
2491          nPossible++;
2492          if (mustRowOut[jColumn]>=0) {
2493            // choose one ???
2494            //printf("Column %d already selected on row %d now on %d\n",
2495            //     jColumn,mustRowOut[jColumn],iRow);
2496            continue;
2497          }
2498          mustRowOut[jColumn]=iRow;
2499          mustRowOut[iRow]=jColumn;
2500        }
2501      }
2502    }
2503  }
2504  if (nPossible) {
2505#if ABC_NORMAL_DEBUG>0
2506    printf("%d possible candidates\n",nPossible);
2507#endif
2508    if ((type&2)==0) {
2509      // triangular
2510      int iVector2=getAvailableArray();
2511      int * COIN_RESTRICT counts = usefulArray_[iVector2].getIndices();
2512      CoinAbcMemset0(counts,numberRows_);
2513      for (int iRow=0;iRow<numberRows_;iRow++) {
2514        int n=0;
2515        for (CoinBigIndex j=rowStart[iRow];j<rowEnd[iRow];j++) {
2516          int iColumn=column[j];
2517          if (getInternalStatus(iColumn)==basic)
2518            n++;
2519        }
2520        counts[iRow]=n;
2521      }
2522      const CoinBigIndex * columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts()
2523        -maximumAbcNumberRows_;
2524      const int * row = abcMatrix_->getPackedMatrix()->getIndices();
2525      for (int iRow=0;iRow<numberRows_;iRow++) {
2526        if (!counts[iRow]) {
2527          int iColumn=mustRowOut[iRow];
2528          if (iColumn>=0) {
2529            setInternalStatus(iRow,isFixed);
2530            solutionSaved_[iRow]=abcLower_[iRow];
2531            setInternalStatus(iColumn,basic);
2532            for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
2533              int jRow=row[j];
2534              counts[jRow]++;
2535            }
2536          }
2537        }
2538      }
2539      setAvailableArray(iVector2);
2540    } else {
2541      // all
2542      for (int iRow=0;iRow<numberRows_;iRow++) {
2543        int iColumn=mustRowOut[iRow];
2544        if (iColumn>=0) {
2545          setInternalStatus(iRow,isFixed);
2546          solutionSaved_[iRow]=abcLower_[iRow];
2547          setInternalStatus(iColumn,basic);
2548        }
2549      }
2550    }
2551    // redo pivot array
2552    int numberBasic=0;
2553    for (int iSequence=0;iSequence<numberTotal_;iSequence++) {
2554      if (getInternalStatus(iSequence)==basic) 
2555        abcPivotVariable_[numberBasic++]=iSequence;
2556    }
2557    assert (numberBasic==numberRows_);
2558  }
2559  setAvailableArray(iVector);
2560}
2561// Computes nonbasic cost and total cost
2562void 
2563AbcSimplex::computeObjective ()
2564{
2565  sumNonBasicCosts_ = 0.0;
2566  rawObjectiveValue_ = 0.0;
2567  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2568    double value = abcSolution_[iSequence]*abcCost_[iSequence];
2569    rawObjectiveValue_ += value;
2570    if (getInternalStatus(iSequence)!=basic)
2571      sumNonBasicCosts_ += value;
2572  }
2573  setClpSimplexObjectiveValue();
2574}
2575// This sets largest infeasibility and most infeasible
2576void
2577AbcSimplex::checkPrimalSolution(bool justBasic)
2578{
2579  rawObjectiveValue_ = CoinAbcInnerProduct(costBasic_,numberRows_,solutionBasic_);
2580  //rawObjectiveValue_ += sumNonBasicCosts_; 
2581  setClpSimplexObjectiveValue();
2582  // now look at primal solution
2583  sumPrimalInfeasibilities_ = 0.0;
2584  numberPrimalInfeasibilities_ = 0;
2585  double primalTolerance = primalTolerance_;
2586  double relaxedTolerance = primalTolerance_;
2587  // we can't really trust infeasibilities if there is primal error
2588  double error = CoinMin(1.0e-2, largestPrimalError_);
2589  // allow tolerance at least slightly bigger than standard
2590  relaxedTolerance = relaxedTolerance +  error;
2591  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2592  const double *  COIN_RESTRICT lowerBasic = lowerBasic_;
2593  const double *  COIN_RESTRICT upperBasic = upperBasic_;
2594  const double *  COIN_RESTRICT solutionBasic = solutionBasic_;
2595  if (justBasic) {
2596    for (int iRow = 0; iRow < numberRows_; iRow++) {
2597      double infeasibility = 0.0;
2598      if (solutionBasic[iRow] > upperBasic[iRow]) {
2599        infeasibility = solutionBasic[iRow] - upperBasic[iRow];
2600      } else if (solutionBasic[iRow] < lowerBasic[iRow]) {
2601        infeasibility = lowerBasic[iRow] - solutionBasic[iRow];
2602      }
2603      if (infeasibility > primalTolerance) {
2604        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2605        if (infeasibility > relaxedTolerance)
2606          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2607        numberPrimalInfeasibilities_ ++;
2608      }
2609    }
2610  } else {
2611    for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2612      double infeasibility = 0.0;
2613      if (abcSolution_[iSequence] > abcUpper_[iSequence]) {
2614        infeasibility = abcSolution_[iSequence] - abcUpper_[iSequence];
2615      } else if (abcSolution_[iSequence] < abcLower_[iSequence]) {
2616        infeasibility = abcLower_[iSequence] - abcSolution_[iSequence];
2617      }
2618      if (infeasibility > primalTolerance) {
2619        //assert (getInternalStatus(iSequence)==basic);
2620        sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2621        if (infeasibility > relaxedTolerance)
2622          sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2623        numberPrimalInfeasibilities_ ++;
2624      }
2625    }
2626  }
2627}
2628void
2629AbcSimplex::checkDualSolution()
2630{
2631 
2632  sumDualInfeasibilities_ = 0.0;
2633  numberDualInfeasibilities_ = 0;
2634  int firstFreePrimal = -1;
2635  int firstFreeDual = -1;
2636  int numberSuperBasicWithDj = 0;
2637  bestPossibleImprovement_ = 0.0;
2638  // we can't really trust infeasibilities if there is dual error
2639  double error = CoinMin(1.0e-2, largestDualError_);
2640  // allow tolerance at least slightly bigger than standard
2641  double relaxedTolerance = currentDualTolerance_ +  error;
2642  // allow bigger tolerance for possible improvement
2643  double possTolerance = 5.0 * relaxedTolerance;
2644  sumOfRelaxedDualInfeasibilities_ = 0.0;
2645  int numberNeedToMove=0;
2646  for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2647    if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
2648      // not basic
2649      double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
2650      double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
2651      double value = abcDj_[iSequence];
2652      if (distanceUp > primalTolerance_) {
2653        // Check if "free"
2654        if (distanceDown <= primalTolerance_) {
2655          // should not be negative
2656          if (value < 0.0) {
2657            value = - value;
2658            if (value > currentDualTolerance_) {
2659              sumDualInfeasibilities_ += value - currentDualTolerance_;
2660              if (value > possTolerance)
2661                bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2662              if (value > relaxedTolerance)
2663                sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2664              numberDualInfeasibilities_ ++;
2665            }
2666          }
2667        } else {
2668          // free
2669          value=fabs(value);
2670          if (value > 1.0 * relaxedTolerance) {
2671            numberSuperBasicWithDj++;
2672            if (firstFreeDual < 0)
2673              firstFreeDual = iSequence;
2674          }
2675          if (value > currentDualTolerance_||getInternalStatus(iSequence)!=AbcSimplex::isFree) {
2676            numberNeedToMove++;
2677            if (firstFreePrimal < 0)
2678              firstFreePrimal = iSequence;
2679            sumDualInfeasibilities_ += value - currentDualTolerance_;
2680            if (value > possTolerance)
2681              bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2682            if (value > relaxedTolerance)
2683              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2684            numberDualInfeasibilities_ ++;
2685          }
2686        }
2687      } else if (distanceDown > primalTolerance_) {
2688        // should not be positive
2689        if (value > 0.0) {
2690          if (value > currentDualTolerance_) {
2691            sumDualInfeasibilities_ += value - currentDualTolerance_;
2692            if (value > possTolerance)
2693              bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2694            if (value > relaxedTolerance)
2695              sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2696            numberDualInfeasibilities_ ++;
2697          }
2698        }
2699      }
2700    }
2701  }
2702  numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_-numberSuperBasicWithDj;
2703  if (algorithm_ < 0 && firstFreeDual >= 0) {
2704    // dual
2705    firstFree_ = firstFreeDual;
2706  } else if (numberSuperBasicWithDj||numberNeedToMove) {
2707    //(abcProgress_.lastIterationNumber(0) <= 0)) {
2708    firstFree_ = firstFreePrimal;
2709    if (!sumDualInfeasibilities_)
2710      sumDualInfeasibilities_=1.0e-8;
2711  }
2712}
2713/* This sets sum and number of infeasibilities (Dual and Primal) */
2714void
2715AbcSimplex::checkBothSolutions()
2716{
2717    // old way
2718  checkPrimalSolution(true);
2719  checkDualSolution();
2720}
2721/*
2722  Unpacks one column of the matrix into indexed array
2723*/
2724void 
2725AbcSimplex::unpack(CoinIndexedVector & rowArray, int sequence) const
2726{ 
2727  abcMatrix_->unpack(rowArray,sequence);
2728}
2729// Sets row pivot choice algorithm in dual
2730void
2731AbcSimplex::setDualRowPivotAlgorithm(AbcDualRowPivot & choice)
2732{
2733  delete abcDualRowPivot_;
2734  abcDualRowPivot_ = choice.clone(true);
2735  abcDualRowPivot_->setModel(this);
2736}
2737// Sets column pivot choice algorithm in primal
2738void
2739AbcSimplex::setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot & choice)
2740{
2741  delete abcPrimalColumnPivot_;
2742  abcPrimalColumnPivot_ = choice.clone(true);
2743  abcPrimalColumnPivot_->setModel(this);
2744}
2745void
2746AbcSimplex::setFactorization( AbcSimplexFactorization & factorization)
2747{
2748  if (abcFactorization_)
2749    abcFactorization_->setFactorization(factorization);
2750  else
2751    abcFactorization_ = new AbcSimplexFactorization(factorization,
2752                                          numberRows_);
2753}
2754
2755// Swaps factorization
2756AbcSimplexFactorization *
2757AbcSimplex::swapFactorization( AbcSimplexFactorization * factorization)
2758{
2759  AbcSimplexFactorization * swap = abcFactorization_;
2760  abcFactorization_ = factorization;
2761  return swap;
2762}
2763
2764/* Tightens primal bounds to make dual faster.  Unless
2765   fixed, bounds are slightly looser than they could be.
2766   This is to make dual go faster and is probably not needed
2767   with a presolve.  Returns non-zero if problem infeasible
2768*/
2769int
2770AbcSimplex::tightenPrimalBounds()
2771{
2772  int tightenType=1;
2773  if (maximumIterations()>=1000000&&maximumIterations()<1000010)
2774    tightenType=maximumIterations()-1000000;
2775  if (!tightenType)
2776    return 0;
2777  if (integerType_) {
2778#if ABC_NORMAL_DEBUG>0
2779    printf("Redo tighten to relax by 1 for integers (and don't be shocked by infeasibility)\n");
2780#endif
2781    return 0;
2782  }
2783  // This needs to work on scaled matrix - then replace   
2784  // Get a row copy in standard format
2785  CoinPackedMatrix copy;
2786  copy.setExtraGap(0.0);
2787  copy.setExtraMajor(0.0);
2788  copy.reverseOrderedCopyOf(*abcMatrix_->matrix());
2789  // get matrix data pointers
2790  const int *  COIN_RESTRICT column = copy.getIndices();
2791  const CoinBigIndex *  COIN_RESTRICT rowStart = copy.getVectorStarts();
2792  const int *  COIN_RESTRICT rowLength = copy.getVectorLengths();
2793  const double *  COIN_RESTRICT element = copy.getElements();
2794  int numberChanged = 1, iPass = 0;
2795  double large = largeValue(); // treat bounds > this as infinite
2796#ifndef NDEBUG
2797  double large2 = 1.0e10 * large;
2798#endif
2799  int numberInfeasible = 0;
2800  int totalTightened = 0;
2801 
2802  double tolerance = primalTolerance();
2803 
2804 
2805  // A copy of bounds is up at top
2806  translate(0); // move down
2807 
2808#define MAXPASS 10
2809 
2810  // Loop round seeing if we can tighten bounds
2811  // Would be faster to have a stack of possible rows
2812  // and we put altered rows back on stack
2813  int numberCheck = -1;
2814  // temp swap signs so can use existing code
2815  double *  COIN_RESTRICT rowLower=abcLower_;
2816  double *  COIN_RESTRICT rowUpper=abcUpper_;
2817  for (int iRow=0;iRow<numberRows_;iRow++) {
2818    double value=-rowLower[iRow];
2819    rowLower[iRow]=-rowUpper[iRow];
2820    rowUpper[iRow]=value;
2821  }
2822  // Keep lower and upper for rows
2823  int whichArray[5];
2824  for (int i=0;i<5;i++)
2825    whichArray[i]=getAvailableArray();
2826  double * COIN_RESTRICT minRhs = usefulArray_[whichArray[0]].denseVector();
2827  double * COIN_RESTRICT maxRhs = usefulArray_[whichArray[1]].denseVector();
2828  int * COIN_RESTRICT changedList = usefulArray_[whichArray[0]].getIndices(); 
2829  int * COIN_RESTRICT changedListNext = usefulArray_[whichArray[1]].getIndices(); 
2830  char * COIN_RESTRICT changed = reinterpret_cast<char*>(usefulArray_[whichArray[2]].getIndices()); 
2831  // -1 no infinite, -2 more than one infinite >=0 column
2832  int * COIN_RESTRICT whichInfiniteDown = usefulArray_[whichArray[3]].getIndices(); 
2833  int * COIN_RESTRICT whichInfiniteUp = usefulArray_[whichArray[3]].getIndices(); 
2834  int numberToDoNext=0;
2835  for (int iRow=0;iRow<numberRows_;iRow++) {
2836    if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
2837      changedListNext[numberToDoNext++]=iRow;
2838      whichInfiniteDown[iRow]=-1;
2839      whichInfiniteUp[iRow]=-1;
2840    } else {
2841      minRhs[iRow]=-COIN_DBL_MAX;
2842      maxRhs[iRow]=COIN_DBL_MAX;
2843      whichInfiniteDown[iRow]=-2;
2844      whichInfiniteUp[iRow]=-2;
2845    }
2846  }
2847  const int * COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
2848  const CoinBigIndex * COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
2849  const double * COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
2850 
2851  double *  COIN_RESTRICT columnLower=abcLower_+maximumAbcNumberRows_;
2852  double *  COIN_RESTRICT columnUpper=abcUpper_+maximumAbcNumberRows_;
2853  while(numberChanged > numberCheck) {
2854    int numberToDo=numberToDoNext;
2855    numberToDoNext=0;
2856    int * COIN_RESTRICT temp=changedList;
2857    changedList=changedListNext;
2858    changedListNext=temp;
2859    CoinAbcMemset0(changed,numberRows_);
2860   
2861    numberChanged = 0; // Bounds tightened this pass
2862   
2863    if (iPass == MAXPASS) break;
2864    iPass++;
2865    for (int k=0;k<numberToDo;k++) {
2866      int iRow=changedList[k];
2867      // possible row
2868      int infiniteUpper = 0;
2869      int infiniteLower = 0;
2870      int firstInfiniteUpper=-1;
2871      int firstInfiniteLower=-1;
2872      double maximumUp = 0.0;
2873      double maximumDown = 0.0;
2874      double newBound;
2875      CoinBigIndex rStart = rowStart[iRow];
2876      CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
2877      CoinBigIndex j;
2878      // Compute possible lower and upper ranges
2879     
2880      for (j = rStart; j < rEnd; ++j) {
2881        double value = element[j];
2882        int iColumn = column[j];
2883        if (value > 0.0) {
2884          if (columnUpper[iColumn] >= large) {
2885            firstInfiniteUpper=iColumn;
2886            ++infiniteUpper;
2887          } else {
2888            maximumUp += columnUpper[iColumn] * value;
2889          }
2890          if (columnLower[iColumn] <= -large) {
2891            firstInfiniteLower=iColumn;
2892            ++infiniteLower;
2893          } else {
2894            maximumDown += columnLower[iColumn] * value;
2895          }
2896        } else if (value < 0.0) {
2897          if (columnUpper[iColumn] >= large) {
2898            firstInfiniteLower=iColumn;
2899            ++infiniteLower;
2900          } else {
2901            maximumDown += columnUpper[iColumn] * value;
2902          }
2903          if (columnLower[iColumn] <= -large) {
2904            firstInfiniteUpper=iColumn;
2905            ++infiniteUpper;
2906          } else {
2907            maximumUp += columnLower[iColumn] * value;
2908          }
2909        }
2910      }
2911      // Build in a margin of error
2912      maximumUp += 1.0e-8 * fabs(maximumUp);
2913      maximumDown -= 1.0e-8 * fabs(maximumDown);
2914      double maxUp = maximumUp + infiniteUpper * 1.0e31;
2915      double maxDown = maximumDown - infiniteLower * 1.0e31;
2916      minRhs[iRow]=infiniteLower ? -COIN_DBL_MAX : maximumDown;
2917      maxRhs[iRow]=infiniteUpper ? COIN_DBL_MAX : maximumUp;
2918      if (infiniteLower==0)
2919        whichInfiniteDown[iRow]=-1;
2920      else if (infiniteLower==1)
2921        whichInfiniteDown[iRow]=firstInfiniteLower;
2922      else
2923        whichInfiniteDown[iRow]=-2;
2924      if (infiniteUpper==0)
2925        whichInfiniteUp[iRow]=-1;
2926      else if (infiniteUpper==1)
2927        whichInfiniteUp[iRow]=firstInfiniteUpper;
2928      else
2929        whichInfiniteUp[iRow]=-2;
2930      if (maxUp <= rowUpper[iRow] + tolerance &&
2931          maxDown >= rowLower[iRow] - tolerance) {
2932       
2933        // Row is redundant - make totally free
2934        // NO - can't do this for postsolve
2935        // rowLower[iRow]=-COIN_DBL_MAX;
2936        // rowUpper[iRow]=COIN_DBL_MAX;
2937        //printf("Redundant row in presolveX %d\n",iRow);
2938       
2939      } else {
2940        if (maxUp < rowLower[iRow] - 100.0 * tolerance ||
2941            maxDown > rowUpper[iRow] + 100.0 * tolerance) {
2942          // problem is infeasible - exit at once
2943          numberInfeasible++;
2944          break;
2945        }
2946        double lower = rowLower[iRow];
2947        double upper = rowUpper[iRow];
2948        for (j = rStart; j < rEnd; ++j) {
2949          double value = element[j];
2950          int iColumn = column[j];
2951          double nowLower = columnLower[iColumn];
2952          double nowUpper = columnUpper[iColumn];
2953          if (value > 0.0) {
2954            // positive value
2955            if (lower > -large) {
2956              if (!infiniteUpper) {
2957                assert(nowUpper < large2);
2958                newBound = nowUpper +
2959                  (lower - maximumUp) / value;
2960                // relax if original was large
2961                if (fabs(maximumUp) > 1.0e8)
2962                  newBound -= 1.0e-12 * fabs(maximumUp);
2963              } else if (infiniteUpper == 1 && nowUpper > large) {
2964                newBound = (lower - maximumUp) / value;
2965                // relax if original was large
2966                if (fabs(maximumUp) > 1.0e8)
2967                  newBound -= 1.0e-12 * fabs(maximumUp);
2968              } else {
2969                newBound = -COIN_DBL_MAX;
2970              }
2971              if (newBound > nowLower + 1.0e-12 && newBound > -large) {
2972                // Tighten the lower bound
2973                numberChanged++;
2974                // check infeasible (relaxed)
2975                if (nowUpper < newBound) {
2976                  if (nowUpper - newBound <
2977                      -100.0 * tolerance)
2978                    numberInfeasible++;
2979                  else
2980                    newBound = nowUpper;
2981                }
2982                columnLower[iColumn] = newBound;
2983                for (CoinBigIndex j1=columnStart[iColumn];j1<columnStart[iColumn+1];j1++) {
2984                  int jRow=row[j1];
2985                  if (!changed[jRow]) {
2986                    changed[jRow]=1;
2987                    changedListNext[numberToDoNext++]=jRow;
2988                  }
2989                }
2990                // adjust
2991                double now;
2992                if (nowLower < -large) {
2993                  now = 0.0;
2994                  infiniteLower--;
2995                } else {
2996                  now = nowLower;
2997                }
2998                maximumDown += (newBound - now) * value;
2999                nowLower = newBound;
3000              }
3001            }
3002            if (upper < large) {
3003              if (!infiniteLower) {
3004                assert(nowLower > - large2);
3005                newBound = nowLower +
3006                  (upper - maximumDown) / value;
3007                // relax if original was large
3008                if (fabs(maximumDown) > 1.0e8)
3009                  newBound += 1.0e-12 * fabs(maximumDown);
3010              } else if (infiniteLower == 1 && nowLower < -large) {
3011                newBound =   (upper - maximumDown) / value;
3012                // relax if original was large
3013                if (fabs(maximumDown) > 1.0e8)
3014                  newBound += 1.0e-12 * fabs(maximumDown);
3015              } else {
3016                newBound = COIN_DBL_MAX;
3017              }
3018              if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3019                // Tighten the upper bound
3020                numberChanged++;
3021                // check infeasible (relaxed)
3022                if (nowLower > newBound) {
3023                  if (newBound - nowLower <
3024                      -100.0 * tolerance)
3025                    numberInfeasible++;
3026                  else
3027                    newBound = nowLower;
3028                }
3029                columnUpper[iColumn] = newBound;
3030                for (CoinBigIndex j1=columnStart[iColumn];j1<columnStart[iColumn+1];j1++) {
3031                  int jRow=row[j1];
3032                  if (!changed[jRow]) {
3033                    changed[jRow]=1;
3034                    changedListNext[numberToDoNext++]=jRow;
3035                  }
3036                }
3037                // adjust
3038                double now;
3039                if (nowUpper > large) {
3040                  now = 0.0;
3041                  infiniteUpper--;
3042                } else {
3043                  now = nowUpper;
3044                }
3045                maximumUp += (newBound - now) * value;
3046                nowUpper = newBound;
3047              }
3048            }
3049          } else {
3050            // negative value
3051            if (lower > -large) {
3052              if (!infiniteUpper) {
3053                assert(nowLower < large2);
3054                newBound = nowLower +
3055                  (lower - maximumUp) / value;
3056                // relax if original was large
3057                if (fabs(maximumUp) > 1.0e8)
3058                  newBound += 1.0e-12 * fabs(maximumUp);
3059              } else if (infiniteUpper == 1 && nowLower < -large) {
3060                newBound = (lower - maximumUp) / value;
3061                // relax if original was large
3062                if (fabs(maximumUp) > 1.0e8)
3063                  newBound += 1.0e-12 * fabs(maximumUp);
3064              } else {
3065                newBound = COIN_DBL_MAX;
3066              }
3067              if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3068                // Tighten the upper bound
3069                numberChanged++;
3070                // check infeasible (relaxed)
3071                if (nowLower > newBound) {
3072                  if (newBound - nowLower <
3073                      -100.0 * tolerance)
3074                    numberInfeasible++;
3075                  else
3076                    newBound = nowLower;
3077                }
3078                columnUpper[iColumn] = newBound;
3079                for (CoinBigIndex j1=columnStart[iColumn];j1<columnStart[iColumn+1];j1++) {
3080                  int jRow=row[j1];
3081                  if (!changed[jRow]) {
3082                    changed[jRow]=1;
3083                    changedListNext[numberToDoNext++]=jRow;
3084                  }
3085                }
3086                // adjust
3087                double now;
3088                if (nowUpper > large) {
3089                  now = 0.0;
3090                  infiniteLower--;
3091                } else {
3092                  now = nowUpper;
3093                }
3094                maximumDown += (newBound - now) * value;
3095                nowUpper = newBound;
3096              }
3097            }
3098            if (upper < large) {
3099              if (!infiniteLower) {
3100                assert(nowUpper < large2);
3101                newBound = nowUpper +
3102                  (upper - maximumDown) / value;
3103                // relax if original was large
3104                if (fabs(maximumDown) > 1.0e8)
3105                  newBound -= 1.0e-12 * fabs(maximumDown);
3106              } else if (infiniteLower == 1 && nowUpper > large) {
3107                newBound =   (upper - maximumDown) / value;
3108                // relax if original was large
3109                if (fabs(maximumDown) > 1.0e8)
3110                  newBound -= 1.0e-12 * fabs(maximumDown);
3111              } else {
3112                newBound = -COIN_DBL_MAX;
3113              }
3114              if (newBound > nowLower + 1.0e-12 && newBound > -large) {
3115                // Tighten the lower bound
3116                numberChanged++;
3117                // check infeasible (relaxed)
3118                if (nowUpper < newBound) {
3119                  if (nowUpper - newBound <
3120                      -100.0 * tolerance)
3121                    numberInfeasible++;
3122                  else
3123                    newBound = nowUpper;
3124                }
3125                columnLower[iColumn] = newBound;
3126                for (CoinBigIndex j1=columnStart[iColumn];j1<columnStart[iColumn+1];j1++) {
3127                  int jRow=row[j1];
3128                  if (!changed[jRow]) {
3129                    changed[jRow]=1;
3130                    changedListNext[numberToDoNext++]=jRow;
3131                  }
3132                }
3133                // adjust
3134                double now;
3135                if (nowLower < -large) {
3136                  now = 0.0;
3137                  infiniteUpper--;
3138                } else {
3139                  now = nowLower;
3140                }
3141                maximumUp += (newBound - now) * value;
3142                nowLower = newBound;
3143              }
3144            }
3145          }
3146        }
3147      }
3148    }
3149#if 0
3150    // see if we can tighten doubletons with infinite bounds
3151    if (iPass==1) {
3152      const double * COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
3153      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
3154        CoinBigIndex start = columnStart[iColumn];
3155        if (start+2==columnStart[iColumn+1]&&
3156          (columnLower[iColumn]<-1.0e30||columnUpper[iColumn])) {
3157          int iRow0=row[start];
3158          int iRow1=row[start+1];
3159          printf("col %d row0 %d el0 %g whichDown0 %d whichUp0 %d row1 %d el1 %g whichDown1 %d whichUp1 %d\n",iColumn,
3160                 iRow0,elementByColumn[start],whichInfiniteDown[iRow0],whichInfiniteUp[iRow0],
3161                 iRow1,elementByColumn[start+1],whichInfiniteDown[iRow1],whichInfiniteUp[iRow1]);
3162        }
3163      }
3164    }
3165#endif
3166    totalTightened += numberChanged;
3167    if (iPass == 1)
3168      numberCheck = numberChanged >> 4;
3169    if (numberInfeasible) break;
3170  }
3171  const double *  COIN_RESTRICT saveLower = abcLower_+maximumNumberTotal_+maximumAbcNumberRows_;
3172  const double *  COIN_RESTRICT saveUpper = abcUpper_+maximumNumberTotal_+maximumAbcNumberRows_;
3173  if (!numberInfeasible) {
3174    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_)
3175      << totalTightened
3176      << CoinMessageEol;
3177    int numberLowerChanged=0;
3178    int numberUpperChanged=0;
3179    if (!integerType_) {
3180    // Set bounds slightly loose
3181    // tighten rows as well
3182      //#define PRINT_TIGHTEN_PROGRESS 0
3183    for (int iRow = 0; iRow < numberRows_; iRow++) {
3184      assert (maxRhs[iRow]>=rowLower[iRow]);
3185      assert (minRhs[iRow]<=rowUpper[iRow]);
3186      rowLower[iRow]=CoinMax(rowLower[iRow],minRhs[iRow]);
3187      rowUpper[iRow]=CoinMin(rowUpper[iRow],maxRhs[iRow]);
3188#if PRINT_TIGHTEN_PROGRESS>3
3189      if (handler_->logLevel()>5)
3190        printf("Row %d min %g max %g lower %g upper %g\n",
3191               iRow,minRhs[iRow],maxRhs[iRow],rowLower[iRow],rowUpper[iRow]);
3192#endif
3193    }
3194#if 0
3195    double useTolerance = 1.0e-4;
3196    double multiplier = 100.0;
3197    // To do this we need to compute min and max JUST for no costs?
3198    const double * COIN_RESTRICT cost = abcCost_+maximumAbcNumberRows_;
3199    for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3200      if (saveUpper[iColumn] > saveLower[iColumn] + useTolerance) {
3201        double awayFromLower=0.0;
3202        double awayFromUpper=0.0;
3203        //double gap=columnUpper[iColumn]-columnLower[iColumn];
3204        for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3205          int iRow=row[j];
3206          double value=elementByColumn[j];
3207#if PRINT_TIGHTEN_PROGRESS>3
3208          if (handler_->logLevel()>5)
3209            printf("row %d element %g\n",iRow,value);
3210#endif
3211          if (value>0.0) {
3212            if (minRhs[iRow]+value*awayFromLower<rowLower[iRow]) {
3213              if (minRhs[iRow]>-1.0e10&&awayFromLower<1.0e10)
3214                awayFromLower = CoinMax(awayFromLower,(rowLower[iRow]-minRhs[iRow])/value);
3215              else
3216                awayFromLower=COIN_DBL_MAX;
3217            }
3218            if (maxRhs[iRow]-value*awayFromUpper>rowUpper[iRow]) {
3219              if (maxRhs[iRow]<1.0e10&&awayFromUpper<1.0e10)
3220                awayFromUpper = CoinMax(awayFromUpper,(maxRhs[iRow]-rowUpper[iRow])/value);
3221              else
3222                awayFromUpper=COIN_DBL_MAX;
3223            }
3224          } else {
3225            if (maxRhs[iRow]+value*awayFromLower>rowUpper[iRow]) {
3226              if (maxRhs[iRow]<1.0e10&&awayFromLower<1.0e10)
3227                awayFromLower = CoinMax(awayFromLower,(rowUpper[iRow]-maxRhs[iRow])/value);
3228              else
3229                awayFromLower=COIN_DBL_MAX;
3230            }
3231            if (minRhs[iRow]-value*awayFromUpper<rowLower[iRow]) {
3232              if (minRhs[iRow]>-1.0e10&&awayFromUpper<1.0e10)
3233                awayFromUpper = CoinMin(awayFromUpper,(minRhs[iRow]-rowLower[iRow])/value);
3234              else
3235                awayFromUpper=COIN_DBL_MAX;
3236            }
3237          }
3238        }
3239        // Might have to go as high as
3240        double upper = CoinMin(columnLower[iColumn]+awayFromLower,columnUpper[iColumn]);
3241        // and as low as
3242        double lower = CoinMax(columnUpper[iColumn]-awayFromUpper,columnLower[iColumn]);
3243        // but be sensible
3244        double gap=0.999*(CoinMin(columnUpper[iColumn]-columnLower[iColumn],1.0e10));
3245        if (awayFromLower>gap||awayFromUpper>gap) {
3246          if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3247#if PRINT_TIGHTEN_PROGRESS>1
3248            printf("YYchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3249                   columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3250#endif
3251            upper=columnUpper[iColumn];
3252          }
3253          if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3254#if PRINT_TIGHTEN_PROGRESS>1
3255            printf("YYchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3256                   columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3257#endif
3258            lower=columnLower[iColumn];
3259          }
3260        }
3261        if (lower>upper) {
3262#if PRINT_TIGHTEN_PROGRESS
3263          printf("XXchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3264                 columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3265          printf("XXchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3266                 columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3267#endif
3268          lower=columnLower[iColumn];
3269          upper=columnUpper[iColumn];
3270        }
3271        //upper=CoinMax(upper,0.0);
3272        //upper=CoinMax(upper,0.0);
3273        if (cost[iColumn]>=0.0&&awayFromLower<1.0e10&&columnLower[iColumn]>-1.0e10) {
3274          // doesn't want to be too positive
3275          if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3276#if PRINT_TIGHTEN_PROGRESS>2
3277            if (handler_->logLevel()>1)
3278              printf("change upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3279                     columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3280#endif
3281            double difference=columnUpper[iColumn]-upper;
3282            for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3283              int iRow=row[j];
3284              double value=elementByColumn[j];
3285              if (value>0.0) {
3286                assert (difference*value>=0.0);
3287                maxRhs[iRow]-=difference*value;
3288              } else {
3289                assert (difference*value<=0.0);
3290                minRhs[iRow]-=difference*value;
3291              }
3292            }
3293            columnUpper[iColumn]=upper;
3294            numberUpperChanged++;
3295          }
3296        }
3297        if (cost[iColumn]<=0.0&&awayFromUpper<1.0e10&&columnUpper[iColumn]<1.0e10) {
3298          // doesn't want to be too negative
3299          if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3300#if PRINT_TIGHTEN_PROGRESS>2
3301            if (handler_->logLevel()>1)
3302              printf("change lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3303                     columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3304#endif
3305            double difference=columnLower[iColumn]-lower;
3306            for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3307              int iRow=row[j];
3308              double value=elementByColumn[j];
3309              if (value>0.0) {
3310                assert (difference*value<=0.0);
3311                minRhs[iRow]-=difference*value;
3312              } else {
3313                assert (difference*value>=0.0);
3314                maxRhs[iRow]-=difference*value;
3315              }
3316            }
3317            columnLower[iColumn]=lower;
3318            numberLowerChanged++;
3319          }
3320        }
3321        // Make large bounds stay infinite
3322        if (saveUpper[iColumn] > 1.0e30 && columnUpper[iColumn] > 1.0e10) {
3323          columnUpper[iColumn] = COIN_DBL_MAX;
3324        }
3325        if (saveLower[iColumn] < -1.0e30 && columnLower[iColumn] < -1.0e10) {
3326          columnLower[iColumn] = -COIN_DBL_MAX;
3327        }
3328        if (columnUpper[iColumn] - columnLower[iColumn] < useTolerance + 1.0e-8) {
3329          // relax enough so will have correct dj
3330          columnLower[iColumn] = CoinMax(saveLower[iColumn],
3331                                          columnLower[iColumn] - multiplier * useTolerance);
3332          columnUpper[iColumn] = CoinMin(saveUpper[iColumn],
3333                                          columnUpper[iColumn] + multiplier * useTolerance);
3334        } else {
3335          if (columnUpper[iColumn] < saveUpper[iColumn]) {
3336            // relax a bit
3337            columnUpper[iColumn] = CoinMin(columnUpper[iColumn] + multiplier * useTolerance,
3338                                            saveUpper[iColumn]);
3339          }
3340          if (columnLower[iColumn] > saveLower[iColumn]) {
3341            // relax a bit
3342            columnLower[iColumn] = CoinMax(columnLower[iColumn] - multiplier * useTolerance,
3343                                            saveLower[iColumn]);
3344          }
3345        }
3346        if (0) {
3347          // debug
3348          // tighten rows as well
3349          for (int iRow = 0; iRow < numberRows_; iRow++) {
3350            if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3351              // possible row
3352              int infiniteUpper = 0;
3353              int infiniteLower = 0;
3354              double maximumUp = 0.0;
3355              double maximumDown = 0.0;
3356              CoinBigIndex rStart = rowStart[iRow];
3357              CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3358              CoinBigIndex j;
3359              // Compute possible lower and upper ranges
3360              for (j = rStart; j < rEnd; ++j) {
3361                double value = element[j];
3362                int iColumn = column[j];
3363                if (value > 0.0) {
3364                  if (columnUpper[iColumn] >= large) {
3365                    ++infiniteUpper;
3366                  } else {
3367                    maximumUp += columnUpper[iColumn] * value;
3368                  }
3369                  if (columnLower[iColumn] <= -large) {
3370                    ++infiniteLower;
3371                  } else {
3372                    maximumDown += columnLower[iColumn] * value;
3373                  }
3374                } else if (value < 0.0) {
3375                  if (columnUpper[iColumn] >= large) {
3376                    ++infiniteLower;
3377                  } else {
3378                    maximumDown += columnUpper[iColumn] * value;
3379                  }
3380                  if (columnLower[iColumn] <= -large) {
3381                    ++infiniteUpper;
3382                  } else {
3383                    maximumUp += columnLower[iColumn] * value;
3384                  }
3385                }
3386              }
3387              // Build in a margin of error
3388              double maxUp = maximumUp+1.0e-8 * fabs(maximumUp);
3389              double maxDown = maximumDown-1.0e-8 * fabs(maximumDown);
3390              double rLower=rowLower[iRow];
3391              double rUpper=rowUpper[iRow];
3392              if (!infiniteUpper&&maxUp < rUpper - tolerance) {
3393                if (rUpper!=COIN_DBL_MAX||maximumUp<1.0e5)
3394                  rUpper=maximumUp;
3395              }
3396              if (!infiniteLower&&maxDown > rLower + tolerance) {
3397                if (rLower!=-COIN_DBL_MAX||maximumDown>-1.0e5)
3398                  rLower=maximumDown;
3399              }
3400              if (infiniteUpper)
3401                maxUp=COIN_DBL_MAX;
3402              if (fabs(maxUp-maxRhs[iRow])>1.0e-3*(1.0+fabs(maxUp)))
3403                printf("bad Maxup row %d maxUp %g maxRhs %g\n",iRow,maxUp,maxRhs[iRow]);
3404              maxRhs[iRow]=maxUp;
3405              if (infiniteLower)
3406                maxDown=-COIN_DBL_MAX;
3407              if (fabs(maxDown-minRhs[iRow])>1.0e-3*(1.0+fabs(maxDown)))
3408                printf("bad Maxdown row %d maxDown %g minRhs %g\n",iRow,maxDown,minRhs[iRow]);
3409              minRhs[iRow]=maxDown;
3410              assert(rLower<=rUpper);
3411            }
3412          }
3413            // end debug
3414        }
3415      }
3416    }
3417    if (tightenType>1) {
3418      // relax
3419      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3420        // relax enough so will have correct dj
3421        double lower=saveLower[iColumn];
3422        double upper=saveUpper[iColumn];
3423        columnLower[iColumn] = CoinMax(lower,columnLower[iColumn] - multiplier);
3424        columnUpper[iColumn] = CoinMin(upper,columnUpper[iColumn] + multiplier);
3425      }
3426    }
3427#endif
3428    } else {
3429#if ABC_NORMAL_DEBUG>0
3430      printf("Redo tighten to relax by 1 for integers\n");
3431#endif
3432    }
3433#if ABC_NORMAL_DEBUG>0
3434    printf("Tighten - phase1 %d bounds changed, phase2 %d lower, %d upper\n",
3435           totalTightened,numberLowerChanged,numberUpperChanged);
3436#endif
3437    if (integerType_) {
3438      const double * columnScale=scaleToExternal_+maximumAbcNumberRows_;
3439      const double * inverseColumnScale=scaleFromExternal_+maximumAbcNumberRows_;
3440      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3441        if (integerType_[iColumn]) {
3442          double value;
3443          double lower = columnLower[iColumn]*inverseColumnScale[iColumn];
3444          double upper = columnUpper[iColumn]*inverseColumnScale[iColumn];
3445          value = floor(lower + 0.5);
3446          if (fabs(value - lower) > primalTolerance_)
3447            value = ceil(lower);
3448          columnLower_[iColumn] = value;
3449          columnLower[iColumn]=value*columnScale[iColumn];
3450          value = floor(upper + 0.5);
3451          if (fabs(value - upper) > primalTolerance_)
3452            value = floor(upper);
3453          columnUpper_[iColumn] = value;
3454          columnUpper[iColumn]=value*columnScale[iColumn];
3455          if (columnLower_[iColumn] > columnUpper_[iColumn])
3456            numberInfeasible++;
3457        }
3458      }
3459      if (numberInfeasible) {
3460        handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3461          << numberInfeasible
3462          << CoinMessageEol;
3463        // restore column bounds
3464        CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3465        CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3466      }
3467    }
3468  } else {
3469    handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3470      << numberInfeasible
3471      << CoinMessageEol;
3472    // restore column bounds
3473    CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3474    CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3475  }
3476  if (!numberInfeasible) {
3477    // tighten rows as well
3478    for (int iRow = 0; iRow < numberRows_; iRow++) {
3479      if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3480        // possible row
3481        int infiniteUpper = 0;
3482        int infiniteLower = 0;
3483        double maximumUp = 0.0;
3484        double maximumDown = 0.0;
3485        CoinBigIndex rStart = rowStart[iRow];
3486        CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3487        CoinBigIndex j;
3488        // Compute possible lower and upper ranges
3489        for (j = rStart; j < rEnd; ++j) {
3490          double value = element[j];
3491          int iColumn = column[j];
3492          if (value > 0.0) {
3493            if (columnUpper[iColumn] >= large) {
3494              ++infiniteUpper;
3495            } else {
3496              maximumUp += columnUpper[iColumn] * value;
3497            }
3498            if (columnLower[iColumn] <= -large) {
3499              ++infiniteLower;
3500            } else {
3501              maximumDown += columnLower[iColumn] * value;
3502            }
3503          } else if (value < 0.0) {
3504            if (columnUpper[iColumn] >= large) {
3505              ++infiniteLower;
3506            } else {
3507              maximumDown += columnUpper[iColumn] * value;
3508            }
3509            if (columnLower[iColumn] <= -large) {
3510              ++infiniteUpper;
3511            } else {
3512              maximumUp += columnLower[iColumn] * value;
3513            }
3514          }
3515        }
3516        // Build in a margin of error
3517        double maxUp = maximumUp+1.0e-8 * fabs(maximumUp);
3518        double maxDown = maximumDown-1.0e-8 * fabs(maximumDown);
3519        if (!infiniteUpper&&maxUp < rowUpper[iRow] - tolerance) {
3520          if (rowUpper[iRow]!=COIN_DBL_MAX||maximumUp<1.0e5) 
3521            rowUpper[iRow]=maximumUp;
3522        }
3523        if (!infiniteLower&&maxDown > rowLower[iRow] + tolerance) {
3524          if (rowLower[iRow]!=-COIN_DBL_MAX||maximumDown>-1.0e5) 
3525            rowLower[iRow]=maximumDown;
3526        }
3527        assert(rowLower[iRow]<=rowUpper[iRow]);
3528        if(rowUpper[iRow]-rowLower[iRow]<1.0e-12) {
3529          if (fabs(rowLower[iRow])>fabs(rowUpper[iRow]))
3530            rowLower[iRow]=rowUpper[iRow];
3531          else
3532            rowUpper[iRow]=rowLower[iRow];
3533        }
3534      }
3535    }
3536    // flip back
3537    for (int iRow=0;iRow<numberRows_;iRow++) {
3538      double value=-rowLower[iRow];
3539      rowLower[iRow]=-rowUpper[iRow];
3540      rowUpper[iRow]=value;
3541    }
3542    // move back - relaxed unless integer
3543    if (!integerType_) {
3544      double * COIN_RESTRICT lowerSaved = abcLower_+maximumNumberTotal_; 
3545      double * COIN_RESTRICT upperSaved = abcUpper_+maximumNumberTotal_; 
3546      double tolerance=2.0*primalTolerance_;
3547      for (int i=0;i<numberTotal_;i++) {
3548        double lower=abcLower_[i];
3549        double upper=abcUpper_[i];
3550        if (lower!=upper) {
3551          lower=CoinMax(lower-tolerance,lowerSaved[i]);
3552          upper=CoinMin(upper+tolerance,upperSaved[i]);
3553        }
3554        abcLower_[i]=lower;
3555        lowerSaved[i]=lower;
3556        abcUpper_[i]=upper;
3557        upperSaved[i]=upper;
3558      }
3559    } else {
3560      CoinAbcMemcpy(abcLower_+maximumNumberTotal_,abcLower_,numberTotal_);
3561      CoinAbcMemcpy(abcUpper_+maximumNumberTotal_,abcUpper_,numberTotal_);
3562    }
3563#if 0
3564    for (int i=0;i<numberTotal_;i++) {
3565      if (abcSolution_[i+maximumNumberTotal_]>abcUpper_[i+maximumNumberTotal_]+1.0e-5)
3566        printf ("above %d %g %g\n",i,
3567                abcSolution_[i+maximumNumberTotal_],abcUpper_[i+maximumNumberTotal_]);
3568      if (abcSolution_[i+maximumNumberTotal_]<abcLower_[i+maximumNumberTotal_]-1.0e-5)
3569        printf ("below %d %g %g\n",i,
3570                abcSolution_[i+maximumNumberTotal_],abcLower_[i+maximumNumberTotal_]);
3571    }
3572#endif
3573    permuteBasis();
3574    // move solution
3575    CoinAbcMemcpy(abcSolution_,solutionSaved_,numberTotal_);
3576  }
3577  CoinAbcMemset0(minRhs,numberRows_);
3578  CoinAbcMemset0(maxRhs,numberRows_);
3579  for (int i=0;i<5;i++)
3580    setAvailableArray(whichArray[i]);
3581  return (numberInfeasible);
3582}
3583// dual
3584#include "AbcSimplexDual.hpp"
3585// Returns 0 if dual can be skipped
3586int 
3587ClpSimplex::doAbcDual()
3588{
3589  if ((abcState_&CLP_ABC_WANTED)==0) {
3590    return 1; // not wanted
3591  } else {
3592    assert (abcSimplex_);
3593    int crashState=(abcState_>>8)&7;
3594    if (crashState&&crashState<4) {
3595      switch (crashState) {
3596      case 1:
3597        crash(1000.0,1);
3598        break;
3599      case 2:
3600        crash(abcSimplex_->dualBound(),0);
3601        break;
3602      case 3:
3603        crash(0.0,3);
3604        break;
3605      }
3606    }
3607    // this and abcSimplex_ are approximately same pointer
3608    if ((abcState_&CLP_ABC_FULL_DONE)==0) {
3609      abcSimplex_->gutsOfResize(numberRows_,numberColumns_);
3610      abcSimplex_->translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
3611      //abcSimplex_->permuteIn();
3612      int maximumPivotsAbc=CoinMin(abcSimplex_->factorization()->maximumPivots(),numberColumns_+200);
3613      abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
3614      if (numberColumns_)
3615        abcSimplex_->tightenPrimalBounds();
3616    }
3617    if ((abcSimplex_->stateOfProblem()&VALUES_PASS2)!=0) {
3618      if (solution_)
3619        crashState=6;
3620    }
3621    if (crashState&&crashState>3) {
3622      abcSimplex_->crash(crashState-2);
3623    }
3624    if (crashState&&crashState<4) {
3625      abcSimplex_->putStuffInBasis(1+2);
3626    }
3627    int returnCode = abcSimplex_->doAbcDual();
3628    //set to force crossover test returnCode=1;
3629    abcSimplex_->permuteOut(ROW_PRIMAL_OK|ROW_DUAL_OK|COLUMN_PRIMAL_OK|COLUMN_DUAL_OK|ALL_STATUS_OK);
3630    if (returnCode) {
3631      // fix as this model is all messed up e.g. scaling
3632      scalingFlag_ = abs(scalingFlag_);
3633      //delete [] rowScale_;
3634      //delete [] columnScale_;
3635      rowScale_ = NULL;
3636      columnScale_ = NULL;
3637#if 0
3638      delete [] savedRowScale_;
3639      delete [] savedColumnScale_;
3640      savedRowScale_ = NULL;
3641      savedColumnScale_ = NULL;
3642      inverseRowScale_ = NULL;
3643      inverseColumnScale_ = NULL;
3644#endif
3645      if (perturbation_==50)
3646        perturbation_=51;
3647      //perturbation_=100;
3648#if ABC_NORMAL_DEBUG>0
3649      printf("Bad exit with status of %d\n",problemStatus_);
3650#endif
3651      //problemStatus_=10;
3652      int status=problemStatus_;
3653      //problemStatus_=-1;
3654      if (status!=10) {
3655        dual(); // Do ClpSimplexDual
3656     } else {
3657        moreSpecialOptions_ |= 256;
3658        primal(1);
3659      }
3660    } else {
3661      // double check objective
3662      //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
3663      perturbation_=100;
3664      moreSpecialOptions_ |= 256;
3665      //primal(1);
3666    }
3667    return returnCode;
3668  }
3669}
3670#include "AbcSimplexPrimal.hpp"
3671// Do dual (return 1 if cleanup needed)
3672int 
3673AbcSimplex::doAbcDual()
3674{
3675  if (objective_) {
3676    objective_->setActivated(0);
3677  } else {
3678    // create dummy stuff
3679    assert (!numberColumns_);
3680    if (!numberRows_)
3681      problemStatus_ = 0; // say optimal
3682    return 0;
3683  }
3684  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
3685      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
3686      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
3687      additional data and have no destructor or (non-default) constructor.
3688     
3689      This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
3690      in primal and dual.
3691     
3692      As far as I can see this is perfectly safe.
3693  */
3694  algorithm_=-1;
3695  baseIteration_=numberIterations_;
3696  if (!abcMatrix_->gotRowCopy())
3697    abcMatrix_->createRowCopy();
3698#ifdef EARLY_FACTORIZE
3699  if (maximumIterations()>1000000&&maximumIterations()<1000999) {
3700    numberEarly_=maximumIterations()-1000000;
3701#if ABC_NORMAL_DEBUG>0
3702    printf("Setting numberEarly_ to %d\n",numberEarly_);
3703#endif
3704    numberEarly_+=numberEarly_<<16;
3705  }
3706#endif
3707  minimumThetaMovement_=1.0e-10;
3708  if (fabs(infeasibilityCost_-1.0e10)<1.1e9
3709      &&fabs(infeasibilityCost_-1.0e10)>1.0e5&&false) {
3710    minimumThetaMovement_=1.0e-3/fabs(infeasibilityCost_-1.0e10);
3711    printf("trying minimum theta movement of %g\n",minimumThetaMovement_);
3712    infeasibilityCost_=1.0e10;
3713  }
3714  int savePerturbation=perturbation_;
3715  static_cast<AbcSimplexDual *> (this)->dual();
3716  minimumThetaMovement_=1.0e-10;
3717  if ((specialOptions_&2048)!=0 && problemStatus_==10 && !numberPrimalInfeasibilities_
3718      && sumDualInfeasibilities_ < 1000.0* dualTolerance_)
3719    problemStatus_ = 0; // ignore
3720  onStopped(); // set secondary status if stopped
3721  if (problemStatus_==1&&numberFlagged_) {
3722    problemStatus_=10;
3723    static_cast<AbcSimplexPrimal *> (this)->unflag();
3724  }
3725  if (perturbation_<101) {
3726    //printf("Perturbed djs?\n");
3727    double * save=abcCost_;
3728    abcCost_=costSaved_;
3729    computeInternalObjectiveValue();
3730    abcCost_=save;
3731  }
3732  int totalNumberIterations=numberIterations_;
3733#ifndef TRY_ABC_GUS
3734  if (problemStatus_==10) {
3735    int savePerturbation=perturbation_;
3736    perturbation_=100;
3737    copyFromSaved(2);
3738    minimumThetaMovement_=1.0e-12;
3739    baseIteration_=numberIterations_;
3740    static_cast<AbcSimplexPrimal *> (this)->primal(0);
3741    totalNumberIterations+=numberIterations_-baseIteration_;
3742    perturbation_=savePerturbation;
3743  }
3744#else
3745  int iPass=0;
3746  while (problemStatus_==10&&minimumThetaMovement_>0.99999e-12) {
3747    iPass++;
3748    if (initialSumInfeasibilities_>=0.0) {
3749      if (savePerturbation>=100) {
3750        perturbation_=100;
3751      } else {
3752        if (savePerturbation==50)
3753          perturbation_=CoinMax(51,HEAVY_PERTURBATION-4-iPass); //smaller
3754        else
3755          perturbation_=CoinMax(51,savePerturbation-1-iPass); //smaller
3756      }
3757    } else {
3758      perturbation_=savePerturbation;
3759      abcFactorization_->setPivots(100000); // force factorization
3760      initialSumInfeasibilities_=1.23456789e-5;
3761      // clean pivots
3762      int numberBasic=0;
3763      int * pivot=pivotVariable();
3764      for (int i=0;i<numberTotal_;i++) {
3765        if (getInternalStatus(i)==basic)
3766          pivot[numberBasic++]=i;
3767      }
3768      assert (numberBasic==numberRows_);
3769    }
3770    copyFromSaved(14);
3771    // Say second call
3772    if (iPass>3)
3773      moreSpecialOptions_ |= 256;
3774    if (iPass>2)
3775      perturbation_=101;
3776    baseIteration_=numberIterations_;
3777    static_cast<AbcSimplexPrimal *> (this)->primal(0);
3778    totalNumberIterations+=numberIterations_-baseIteration_;
3779    if (problemStatus_==10) {
3780      // reduce
3781      minimumThetaMovement_*=1.0e-1;
3782      copyFromSaved(14);
3783      baseIteration_=numberIterations_;
3784      static_cast<AbcSimplexDual *> (this)->dual();
3785      totalNumberIterations+=numberIterations_-baseIteration_;
3786    }
3787    perturbation_=savePerturbation;
3788  }
3789#endif
3790  numberIterations_=totalNumberIterations;
3791  return (problemStatus_<0||problemStatus_>3) ? 1 : 0;
3792}
3793int AbcSimplex::dual ()
3794{
3795  if (!abcMatrix_->gotRowCopy())
3796    abcMatrix_->createRowCopy();
3797  assert (!numberFlagged_);
3798  baseIteration_=numberIterations_;
3799  //double savedPivotTolerance = abcFactorization_->pivotTolerance();
3800  if (objective_) {
3801    objective_->setActivated(0);
3802  } else {
3803    // create dummy stuff
3804    assert (!numberColumns_);
3805    if (!numberRows_)
3806      problemStatus_ = 0; // say optimal
3807    return 0;
3808  }
3809  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
3810      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
3811      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
3812      additional data and have no destructor or (non-default) constructor.
3813     
3814      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3815      in primal and dual.
3816     
3817      As far as I can see this is perfectly safe.
3818  */
3819  minimumThetaMovement_=1.0e-12;
3820  int returnCode = static_cast<AbcSimplexDual *> (this)->dual();
3821  if ((specialOptions_ & 2048) != 0 && problemStatus_ == 10 && !numberPrimalInfeasibilities_
3822      && sumDualInfeasibilities_ < 1000.0 * dualTolerance_ && perturbation_ >= 100)
3823    problemStatus_ = 0; // ignore
3824  if (problemStatus_ == 10) {
3825#if 1 //ABC_NORMAL_DEBUG>0
3826    printf("return and use ClpSimplexPrimal\n");
3827    abort();
3828#else
3829    //printf("Cleaning up with primal\n");
3830    //lastAlgorithm=1;
3831    int savePerturbation = perturbation_;
3832    int saveLog = handler_->logLevel();
3833    //handler_->setLogLevel(63);
3834    perturbation_ = 101;
3835    bool denseFactorization = initialDenseFactorization();
3836    // It will be safe to allow dense
3837    setInitialDenseFactorization(true);
3838    // Allow for catastrophe
3839    int saveMax = intParam_[ClpMaxNumIteration];
3840    if (numberIterations_) {
3841      // normal
3842      if (intParam_[ClpMaxNumIteration] > 100000 + numberIterations_)
3843        intParam_[ClpMaxNumIteration]
3844          = numberIterations_ + 1000 + 2 * numberRows_ + numberColumns_;
3845    } else {
3846      // Not normal allow more
3847      baseIteration_ += 2 * (numberRows_ + numberColumns_);
3848    }
3849    // check which algorithms allowed
3850    int dummy;
3851    if (problemStatus_ == 10 && saveObjective == objective_)
3852      startFinishOptions |= 2;
3853    baseIteration_ = numberIterations_;
3854    // Say second call
3855    moreSpecialOptions_ |= 256;
3856    minimumThetaMovement_=1.0e-12;
3857    returnCode = static_cast<AbcSimplexPrimal *> (this)->primal(1, startFinishOptions);
3858    // Say not second call
3859    moreSpecialOptions_ &= ~256;
3860    baseIteration_ = 0;
3861    if (saveObjective != objective_) {
3862      // We changed objective to see if infeasible
3863      delete objective_;
3864      objective_ = saveObjective;
3865      if (!problemStatus_) {
3866        // carry on
3867        returnCode = static_cast<AbcSimplexPrimal *> (this)->primal(1, startFinishOptions);
3868      }
3869    }
3870    if (problemStatus_ == 3 && numberIterations_ < saveMax) {
3871      // flatten solution and try again
3872      for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
3873        if (getInternalStatus(iSequence) != basic) {
3874          setInternalStatus(iSequence, superBasic);
3875          // but put to bound if close
3876          if (fabs(rowActivity_[iSequence] - rowLower_[iSequence])
3877              <= primalTolerance_) {
3878            abcSolution_[iSequence] = abcLower_[iSequence];
3879            setInternalStatus(iSequence, atLowerBound);
3880          } else if (fabs(rowActivity_[iSequence] - rowUpper_[iSequence])
3881                     <= primalTolerance_) {
3882            abcSolution_[iSequence] = abcUpper_[iSequence];
3883            setInternalStatus(iSequence, atUpperBound);
3884          }
3885        }
3886      }
3887      problemStatus_ = -1;
3888      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 +
3889                                              2 * numberRows_ + numberColumns_, saveMax);
3890      perturbation_ = savePerturbation;
3891      baseIteration_ = numberIterations_;
3892      // Say second call
3893      moreSpecialOptions_ |= 256;
3894      returnCode = static_cast<AbcSimplexPrimal *> (this)->primal(0, startFinishOptions);
3895      // Say not second call
3896      moreSpecialOptions_ &= ~256;
3897      baseIteration_ = 0;
3898      computeObjectiveValue();
3899      // can't rely on djs either
3900      memset(reducedCost_, 0, numberColumns_ * sizeof(double));
3901    }
3902    intParam_[ClpMaxNumIteration] = saveMax;
3903   
3904    setInitialDenseFactorization(denseFactorization);
3905    perturbation_ = savePerturbation;
3906    if (problemStatus_ == 10) {
3907      if (!numberPrimalInfeasibilities_)
3908        problemStatus_ = 0;
3909      else
3910        problemStatus_ = 4;
3911    }
3912    handler_->setLogLevel(saveLog);
3913#endif
3914  }
3915  onStopped(); // set secondary status if stopped
3916  return returnCode;
3917}
3918// primal
3919// Returns 0 if primal can be skipped
3920int 
3921ClpSimplex::doAbcPrimal(int ifValuesPass)
3922{
3923  if ((abcState_&CLP_ABC_WANTED)==0) {
3924    return 1; // not wanted
3925  } else {
3926    assert (abcSimplex_);
3927    if (ifValuesPass) 
3928      abcSimplex_->setStateOfProblem(abcSimplex_->stateOfProblem()|VALUES_PASS);
3929    int crashState=(abcState_>>8)&7;
3930    if (crashState&&crashState<4) {
3931      switch (crashState) {
3932      case 1:
3933        crash(1000.0,1);
3934        break;
3935      case 2:
3936        crash(abcSimplex_->dualBound(),0);
3937        break;
3938      case 3:
3939        crash(0.0,3);
3940        break;
3941      }
3942    }
3943    // this and abcSimplex_ are approximately same pointer
3944    if ((abcState_&CLP_ABC_FULL_DONE)==0) {
3945      abcSimplex_->gutsOfResize(numberRows_,numberColumns_);
3946      abcSimplex_->translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
3947      //abcSimplex_->permuteIn();
3948      int maximumPivotsAbc=CoinMin(abcSimplex_->factorization()->maximumPivots(),numberColumns_+200);
3949      abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
3950      abcSimplex_->copyFromSaved(1);
3951    }
3952    if ((abcSimplex_->stateOfProblem()&VALUES_PASS2)!=0) {
3953      if (solution_)
3954        crashState=6;
3955    }
3956    if (crashState&&crashState>3) {
3957      abcSimplex_->crash(crashState-2);
3958    }
3959    if (crashState&&crashState<4) {
3960      abcSimplex_->putStuffInBasis(1+2);
3961    }
3962    int returnCode = abcSimplex_->doAbcPrimal(ifValuesPass);
3963    //set to force crossover test returnCode=1;
3964    abcSimplex_->permuteOut(ROW_PRIMAL_OK|ROW_DUAL_OK|COLUMN_PRIMAL_OK|COLUMN_DUAL_OK|ALL_STATUS_OK);
3965    if (returnCode) {
3966      // fix as this model is all messed up e.g. scaling
3967     scalingFlag_ = abs(scalingFlag_);
3968     //delete [] rowScale_;
3969     //delete [] columnScale_;
3970     rowScale_ = NULL;
3971     columnScale_ = NULL;
3972#if 0
3973     delete [] savedRowScale_;
3974     delete [] savedColumnScale_;
3975     savedRowScale_ = NULL;
3976     savedColumnScale_ = NULL;
3977     inverseRowScale_ = NULL;
3978     inverseColumnScale_ = NULL;
3979#endif
3980     if (perturbation_==50)
3981       perturbation_=51;
3982     //perturbation_=100;
3983#if ABC_NORMAL_DEBUG>0
3984     if (problemStatus_)
3985       printf("Bad exit with status of %d\n",problemStatus_);
3986#endif
3987     //problemStatus_=10;
3988     int status=problemStatus_;
3989     // should not be using
3990     for (int i = 0; i < 6; i++) {
3991       if (rowArray_[i])
3992         rowArray_[i]->clear();
3993       if (columnArray_[i])
3994         columnArray_[i]->clear();
3995     }
3996     //problemStatus_=-1;
3997     if (status!=3&&status!=0) {
3998       if (status!=10) {
3999         primal(); // Do ClpSimplexPrimal
4000       } else {
4001         moreSpecialOptions_ |= 256;
4002         dual();
4003       }
4004     }
4005    } else {
4006      // double check objective
4007      //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
4008    }
4009    numberIterations_=abcSimplex_->numberIterations();
4010    return returnCode;
4011  }
4012}
4013// Do primal (return 1 if cleanup needed)
4014int 
4015AbcSimplex::doAbcPrimal(int ifValuesPass)
4016{
4017  if (objective_) {
4018    objective_->setActivated(0);
4019  } else {
4020    // create dummy stuff
4021    assert (!numberColumns_);
4022    if (!numberRows_)
4023      problemStatus_ = 0; // say optimal
4024    return 0;
4025  }
4026  /*  Note use of "down casting".  The only class the user sees is AbcSimplex.
4027      Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
4028      and AbcSimplexOther all exist and inherit from AbcSimplex but have no
4029      additional data and have no destructor or (non-default) constructor.
4030     
4031      This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
4032      in primal and dual.
4033     
4034      As far as I can see this is perfectly safe.
4035  */
4036  algorithm_=-1;
4037  int savePerturbation=perturbation_;
4038  baseIteration_=numberIterations_;
4039  if (!abcMatrix_->gotRowCopy())
4040    abcMatrix_->createRowCopy();
4041#ifdef EARLY_FACTORIZE
4042  if (maximumIterations()>1000000&&maximumIterations()<1000999) {
4043    numberEarly_=maximumIterations()-1000000;
4044#if ABC_NORMAL_DEBUG>0
4045    printf("Setting numberEarly_ to %d\n",numberEarly_);
4046#endif
4047    numberEarly_+=numberEarly_<<16;
4048  }
4049#endif
4050  // add values pass options
4051  minimumThetaMovement_=1.0e-12;
4052  if ((specialOptions_&8192)!=0)
4053    minimumThetaMovement_=1.0e-15;
4054  int returnCode=static_cast<AbcSimplexPrimal *> (this)->primal(ifValuesPass);
4055  stateOfProblem_ &= ~VALUES_PASS;
4056#ifndef TRY_ABC_GUS
4057  if (returnCode==10) {
4058    copyFromSaved(14);
4059    int savePerturbation=perturbation_;
4060    perturbation_=51;
4061    minimumThetaMovement_=1.0e-12;
4062    returnCode=static_cast<AbcSimplexDual *> (this)->dual();
4063    perturbation_=savePerturbation;
4064  }
4065#else
4066  minimumThetaMovement_=1.0e-12;
4067  int iPass=0;
4068  while (problemStatus_==10&&minimumThetaMovement_>1.0e-15) {
4069    iPass++;
4070    if (minimumThetaMovement_==1.0e-12)
4071      perturbation_=CoinMin(savePerturbation,55);
4072    else
4073      perturbation_=100;
4074    copyFromSaved(14);
4075    // reduce ?
4076    // Say second call
4077    // Say second call
4078    if (iPass>3)
4079      moreSpecialOptions_ |= 256;
4080    if (iPass>2)
4081      perturbation_=101;
4082    baseIteration_=numberIterations_;
4083    if ((specialOptions_&8192)==0)
4084      static_cast<AbcSimplexDual *> (this)->dual();
4085    baseIteration_=numberIterations_;
4086    if (problemStatus_==10) {
4087      if (initialSumInfeasibilities_>=0.0) {
4088        if (savePerturbation>=100) {
4089          perturbation_=100;
4090        } else {
4091          if (savePerturbation==50)
4092            perturbation_=CoinMax(51,HEAVY_PERTURBATION-4-iPass); //smaller
4093          else
4094            perturbation_=CoinMax(51,savePerturbation-1-iPass); //smaller
4095        }
4096      } else {
4097        specialOptions_ |= 8192; // stop going to dual
4098        perturbation_=savePerturbation;
4099        abcFactorization_->setPivots(100000); // force factorization
4100        initialSumInfeasibilities_=1.23456789e-5;
4101        // clean pivots
4102        int numberBasic=0;
4103        int * pivot=pivotVariable();
4104        for (int i=0;i<numberTotal_;i++) {
4105          if (getInternalStatus(i)==basic)
4106            pivot[numberBasic++]=i;
4107        }
4108        assert (numberBasic==numberRows_);
4109      }
4110      if (iPass>2)
4111        perturbation_=100;
4112      copyFromSaved(14);
4113      baseIteration_=numberIterations_;
4114      static_cast<AbcSimplexPrimal *> (this)->primal(0);
4115    }
4116    minimumThetaMovement_*=0.5;
4117    perturbation_=savePerturbation;
4118  }
4119#endif
4120  return returnCode;
4121}
4122/* Sets up all slack basis and resets solution to
4123   as it was after initial load or readMps */
4124void AbcSimplex::allSlackBasis()
4125{
4126  // set column status to one nearest zero
4127  CoinFillN(internalStatus_,numberRows_,static_cast<unsigned char>(basic));
4128  const double *  COIN_RESTRICT lower = abcLower_;
4129  const double *  COIN_RESTRICT upper = abcUpper_;
4130  double *  COIN_RESTRICT solution = abcSolution_;
4131  // But set value to zero if lb <0.0 and ub>0.0
4132  for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
4133    if (lower[i] >= 0.0) {
4134      solution[i] = lower[i];
4135      setInternalStatus(i, atLowerBound);
4136    } else if (upper[i] <= 0.0) {
4137      solution[i] = upper[i];
4138      setInternalStatus(i, atUpperBound);
4139    } else if (lower[i] < -1.0e20 && upper[i] > 1.0e20) {
4140      // free
4141      solution[i] = 0.0;
4142      setInternalStatus(i, isFree);
4143    } else if (fabs(lower[i]) < fabs(upper[i])) {
4144      solution[i] = 0.0;
4145      setInternalStatus(i, atLowerBound);
4146    } else {
4147      solution[i] = 0.0;
4148      setInternalStatus(i, atUpperBound);
4149    }
4150  }
4151}
4152#if 0 //ndef SLIM_NOIO
4153// Read an mps file from the given filename
4154int
4155AbcSimplex::readMps(const char *filename,
4156                    bool keepNames,
4157                    bool ignoreErrors)
4158{
4159  int status = ClpSimplex::readMps(filename, keepNames, ignoreErrors);
4160  ClpSimplex * thisModel=this;
4161  gutsOfResize(thisModel->numberRows(),thisModel->numberColumns());
4162  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4163  return status;
4164}
4165// Read GMPL files from the given filenames
4166int
4167AbcSimplex::readGMPL(const char *filename, const char * dataName,
4168                     bool keepNames)
4169{
4170  int status = ClpSimplex::readGMPL(filename, dataName, keepNames);
4171  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4172  return status;
4173}
4174// Read file in LP format from file with name filename.
4175int
4176AbcSimplex::readLp(const char *filename, const double epsilon )
4177{
4178  FILE *fp = fopen(filename, "r");
4179 
4180  if(!fp) {
4181    printf("### ERROR: AbcSimplex::readLp():  Unable to open file %s for reading\n",
4182           filename);
4183    return(1);
4184  }
4185  CoinLpIO m;
4186  m.readLp(fp, epsilon);
4187  fclose(fp);
4188 
4189  // set problem name
4190  setStrParam(ClpProbName, m.getProblemName());
4191  // no errors
4192  loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
4193              m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
4194 
4195  if (m.integerColumns()) {
4196    integerType_ = new char[numberColumns_];
4197    CoinMemcpyN(m.integerColumns(), numberColumns_, integerType_);
4198  } else {
4199    integerType_ = NULL;
4200  }
4201  translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4202  unsigned int maxLength = 0;
4203  int iRow;
4204  rowNames_ = std::vector<std::string> ();
4205  columnNames_ = std::vector<std::string> ();
4206  rowNames_.reserve(numberRows_);
4207  for (iRow = 0; iRow < numberRows_; iRow++) {
4208    const char * name = m.rowName(iRow);
4209    if (name) {
4210      maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4211      rowNames_.push_back(name);
4212    } else {
4213      rowNames_.push_back("");
4214    }
4215  }
4216 
4217  int iColumn;
4218  columnNames_.reserve(numberColumns_);
4219  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4220    const char * name = m.columnName(iColumn);
4221    if (name) {
4222      maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4223      columnNames_.push_back(name);
4224    } else {
4225      columnNames_.push_back("");
4226    }
4227  }
4228  lengthNames_ = static_cast<int> (maxLength);
4229 
4230  return 0;
4231}
4232#endif
4233// Factorization frequency
4234int
4235AbcSimplex::factorizationFrequency() const
4236{
4237  if (abcFactorization_)
4238    return abcFactorization_->maximumPivots();
4239  else
4240    return -1;
4241}
4242void
4243AbcSimplex::setFactorizationFrequency(int value)
4244{
4245  if (abcFactorization_)
4246    abcFactorization_->maximumPivots(value);
4247}
4248
4249/* Get a clean factorization - i.e. throw out singularities
4250   may do more later */
4251int
4252AbcSimplex::cleanFactorization(int ifValuesPass)
4253{
4254  int status = internalFactorize(ifValuesPass ? 10 : 0);
4255  if (status < 0) {
4256    return 1; // some error
4257  } else {
4258    firstFree_=0;
4259    return 0;
4260  }
4261}
4262// Save data
4263ClpDataSave
4264AbcSimplex::saveData()
4265{
4266  ClpDataSave saved;
4267  saved.dualBound_ = dualBound_;
4268  saved.infeasibilityCost_ = infeasibilityCost_;
4269  saved.pivotTolerance_ = abcFactorization_->pivotTolerance();
4270  saved.zeroFactorizationTolerance_ = abcFactorization_->zeroTolerance();
4271  saved.zeroSimplexTolerance_ = zeroTolerance_;
4272  saved.perturbation_ = perturbation_;
4273  saved.forceFactorization_ = forceFactorization_;
4274  saved.acceptablePivot_ = acceptablePivot_;
4275  saved.objectiveScale_ = objectiveScale_;
4276  // Progress indicator
4277  abcProgress_.fillFromModel (this);
4278  return saved;
4279}
4280// Restore data
4281void
4282AbcSimplex::restoreData(ClpDataSave saved)
4283{
4284  //abcFactorization_->sparseThreshold(saved.sparseThreshold_);
4285  abcFactorization_->pivotTolerance(saved.pivotTolerance_);
4286  abcFactorization_->zeroTolerance(saved.zeroFactorizationTolerance_);
4287  zeroTolerance_ = saved.zeroSimplexTolerance_;
4288  perturbation_ = saved.perturbation_;
4289  infeasibilityCost_ = saved.infeasibilityCost_;
4290  dualBound_ = saved.dualBound_;
4291  forceFactorization_ = saved.forceFactorization_;
4292  objectiveScale_ = saved.objectiveScale_;
4293  acceptablePivot_ = saved.acceptablePivot_;
4294}
4295// To flag a variable
4296void
4297AbcSimplex::setFlagged( int sequence)
4298{
4299  assert (sequence>=0&&sequence<numberTotal_);
4300  if (sequence>=0) {
4301    internalStatus_[sequence] |= 64;
4302    // use for real disaster lastFlaggedIteration_ = numberIterations_;
4303    numberFlagged_++;
4304  }
4305}
4306#ifndef NDEBUG
4307// For errors to make sure print to screen
4308// only called in debug mode
4309static void indexError(int index,
4310                       std::string methodName)
4311{
4312  std::cerr << "Illegal index " << index << " in AbcSimplex::" << methodName << std::endl;
4313  throw CoinError("Illegal index", methodName, "AbcSimplex");
4314}
4315#endif
4316// After modifying first copy refreshes second copy and marks as updated
4317void 
4318AbcSimplex::refreshCosts()
4319{
4320  whatsChanged_ &= ~OBJECTIVE_SAME;
4321  CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);
4322}
4323void 
4324AbcSimplex::refreshLower(unsigned int type)
4325{
4326  if (type)
4327    whatsChanged_ &= type;
4328  CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_);
4329}
4330void 
4331AbcSimplex::refreshUpper(unsigned int type)
4332{
4333  if (type)
4334    whatsChanged_ &= type;
4335  CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_);
4336}
4337/* Set an objective function coefficient */
4338void
4339AbcSimplex::setObjectiveCoefficient( int elementIndex, double elementValue )
4340{
4341#ifndef NDEBUG
4342  if (elementIndex < 0 || elementIndex >= numberColumns_) {
4343    indexError(elementIndex, "setObjectiveCoefficient");
4344  }
4345#endif
4346  if (objective()[elementIndex] != elementValue) {
4347    objective()[elementIndex] = elementValue;
4348    // work arrays exist - update as well
4349    whatsChanged_ &= ~OBJECTIVE_SAME;
4350    double direction = optimizationDirection_ * objectiveScale_;
4351    costSaved_[elementIndex+maximumAbcNumberRows_] = direction * elementValue
4352      * columnUseScale_[elementIndex+maximumAbcNumberRows_];
4353  }
4354}
4355/* Set a single row lower bound<br>
4356   Use -DBL_MAX for -infinity. */
4357void
4358AbcSimplex::setRowLower( int elementIndex, double elementValue )
4359{
4360#ifndef NDEBUG
4361  int n = numberRows_;
4362  if (elementIndex < 0 || elementIndex >= n) {
4363    indexError(elementIndex, "setRowLower");
4364  }
4365#endif
4366  if (elementValue < -1.0e27)
4367    elementValue = -COIN_DBL_MAX;
4368  if (rowLower_[elementIndex] != elementValue) {
4369    rowLower_[elementIndex] = elementValue;
4370    // work arrays exist - update as well
4371    whatsChanged_ &= ~ROW_LOWER_SAME;
4372    if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4373      lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4374    } else {
4375      lowerSaved_[elementIndex] = elementValue * rhsScale_
4376        * rowUseScale_[elementIndex];
4377    }
4378  }
4379}
4380
4381/* Set a single row upper bound<br>
4382   Use DBL_MAX for infinity. */
4383void
4384AbcSimplex::setRowUpper( int elementIndex, double elementValue )
4385{
4386#ifndef NDEBUG
4387  int n = numberRows_;
4388  if (elementIndex < 0 || elementIndex >= n) {
4389    indexError(elementIndex, "setRowUpper");
4390  }
4391#endif
4392  if (elementValue > 1.0e27)
4393    elementValue = COIN_DBL_MAX;
4394  if (rowUpper_[elementIndex] != elementValue) {
4395    rowUpper_[elementIndex] = elementValue;
4396    // work arrays exist - update as well
4397    whatsChanged_ &= ~ROW_UPPER_SAME;
4398    if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4399      upperSaved_[elementIndex] = COIN_DBL_MAX;
4400    } else {
4401      upperSaved_[elementIndex] = elementValue * rhsScale_
4402        * rowUseScale_[elementIndex];
4403    }
4404  }
4405}
4406
4407/* Set a single row lower and upper bound */
4408void
4409AbcSimplex::setRowBounds( int elementIndex,
4410                          double lowerValue, double upperValue )
4411{
4412#ifndef NDEBUG
4413  int n = numberRows_;
4414  if (elementIndex < 0 || elementIndex >= n) {
4415    indexError(elementIndex, "setRowBounds");
4416  }
4417#endif
4418  if (lowerValue < -1.0e27)
4419    lowerValue = -COIN_DBL_MAX;
4420  if (upperValue > 1.0e27)
4421    upperValue = COIN_DBL_MAX;
4422  //CoinAssert (upperValue>=lowerValue);
4423  if (rowLower_[elementIndex] != lowerValue) {
4424    rowLower_[elementIndex] = lowerValue;
4425    // work arrays exist - update as well
4426    whatsChanged_ &= ~ROW_LOWER_SAME;
4427    if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4428      lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4429    } else {
4430      lowerSaved_[elementIndex] = lowerValue * rhsScale_
4431        * rowUseScale_[elementIndex];
4432    }
4433  }
4434  if (rowUpper_[elementIndex] != upperValue) {
4435    rowUpper_[elementIndex] = upperValue;
4436    // work arrays exist - update as well
4437    whatsChanged_ &= ~ROW_UPPER_SAME;
4438    if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4439      upperSaved_[elementIndex] = COIN_DBL_MAX;
4440    } else {
4441      upperSaved_[elementIndex] = upperValue * rhsScale_
4442        * rowUseScale_[elementIndex];
4443    }
4444  }
4445}
4446void AbcSimplex::setRowSetBounds(const int* indexFirst,
4447                                 const int* indexLast,
4448                                 const double* boundList)
4449{
4450#ifndef NDEBUG
4451  int n = numberRows_;
4452#endif
4453  int numberChanged = 0;
4454  const int * saveFirst = indexFirst;
4455  while (indexFirst != indexLast) {
4456    const int iRow = *indexFirst++;
4457#ifndef NDEBUG
4458    if (iRow < 0 || iRow >= n) {
4459      indexError(iRow, "setRowSetBounds");
4460    }
4461#endif
4462    double lowerValue = *boundList++;
4463    double upperValue = *boundList++;
4464    if (lowerValue < -1.0e27)
4465      lowerValue = -COIN_DBL_MAX;
4466    if (upperValue > 1.0e27)
4467      upperValue = COIN_DBL_MAX;
4468    //CoinAssert (upperValue>=lowerValue);
4469    if (rowLower_[iRow] != lowerValue) {
4470      rowLower_[iRow] = lowerValue;
4471      whatsChanged_ &= ~ROW_LOWER_SAME;
4472      numberChanged++;
4473    }
4474    if (rowUpper_[iRow] != upperValue) {
4475      rowUpper_[iRow] = upperValue;
4476      whatsChanged_ &= ~ROW_UPPER_SAME;
4477      numberChanged++;
4478    }
4479  }
4480  if (numberChanged) {
4481    indexFirst = saveFirst;
4482    while (indexFirst != indexLast) {
4483      const int iRow = *indexFirst++;
4484      if (rowLower_[iRow] == -COIN_DBL_MAX) {
4485        lowerSaved_[iRow] = -COIN_DBL_MAX;
4486      } else {
4487        lowerSaved_[iRow] = rowLower_[iRow] * rhsScale_
4488          * rowUseScale_[iRow];
4489      }
4490      if (rowUpper_[iRow] == COIN_DBL_MAX) {
4491        upperSaved_[iRow] = COIN_DBL_MAX;
4492      } else {
4493        upperSaved_[iRow] = rowUpper_[iRow] * rhsScale_
4494          * rowUseScale_[iRow];
4495      }
4496    }
4497  }
4498}
4499//-----------------------------------------------------------------------------
4500/* Set a single column lower bound<br>
4501   Use -DBL_MAX for -infinity. */
4502void
4503AbcSimplex::setColumnLower( int elementIndex, double elementValue )
4504{
4505#ifndef NDEBUG
4506  int n = numberColumns_;
4507  if (elementIndex < 0 || elementIndex >= n) {
4508    indexError(elementIndex, "setColumnLower");
4509  }
4510#endif
4511  if (elementValue < -1.0e27)
4512    elementValue = -COIN_DBL_MAX;
4513  if (columnLower_[elementIndex] != elementValue) {
4514    columnLower_[elementIndex] = elementValue;
4515    // work arrays exist - update as well
4516    whatsChanged_ &= ~COLUMN_LOWER_SAME;
4517    double value;
4518    if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4519      value = -COIN_DBL_MAX;
4520    } else {
4521      value = elementValue * rhsScale_
4522        * inverseColumnUseScale_[elementIndex];
4523    }
4524    lowerSaved_[elementIndex+maximumAbcNumberRows_] = value;
4525  }
4526}
4527
4528/* Set a single column upper bound<br>
4529   Use DBL_MAX for infinity. */
4530void
4531AbcSimplex::setColumnUpper( int elementIndex, double elementValue )
4532{
4533#ifndef NDEBUG
4534  int n = numberColumns_;
4535  if (elementIndex < 0 || elementIndex >= n) {
4536    indexError(elementIndex, "setColumnUpper");
4537  }
4538#endif
4539  if (elementValue > 1.0e27)
4540    elementValue = COIN_DBL_MAX;
4541  if (columnUpper_[elementIndex] != elementValue) {
4542    columnUpper_[elementIndex] = elementValue;
4543    // work arrays exist - update as well
4544    whatsChanged_ &= ~COLUMN_UPPER_SAME;
4545    double value;
4546    if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4547      value = COIN_DBL_MAX;
4548    } else {
4549      value = elementValue * rhsScale_
4550        * inverseColumnUseScale_[elementIndex];
4551    }
4552    upperSaved_[elementIndex+maximumAbcNumberRows_] = value;
4553  }
4554}
4555
4556/* Set a single column lower and upper bound */
4557void
4558AbcSimplex::setColumnBounds( int elementIndex,
4559                             double lowerValue, double upperValue )
4560{
4561#ifndef NDEBUG
4562  int n = numberColumns_;
4563  if (elementIndex < 0 || elementIndex >= n) {
4564    indexError(elementIndex, "setColumnBounds");
4565  }
4566#endif
4567  if (lowerValue < -1.0e27)
4568    lowerValue = -COIN_DBL_MAX;
4569  if (columnLower_[elementIndex] != lowerValue) {
4570    columnLower_[elementIndex] = lowerValue;
4571    // work arrays exist - update as well
4572    whatsChanged_ &= ~COLUMN_LOWER_SAME;
4573    if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4574      lowerSaved_[elementIndex+maximumAbcNumberRows_] = -COIN_DBL_MAX;
4575    } else {
4576      lowerSaved_[elementIndex+maximumAbcNumberRows_] = lowerValue * rhsScale_
4577        * inverseColumnUseScale_[elementIndex];
4578    }
4579  }
4580  if (upperValue > 1.0e27)
4581    upperValue = COIN_DBL_MAX;
4582  if (columnUpper_[elementIndex] != upperValue) {
4583    columnUpper_[elementIndex] = upperValue;
4584    // work arrays exist - update as well
4585    whatsChanged_ &= ~COLUMN_UPPER_SAME;
4586    if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4587      upperSaved_[elementIndex+maximumAbcNumberRows_] = COIN_DBL_MAX;
4588    } else {
4589      upperSaved_[elementIndex+maximumAbcNumberRows_] = upperValue * rhsScale_
4590        * inverseColumnUseScale_[elementIndex];
4591    }
4592  }
4593}
4594void AbcSimplex::setColumnSetBounds(const int* indexFirst,
4595                                    const int* indexLast,
4596                                    const double* boundList)
4597{
4598#ifndef NDEBUG
4599  int n = numberColumns_;
4600#endif
4601  int numberChanged = 0;
4602  const int * saveFirst = indexFirst;
4603  while (indexFirst != indexLast) {
4604    const int iColumn = *indexFirst++;
4605#ifndef NDEBUG
4606    if (iColumn < 0 || iColumn >= n) {
4607      indexError(iColumn, "setColumnSetBounds");
4608    }
4609#endif
4610    double lowerValue = *boundList++;
4611    double upperValue = *boundList++;
4612    if (lowerValue < -1.0e27)
4613      lowerValue = -COIN_DBL_MAX;
4614    if (upperValue > 1.0e27)
4615      upperValue = COIN_DBL_MAX;
4616    //CoinAssert (upperValue>=lowerValue);
4617    if (columnLower_[iColumn] != lowerValue) {
4618      columnLower_[iColumn] = lowerValue;
4619      whatsChanged_ &= ~COLUMN_LOWER_SAME;
4620      numberChanged++;
4621    }
4622    if (columnUpper_[iColumn] != upperValue) {
4623      columnUpper_[iColumn] = upperValue;
4624      whatsChanged_ &= ~COLUMN_UPPER_SAME;
4625      numberChanged++;
4626    }
4627  }
4628  if (numberChanged) {
4629    indexFirst = saveFirst;
4630    while (indexFirst != indexLast) {
4631      const int iColumn = *indexFirst++;
4632      if (columnLower_[iColumn] == -COIN_DBL_MAX) {
4633        lowerSaved_[iColumn+maximumAbcNumberRows_] = -COIN_DBL_MAX;
4634      } else {
4635        lowerSaved_[iColumn+maximumAbcNumberRows_] = columnLower_[iColumn] * rhsScale_
4636          * inverseColumnUseScale_[iColumn];
4637      }
4638      if (columnUpper_[iColumn] == COIN_DBL_MAX) {
4639        upperSaved_[iColumn+maximumAbcNumberRows_] = COIN_DBL_MAX;
4640      } else {
4641        upperSaved_[iColumn+maximumAbcNumberRows_] = columnUpper_[iColumn] * rhsScale_
4642          * inverseColumnUseScale_[iColumn];
4643      }
4644    }
4645  }
4646}
4647#ifndef SLIM_CLP
4648#include "CoinWarmStartBasis.hpp"
4649
4650// Returns a basis (to be deleted by user)
4651CoinWarmStartBasis *
4652AbcSimplex::getBasis() const
4653{
4654  CoinWarmStartBasis * basis = new CoinWarmStartBasis();
4655  basis->setSize(numberColumns_, numberRows_);
4656 
4657  unsigned char lookup[8]={2,3,255,255,0,0,1,3};
4658  for (int iRow = 0; iRow < numberRows_; iRow++) {
4659    int iStatus = getInternalStatus(iRow);
4660    iStatus = lookup[iStatus];
4661    basis->setArtifStatus(iRow, static_cast<CoinWarmStartBasis::Status> (iStatus));
4662  }
4663  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4664    int iStatus = getInternalStatus(iColumn+maximumAbcNumberRows_);
4665    iStatus = lookup[iStatus];
4666    basis->setStructStatus(iColumn, static_cast<CoinWarmStartBasis::Status> (iStatus));
4667  }
4668  return basis;
4669}
4670#endif
4671// Compute objective value from solution
4672void
4673AbcSimplex::computeObjectiveValue(bool useInternalArrays)
4674{
4675  const double * obj = objective();
4676  if (!useInternalArrays) {
4677    objectiveValue_ = 0.0;
4678    for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
4679      double value = columnActivity_[iSequence];
4680      objectiveValue_ += value * obj[iSequence];
4681    }
4682    // But remember direction as we are using external objective
4683    objectiveValue_ *= optimizationDirection_;
4684  } else {
4685    rawObjectiveValue_ = 0.0;
4686    for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4687      double value = abcSolution_[iSequence];
4688      rawObjectiveValue_ += value * abcCost_[iSequence];
4689    }
4690    setClpSimplexObjectiveValue();
4691  }
4692}
4693// Compute minimization objective value from internal solution
4694double
4695AbcSimplex::computeInternalObjectiveValue()
4696{
4697  rawObjectiveValue_ = 0.0;
4698  for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4699    double value = abcSolution_[iSequence];
4700    rawObjectiveValue_ += value * abcCost_[iSequence];
4701  }
4702  setClpSimplexObjectiveValue();
4703  return objectiveValue_-dblParam_[ClpObjOffset];
4704}
4705// If user left factorization frequency then compute
4706void
4707AbcSimplex::defaultFactorizationFrequency()
4708{
4709  if (factorizationFrequency() == 200) {
4710    // User did not touch preset
4711    const int cutoff1 = 10000;
4712    const int cutoff2 = 100000;
4713    const int base = 75;
4714    const int freq0 = 50;
4715    const int freq1 = 150;
4716    const int maximum = 10000;
4717    int frequency;
4718    if (numberRows_ < cutoff1)
4719      frequency = base + numberRows_ / freq0;
4720    else 
4721      frequency = base + cutoff1 / freq0 + (numberRows_ - cutoff1) / freq1;
4722    setFactorizationFrequency(CoinMin(maximum, frequency));
4723  }
4724}
4725// Gets clean and emptyish factorization
4726AbcSimplexFactorization *
4727AbcSimplex::getEmptyFactorization()
4728{
4729  if ((specialOptions_ & 65536) == 0) {
4730    assert (!abcFactorization_);
4731    abcFactorization_ = new AbcSimplexFactorization();
4732  } else if (!abcFactorization_) {
4733    abcFactorization_ = new AbcSimplexFactorization();
4734  }
4735  return abcFactorization_;
4736}
4737// Resizes rim part of model
4738void
4739AbcSimplex::resize (int newNumberRows, int newNumberColumns)
4740{
4741  assert (maximumAbcNumberRows_>=0);
4742  // save
4743  int numberRows=numberRows_;
4744  int numberColumns=numberColumns_;
4745  ClpSimplex::resize(newNumberRows, newNumberColumns);
4746  // back out
4747  numberRows_=numberRows;
4748  numberColumns_=numberColumns;
4749  gutsOfResize(newNumberRows,newNumberColumns);
4750}
4751// Return true if the objective limit test can be relied upon
4752bool
4753AbcSimplex::isObjectiveLimitTestValid() const
4754{
4755  if (problemStatus_ == 0) {
4756    return true;
4757  } else if (problemStatus_ == 1) {
4758    // ok if dual
4759    return (algorithm_ < 0);
4760  } else if (problemStatus_ == 2) {
4761    // ok if primal
4762    return (algorithm_ > 0);
4763  } else {
4764    return false;
4765  }
4766}
4767/*
4768  Permutes in from ClpModel data - assumes scale factors done
4769  and AbcMatrix exists but is in original order (including slacks)
4770  For now just add basicArray at end
4771  ==
4772  But could partition into
4773  normal (i.e. reasonable lower/upper)
4774  abnormal - free, odd bounds
4775  fixed
4776  ==
4777  sets a valid pivotVariable
4778  Slacks always shifted by offset
4779  Fixed variables always shifted by offset
4780  Recode to allow row objective so can use pi from idiot etc
4781*/
4782void 
4783AbcSimplex::permuteIn()
4784{
4785  // for now only if maximumAbcNumberRows_==numberRows_
4786  //assert(maximumAbcNumberRows_==numberRows_);
4787  numberTotal_=maximumAbcNumberRows_+numberColumns_;
4788  double direction = optimizationDirection_;
4789  // all this could use cilk
4790  // move primal stuff to end
4791  const double *  COIN_RESTRICT rowScale=scaleFromExternal_;
4792  const double *  COIN_RESTRICT inverseRowScale=scaleToExternal_;
4793  const double *  COIN_RESTRICT columnScale=scaleToExternal_+maximumAbcNumberRows_;
4794  const double *  COIN_RESTRICT inverseColumnScale=scaleFromExternal_+maximumAbcNumberRows_;
4795  double *  COIN_RESTRICT offsetRhs=offsetRhs_;
4796  double *  COIN_RESTRICT saveLower=lowerSaved_+maximumAbcNumberRows_;
4797  double *  COIN_RESTRICT saveUpper=upperSaved_+maximumAbcNumberRows_;
4798  double *  COIN_RESTRICT saveCost=costSaved_+maximumAbcNumberRows_;
4799  double *  COIN_RESTRICT saveSolution=solutionSaved_+maximumAbcNumberRows_;
4800  CoinAbcMemset0(offsetRhs,maximumAbcNumberRows_);
4801  const double *  COIN_RESTRICT objective=this->objective();
4802  objectiveOffset_=0.0;
4803  double *  COIN_RESTRICT offset=offset_+maximumAbcNumberRows_;
4804  const int * COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
4805  const CoinBigIndex * COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
4806  const double * COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
4807  largestGap_=1.0e-12;
4808  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
4809    double scale=inverseColumnScale[iColumn];
4810    double lowerValue=columnLower_[iColumn];
4811    double upperValue=columnUpper_[iColumn];
4812    double thisOffset=0.0;
4813    double lowerValue2 = fabs(lowerValue);
4814    double upperValue2 = fabs(upperValue);
4815    if (CoinMin(lowerValue2,upperValue2)<1000.0) {
4816      // move to zero
4817      if (lowerValue2>upperValue2) 
4818        thisOffset=upperValue;
4819      else
4820        thisOffset=lowerValue;
4821    }
4822    offset[iColumn]=thisOffset;
4823    objectiveOffset_ += thisOffset*objective[iColumn]*optimizationDirection_;
4824    double scaledOffset = thisOffset*scale;
4825    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
4826      int iRow=row[j];
4827      offsetRhs[iRow]+= scaledOffset*elementByColumn[j];
4828    }
4829    lowerValue-=thisOffset;
4830    if (lowerValue>-1.0e30) 
4831      lowerValue*=scale;
4832    saveLower[iColumn]=lowerValue;
4833    upperValue-=thisOffset;
4834    if (upperValue<1.0e30) 
4835      upperValue*=scale;
4836    saveUpper[iColumn]=upperValue;
4837    largestGap_=CoinMax(largestGap_,upperValue-lowerValue);
4838    saveSolution[iColumn]=scale*(columnActivity_[iColumn]-thisOffset);
4839    saveCost[iColumn]=objective[iColumn]*direction*columnScale[iColumn];
4840  }
4841  CoinAbcMemset0(offset_,maximumAbcNumberRows_);
4842  saveLower-=maximumAbcNumberRows_;
4843  saveUpper-=maximumAbcNumberRows_;
4844  saveCost-=maximumAbcNumberRows_;
4845  saveSolution-=maximumAbcNumberRows_;
4846  for (int iRow=0;iRow<numberRows_;iRow++) {
4847    // note switch of lower to upper
4848    double upperValue=-rowLower_[iRow];
4849    double lowerValue=-rowUpper_[iRow];
4850    double thisOffset=offsetRhs_[iRow];
4851    double scale=rowScale[iRow];
4852    if (lowerValue>-1.0e30) 
4853      lowerValue=lowerValue*scale+thisOffset;
4854    saveLower[iRow]=lowerValue;
4855    if (upperValue<1.0e30) 
4856      upperValue=upperValue*scale+thisOffset;
4857    saveUpper[iRow]=upperValue;
4858    largestGap_=CoinMax(largestGap_,upperValue-lowerValue);
4859    saveCost[iRow]=0.0;
4860    dual_[iRow]*=direction*inverseRowScale[iRow];
4861    saveSolution[iRow]=0.0; // not necessary
4862  }
4863  dualBound_=CoinMin(dualBound_,largestGap_);
4864  // Compute rhsScale_ and objectiveScale_
4865  double minValue=COIN_DBL_MAX;
4866  double maxValue=0.0;
4867  CoinAbcMinMaxAbsNormalValues(costSaved_+maximumAbcNumberRows_,numberTotal_-maximumAbcNumberRows_,minValue,maxValue);
4868  // scale to 1000.0 ?
4869  if (minValue&&false) {
4870    objectiveScale_= 1000.0/sqrt(minValue*maxValue);
4871    objectiveScale_=CoinMin(1.0,1000.0/maxValue);
4872#ifndef NDEBUG
4873    double smallestNormal=COIN_DBL_MAX;
4874    double smallestAny=COIN_DBL_MAX;
4875    double largestAny=0.0;
4876    for (int i=0;i<numberTotal_;i++) {
4877      double value=fabs(costSaved_[i]);
4878      if (value) {
4879        if (value>1.0e-8)
4880          smallestNormal=CoinMin(smallestNormal,value);
4881        smallestAny=CoinMin(smallestAny,value);
4882        largestAny=CoinMax(largestAny,value);
4883      }
4884    }
4885    printf("objectiveScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
4886           objectiveScale_,minValue,smallestNormal,smallestAny,maxValue,largestAny);
4887#endif
4888  } else {
4889    //objectiveScale_=1.0;
4890  }
4891  CoinAbcScale(costSaved_,objectiveScale_,numberTotal_);
4892  minValue=COIN_DBL_MAX;
4893  maxValue=0.0;
4894  CoinAbcMinMaxAbsNormalValues(lowerSaved_,numberTotal_,minValue,maxValue);
4895  CoinAbcMinMaxAbsNormalValues(upperSaved_,numberTotal_,minValue,maxValue);
4896  // scale to 100.0 ?
4897  if (minValue&&false) {
4898    rhsScale_= 100.0/sqrt(minValue*maxValue);
4899#ifndef NDEBUG
4900    double smallestNormal=COIN_DBL_MAX;
4901    double smallestAny=COIN_DBL_MAX;
4902    double largestAny=0.0;
4903    for (int i=0;i<numberTotal_;i++) {
4904      double value=fabs(lowerSaved_[i]);
4905      if (value&&value!=COIN_DBL_MAX) {
4906        if (value>1.0e-8)
4907          smallestNormal=CoinMin(smallestNormal,value);
4908        smallestAny=CoinMin(smallestAny,value);
4909        largestAny=CoinMax(largestAny,value);
4910      }
4911    }
4912    for (int i=0;i<numberTotal_;i++) {
4913      double value=fabs(upperSaved_[i]);
4914      if (value&&value!=COIN_DBL_MAX) {
4915        if (value>1.0e-8)
4916          smallestNormal=CoinMin(smallestNormal,value);
4917        smallestAny=CoinMin(smallestAny,value);
4918        largestAny=CoinMax(largestAny,value);
4919      }
4920    }
4921    printf("rhsScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
4922           rhsScale_,minValue,smallestNormal,smallestAny,maxValue,largestAny);
4923#endif
4924  } else {
4925    rhsScale_=1.0;
4926  }
4927  CoinAbcScaleNormalValues(lowerSaved_,rhsScale_,1.0e-13,numberTotal_);
4928  CoinAbcScaleNormalValues(upperSaved_,rhsScale_,1.0e-13,numberTotal_);
4929  // copy
4930  CoinAbcMemcpy(abcLower_,abcLower_+maximumNumberTotal_,numberTotal_);
4931  CoinAbcMemcpy(abcUpper_,abcUpper_+maximumNumberTotal_,numberTotal_);
4932  CoinAbcMemcpy(abcSolution_,abcSolution_+maximumNumberTotal_,numberTotal_);
4933  CoinAbcMemcpy(abcCost_,abcCost_+maximumNumberTotal_,numberTotal_);
4934}
4935void 
4936AbcSimplex::permuteBasis()
4937{
4938  assert (abcPivotVariable_||(!numberRows_&&!numberColumns_));
4939  int numberBasic=0;
4940  // from Clp enum to Abc enum (and bound flip)
4941  unsigned char lookupToAbcSlack[6]={4,6,0/*1*/,1/*0*/,5,7};
4942  unsigned char *  COIN_RESTRICT getStatus = status_+numberColumns_;
4943  double *  COIN_RESTRICT solutionSaved=solutionSaved_;
4944  double *  COIN_RESTRICT lowerSaved=lowerSaved_;
4945  double *  COIN_RESTRICT upperSaved=upperSaved_;
4946  bool ordinaryVariables=true;
4947  bool valuesPass=(stateOfProblem_&VALUES_PASS)!=0;
4948  if (valuesPass) {
4949    // get solution
4950    CoinAbcMemset0(abcSolution_,numberRows_);
4951    abcMatrix_->timesIncludingSlacks(-1.0,abcSolution_,abcSolution_);
4952    //double * temp = new double[numberRows_];
4953    //memset(temp,0,numberRows_*sizeof(double));
4954    //matrix_->times(1.0,columnActivity_,temp);
4955    CoinAbcMemcpy(solutionSaved_,abcSolution_,numberRows_);
4956    int n=0;
4957    for (int i=0;i<numberTotal_;i++) {
4958      if (solutionSaved_[i]>upperSaved_[i]+1.0e-5)
4959        n++;
4960      else if (solutionSaved_[i]<lowerSaved_[i]-1.0e-5)
4961        n++;
4962    }
4963#if ABC_NORMAL_DEBUG>0
4964    if (n)
4965      printf("%d infeasibilities\n",n);
4966#endif
4967  }
4968  // dual at present does not like superbasic
4969  for (int iRow=0;iRow<numberRows_;iRow++) {
4970    unsigned char status=getStatus[iRow]&7;
4971    AbcSimplex::Status abcStatus=static_cast<AbcSimplex::Status>(lookupToAbcSlack[status]);
4972    if (status!=ClpSimplex::basic) {
4973      double lowerValue=lowerSaved[iRow];
4974      double upperValue=upperSaved[iRow];
4975      if (lowerValue==-COIN_DBL_MAX) {
4976        if(upperValue==COIN_DBL_MAX) {
4977          // free
4978          abcStatus=isFree;
4979          ordinaryVariables=false;
4980        } else {
4981          abcStatus=atUpperBound;
4982        }
4983      } else if (upperValue==COIN_DBL_MAX) {
4984        abcStatus=atLowerBound;
4985      } else if (lowerValue==upperValue) {
4986        abcStatus=isFixed;
4987      } else if (abcStatus==isFixed) {
4988        double value=solutionSaved[iRow];
4989        if (value-lowerValue<upperValue-value)
4990          abcStatus=atLowerBound;
4991        else
4992          abcStatus=atUpperBound;
4993      }
4994      if (valuesPass) {
4995        double value=solutionSaved[iRow];
4996        if (value>lowerValue+primalTolerance_&&value<upperValue-primalTolerance_&&
4997            (abcStatus==atLowerBound||abcStatus==atUpperBound))
4998            abcStatus=superBasic;
4999      }
5000      switch (abcStatus) {
5001      case isFixed:
5002      case atLowerBound:
5003        solutionSaved[iRow]=lowerValue;
5004        break;
5005      case atUpperBound:
5006        solutionSaved[iRow]=upperValue;
5007        break;
5008      default:
5009        ordinaryVariables=false;
5010        break;
5011      }
5012    } else {
5013      // basic
5014      abcPivotVariable_[numberBasic++]=iRow;
5015    }
5016    internalStatus_[iRow]=abcStatus;
5017  }
5018  // from Clp enum to Abc enum
5019  unsigned char lookupToAbc[6]={4,6,1,0,5,7};
5020  unsigned char *  COIN_RESTRICT putStatus=internalStatus_+maximumAbcNumberRows_;
5021  getStatus = status_;
5022  solutionSaved+= maximumAbcNumberRows_;
5023  lowerSaved+= maximumAbcNumberRows_;
5024  upperSaved+= maximumAbcNumberRows_;
5025  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
5026    unsigned char status=getStatus[iColumn]&7;
5027    AbcSimplex::Status abcStatus=static_cast<AbcSimplex::Status>(lookupToAbc[status]);
5028    if (status!=ClpSimplex::basic) {
5029      double lowerValue=lowerSaved[iColumn];
5030      double upperValue=upperSaved[iColumn];
5031      if (lowerValue==-COIN_DBL_MAX) {
5032        if(upperValue==COIN_DBL_MAX) {
5033          // free
5034          abcStatus=isFree;
5035          ordinaryVariables=false;
5036        } else {
5037          abcStatus=atUpperBound;
5038        }
5039      } else if (upperValue==COIN_DBL_MAX) {
5040        abcStatus=atLowerBound;
5041      } else if (lowerValue==upperValue) {
5042        abcStatus=isFixed;
5043      } else if (abcStatus==isFixed) {
5044        double value=solutionSaved[iColumn];
5045        if (value-lowerValue<upperValue-value)
5046          abcStatus=atLowerBound;
5047        else
5048          abcStatus=atUpperBound;
5049      } else if (abcStatus==isFree) {
5050        abcStatus=superBasic;
5051      }
5052      if (valuesPass&&(abcStatus==atLowerBound||abcStatus==atUpperBound)) {
5053        double value=solutionSaved[iColumn];
5054        if (value>lowerValue+primalTolerance_) {
5055          if(value<upperValue-primalTolerance_) {
5056            abcStatus=superBasic;
5057          } else {
5058            abcStatus=atUpperBound;
5059          }
5060        } else {
5061          abcStatus=atLowerBound;
5062        }
5063      }
5064      switch (abcStatus) {
5065      case isFixed:
5066      case atLowerBound:
5067        solutionSaved[iColumn]=lowerValue;
5068        break;
5069      case atUpperBound:
5070        solutionSaved[iColumn]=upperValue;
5071        break;
5072      default:
5073        ordinaryVariables=false;
5074        break;
5075      }
5076    } else {
5077      // basic
5078      if (numberBasic<numberRows_) 
5079        abcPivotVariable_[numberBasic++]=iColumn+maximumAbcNumberRows_;
5080      else
5081        abcStatus=superBasic;
5082    }
5083    putStatus[iColumn]=abcStatus;
5084  }
5085  ordinaryVariables_ =  ordinaryVariables ? 1 : 0;
5086  if (numberBasic<numberRows_) {
5087    for (int iRow=0;iRow<numberRows_;iRow++) {
5088      AbcSimplex::Status status=getInternalStatus(iRow);
5089      if (status!=AbcSimplex::basic) {
5090        abcPivotVariable_[numberBasic++]=iRow;
5091        setInternalStatus(iRow,basic);
5092        if (numberBasic==numberRows_)
5093          break;
5094      }
5095    }
5096  }
5097  // copy
5098  CoinAbcMemcpy(internalStatus_+maximumNumberTotal_,internalStatus_,numberTotal_);
5099}
5100// Permutes out - bit settings same as stateOfProblem
5101void 
5102AbcSimplex::permuteOut(int whatsWanted)
5103{
5104  assert (whatsWanted);
5105  if ((whatsWanted&ALL_STATUS_OK)!=0&&(stateOfProblem_&ALL_STATUS_OK)==0) {
5106    stateOfProblem_ |= ALL_STATUS_OK;
5107    // from Abc enum to Clp enum (and bound flip)
5108    unsigned char lookupToClpSlack[8]={2,3,255,255,0,0,1,5};
5109    unsigned char *  COIN_RESTRICT putStatus = status_+numberColumns_;
5110    const unsigned char *  COIN_RESTRICT getStatus=internalStatus_;
5111    for (int iRow=0;iRow<numberRows_;iRow++) {
5112      putStatus[iRow]=lookupToClpSlack[getStatus[iRow]&7];
5113    }
5114    // from Abc enum to Clp enum
5115    unsigned char lookupToClp[8]={3,2,255,255,0,0,1,5};
5116    putStatus = status_;
5117    getStatus=internalStatus_+maximumAbcNumberRows_;
5118    for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
5119      putStatus[iColumn]=lookupToClp[getStatus[iColumn]&7];
5120    }
5121  }
5122  const double *  COIN_RESTRICT rowScale=scaleFromExternal_;
5123  const double *  COIN_RESTRICT inverseRowScale=scaleToExternal_;
5124  const double *  COIN_RESTRICT columnScale=scaleToExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5125  const double *  COIN_RESTRICT inverseColumnScale=scaleFromExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5126  double scaleC = 1.0 / objectiveScale_;
5127  double scaleR = 1.0 / rhsScale_;
5128  int numberPrimalScaled = 0;
5129  int numberPrimalUnscaled = 0;
5130  int numberDualScaled = 0;
5131  int numberDualUnscaled = 0;
5132  bool filledInSolution=false;
5133  if ((whatsWanted&ROW_PRIMAL_OK)!=0&&(stateOfProblem_&ROW_PRIMAL_OK)==0) {
5134    stateOfProblem_ |= ROW_PRIMAL_OK;
5135    if (!filledInSolution) {
5136      filledInSolution=true;
5137      CoinAbcScatterTo(solutionBasic_,abcSolution_,abcPivotVariable_,numberRows_);
5138      CoinAbcScatterZeroTo(abcDj_,abcPivotVariable_,numberRows_);
5139    }
5140    // Collect infeasibilities
5141    double *  COIN_RESTRICT putSolution = rowActivity_;
5142    const double *  COIN_RESTRICT rowLower =rowLower_;
5143    const double *  COIN_RESTRICT rowUpper =rowUpper_;
5144    const double *  COIN_RESTRICT abcLower =abcLower_;
5145    const double *  COIN_RESTRICT abcUpper =abcUpper_;
5146    const double *  COIN_RESTRICT abcSolution =abcSolution_;
5147    const double *  COIN_RESTRICT offsetRhs =offsetRhs_;
5148    for (int i = 0; i < numberRows_; i++) {
5149      double scaleFactor = inverseRowScale[i];
5150      double valueScaled = abcSolution[i];
5151      double lowerScaled = abcLower[i];
5152      double upperScaled = abcUpper[i];
5153      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5154        if (valueScaled < lowerScaled - primalTolerance_ ||
5155            valueScaled > upperScaled + primalTolerance_)
5156          numberPrimalScaled++;
5157        else
5158          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5159      }
5160      double value = (offsetRhs[i]-valueScaled * scaleR) * scaleFactor;
5161      putSolution[i]=value;
5162      if (value < rowLower[i] - primalTolerance_)
5163        numberPrimalUnscaled++;
5164      else if (value > rowUpper[i] + primalTolerance_)
5165        numberPrimalUnscaled++;
5166    }
5167  }
5168  if ((whatsWanted&COLUMN_PRIMAL_OK)!=0&&(stateOfProblem_&COLUMN_PRIMAL_OK)==0) {
5169    stateOfProblem_ |= COLUMN_PRIMAL_OK;
5170    // Collect infeasibilities
5171    if (!filledInSolution) {
5172      filledInSolution=true;
5173      CoinAbcScatterTo(solutionBasic_,abcSolution_,abcPivotVariable_,numberRows_);
5174      CoinAbcScatterZeroTo(abcDj_,abcPivotVariable_,numberRows_);
5175    }
5176    // Collect infeasibilities
5177    double *  COIN_RESTRICT putSolution = columnActivity_-maximumAbcNumberRows_;
5178    const double *  COIN_RESTRICT columnLower=columnLower_-maximumAbcNumberRows_;
5179    const double *  COIN_RESTRICT columnUpper=columnUpper_-maximumAbcNumberRows_;
5180    for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5181      double scaleFactor = columnScale[i];
5182      double valueScaled = abcSolution_[i];
5183      double lowerScaled = abcLower_[i];
5184      double upperScaled = abcUpper_[i];
5185      if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5186        if (valueScaled < lowerScaled - primalTolerance_ ||
5187            valueScaled > upperScaled + primalTolerance_)
5188          numberPrimalScaled++;
5189        else
5190          upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5191      }
5192      double value = (valueScaled * scaleR) * scaleFactor+offset_[i];
5193      putSolution[i]=value;
5194      if (value < columnLower[i] - primalTolerance_)
5195        numberPrimalUnscaled++;
5196      else if (value > columnUpper[i] + primalTolerance_)
5197        numberPrimalUnscaled++;
5198    }
5199  }
5200  if ((whatsWanted&COLUMN_DUAL_OK)!=0&&(stateOfProblem_&COLUMN_DUAL_OK)==0) {
5201    // Get fixed djs
5202    CoinAbcMemcpy(abcDj_,abcCost_,numberTotal_);
5203    abcMatrix_->transposeTimesAll( dual_,abcDj_);
5204    stateOfProblem_ |= COLUMN_DUAL_OK;
5205    // Collect infeasibilities
5206    double *  COIN_RESTRICT putDj=reducedCost_-maximumAbcNumberRows_;
5207    for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5208      double scaleFactor = inverseColumnScale[i];
5209      double valueDual = abcDj_[i];
5210      double value=abcSolution_[i];
5211      bool aboveLower= value>abcLower_[i]+primalTolerance_;
5212      bool belowUpper= value<abcUpper_[i]-primalTolerance_;
5213      if (aboveLower&& valueDual > dualTolerance_)
5214        numberDualScaled++;
5215      if (belowUpper && valueDual < -dualTolerance_)
5216        numberDualScaled++;
5217      valueDual *= scaleFactor * scaleC;
5218      putDj[i]=valueDual;
5219      if (aboveLower&& valueDual > dualTolerance_)
5220        numberDualUnscaled++;
5221      if (belowUpper && valueDual < -dualTolerance_)
5222        numberDualUnscaled++;
5223    }
5224  }
5225  if ((whatsWanted&ROW_DUAL_OK)!=0&&(stateOfProblem_&ROW_DUAL_OK)==0) {
5226    stateOfProblem_ |= ROW_DUAL_OK;
5227    // Collect infeasibilities
5228    for (int i = 0; i < numberRows_; i++) {
5229      double scaleFactor = rowScale[i];
5230      double valueDual = abcDj_[i]; // ? +- and direction
5231      double value=abcSolution_[i];
5232      bool aboveLower= value>abcLower_[i]+primalTolerance_;
5233      bool belowUpper= value<abcUpper_[i]-primalTolerance_;
5234      if (aboveLower&& valueDual > dualTolerance_)
5235        numberDualScaled++;
5236      if (belowUpper && valueDual < -dualTolerance_)
5237        numberDualScaled++;
5238      valueDual *= scaleFactor * scaleC;
5239      dual_[i]=-(valueDual-abcCost_[i]); // sign
5240      if (aboveLower&& valueDual > dualTolerance_)
5241        numberDualUnscaled++;
5242      if (belowUpper && valueDual < -dualTolerance_)
5243        numberDualUnscaled++;
5244    }
5245  }
5246  // do after djs
5247  if (!problemStatus_ && (!secondaryStatus_||secondaryStatus_==2||secondaryStatus_==3)) {
5248    // See if we need to set secondary status
5249    if (numberPrimalUnscaled) {
5250      if (numberDualUnscaled||secondaryStatus_==3)
5251        secondaryStatus_ = 4;
5252      else
5253        secondaryStatus_ = 2;
5254    } else if (numberDualUnscaled) {
5255      if (secondaryStatus_==0)
5256        secondaryStatus_ = 3;
5257      else
5258        secondaryStatus_ = 4;
5259    }
5260  }
5261  if (problemStatus_ == 2) {
5262    for (int i = 0; i < numberColumns_; i++) {
5263      ray_[i] *= columnScale[i];
5264    }
5265  } else if (problemStatus_ == 1 && ray_) {
5266    for (int i = 0; i < numberRows_; i++) {
5267      ray_[i] *= rowScale[i];
5268    }
5269  }
5270}
5271#if ABC_DEBUG>1
5272// For debug - moves solution back to external and computes stuff
5273void 
5274AbcSimplex::checkMoveBack(bool checkDuals)
5275{
5276  stateOfProblem_ &= ~(ROW_PRIMAL_OK|ROW_DUAL_OK|COLUMN_PRIMAL_OK|COLUMN_DUAL_OK|ALL_STATUS_OK);
5277  permuteOut(ROW_PRIMAL_OK|ROW_DUAL_OK|COLUMN_PRIMAL_OK|COLUMN_DUAL_OK|ALL_STATUS_OK);
5278  ClpSimplex::computeObjectiveValue(false);
5279#if ABC_NORMAL_DEBUG>0
5280  printf("Check objective %g\n",objectiveValue()-objectiveOffset_);
5281#endif
5282  double * region = new double [numberRows_+numberColumns_];
5283  CoinAbcMemset0(region,numberRows_);
5284  ClpModel::times(1.0,columnActivity_,region);
5285  int numberInf;
5286  double sumInf;
5287  numberInf=0;
5288  sumInf=0.0;
5289  for (int i=0;i<numberColumns_;i++) {
5290    if (columnActivity_[i]>columnUpper_[i]+1.0e-7) {
5291      numberInf++;
5292      sumInf+=columnActivity_[i]-columnUpper_[i];
5293    } else if (columnActivity_[i]<columnLower_[i]-1.0e-7) {
5294      numberInf++;
5295      sumInf-=columnActivity_[i]-columnLower_[i];
5296    }
5297  }
5298#if ABC_NORMAL_DEBUG>0
5299  if (numberInf)
5300    printf("Check column infeasibilities %d sum %g\n",numberInf,sumInf);
5301#endif
5302  numberInf=0;
5303  sumInf=0.0;
5304  for (int i=0;i<numberRows_;i++) {
5305    if (rowActivity_[i]>rowUpper_[i]+1.0e-7) {
5306      numberInf++;
5307      sumInf+=rowActivity_[i]-rowUpper_[i];
5308    } else if (rowActivity_[i]<rowLower_[i]-1.0e-7) {
5309      numberInf++;
5310      sumInf-=rowActivity_[i]-rowLower_[i];
5311    }
5312  }
5313#if ABC_NORMAL_DEBUG>0
5314  if (numberInf)
5315    printf("Check row infeasibilities %d sum %g\n",numberInf,sumInf);
5316#endif
5317  CoinAbcMemcpy(region,objective(),numberColumns_);
5318  ClpModel::transposeTimes(-1.0,dual_,region);
5319  numberInf=0;
5320  sumInf=0.0;
5321  for (int i=0;i<numberColumns_;i++) {
5322    if (columnUpper_[i]>columnLower_[i]) {
5323      if (getColumnStatus(i)==ClpSimplex::atLowerBound) {
5324        if (region[i]<-1.0e-7) {
5325          numberInf++;
5326          sumInf-=region[i];
5327        }
5328      } else if (getColumnStatus(i)==ClpSimplex::atUpperBound) {
5329        if (region[i]>1.0e-7) {
5330          numberInf++;
5331          sumInf+=region[i];
5332        }
5333      }
5334    }
5335  }
5336#if ABC_NORMAL_DEBUG>0
5337  if (numberInf)
5338    printf("Check column dj infeasibilities %d sum %g\n",numberInf,sumInf);
5339#endif
5340  if (checkDuals) {
5341    numberInf=0;
5342    sumInf=0.0;
5343    for (int i=0;i<numberRows_;i++) {
5344      if (rowUpper_[i]>rowLower_[i]) {
5345        if (getRowStatus(i)==ClpSimplex::atLowerBound) {
5346          if (dual_[i]<-1.0e-7) {
5347            numberInf++;
5348            sumInf-=region[i];
5349          }
5350        } else if (getRowStatus(i)==ClpSimplex::atUpperBound) {
5351          if (dual_[i]>1.0e-7) {
5352            numberInf++;
5353            sumInf+=region[i];
5354          }
5355        }
5356      }
5357    }
5358#if ABC_NORMAL_DEBUG>0
5359    if (numberInf)
5360      printf("Check row dual infeasibilities %d sum %g\n",numberInf,sumInf);
5361#endif
5362  }
5363  delete [] region;
5364}
5365#if 0
5366void xxxxxx(const char * where) {
5367  printf("xxxxx %s\n",where);
5368}
5369#endif
5370// For debug - checks solutionBasic
5371void 
5372AbcSimplex::checkSolutionBasic() const
5373{
5374  //work space
5375  int whichArray[2];
5376  for (int i=0;i<2;i++)
5377    whichArray[i]=getAvailableArray();
5378  CoinIndexedVector * arrayVector = &usefulArray_[whichArray[0]];
5379  double * solution = usefulArray_[whichArray[1]].denseVector();
5380  CoinAbcMemcpy(solution,abcSolution_,numberTotal_);
5381  // accumulate non basic stuff
5382  double * array = arrayVector->denseVector();
5383  CoinAbcScatterZeroTo(solution,abcPivotVariable_,numberRows_);
5384  abcMatrix_->timesIncludingSlacks(-1.0, solution, array);
5385  //arrayVector->scan(0,numberRows_,zeroTolerance_);
5386  // Ftran adjusted RHS
5387  //if (arrayVector->getNumElements())
5388  abcFactorization_->updateFullColumn(*arrayVector);
5389  CoinAbcScatterTo(array,solution,abcPivotVariable_,numberRows_);
5390  double largestDifference=0.0;
5391  int whichDifference=-1;
5392  for (int i=0;i<numberRows_;i++) {
5393    double difference = fabs(solutionBasic_[i]-array[i]);
5394    if (difference>1.0e-5&&numberRows_<100)
5395      printf("solutionBasic difference is %g on row %d solutionBasic_ %g computed %g\n",
5396           difference,i,solutionBasic_[i],array[i]);
5397    if (difference>largestDifference) {
5398      largestDifference=difference;
5399      whichDifference=i;
5400    }
5401  }
5402  if (largestDifference>1.0e-9)
5403    printf("Debug largest solutionBasic difference is %g on row %d solutionBasic_ %g computed %g\n",
5404           largestDifference,whichDifference,solutionBasic_[whichDifference],array[whichDifference]);
5405    arrayVector->clear();
5406  CoinAbcMemset0(solution,numberTotal_);
5407  for (int i=0;i<2;i++)
5408    setAvailableArray(whichArray[i]);
5409}
5410
5411// For debug - summarizes dj situation
5412void 
5413AbcSimplex::checkDjs(int type) const
5414{
5415  if (type) {
5416    //work space
5417    int whichArrays[2];
5418    for (int i=0;i<2;i++)
5419      whichArrays[i]=getAvailableArray();
5420   
5421    CoinIndexedVector * arrayVector = &usefulArray_[whichArrays[0]];
5422    double * array = arrayVector->denseVector();
5423    int * index = arrayVector->getIndices();
5424    int number = 0;
5425    for (int iRow = 0; iRow < numberRows_; iRow++) {
5426      double value = costBasic_[iRow];
5427      if (value) {
5428        array[iRow] = value;
5429        index[number++] = iRow;
5430      }
5431    }
5432    arrayVector->setNumElements(number);
5433    // Btran basic costs
5434    abcFactorization_->updateFullColumnTranspose(*arrayVector);
5435    double largestDifference=0.0;
5436    int whichDifference=-1;
5437    if (type==2) {
5438      for (int i=0;i<numberRows_;i++) {
5439        double difference = fabs(dual_[i]-array[i]);
5440        if (difference>largestDifference) {
5441          largestDifference=difference;
5442          whichDifference=i;
5443        }
5444      }
5445      if (largestDifference>1.0e-9)
5446        printf("Debug largest dual difference is %g on row %d dual_ %g computed %g\n",
5447               largestDifference,whichDifference,dual_[whichDifference],array[whichDifference]);
5448    }
5449    // now look at djs
5450    double * djs=usefulArray_[whichArrays[1]].denseVector();
5451    CoinAbcMemcpy(djs,abcCost_,numberTotal_);
5452    abcMatrix_->transposeTimesNonBasic(-1.0, array,djs);
5453    largestDifference=0.0;
5454    whichDifference=-1;
5455    for (int i=0;i<numberTotal_;i++) {
5456      Status status = getInternalStatus(i);
5457      double difference = fabs(abcDj_[i]-djs[i]);
5458      switch(status) {
5459      case basic:
5460        // probably not an error
5461        break;
5462      case AbcSimplex::isFixed:
5463        break;
5464      case isFree:
5465      case superBasic:
5466      case atUpperBound:
5467      case atLowerBound:
5468        if (difference>1.0e-5&&numberTotal_<200) 
5469          printf("Debug dj difference is %g on column %d abcDj_ %g computed %g\n",
5470                 difference,i,abcDj_[i],djs[i]);
5471        if (difference>largestDifference) {
5472          largestDifference=difference;
5473          whichDifference=i;
5474        }
5475        break;
5476      }
5477    }
5478    if (largestDifference>1.0e-9)
5479      printf("Debug largest dj difference is %g on column %d abcDj_ %g computed %g\n",
5480             largestDifference,whichDifference,abcDj_[whichDifference],djs[whichDifference]);
5481    CoinAbcMemset0(djs,numberTotal_);
5482    arrayVector->clear();
5483    for (int i=0;i<2;i++)
5484      setAvailableArray(whichArrays[i]);
5485  }
5486  int state[8]={0,0,0,0,0,0,0,0};
5487  int badT[8]={0,0,0,0,0,0,0,0};
5488  //int badP[8]={0,0,0,0,0,0,0,0};
5489  int numberInfeasibilities=0;
5490  for (int i=0;i<numberTotal_;i++) {
5491    int iStatus = internalStatus_[i] & 7;
5492    state[iStatus]++;
5493    Status status = getInternalStatus(i);
5494    double djValue=abcDj_[i];
5495    double value=abcSolution_[i];
5496    double lowerValue=abcLower_[i];
5497    double upperValue=abcUpper_[i];
5498    //double perturbationValue=abcPerturbation_[i];
5499    switch(status) {
5500    case basic:
5501      // probably not an error
5502      if(fabs(djValue)>currentDualTolerance_)
5503        badT[iStatus]++;
5504      //if(fabs(djValue)>perturbationValue)
5505      //badP[iStatus]++;
5506      break;
5507    case AbcSimplex::isFixed:
5508      break;
5509    case isFree:
5510    case superBasic:
5511      if(fabs(djValue)>currentDualTolerance_)
5512        badT[iStatus]++;
5513      //if(fabs(djValue)>perturbationValue)
5514      //badP[iStatus]++;
5515      break;
5516    case atUpperBound:
5517      if (fabs(value - upperValue) > primalTolerance_)
5518        numberInfeasibilities++;
5519      if(djValue>currentDualTolerance_)
5520        badT[iStatus]++;
5521      //if(djValue>perturbationValue)
5522      //badP[iStatus]++;
5523      break;
5524    case atLowerBound:
5525      if (fabs(value - lowerValue) > primalTolerance_)
5526        numberInfeasibilities++;
5527      if(-djValue>currentDualTolerance_)
5528        badT[iStatus]++;
5529      //if(-djValue>perturbationValue)
5530      //badP[iStatus]++;
5531      break;
5532    }
5533  }
5534  if (numberInfeasibilities)
5535    printf("Debug %d variables away from bound when should be\n",numberInfeasibilities);
5536  int numberBad=0;
5537  for (int i=0;i<8;i++) {
5538    numberBad += badT[i]/*+badP[i]*/;
5539  }
5540  if (numberBad) {
5541    const char * type[]={"atLowerBound",
5542                         "atUpperBound",
5543                         "??",
5544                         "??",
5545                         "isFree",
5546                         "superBasic",
5547                         "basic",
5548                         "isFixed"};
5549    for (int i=0;i<8;i++) {
5550      if (state[i])
5551        printf("Debug - %s %d total %d bad with tolerance %d bad with perturbation\n",
5552               type[i],state[i],badT[i],0/*badP[i]*/);
5553    }
5554  }
5555}
5556#else
5557// For debug - moves solution back to external and computes stuff
5558void 
5559AbcSimplex::checkMoveBack(bool )
5560{
5561}
5562// For debug - checks solutionBasic
5563void 
5564AbcSimplex::checkSolutionBasic() const
5565{
5566}
5567
5568// For debug - summarizes dj situation
5569void 
5570AbcSimplex::checkDjs(int) const
5571{
5572}
5573#endif
5574#ifndef EARLY_FACTORIZE
5575#define ABC_NUMBER_USEFUL_NORMAL ABC_NUMBER_USEFUL
5576#else
5577#define ABC_NUMBER_USEFUL_NORMAL ABC_NUMBER_USEFUL-1
5578#endif
5579static double * elAddress[ABC_NUMBER_USEFUL_NORMAL];
5580// For debug - prints summary of arrays which are out of kilter
5581void 
5582AbcSimplex::checkArrays(int ignoreEmpty) const
5583{
5584  if (!numberIterations_||!elAddress[0]) {
5585    for (int i=0;i<ABC_NUMBER_USEFUL_NORMAL;i++) 
5586      elAddress[i]=usefulArray_[i].denseVector();
5587  } else {
5588    for (int i=0;i<ABC_NUMBER_USEFUL_NORMAL;i++) 
5589      assert(elAddress[i]==usefulArray_[i].denseVector());
5590  }
5591  for (int i=0;i<ABC_NUMBER_USEFUL_NORMAL;i++) {
5592    int check=1<<i;
5593    if (usefulArray_[i].getNumElements()) {
5594      if ((stateOfProblem_&check)==0)
5595        printf("Array %d has %d elements but says it is available!\n",
5596               i,usefulArray_[i].getNumElements());
5597      if (usefulArray_[i].getNumPartitions())
5598        usefulArray_[i].checkClean();
5599      else
5600        usefulArray_[i].CoinIndexedVector::checkClean();
5601    } else {
5602      if ((stateOfProblem_&check)!=0&&(ignoreEmpty&check)==0)
5603        printf("Array %d has no elements but says it is not available\n",i);
5604      if (usefulArray_[i].getNumPartitions())
5605        usefulArray_[i].checkClear();
5606      else
5607        usefulArray_[i].CoinIndexedVector::checkClear();
5608    }
5609  }
5610}
5611#include "CoinAbcCommon.hpp"
5612#if 0
5613#include "ClpFactorization.hpp"
5614// Loads tolerances etc
5615void
5616ClpSimplex::loadTolerancesEtc(const AbcTolerancesEtc & data)
5617{
5618  zeroTolerance_ = data.zeroTolerance_;
5619  primalToleranceToGetOptimal_ = data.primalToleranceToGetOptimal_;
5620  largeValue_ = data.largeValue_;
5621  alphaAccuracy_ = data.alphaAccuracy_;
5622  dualBound_ = data.dualBound_;
5623  dualTolerance_ = data.dualTolerance_;
5624  primalTolerance_ = data.primalTolerance_;
5625  infeasibilityCost_ = data.infeasibilityCost_;
5626  incomingInfeasibility_ = data.incomingInfeasibility_;
5627  allowedInfeasibility_ = data.allowedInfeasibility_;
5628  baseIteration_ = data.baseIteration_;
5629  numberRefinements_ = data.numberRefinements_;
5630  forceFactorization_ = data.forceFactorization_;
5631  perturbation_ = data.perturbation_;
5632  dontFactorizePivots_ = data.dontFactorizePivots_;
5633  /// For factorization
5634  abcFactorization_->maximumPivots(data.maximumPivots_);
5635}
5636// Loads tolerances etc
5637void
5638ClpSimplex::unloadTolerancesEtc(AbcTolerancesEtc & data)
5639{
5640  data.zeroTolerance_ = zeroTolerance_;
5641  data.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5642  data.largeValue_ = largeValue_;
5643  data.alphaAccuracy_ = alphaAccuracy_;
5644  data.dualBound_ = dualBound_;
5645  data.dualTolerance_ = dualTolerance_;
5646  data.primalTolerance_ = primalTolerance_;
5647  data.infeasibilityCost_ = infeasibilityCost_;
5648  data.incomingInfeasibility_ = incomingInfeasibility_;
5649  data.allowedInfeasibility_ = allowedInfeasibility_;
5650  data.baseIteration_ = baseIteration