source: trunk/Clp/src/ClpSimplex.cpp @ 1060

Last change on this file since 1060 was 1060, checked in by forrest, 12 years ago

for secondary status on time

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