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

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