source: stable/1.15/Clp/src/AbcSimplex.cpp @ 1989

Last change on this file since 1989 was 1989, checked in by forrest, 6 years ago

Minor stability fixes and an option to allow perturbation after presolve

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