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

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

fix some ampl stuff, add ClpSolver? and a few fixes

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