source: branches/devel/Clp/src/ClpSimplex.cpp @ 990

Last change on this file since 990 was 990, checked in by forrest, 14 years ago

add lpio and correct tightenbounds

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 285.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 {
1585    if (forceFactorization_>0&&
1586        factorization_->pivots()==forceFactorization_) {
1587      // relax
1588      forceFactorization_ = (3+5*forceFactorization_)/4;
1589      if (forceFactorization_>factorization_->maximumPivots())
1590        forceFactorization_ = -1; //off
1591      return 1;
1592    } else {
1593      // carry on iterating
1594      return 0;
1595    }
1596  }
1597}
1598// Copy constructor.
1599ClpSimplex::ClpSimplex(const ClpSimplex &rhs,int scalingMode) :
1600  ClpModel(rhs,scalingMode),
1601  columnPrimalInfeasibility_(0.0),
1602  rowPrimalInfeasibility_(0.0),
1603  columnPrimalSequence_(-2),
1604  rowPrimalSequence_(-2), 
1605  columnDualInfeasibility_(0.0),
1606  rowDualInfeasibility_(0.0),
1607  columnDualSequence_(-2),
1608  rowDualSequence_(-2),
1609  primalToleranceToGetOptimal_(-1.0),
1610  remainingDualInfeasibility_(0.0),
1611  largeValue_(1.0e15),
1612  largestPrimalError_(0.0),
1613  largestDualError_(0.0),
1614  alphaAccuracy_(-1.0),
1615  dualBound_(1.0e10),
1616  alpha_(0.0),
1617  theta_(0.0),
1618  lowerIn_(0.0),
1619  valueIn_(0.0),
1620  upperIn_(-COIN_DBL_MAX),
1621  dualIn_(0.0),
1622  lowerOut_(-1),
1623  valueOut_(-1),
1624  upperOut_(-1),
1625  dualOut_(-1),
1626  dualTolerance_(0.0),
1627  primalTolerance_(0.0),
1628  sumDualInfeasibilities_(0.0),
1629  sumPrimalInfeasibilities_(0.0),
1630  infeasibilityCost_(1.0e10),
1631  sumOfRelaxedDualInfeasibilities_(0.0),
1632  sumOfRelaxedPrimalInfeasibilities_(0.0),
1633  acceptablePivot_(1.0e-8),
1634  lower_(NULL),
1635  rowLowerWork_(NULL),
1636  columnLowerWork_(NULL),
1637  upper_(NULL),
1638  rowUpperWork_(NULL),
1639  columnUpperWork_(NULL),
1640  cost_(NULL),
1641  rowObjectiveWork_(NULL),
1642  objectiveWork_(NULL),
1643  sequenceIn_(-1),
1644  directionIn_(-1),
1645  sequenceOut_(-1),
1646  directionOut_(-1),
1647  pivotRow_(-1),
1648  lastGoodIteration_(-100),
1649  dj_(NULL),
1650  rowReducedCost_(NULL),
1651  reducedCostWork_(NULL),
1652  solution_(NULL),
1653  rowActivityWork_(NULL),
1654  columnActivityWork_(NULL),
1655  auxiliaryModel_(NULL),
1656  numberDualInfeasibilities_(0),
1657  numberDualInfeasibilitiesWithoutFree_(0),
1658  numberPrimalInfeasibilities_(100),
1659  numberRefinements_(0),
1660  pivotVariable_(NULL),
1661  factorization_(NULL),
1662  savedSolution_(NULL),
1663  numberTimesOptimal_(0),
1664  disasterArea_(NULL),
1665  changeMade_(1),
1666  algorithm_(0),
1667  forceFactorization_(-1),
1668  perturbation_(100),
1669  nonLinearCost_(NULL),
1670  lastBadIteration_(-999999),
1671  lastFlaggedIteration_(-999999),
1672  numberFake_(0),
1673  numberChanged_(0),
1674  progressFlag_(0),
1675  firstFree_(-1),
1676  numberExtraRows_(0),
1677  maximumBasic_(0),
1678  incomingInfeasibility_(1.0),
1679  allowedInfeasibility_(10.0),
1680  automaticScale_(0),
1681  progress_(NULL)
1682{
1683  int i;
1684  for (i=0;i<6;i++) {
1685    rowArray_[i]=NULL;
1686    columnArray_[i]=NULL;
1687  }
1688  for (i=0;i<4;i++) {
1689    spareIntArray_[i]=0;
1690    spareDoubleArray_[i]=0.0;
1691  }
1692  saveStatus_=NULL;
1693  factorization_ = NULL;
1694  dualRowPivot_ = NULL;
1695  primalColumnPivot_ = NULL;
1696  gutsOfDelete(0);
1697  delete nonLinearCost_;
1698  nonLinearCost_ = NULL;
1699  gutsOfCopy(rhs);
1700  solveType_=1; // say simplex based life form
1701}
1702// Copy constructor from model
1703ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
1704  ClpModel(rhs,scalingMode),
1705  columnPrimalInfeasibility_(0.0),
1706  rowPrimalInfeasibility_(0.0),
1707  columnPrimalSequence_(-2),
1708  rowPrimalSequence_(-2), 
1709  columnDualInfeasibility_(0.0),
1710  rowDualInfeasibility_(0.0),
1711  columnDualSequence_(-2),
1712  rowDualSequence_(-2),
1713  primalToleranceToGetOptimal_(-1.0),
1714  remainingDualInfeasibility_(0.0),
1715  largeValue_(1.0e15),
1716  largestPrimalError_(0.0),
1717  largestDualError_(0.0),
1718  alphaAccuracy_(-1.0),
1719  dualBound_(1.0e10),
1720  alpha_(0.0),
1721  theta_(0.0),
1722  lowerIn_(0.0),
1723  valueIn_(0.0),
1724  upperIn_(-COIN_DBL_MAX),
1725  dualIn_(0.0),
1726  lowerOut_(-1),
1727  valueOut_(-1),
1728  upperOut_(-1),
1729  dualOut_(-1),
1730  dualTolerance_(0.0),
1731  primalTolerance_(0.0),
1732  sumDualInfeasibilities_(0.0),
1733  sumPrimalInfeasibilities_(0.0),
1734  infeasibilityCost_(1.0e10),
1735  sumOfRelaxedDualInfeasibilities_(0.0),
1736  sumOfRelaxedPrimalInfeasibilities_(0.0),
1737  acceptablePivot_(1.0e-8),
1738  lower_(NULL),
1739  rowLowerWork_(NULL),
1740  columnLowerWork_(NULL),
1741  upper_(NULL),
1742  rowUpperWork_(NULL),
1743  columnUpperWork_(NULL),
1744  cost_(NULL),
1745  rowObjectiveWork_(NULL),
1746  objectiveWork_(NULL),
1747  sequenceIn_(-1),
1748  directionIn_(-1),
1749  sequenceOut_(-1),
1750  directionOut_(-1),
1751  pivotRow_(-1),
1752  lastGoodIteration_(-100),
1753  dj_(NULL),
1754  rowReducedCost_(NULL),
1755  reducedCostWork_(NULL),
1756  solution_(NULL),
1757  rowActivityWork_(NULL),
1758  columnActivityWork_(NULL),
1759  auxiliaryModel_(NULL),
1760  numberDualInfeasibilities_(0),
1761  numberDualInfeasibilitiesWithoutFree_(0),
1762  numberPrimalInfeasibilities_(100),
1763  numberRefinements_(0),
1764  pivotVariable_(NULL),
1765  factorization_(NULL),
1766  savedSolution_(NULL),
1767  numberTimesOptimal_(0),
1768  disasterArea_(NULL),
1769  changeMade_(1),
1770  algorithm_(0),
1771  forceFactorization_(-1),
1772  perturbation_(100),
1773  nonLinearCost_(NULL),
1774  lastBadIteration_(-999999),
1775  lastFlaggedIteration_(-999999),
1776  numberFake_(0),
1777  numberChanged_(0),
1778  progressFlag_(0),
1779  firstFree_(-1),
1780  numberExtraRows_(0),
1781  maximumBasic_(0),
1782  incomingInfeasibility_(1.0),
1783  allowedInfeasibility_(10.0),
1784  automaticScale_(0),
1785  progress_(NULL)
1786{
1787  int i;
1788  for (i=0;i<6;i++) {
1789    rowArray_[i]=NULL;
1790    columnArray_[i]=NULL;
1791  }
1792  for (i=0;i<4;i++) {
1793    spareIntArray_[i]=0;
1794    spareDoubleArray_[i]=0.0;
1795  }
1796  saveStatus_=NULL;
1797  // get an empty factorization so we can set tolerances etc
1798  getEmptyFactorization();
1799  // say Steepest pricing
1800  dualRowPivot_ = new ClpDualRowSteepest();
1801  // say Steepest pricing
1802  primalColumnPivot_ = new ClpPrimalColumnSteepest();
1803  solveType_=1; // say simplex based life form
1804 
1805}
1806// Assignment operator. This copies the data
1807ClpSimplex & 
1808ClpSimplex::operator=(const ClpSimplex & rhs)
1809{
1810  if (this != &rhs) {
1811    gutsOfDelete(0);
1812    delete nonLinearCost_;
1813    nonLinearCost_ = NULL;
1814    ClpModel::operator=(rhs);
1815    gutsOfCopy(rhs);
1816  }
1817  return *this;
1818}
1819void 
1820ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1821{
1822  assert (numberRows_==rhs.numberRows_);
1823  assert (numberColumns_==rhs.numberColumns_);
1824  numberExtraRows_ = rhs.numberExtraRows_;
1825  maximumBasic_ = rhs.maximumBasic_;
1826  int numberRows2 = numberRows_+numberExtraRows_;
1827  auxiliaryModel_ = NULL;
1828  if ((whatsChanged_&1)!=0) {
1829    lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows2);
1830    rowLowerWork_ = lower_+numberColumns_;
1831    columnLowerWork_ = lower_;
1832    upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows2);
1833    rowUpperWork_ = upper_+numberColumns_;
1834    columnUpperWork_ = upper_;
1835    //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
1836    cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows2));
1837    objectiveWork_ = cost_;
1838    rowObjectiveWork_ = cost_+numberColumns_;
1839    dj_ = ClpCopyOfArray(rhs.dj_,numberRows2+numberColumns_);
1840    if (dj_) {
1841      reducedCostWork_ = dj_;
1842      rowReducedCost_ = dj_+numberColumns_;
1843    }
1844    solution_ = ClpCopyOfArray(rhs.solution_,numberRows2+numberColumns_);
1845    if (solution_) {
1846      columnActivityWork_ = solution_;
1847      rowActivityWork_ = solution_+numberColumns_;
1848    }
1849    if (rhs.pivotVariable_) {
1850      pivotVariable_ = new int[numberRows2];
1851      CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
1852    } else {
1853      pivotVariable_=NULL;
1854    }
1855    savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1856    int i;
1857    for (i=0;i<6;i++) {
1858      rowArray_[i]=NULL;
1859      if (rhs.rowArray_[i]) 
1860        rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
1861      columnArray_[i]=NULL;
1862      if (rhs.columnArray_[i]) 
1863        columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
1864    }
1865    if (rhs.saveStatus_) {
1866      saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1867    }
1868  } else {
1869    lower_ = NULL;
1870    rowLowerWork_ = NULL;
1871    columnLowerWork_ = NULL; 
1872    upper_ = NULL; 
1873    rowUpperWork_ = NULL; 
1874    columnUpperWork_ = NULL;
1875    cost_ = NULL; 
1876    objectiveWork_ = NULL; 
1877    rowObjectiveWork_ = NULL; 
1878    dj_ = NULL; 
1879    reducedCostWork_ = NULL;
1880    rowReducedCost_ = NULL;
1881    solution_ = NULL; 
1882    columnActivityWork_ = NULL; 
1883    rowActivityWork_ = NULL; 
1884    pivotVariable_=NULL;
1885    savedSolution_ = NULL;
1886    int i;
1887    for (i=0;i<6;i++) {
1888      rowArray_[i]=NULL;
1889      columnArray_[i]=NULL;
1890    }
1891    saveStatus_ = NULL;
1892  }
1893  if (rhs.factorization_) {
1894    delete factorization_;
1895    factorization_ = new ClpFactorization(*rhs.factorization_);
1896  } else {
1897    factorization_=NULL;
1898  }
1899  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1900  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1901  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1902  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1903  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1904  columnDualSequence_ = rhs.columnDualSequence_;
1905  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1906  rowDualSequence_ = rhs.rowDualSequence_;
1907  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1908  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1909  largeValue_ = rhs.largeValue_;
1910  largestPrimalError_ = rhs.largestPrimalError_;
1911  largestDualError_ = rhs.largestDualError_;
1912  alphaAccuracy_ = rhs.alphaAccuracy_;
1913  dualBound_ = rhs.dualBound_;
1914  alpha_ = rhs.alpha_;
1915  theta_ = rhs.theta_;
1916  lowerIn_ = rhs.lowerIn_;
1917  valueIn_ = rhs.valueIn_;
1918  upperIn_ = rhs.upperIn_;
1919  dualIn_ = rhs.dualIn_;
1920  sequenceIn_ = rhs.sequenceIn_;
1921  directionIn_ = rhs.directionIn_;
1922  lowerOut_ = rhs.lowerOut_;
1923  valueOut_ = rhs.valueOut_;
1924  upperOut_ = rhs.upperOut_;
1925  dualOut_ = rhs.dualOut_;
1926  sequenceOut_ = rhs.sequenceOut_;
1927  directionOut_ = rhs.directionOut_;
1928  pivotRow_ = rhs.pivotRow_;
1929  lastGoodIteration_ = rhs.lastGoodIteration_;
1930  numberRefinements_ = rhs.numberRefinements_;
1931  dualTolerance_ = rhs.dualTolerance_;
1932  primalTolerance_ = rhs.primalTolerance_;
1933  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1934  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1935  numberDualInfeasibilitiesWithoutFree_ = 
1936    rhs.numberDualInfeasibilitiesWithoutFree_;
1937  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1938  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1939  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1940  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1941  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1942  disasterArea_ = NULL;
1943  changeMade_ = rhs.changeMade_;
1944  algorithm_ = rhs.algorithm_;
1945  forceFactorization_ = rhs.forceFactorization_;
1946  perturbation_ = rhs.perturbation_;
1947  infeasibilityCost_ = rhs.infeasibilityCost_;
1948  lastBadIteration_ = rhs.lastBadIteration_;
1949  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
1950  numberFake_ = rhs.numberFake_;
1951  numberChanged_ = rhs.numberChanged_;
1952  progressFlag_ = rhs.progressFlag_;
1953  firstFree_ = rhs.firstFree_;
1954  incomingInfeasibility_ = rhs.incomingInfeasibility_;
1955  allowedInfeasibility_ = rhs.allowedInfeasibility_;
1956  automaticScale_ = rhs.automaticScale_;
1957  if (rhs.progress_)
1958    progress_ = new ClpSimplexProgress(*rhs.progress_);
1959  else
1960    progress_=NULL;
1961  for (int i=0;i<4;i++) {
1962    spareIntArray_[i]=rhs.spareIntArray_[i];
1963    spareDoubleArray_[i]=rhs.spareDoubleArray_[i];
1964  }
1965  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
1966  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
1967  acceptablePivot_ = rhs.acceptablePivot_;
1968  if (rhs.nonLinearCost_!=NULL)
1969    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1970  else
1971    nonLinearCost_=NULL;
1972  solveType_=rhs.solveType_;
1973}
1974// type == 0 do everything, most + pivot data, 2 factorization data as well
1975void 
1976ClpSimplex::gutsOfDelete(int type)
1977{
1978  delete [] lower_;
1979  lower_=NULL;
1980  rowLowerWork_=NULL;
1981  columnLowerWork_=NULL;
1982  delete [] upper_;
1983  upper_=NULL;
1984  rowUpperWork_=NULL;
1985  columnUpperWork_=NULL;
1986  delete [] cost_;
1987  cost_=NULL;
1988  objectiveWork_=NULL;
1989  rowObjectiveWork_=NULL;
1990  delete [] dj_;
1991  dj_=NULL;
1992  reducedCostWork_=NULL;
1993  rowReducedCost_=NULL;
1994  delete [] solution_;
1995  solution_=NULL;
1996  rowActivityWork_=NULL;
1997  columnActivityWork_=NULL;
1998  delete [] savedSolution_;
1999  savedSolution_ = NULL;
2000  if ((specialOptions_&2)==0) {
2001    delete nonLinearCost_;
2002    nonLinearCost_ = NULL;
2003  }
2004  int i;
2005  if ((specialOptions_&65536)==0) {
2006    for (i=0;i<6;i++) {
2007      delete rowArray_[i];
2008      rowArray_[i]=NULL;
2009      delete columnArray_[i];
2010      columnArray_[i]=NULL;
2011    }
2012  }
2013  delete rowCopy_;
2014  rowCopy_=NULL;
2015  delete [] saveStatus_;
2016  saveStatus_=NULL;
2017  if (!type) {
2018    // delete everything
2019    delete auxiliaryModel_;
2020    auxiliaryModel_ = NULL;
2021    setEmptyFactorization();
2022    delete [] pivotVariable_;
2023    pivotVariable_=NULL;
2024    delete dualRowPivot_;
2025    dualRowPivot_ = NULL;
2026    delete primalColumnPivot_;
2027    primalColumnPivot_ = NULL;
2028    delete progress_;
2029    progress_=NULL;
2030  } else {
2031    // delete any size information in methods
2032    if (type>1) {
2033      factorization_->clearArrays();
2034      delete [] pivotVariable_;
2035      pivotVariable_=NULL;
2036    }
2037    dualRowPivot_->clearArrays();
2038    primalColumnPivot_->clearArrays();
2039  }
2040}
2041// This sets largest infeasibility and most infeasible
2042void 
2043ClpSimplex::checkPrimalSolution(const double * rowActivities,
2044                                        const double * columnActivities)
2045{
2046  double * solution;
2047  int iRow,iColumn;
2048
2049  objectiveValue_ = 0.0;
2050  // now look at primal solution
2051  solution = rowActivityWork_;
2052  sumPrimalInfeasibilities_=0.0;
2053  numberPrimalInfeasibilities_=0;
2054  double primalTolerance = primalTolerance_;
2055  double relaxedTolerance=primalTolerance_;
2056  // we can't really trust infeasibilities if there is primal error
2057  double error = CoinMin(1.0e-2,largestPrimalError_);
2058  // allow tolerance at least slightly bigger than standard
2059  relaxedTolerance = relaxedTolerance +  error;
2060  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2061  for (iRow=0;iRow<numberRows_;iRow++) {
2062    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
2063    double infeasibility=0.0;
2064    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
2065    if (solution[iRow]>rowUpperWork_[iRow]) {
2066      infeasibility=solution[iRow]-rowUpperWork_[iRow];
2067    } else if (solution[iRow]<rowLowerWork_[iRow]) {
2068      infeasibility=rowLowerWork_[iRow]-solution[iRow];
2069    }
2070    if (infeasibility>primalTolerance) {
2071      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2072      if (infeasibility>relaxedTolerance) 
2073        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2074      numberPrimalInfeasibilities_ ++;
2075    }
2076    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
2077  }
2078  // Check any infeasibilities from dynamic rows
2079  matrix_->primalExpanded(this,2);
2080  solution = columnActivityWork_;
2081  if (!matrix_->rhsOffset(this)) {
2082    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2083      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2084      double infeasibility=0.0;
2085      objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
2086      if (solution[iColumn]>columnUpperWork_[iColumn]) {
2087        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
2088      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
2089        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
2090      }
2091      if (infeasibility>primalTolerance) {
2092        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2093        if (infeasibility>relaxedTolerance) 
2094          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2095        numberPrimalInfeasibilities_ ++;
2096      }
2097      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
2098    }
2099  } else {
2100    // as we are using effective rhs we only check basics
2101    // But we do need to get objective
2102    objectiveValue_ += innerProduct(objectiveWork_,numberColumns_,solution);
2103    for (int j=0;j<numberRows_;j++) {
2104      int iColumn = pivotVariable_[j];
2105      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2106      double infeasibility=0.0;
2107      if (solution[iColumn]>columnUpperWork_[iColumn]) {
2108        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
2109      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
2110        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
2111      }
2112      if (infeasibility>primalTolerance) {
2113        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2114        if (infeasibility>relaxedTolerance) 
2115          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2116        numberPrimalInfeasibilities_ ++;
2117      }
2118      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
2119    }
2120  }
2121  objectiveValue_ += objective_->nonlinearOffset();
2122  objectiveValue_ /= (objectiveScale_*rhsScale_);
2123}
2124void 
2125ClpSimplex::checkDualSolution()
2126{
2127
2128  int iRow,iColumn;
2129  sumDualInfeasibilities_=0.0;
2130  numberDualInfeasibilities_=0;
2131  numberDualInfeasibilitiesWithoutFree_=0;
2132  if (matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2) {
2133    // pretend we found dual infeasibilities
2134    sumOfRelaxedDualInfeasibilities_ = 1.0;
2135    sumDualInfeasibilities_=1.0;
2136    numberDualInfeasibilities_=1;
2137    return;
2138  }
2139  int firstFreePrimal = -1;
2140  int firstFreeDual = -1;
2141  int numberSuperBasicWithDj=0;
2142  double relaxedTolerance=dualTolerance_;
2143  // we can't really trust infeasibilities if there is dual error
2144  double error = CoinMin(1.0e-2,largestDualError_);
2145  // allow tolerance at least slightly bigger than standard
2146  relaxedTolerance = relaxedTolerance +  error;
2147  sumOfRelaxedDualInfeasibilities_ = 0.0;
2148
2149  // Check any djs from dynamic rows
2150  matrix_->dualExpanded(this,NULL,NULL,3);
2151  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
2152  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2153    if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
2154      // not basic
2155      double distanceUp = columnUpperWork_[iColumn]-
2156        columnActivityWork_[iColumn];
2157      double distanceDown = columnActivityWork_[iColumn] -
2158        columnLowerWork_[iColumn];
2159      if (distanceUp>primalTolerance_) {
2160        double value = reducedCostWork_[iColumn];
2161        // Check if "free"
2162        if (distanceDown>primalTolerance_) {
2163          if (fabs(value)>1.0e2*relaxedTolerance) {
2164            numberSuperBasicWithDj++;
2165            if (firstFreeDual<0)
2166              firstFreeDual = iColumn;
2167          }
2168          if (firstFreePrimal<0)
2169            firstFreePrimal = iColumn;
2170        }
2171        // should not be negative
2172        if (value<0.0) {
2173          value = - value;
2174          if (value>dualTolerance_) {
2175            if (getColumnStatus(iColumn) != isFree) {
2176              numberDualInfeasibilitiesWithoutFree_ ++;
2177              sumDualInfeasibilities_ += value-dualTolerance_;
2178              if (value>relaxedTolerance) 
2179                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2180              numberDualInfeasibilities_ ++;
2181            } else {
2182              // free so relax a lot
2183              value *= 0.01;
2184              if (value>dualTolerance_) {
2185                sumDualInfeasibilities_ += value-dualTolerance_;
2186                if (value>relaxedTolerance) 
2187                  sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2188                numberDualInfeasibilities_ ++;
2189              }
2190            }
2191          }
2192        }
2193      }
2194      if (distanceDown>primalTolerance_) {
2195        double value = reducedCostWork_[iColumn];
2196        // should not be positive
2197        if (value>0.0) {
2198          if (value>dualTolerance_) {
2199            sumDualInfeasibilities_ += value-dualTolerance_;
2200            if (value>relaxedTolerance) 
2201              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2202            numberDualInfeasibilities_ ++;
2203            if (getColumnStatus(iColumn) != isFree) 
2204              numberDualInfeasibilitiesWithoutFree_ ++;
2205            // maybe we can make feasible by increasing tolerance
2206          }
2207        }
2208      }
2209    }
2210  }
2211  for (iRow=0;iRow<numberRows_;iRow++) {
2212    if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
2213      // not basic
2214      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
2215      double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
2216      if (distanceUp>primalTolerance_) {
2217        double value = rowReducedCost_[iRow];
2218        // Check if "free"
2219        if (distanceDown>primalTolerance_) {
2220          if (fabs(value)>1.0e2*relaxedTolerance) {
2221            numberSuperBasicWithDj++;
2222            if (firstFreeDual<0)
2223              firstFreeDual = iRow+numberColumns_;
2224          }
2225          if (firstFreePrimal<0)
2226            firstFreePrimal = iRow+numberColumns_;
2227        }
2228        // should not be negative
2229        if (value<0.0) {
2230          value = - value;
2231          if (value>dualTolerance_) {
2232            sumDualInfeasibilities_ += value-dualTolerance_;
2233            if (value>relaxedTolerance) 
2234              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2235            numberDualInfeasibilities_ ++;
2236            if (getRowStatus(iRow) != isFree) 
2237              numberDualInfeasibilitiesWithoutFree_ ++;
2238          }
2239        }
2240      }
2241      if (distanceDown>primalTolerance_) {
2242        double value = rowReducedCost_[iRow];
2243        // should not be positive
2244        if (value>0.0) {
2245          if (value>dualTolerance_) {
2246            sumDualInfeasibilities_ += value-dualTolerance_;
2247            if (value>relaxedTolerance) 
2248              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2249            numberDualInfeasibilities_ ++;
2250            if (getRowStatus(iRow) != isFree) 
2251              numberDualInfeasibilitiesWithoutFree_ ++;
2252            // maybe we can make feasible by increasing tolerance
2253          }
2254        }
2255      }
2256    }
2257  }
2258  if (algorithm_<0&&firstFreeDual>=0) {
2259    // dual
2260    firstFree_ = firstFreeDual;
2261  } else if (numberSuperBasicWithDj||
2262             (progress_&&progress_->lastIterationNumber(0)<=0)) {
2263    firstFree_=firstFreePrimal;
2264  }
2265}
2266/* This sets sum and number of infeasibilities (Dual and Primal) */
2267void 
2268ClpSimplex::checkBothSolutions()
2269{
2270  if ((matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2)||
2271      matrix_->rhsOffset(this)) {
2272    // old way
2273    checkPrimalSolution(rowActivityWork_,columnActivityWork_);
2274    checkDualSolution();
2275    return;
2276  }
2277  int iSequence;
2278
2279  objectiveValue_ = 0.0;
2280  sumPrimalInfeasibilities_=0.0;
2281  numberPrimalInfeasibilities_=0;
2282  double primalTolerance = primalTolerance_;
2283  double relaxedToleranceP=primalTolerance_;
2284  // we can't really trust infeasibilities if there is primal error
2285  double error = CoinMin(1.0e-2,largestPrimalError_);
2286  // allow tolerance at least slightly bigger than standard
2287  relaxedToleranceP = relaxedToleranceP +  error;
2288  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2289  sumDualInfeasibilities_=0.0;
2290  numberDualInfeasibilities_=0;
2291  double dualTolerance=dualTolerance_;
2292  double relaxedToleranceD=dualTolerance;
2293  // we can't really trust infeasibilities if there is dual error
2294  error = CoinMin(1.0e-2,largestDualError_);
2295  // allow tolerance at least slightly bigger than standard
2296  relaxedToleranceD = relaxedToleranceD +  error;
2297  sumOfRelaxedDualInfeasibilities_ = 0.0;
2298
2299  // Check any infeasibilities from dynamic rows
2300  matrix_->primalExpanded(this,2);
2301  // Check any djs from dynamic rows
2302  matrix_->dualExpanded(this,NULL,NULL,3);
2303  int numberDualInfeasibilitiesFree= 0;
2304  int firstFreePrimal = -1;
2305  int firstFreeDual = -1;
2306  int numberSuperBasicWithDj=0;
2307
2308  int numberTotal = numberRows_ + numberColumns_;
2309  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2310    double value = solution_[iSequence];
2311#ifdef COIN_DEBUG
2312    if (fabs(value)>1.0e20)
2313      printf("%d values %g %g %g - status %d\n",iSequence,lower_[iSequence],
2314             solution_[iSequence],upper_[iSequence],status_[iSequence]);
2315#endif
2316    objectiveValue_ += value*cost_[iSequence];
2317    double distanceUp =upper_[iSequence]-value;
2318    double distanceDown =value-lower_[iSequence];
2319    if (distanceUp<-primalTolerance) {
2320      double infeasibility=-distanceUp;
2321      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2322      if (infeasibility>relaxedToleranceP) 
2323        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedToleranceP;
2324      numberPrimalInfeasibilities_ ++;
2325    } else if (distanceDown<-primalTolerance) {
2326      double infeasibility=-distanceDown;
2327      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2328      if (infeasibility>relaxedToleranceP) 
2329        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedToleranceP;
2330      numberPrimalInfeasibilities_ ++;
2331    } else {
2332      // feasible (so could be free)
2333      if (getStatus(iSequence) != basic&&!flagged(iSequence)) {
2334        // not basic
2335        double djValue = dj_[iSequence];
2336        if (distanceDown<primalTolerance) {
2337          if (distanceUp>primalTolerance&&djValue<-dualTolerance) {
2338            sumDualInfeasibilities_ -= djValue+dualTolerance;
2339            if (djValue<-relaxedToleranceD) 
2340              sumOfRelaxedDualInfeasibilities_ -= djValue+relaxedToleranceD;
2341            numberDualInfeasibilities_ ++;
2342          } 
2343        } else if (distanceUp<primalTolerance) {
2344          if (djValue>dualTolerance) {
2345            sumDualInfeasibilities_ += djValue-dualTolerance;
2346            if (djValue>relaxedToleranceD) 
2347              sumOfRelaxedDualInfeasibilities_ += djValue-relaxedToleranceD;
2348            numberDualInfeasibilities_ ++;
2349          }
2350        } else {
2351          // may be free
2352          djValue *= 0.01;
2353          if (fabs(djValue)>dualTolerance) {
2354            if (getStatus(iSequence) == isFree) 
2355              numberDualInfeasibilitiesFree++;
2356            sumDualInfeasibilities_ += fabs(djValue)-dualTolerance;
2357            numberDualInfeasibilities_ ++;
2358            if (fabs(djValue)>relaxedToleranceD) {
2359              sumOfRelaxedDualInfeasibilities_ += value-relaxedToleranceD;
2360              numberSuperBasicWithDj++;
2361              if (firstFreeDual<0)
2362                firstFreeDual = iSequence;
2363            }
2364          }
2365          if (firstFreePrimal<0)
2366            firstFreePrimal = iSequence;
2367        }
2368      }
2369    }
2370  }
2371  objectiveValue_ += objective_->nonlinearOffset();
2372  objectiveValue_ /= (objectiveScale_*rhsScale_);
2373  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_-
2374    numberDualInfeasibilitiesFree;
2375  if (algorithm_<0&&firstFreeDual>=0) {
2376    // dual
2377    firstFree_ = firstFreeDual;
2378  } else if (numberSuperBasicWithDj||
2379             (progress_&&progress_->lastIterationNumber(0)<=0)) {
2380    firstFree_=firstFreePrimal;
2381  }
2382}
2383/* Adds multiple of a column into an array */
2384void 
2385ClpSimplex::add(double * array,
2386                int sequence, double multiplier) const
2387{
2388  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2389    //slack
2390    array [sequence-numberColumns_] -= multiplier;
2391  } else {
2392    // column
2393    matrix_->add(this,array,sequence,multiplier);
2394  }
2395}
2396/*
2397  Unpacks one column of the matrix into indexed array
2398*/
2399void 
2400ClpSimplex::unpack(CoinIndexedVector * rowArray) const
2401{
2402  rowArray->clear();
2403  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2404    //slack
2405    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
2406  } else {
2407    // column
2408    matrix_->unpack(this,rowArray,sequenceIn_);
2409  }
2410}
2411void 
2412ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence) const
2413{
2414  rowArray->clear();
2415  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2416    //slack
2417    rowArray->insert(sequence-numberColumns_,-1.0);
2418  } else {
2419    // column
2420    matrix_->unpack(this,rowArray,sequence);
2421  }
2422}
2423/*
2424  Unpacks one column of the matrix into indexed array
2425*/
2426void 
2427ClpSimplex::unpackPacked(CoinIndexedVector * rowArray) 
2428{
2429  rowArray->clear();
2430  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2431    //slack
2432    int * index = rowArray->getIndices();
2433    double * array = rowArray->denseVector();
2434    array[0]=-1.0;
2435    index[0]=sequenceIn_-numberColumns_;
2436    rowArray->setNumElements(1);
2437    rowArray->setPackedMode(true);
2438  } else {
2439    // column
2440    matrix_->unpackPacked(this,rowArray,sequenceIn_);
2441  }
2442}
2443void 
2444ClpSimplex::unpackPacked(CoinIndexedVector * rowArray,int sequence)
2445{
2446  rowArray->clear();
2447  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2448    //slack
2449    int * index = rowArray->getIndices();
2450    double * array = rowArray->denseVector();
2451    array[0]=-1.0;
2452    index[0]=sequence-numberColumns_;
2453    rowArray->setNumElements(1);
2454    rowArray->setPackedMode(true);
2455  } else {
2456    // column
2457    matrix_->unpackPacked(this,rowArray,sequence);
2458  }
2459}
2460//static int scale_times[]={0,0,0,0};
2461bool
2462ClpSimplex::createRim(int what,bool makeRowCopy, int startFinishOptions)
2463{
2464  bool goodMatrix=true;
2465  spareIntArray_[0]=0;
2466  if (!matrix_->canGetRowCopy())
2467    makeRowCopy=false; // switch off row copy if can't produce
2468  int saveLevel=handler_->logLevel();
2469  // Arrays will be there and correct size unless what is 63
2470  bool newArrays = (what==63);
2471  // We may be restarting with same size
2472  bool keepPivots=false;
2473  if (startFinishOptions==-1) {
2474    startFinishOptions=0;
2475    keepPivots=true;
2476  }
2477  bool oldMatrix = ((startFinishOptions&4)!=0&&(whatsChanged_&1)!=0);
2478  if (auxiliaryModel_) {
2479    if (auxiliaryModel_->numberRows_==numberRows_&&
2480        auxiliaryModel_->numberColumns_==numberColumns_&&
2481        (whatsChanged_&511)==511) {
2482      oldMatrix=true; 
2483    } else {
2484      deleteAuxiliaryModel();
2485      oldMatrix=false;
2486    }
2487  }
2488  if (what==63) {
2489    if (!status_)
2490      createStatus();
2491    if (oldMatrix)
2492      newArrays=false;
2493    if (problemStatus_==10) {
2494      handler_->setLogLevel(0); // switch off messages
2495      if (rowArray_[0]) {
2496        // stuff is still there
2497        oldMatrix=true;
2498        newArrays=false;
2499        keepPivots=true;
2500        for (int iRow=0;iRow<4;iRow++) {
2501          rowArray_[iRow]->clear();
2502        }
2503        for (int iColumn=0;iColumn<2;iColumn++) {
2504          columnArray_[iColumn]->clear();
2505        }
2506      }
2507    } else if (factorization_) {
2508      // match up factorization messages
2509      if (handler_->logLevel()<3)
2510        factorization_->messageLevel(0);
2511      else
2512        factorization_->messageLevel(CoinMax(3,factorization_->messageLevel()));
2513      /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
2514         i.e. oldMatrix false but okay as long as same number rows and status array exists
2515      */
2516      if ((startFinishOptions&2)!=0&&factorization_->numberRows()==numberRows_&&status_)
2517        keepPivots=true;
2518    }
2519    numberExtraRows_ = matrix_->generalExpanded(this,2,maximumBasic_);
2520    if (numberExtraRows_&&newArrays) {
2521      // make sure status array large enough
2522      assert (status_);
2523      int numberOld = numberRows_+numberColumns_;
2524      int numberNew = numberRows_+numberColumns_+numberExtraRows_;
2525      unsigned char * newStatus = new unsigned char [numberNew];
2526      memset(newStatus+numberOld,0,numberExtraRows_);
2527      memcpy(newStatus,status_,numberOld);
2528      delete [] status_;
2529      status_=newStatus;
2530    }
2531  }
2532  int numberRows2 = numberRows_+numberExtraRows_;
2533  int numberTotal = numberRows2+numberColumns_;
2534  int i;
2535  bool doSanityCheck=true;
2536  if (what==63&&!auxiliaryModel_) {
2537    // We may want to switch stuff off for speed
2538    if ((specialOptions_&256)!=0)
2539      makeRowCopy=false; // no row copy
2540    if ((specialOptions_&128)!=0)
2541      doSanityCheck=false; // no sanity check
2542    //check matrix
2543    if (!matrix_)
2544      matrix_=new ClpPackedMatrix();
2545    int checkType=(doSanityCheck) ? 15 : 14;
2546    if (oldMatrix)
2547      checkType = 14;
2548    if ((specialOptions_&0x1000000)!=0)
2549      checkType -= 4; // don't check for duplicates
2550    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20,checkType)) {
2551      problemStatus_=4;
2552      secondaryStatus_=8;
2553      //goodMatrix= false;
2554      return false;
2555    }
2556    if (makeRowCopy&&!oldMatrix) {
2557      delete rowCopy_;
2558      // may return NULL if can't give row copy
2559      rowCopy_ = matrix_->reverseOrderedCopy();
2560    }
2561#if 0
2562    if (what==63&&(specialOptions_&131072)!=0) {
2563      int k=rowScale_ ? 1 : 0;
2564      if (oldMatrix)
2565        k+=2;
2566      scale_times[k]++;
2567      if ((scale_times[0]+scale_times[1]+scale_times[2]+scale_times[3])%1000==0)
2568        printf("scale counts %d %d %d %d\n",
2569               scale_times[0],scale_times[1],scale_times[2],scale_times[3]);
2570    }
2571#endif
2572    // do scaling if needed
2573    if (!oldMatrix&&scalingFlag_<0) {
2574      if (scalingFlag_<0&&rowScale_) {
2575        //if (handler_->logLevel()>0)
2576          printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
2577                 scalingFlag_);
2578        scalingFlag_=0;
2579      }
2580      delete [] rowScale_;
2581      delete [] columnScale_;
2582      rowScale_=NULL;
2583      columnScale_=NULL;
2584    }
2585    if (scalingFlag_>0&&!rowScale_) {
2586      if (matrix_->scale(this))
2587        scalingFlag_=-scalingFlag_; // not scaled after all
2588      if (rowScale_&&automaticScale_) {
2589        // try automatic scaling
2590        double smallestObj=1.0e100;
2591        double largestObj=0.0;
2592        double largestRhs=0.0;
2593        const double * obj = objective();
2594        for (i=0;i<numberColumns_;i++) {
2595          double value = fabs(obj[i]);
2596          value *= columnScale_[i];
2597          if (value&&columnLower_[i]!=columnUpper_[i]) {
2598            smallestObj = CoinMin(smallestObj,value);
2599            largestObj = CoinMax(largestObj,value);
2600          }
2601          if (columnLower_[i]>0.0||columnUpper_[i]<0.0) {
2602            double scale = 1.0/columnScale_[i];
2603            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2604            if (columnLower_[i]>0)
2605              largestRhs=CoinMax(largestRhs,columnLower_[i]*scale);
2606            if (columnUpper_[i]<0.0)
2607              largestRhs=CoinMax(largestRhs,-columnUpper_[i]*scale);
2608          }
2609        }
2610        for (i=0;i<numberRows_;i++) {
2611          if (rowLower_[i]>0.0||rowUpper_[i]<0.0) {
2612            double scale = rowScale_[i];
2613            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2614            if (rowLower_[i]>0)
2615              largestRhs=CoinMax(largestRhs,rowLower_[i]*scale);
2616            if (rowUpper_[i]<0.0)
2617              largestRhs=CoinMax(largestRhs,-rowUpper_[i]*scale);
2618          }
2619        }
2620        printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2621        bool scalingDone=false;
2622        // look at element range
2623        double smallestNegative;
2624        double largestNegative;
2625        double smallestPositive;
2626        double largestPositive;
2627        matrix_->rangeOfElements(smallestNegative, largestNegative,
2628                                 smallestPositive, largestPositive);
2629        smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2630        largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2631        if (largestObj) {
2632          double ratio = largestObj/smallestObj;
2633          double scale=1.0;
2634          if (ratio<1.0e8) {
2635            // reasonable
2636            if (smallestObj<1.0e-4) {
2637              // may as well scale up
2638              scalingDone=true;
2639              scale=1.0e-3/smallestObj;
2640            } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2641              //done=true;
2642            } else {
2643              scalingDone=true;
2644              if (algorithm_<0) {
2645                scale = 1.0e6/largestObj;
2646              } else {
2647                scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2648              }
2649            }
2650          } else if (ratio<1.0e12) {
2651            // not so good
2652            if (smallestObj<1.0e-7) {
2653              // may as well scale up
2654              scalingDone=true;
2655              scale=1.0e-6/smallestObj;
2656            } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2657              //done=true;
2658            } else {
2659              scalingDone=true;
2660              if (algorithm_<0) {
2661                scale = 1.0e7/largestObj;
2662              } else {
2663                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2664              }
2665            }
2666          } else {
2667            // Really nasty problem
2668            if (smallestObj<1.0e-8) {
2669              // may as well scale up
2670              scalingDone=true;
2671              scale=1.0e-7/smallestObj;
2672              largestObj *= scale;
2673            } 
2674            if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2675              //done=true;
2676            } else {
2677              scalingDone=true;
2678              if (algorithm_<0) {
2679                scale = 1.0e7/largestObj;
2680              } else {
2681                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2682              }
2683            }
2684          }
2685          objectiveScale_=scale;
2686        }
2687        if (largestRhs>1.0e12) {
2688          scalingDone=true;
2689          rhsScale_=1.0e9/largestRhs;
2690        } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
2691          scalingDone=true;
2692          rhsScale_=1.0e6/largestRhs;
2693        } else {
2694          rhsScale_=1.0;
2695        }
2696        if (scalingDone) {
2697          handler_->message(CLP_RIM_SCALE,messages_)
2698            <<objectiveScale_<<rhsScale_
2699            <<CoinMessageEol;
2700        }
2701      }
2702    } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
2703      matrix_->scaleRowCopy(this);
2704    }
2705    // See if we can try for faster row copy
2706    if (makeRowCopy&&!oldMatrix) {
2707      ClpPackedMatrix* clpMatrix =
2708        dynamic_cast< ClpPackedMatrix*>(matrix_);
2709      if (clpMatrix&&numberThreads_) 
2710        clpMatrix->specialRowCopy(this,rowCopy_);
2711      if (clpMatrix)
2712        clpMatrix->specialColumnCopy(this);
2713    }
2714  }
2715  if (what==63) {
2716    if (newArrays) {
2717      delete [] cost_;
2718      // extra copy with original costs
2719      //cost_ = new double[2*numberTotal];
2720      cost_ = new double[numberTotal];
2721      delete [] lower_;
2722      delete [] upper_;
2723      lower_ = new double[numberColumns_+numberRows2];
2724      upper_ = new double[numberColumns_+numberRows2];
2725      delete [] dj_;
2726      dj_ = new double[numberRows2+numberColumns_];
2727      delete [] solution_;
2728      solution_ = new double[numberRows2+numberColumns_];
2729    } else if (auxiliaryModel_) {
2730      assert (!cost_);
2731      cost_=auxiliaryModel_->cost_;
2732      assert (!lower_);
2733      lower_=auxiliaryModel_->lower_;
2734      assert (!upper_);
2735      upper_=auxiliaryModel_->upper_;
2736      assert (!dj_);
2737      dj_=auxiliaryModel_->dj_;
2738      assert (!solution_);
2739      solution_=auxiliaryModel_->solution_;
2740      assert (!rowScale_);
2741      assert (!columnScale_);
2742    }
2743    reducedCostWork_ = dj_;
2744    rowReducedCost_ = dj_+numberColumns_;
2745    columnActivityWork_ = solution_;
2746    rowActivityWork_ = solution_+numberColumns_;
2747    objectiveWork_ = cost_;
2748    rowObjectiveWork_ = cost_+numberColumns_;
2749    rowLowerWork_ = lower_+numberColumns_;
2750    columnLowerWork_ = lower_;
2751    rowUpperWork_ = upper_+numberColumns_;
2752    columnUpperWork_ = upper_;
2753  }
2754  if ((what&4)!=0) {
2755    if (!auxiliaryModel_||(what==63&&(auxiliaryModel_->specialOptions_&4)==0)) {
2756      double direction = optimizationDirection_*objectiveScale_;
2757      const double * obj = objective();
2758      const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
2759      const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
2760      // and also scale by scale factors
2761      if (rowScale) {
2762        if (rowObjective_) {
2763          for (i=0;i<numberRows_;i++)
2764            rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
2765        } else {
2766          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2767        }
2768        // If scaled then do all columns later in one loop
2769        if (what!=63||auxiliaryModel_) {
2770          for (i=0;i<numberColumns_;i++) {
2771            CoinAssert(fabs(obj[i])<1.0e25);
2772            objectiveWork_[i] = obj[i]*direction*columnScale[i];
2773          }
2774        }
2775      } else {
2776        if (rowObjective_) {
2777          for (i=0;i<numberRows_;i++)
2778            rowObjectiveWork_[i] = rowObjective_[i]*direction;
2779        } else {
2780          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2781        }
2782        for (i=0;i<numberColumns_;i++) {
2783          CoinAssert(fabs(obj[i])<1.0e25);
2784          objectiveWork_[i] = obj[i]*direction;
2785        }
2786      }
2787      if (auxiliaryModel_) {
2788        // save costs
2789        memcpy(auxiliaryModel_->cost_+numberTotal,cost_,numberTotal*sizeof(double));
2790      }
2791    } else {
2792      // just copy
2793      memcpy(cost_,auxiliaryModel_->cost_+numberTotal,numberTotal*sizeof(double));
2794    }
2795  }
2796  if ((what&1)!=0) {
2797    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
2798    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
2799    // clean up any mismatches on infinity
2800    // and fix any variables with tiny gaps
2801    double primalTolerance=dblParam_[ClpPrimalTolerance];
2802    if (!auxiliaryModel_) {
2803      if(rowScale) {
2804        // If scaled then do all columns later in one loop
2805        if (what!=63) {
2806          if ((specialOptions_&131072)==0) {
2807            for (i=0;i<numberColumns_;i++) {
2808              double multiplier = rhsScale_/columnScale[i];
2809              double lowerValue = columnLower_[i];
2810              double upperValue = columnUpper_[i];
2811              if (lowerValue>-1.0e20) {
2812                columnLowerWork_[i]=lowerValue*multiplier;
2813                if (upperValue>=1.0e20) {
2814                  columnUpperWork_[i]=COIN_DBL_MAX;
2815                } else {
2816                  columnUpperWork_[i]=upperValue*multiplier;
2817                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
2818                    if (columnLowerWork_[i]>=0.0) {
2819                      columnUpperWork_[i] = columnLowerWork_[i];
2820                    } else if (columnUpperWork_[i]<=0.0) {
2821                      columnLowerWork_[i] = columnUpperWork_[i];
2822                    } else {
2823                      columnUpperWork_[i] = 0.0;
2824                      columnLowerWork_[i] = 0.0;
2825                    }
2826                  }
2827                }
2828              } else if (upperValue<1.0e20) {
2829                columnLowerWork_[i]=-COIN_DBL_MAX;
2830                columnUpperWork_[i]=upperValue*multiplier;
2831              } else {
2832                // free
2833                columnLowerWork_[i]=-COIN_DBL_MAX;
2834                columnUpperWork_[i]=COIN_DBL_MAX;
2835              }
2836            }
2837          } else {
2838            const double * inverseScale = columnScale_+numberColumns_;
2839            for (i=0;i<numberColumns_;i++) {
2840              double multiplier = rhsScale_*inverseScale[i];
2841              double lowerValue = columnLower_[i];
2842              double upperValue = columnUpper_[i];
2843              if (lowerValue>-1.0e20) {
2844                columnLowerWork_[i]=lowerValue*multiplier;
2845                if (upperValue>=1.0e20) {
2846                  columnUpperWork_[i]=COIN_DBL_MAX;
2847                } else {
2848                  columnUpperWork_[i]=upperValue*multiplier;
2849                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
2850                    if (columnLowerWork_[i]>=0.0) {
2851                      columnUpperWork_[i] = columnLowerWork_[i];
2852                    } else if (columnUpperWork_[i]<=0.0) {
2853                      columnLowerWork_[i] = columnUpperWork_[i];
2854                    } else {
2855                      columnUpperWork_[i] = 0.0;
2856                      columnLowerWork_[i] = 0.0;
2857                    }
2858                  }
2859                }
2860              } else if (upperValue<1.0e20) {
2861                columnLowerWork_[i]=-COIN_DBL_MAX;
2862                columnUpperWork_[i]=upperValue*multiplier;
2863              } else {
2864                // free
2865                columnLowerWork_[i]=-COIN_DBL_MAX;
2866                columnUpperWork_[i]=COIN_DBL_MAX;
2867              }
2868            }
2869          }
2870        }
2871        for (i=0;i<numberRows_;i++) {
2872          double multiplier = rhsScale_*rowScale[i];
2873          double lowerValue = rowLower_[i];
2874          double upperValue = rowUpper_[i];
2875          if (lowerValue>-1.0e20) {
2876            rowLowerWork_[i]=lowerValue*multiplier;
2877            if (upperValue>=1.0e20) {
2878              rowUpperWork_[i]=COIN_DBL_MAX;
2879            } else {
2880              rowUpperWork_[i]=upperValue*multiplier;
2881              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
2882                if (rowLowerWork_[i]>=0.0) {
2883                  rowUpperWork_[i] = rowLowerWork_[i];
2884                } else if (rowUpperWork_[i]<=0.0) {
2885                  rowLowerWork_[i] = rowUpperWork_[i];
2886                } else {
2887                  rowUpperWork_[i] = 0.0;
2888                  rowLowerWork_[i] = 0.0;
2889                }
2890              }
2891            }
2892          } else if (upperValue<1.0e20) {
2893            rowLowerWork_[i]=-COIN_DBL_MAX;
2894            rowUpperWork_[i]=upperValue*multiplier;
2895          } else {
2896            // free
2897            rowLowerWork_[i]=-COIN_DBL_MAX;
2898            rowUpperWork_[i]=COIN_DBL_MAX;
2899          }
2900        }
2901      } else if (rhsScale_!=1.0) {
2902        for (i=0;i<numberColumns_;i++) {
2903          double lowerValue = columnLower_[i];
2904          double upperValue = columnUpper_[i];
2905          if (lowerValue>-1.0e20) {
2906            columnLowerWork_[i]=lowerValue*rhsScale_;
2907            if (upperValue>=1.0e20) {
2908              columnUpperWork_[i]=COIN_DBL_MAX;
2909            } else {
2910              columnUpperWork_[i]=upperValue*rhsScale_;
2911              if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
2912                if (columnLowerWork_[i]>=0.0) {
2913                  columnUpperWork_[i] = columnLowerWork_[i];
2914                } else if (columnUpperWork_[i]<=0.0) {
2915                  columnLowerWork_[i] = columnUpperWork_[i];
2916                } else {
2917                  columnUpperWork_[i] = 0.0;
2918                  columnLowerWork_[i] = 0.0;
2919                }
2920              }
2921            }
2922          } else if (upperValue<1.0e20) {
2923            columnLowerWork_[i]=-COIN_DBL_MAX;
2924            columnUpperWork_[i]=upperValue*rhsScale_;
2925          } else {
2926            // free
2927            columnLowerWork_[i]=-COIN_DBL_MAX;
2928            columnUpperWork_[i]=COIN_DBL_MAX;
2929          }
2930        }
2931        for (i=0;i<numberRows_;i++) {
2932          double lowerValue = rowLower_[i];
2933          double upperValue = rowUpper_[i];
2934          if (lowerValue>-1.0e20) {
2935            rowLowerWork_[i]=lowerValue*rhsScale_;
2936            if (upperValue>=1.0e20) {
2937              rowUpperWork_[i]=COIN_DBL_MAX;
2938            } else {
2939              rowUpperWork_[i]=upperValue*rhsScale_;
2940              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
2941                if (rowLowerWork_[i]>=0.0) {
2942                  rowUpperWork_[i] = rowLowerWork_[i];
2943                } else if (rowUpperWork_[i]<=0.0) {
2944                  rowLowerWork_[i] = rowUpperWork_[i];
2945                } else {
2946                  rowUpperWork_[i] = 0.0;
2947                  rowLowerWork_[i] = 0.0;
2948                }
2949              }
2950            }
2951          } else if (upperValue<1.0e20) {
2952            rowLowerWork_[i]=-COIN_DBL_MAX;
2953            rowUpperWork_[i]=upperValue*rhsScale_;
2954          } else {
2955            // free
2956            rowLowerWork_[i]=-COIN_DBL_MAX;
2957            rowUpperWork_[i]=COIN_DBL_MAX;
2958          }
2959        }
2960      } else {
2961        for (i=0;i<numberColumns_;i++) {
2962          double lowerValue = columnLower_[i];
2963          double upperValue = columnUpper_[i];
2964          if (lowerValue>-1.0e20) {
2965            columnLowerWork_[i]=lowerValue;
2966            if (upperValue>=1.0e20) {
2967              columnUpperWork_[i]=COIN_DBL_MAX;
2968            } else {
2969              columnUpperWork_[i]=upperValue;
2970              if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
2971                if (columnLowerWork_[i]>=0.0) {
2972                  columnUpperWork_[i] = columnLowerWork_[i];
2973                } else if (columnUpperWork_[i]<=0.0) {
2974                  columnLowerWork_[i] = columnUpperWork_[i];
2975                } else {
2976                  columnUpperWork_[i] = 0.0;
2977                  columnLowerWork_[i] = 0.0;
2978                }
2979              }
2980            }
2981          } else if (upperValue<1.0e20) {
2982            columnLowerWork_[i]=-COIN_DBL_MAX;
2983            columnUpperWork_[i]=upperValue;
2984          } else {
2985            // free
2986            columnLowerWork_[i]=-COIN_DBL_MAX;
2987            columnUpperWork_[i]=COIN_DBL_MAX;
2988          }
2989        }
2990        for (i=0;i<numberRows_;i++) {
2991          double lowerValue = rowLower_[i];
2992          double upperValue = rowUpper_[i];
2993          if (lowerValue>-1.0e20) {
2994            rowLowerWork_[i]=lowerValue;
2995            if (upperValue>=1.0e20) {
2996              rowUpperWork_[i]=COIN_DBL_MAX;
2997            } else {
2998              rowUpperWork_[i]=upperValue;
2999              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
3000                if (rowLowerWork_[i]>=0.0) {
3001                  rowUpperWork_[i] = rowLowerWork_[i];
3002                } else if (rowUpperWork_[i]<=0.0) {
3003                  rowLowerWork_[i] = rowUpperWork_[i];
3004                } else {
3005                  rowUpperWork_[i] = 0.0;
3006                  rowLowerWork_[i] = 0.0;
3007                }
3008              }
3009            }
3010          } else if (upperValue<1.0e20) {
3011            rowLowerWork_[i]=-COIN_DBL_MAX;
3012            rowUpperWork_[i]=upperValue;
3013          } else {
3014            // free
3015            rowLowerWork_[i]=-COIN_DBL_MAX;
3016            rowUpperWork_[i]=COIN_DBL_MAX;
3017          }
3018        }
3019      }
3020    } else {
3021      // auxiliary model
3022      if (what!=63) {
3023        // just copy
3024        memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberTotal*sizeof(double));
3025        memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberTotal*sizeof(double));
3026      } else {
3027        assert (rhsScale_==1.0);
3028        assert (objectiveScale_==1.0);
3029        if ((auxiliaryModel_->specialOptions_&2)==0) {
3030          // need to do column bounds
3031          if (columnScale) {
3032            const double * inverseScale = columnScale+numberColumns_;
3033            // scaled
3034            for (i=0;i<numberColumns_;i++) {
3035              double value;
3036              value = columnLower_[i];
3037              if (value>-1.0e20) 
3038                value *= inverseScale[i];
3039              lower_[i] = value;
3040              value = columnUpper_[i];
3041              if (value<1.0e20) 
3042                value *= inverseScale[i];
3043              upper_[i] = value;
3044            }
3045          } else {
3046            for (i=0;i<numberColumns_;i++) {
3047              lower_[i] = columnLower_[i];
3048              upper_[i] = columnUpper_[i];
3049            }
3050          }
3051          // save
3052          memcpy(auxiliaryModel_->lower_+numberTotal,lower_,numberColumns_*sizeof(double));
3053          memcpy(auxiliaryModel_->upper_+numberTotal,upper_,numberColumns_*sizeof(double));
3054        } else {
3055          // just copy
3056          memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberColumns_*sizeof(double));
3057          memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberColumns_*sizeof(double));
3058        }
3059        if ((auxiliaryModel_->specialOptions_&1)==0) {
3060          // need to do row bounds
3061          if (rowScale) {
3062            // scaled
3063            for (i=0;i<numberRows_;i++) {
3064              double value;
3065              value = rowLower_[i];
3066              if (value>-1.0e20) 
3067                value *= rowScale[i];
3068              lower_[i+numberColumns_] = value;
3069              value = rowUpper_[i];
3070              if (value<1.0e20) 
3071                value *= rowScale[i];
3072              upper_[i+numberColumns_] = value;
3073            }
3074          } else {
3075            for (i=0;i<numberRows_;i++) {
3076              lower_[i+numberColumns_] = rowLower_[i];
3077              upper_[i+numberColumns_] = rowUpper_[i];
3078            }
3079          }
3080          // save
3081          memcpy(auxiliaryModel_->lower_+numberTotal+numberColumns_,
3082                 lower_+numberColumns_,numberRows_*sizeof(double));
3083          memcpy(auxiliaryModel_->upper_+numberTotal+numberColumns_,
3084                 upper_+numberColumns_,numberRows_*sizeof(double));
3085        } else {
3086          // just copy
3087          memcpy(lower_+numberColumns_,
3088                 auxiliaryModel_->lower_+numberTotal+numberColumns_,numberRows_*sizeof(double));
3089          memcpy(upper_+numberColumns_,
3090                 auxiliaryModel_->upper_+numberTotal+numberColumns_,numberRows_*sizeof(double));
3091        }
3092      }
3093      if (what==63) {
3094        // do basis
3095        assert ((auxiliaryModel_->specialOptions_&8)!=0);
3096        // clean solution trusting basis
3097        for ( i=0;i<numberTotal;i++) {
3098          Status status =getStatus(i);
3099          double value=solution_[i]; // may as well leave if basic (values pass)
3100          if (status!=basic) {
3101            if (upper_[i]==lower_[i]) {
3102              setStatus(i,isFixed);
3103              value=lower_[i];
3104            } else if (status==atLowerBound) {
3105              if (lower_[i]>-1.0e20) {
3106                value=lower_[i];
3107              } else {
3108                if (upper_[i]<1.0e20) {
3109                  value=upper_[i];
3110                  setStatus(i,atUpperBound);
3111                } else {
3112                  setStatus(i,isFree);
3113                }
3114              }
3115            } else if (status==atUpperBound) {
3116              if (upper_[i]<1.0e20) {
3117              value=upper_[i];
3118              } else {
3119                if (lower_[i]>-1.0e20) {
3120                  value=lower_[i];
3121                  setStatus(i,atLowerBound);
3122                } else {
3123                  setStatus(i,isFree);
3124                }
3125              }
3126            }
3127          }
3128          solution_[i]=value;
3129        }
3130      }
3131    }
3132  }
3133  if (what==63) {
3134    // move information to work arrays
3135    double direction = optimizationDirection_;
3136    // direction is actually scale out not scale in
3137    if (direction)
3138      direction = 1.0/direction;
3139    if (direction!=1.0&&(!auxiliaryModel_||(auxiliaryModel_->specialOptions_&8)==0)) {
3140      // reverse all dual signs
3141      for (i=0;i<numberColumns_;i++) 
3142        reducedCost_[i] *= direction;
3143      for (i=0;i<numberRows_;i++) 
3144        dual_[i] *= direction;
3145    }
3146    for (i=0;i<numberRows_+numberColumns_;i++) {
3147      setFakeBound(i,noFake);
3148    }
3149    if (!auxiliaryModel_) {
3150      if (rowScale_) {
3151        const double * obj = objective();
3152        double direction = optimizationDirection_*objectiveScale_;
3153        // clean up any mismatches on infinity
3154        // and fix any variables with tiny gaps
3155        double primalTolerance=dblParam_[ClpPrimalTolerance];
3156        // on entry
3157        if ((specialOptions_&131072)==0) {
3158          for (i=0;i<numberColumns_;i++) {
3159            CoinAssert(fabs(obj[i])<1.0e25);
3160            double scaleFactor = columnScale_[i];
3161            double multiplier = rhsScale_/scaleFactor;
3162            scaleFactor *= direction;
3163            objectiveWork_[i] = obj[i]*scaleFactor;
3164            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
3165            double lowerValue = columnLower_[i];
3166            double upperValue = columnUpper_[i];
3167            if (lowerValue>-1.0e20) {
3168              columnLowerWork_[i]=lowerValue*multiplier;
3169              if (upperValue>=1.0e20) {
3170                columnUpperWork_[i]=COIN_DBL_MAX;
3171              } else {
3172                columnUpperWork_[i]=upperValue*multiplier;
3173                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3174                  if (columnLowerWork_[i]>=0.0) {
3175                    columnUpperWork_[i] = columnLowerWork_[i];
3176                  } else if (columnUpperWork_[i]<=0.0) {
3177                    columnLowerWork_[i] = columnUpperWork_[i];
3178                  } else {
3179                    columnUpperWork_[i] = 0.0;
3180                    columnLowerWork_[i] = 0.0;
3181                  }
3182                }
3183              }
3184            } else if (upperValue<1.0e20) {
3185              columnLowerWork_[i]=-COIN_DBL_MAX;
3186              columnUpperWork_[i]=upperValue*multiplier;
3187            } else {
3188              // free
3189              columnLowerWork_[i]=-COIN_DBL_MAX;
3190              columnUpperWork_[i]=COIN_DBL_MAX;
3191            }
3192            double value = columnActivity_[i] * multiplier;
3193            if (fabs(value)>1.0e20) {
3194              //printf("bad value of %g for column %d\n",value,i);
3195              setColumnStatus(i,superBasic);
3196              if (columnUpperWork_[i]<0.0) {
3197                value=columnUpperWork_[i];
3198              } else if (columnLowerWork_[i]>0.0) {
3199                value=columnLowerWork_[i];
3200              } else {
3201                value=0.0;
3202              }
3203            }
3204            columnActivityWork_[i] = value;
3205          }
3206          for (i=0;i<numberRows_;i++) {
3207            dual_[i] /= rowScale_[i];
3208            dual_[i] *= objectiveScale_;
3209            rowReducedCost_[i] = dual_[i];
3210            double multiplier = rhsScale_*rowScale_[i];
3211            double value = rowActivity_[i] * multiplier;
3212            if (fabs(value)>1.0e20) {
3213              //printf("bad value of %g for row %d\n",value,i);
3214              setRowStatus(i,superBasic);
3215              if (rowUpperWork_[i]<0.0) {
3216                value=rowUpperWork_[i];
3217              } else if (rowLowerWork_[i]>0.0) {
3218                value=rowLowerWork_[i];
3219              } else {
3220                value=0.0;
3221              }
3222            }
3223            rowActivityWork_[i] = value;
3224          }
3225        } else {
3226          const double * inverseScale = columnScale_+numberColumns_;
3227          for (i=0;i<numberColumns_;i++) {
3228            CoinAssert(fabs(obj[i])<1.0e25);
3229            double scaleFactor = columnScale_[i];
3230            double multiplier = rhsScale_*inverseScale[i];
3231            scaleFactor *= direction;
3232            objectiveWork_[i] = obj[i]*scaleFactor;
3233            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
3234            double lowerValue = columnLower_[i];
3235            double upperValue = columnUpper_[i];
3236            if (lowerValue>-1.0e20) {
3237              columnLowerWork_[i]=lowerValue*multiplier;
3238              if (upperValue>=1.0e20) {
3239                columnUpperWork_[i]=COIN_DBL_MAX;
3240              } else {
3241                columnUpperWork_[i]=upperValue*multiplier;
3242                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3243                  if (columnLowerWork_[i]>=0.0) {
3244                    columnUpperWork_[i] = columnLowerWork_[i];
3245                  } else if (columnUpperWork_[i]<=0.0) {
3246                    columnLowerWork_[i] = columnUpperWork_[i];
3247                  } else {
3248                    columnUpperWork_[i] = 0.0;
3249                    columnLowerWork_[i] = 0.0;
3250                  }
3251                }
3252              }
3253            } else if (upperValue<1.0e20) {
3254              columnLowerWork_[i]=-COIN_DBL_MAX;
3255              columnUpperWork_[i]=upperValue*multiplier;
3256            } else {
3257              // free
3258              columnLowerWork_[i]=-COIN_DBL_MAX;
3259              columnUpperWork_[i]=COIN_DBL_MAX;
3260            }
3261            double value = columnActivity_[i] * multiplier;
3262            if (fabs(value)>1.0e20) {
3263              //printf("bad value of %g for column %d\n",value,i);
3264              setColumnStatus(i,superBasic);
3265              if (columnUpperWork_[i]<0.0) {
3266                value=columnUpperWork_[i];
3267              } else if (columnLowerWork_[i]>0.0) {
3268                value=columnLowerWork_[i];
3269              } else {
3270                value=0.0;
3271              }
3272            }
3273            columnActivityWork_[i] = value;
3274          }
3275          inverseScale = rowScale_+numberRows_;
3276          for (i=0;i<numberRows_;i++) {
3277            dual_[i] *= inverseScale[i];
3278            dual_[i] *= objectiveScale_;
3279            rowReducedCost_[i] = dual_[i];
3280            double multiplier = rhsScale_*rowScale_[i];
3281            double value = rowActivity_[i] * multiplier;
3282            if (fabs(value)>1.0e20) {
3283              //printf("bad value of %g for row %d\n",value,i);
3284              setRowStatus(i,superBasic);
3285              if (rowUpperWork_[i]<0.0) {
3286                value=rowUpperWork_[i];
3287              } else if (rowLowerWork_[i]>0.0) {
3288                value=rowLowerWork_[i];
3289              } else {
3290                value=0.0;
3291              }
3292            }
3293            rowActivityWork_[i] = value;
3294          }
3295        }
3296      } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
3297        // on entry
3298        for (i=0;i<numberColumns_;i++) {
3299          double value = columnActivity_[i];
3300          value *= rhsScale_;
3301          if (fabs(value)>1.0e20) {
3302            //printf("bad value of %g for column %d\n",value,i);
3303            setColumnStatus(i,superBasic);
3304            if (columnUpperWork_[i]<0.0) {
3305              value=columnUpperWork_[i];
3306            } else if (columnLowerWork_[i]>0.0) {
3307              value=columnLowerWork_[i];
3308            } else {
3309              value=0.0;
3310            }
3311          }
3312          columnActivityWork_[i] = value;
3313          reducedCostWork_[i] = reducedCost_[i]*objectiveScale_;
3314        }
3315        for (i=0;i<numberRows_;i++) {
3316          double value = rowActivity_[i];
3317          value *= rhsScale_;
3318          if (fabs(value)>1.0e20) {
3319            //printf("bad value of %g for row %d\n",value,i);
3320            setRowStatus(i,superBasic);
3321            if (rowUpperWork_[i]<0.0) {
3322              value=rowUpperWork_[i];
3323            } else if (rowLowerWork_[i]>0.0) {
3324              value=rowLowerWork_[i];
3325            } else {
3326              value=0.0;
3327            }
3328          }
3329          rowActivityWork_[i] = value;
3330          dual_[i] *= objectiveScale_;
3331          rowReducedCost_[i] = dual_[i];
3332        }
3333      } else {
3334        // on entry
3335        for (i=0;i<numberColumns_;i++) {
3336          double value = columnActivity_[i];
3337          if (fabs(value)>1.0e20) {
3338            //printf("bad value of %g for column %d\n",value,i);
3339            setColumnStatus(i,superBasic);
3340            if (columnUpperWork_[i]<0.0) {
3341              value=columnUpperWork_[i];
3342            } else if (columnLowerWork_[i]>0.0) {
3343              value=columnLowerWork_[i];
3344            } else {
3345              value=0.0;
3346            }
3347          }
3348          columnActivityWork_[i] = value;
3349          reducedCostWork_[i] = reducedCost_[i];
3350        }
3351        for (i=0;i<numberRows_;i++) {
3352          double value = rowActivity_[i];
3353          if (fabs(value)>1.0e20) {
3354            //printf("bad value of %g for row %d\n",value,i);
3355            setRowStatus(i,superBasic);
3356            if (rowUpperWork_[i]<0.0) {
3357              value=rowUpperWork_[i];
3358            } else if (rowLowerWork_[i]>0.0) {
3359              value=rowLowerWork_[i];
3360            } else {
3361              value=0.0;
3362            }
3363          }
3364          rowActivityWork_[i] = value;
3365          rowReducedCost_[i] = dual_[i];
3366        }
3367      }
3368    }
3369  }
3370 
3371  if (what==63&&doSanityCheck&&!auxiliaryModel_) {
3372    // check rim of problem okay
3373    if (!sanityCheck())
3374      goodMatrix=false;
3375  } 
3376  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
3377  // maybe we need to move scales to SimplexModel for factorization?
3378  if ((what==63&&!pivotVariable_)||(newArrays&&!keepPivots)) {
3379    delete [] pivotVariable_;
3380    pivotVariable_=new int[numberRows2];
3381    for (int i=0;i<numberRows2;i++)
3382      pivotVariable_[i]=-1;
3383  } else if (what==63&&!keepPivots) {
3384    // just reset
3385    for (int i=0;i<numberRows2;i++)
3386      pivotVariable_[i]=-1;
3387  } else if (what==63) {
3388    // check pivots
3389    for (int i=0;i<numberRows2;i++) {
3390      int iSequence = pivotVariable_[i];
3391      if (iSequence<numberRows_+numberColumns_&&
3392          getStatus(iSequence)!=basic) {
3393        keepPivots =false;
3394        break;
3395      }
3396    }
3397    if (!keepPivots) {
3398      // reset
3399      for (int i=0;i<numberRows2;i++)
3400        pivotVariable_[i]=-1;
3401    } else {
3402      // clean
3403      for (int i=0;i<numberColumns_+numberRows_;i++) {
3404        Status status =getStatus(i);
3405        if (status!=basic) {
3406          if (upper_[i]==lower_[i]) {
3407            setStatus(i,isFixed);
3408            solution_[i]=lower_[i];
3409          } else if (status==atLowerBound) {
3410            if (lower_[i]>-1.0e20) {
3411              solution_[i]=lower_[i];
3412            } else {
3413              //printf("seq %d at lower of %g\n",i,lower_[i]);
3414              if (upper_[i]<1.0e20) {
3415                solution_[i]=upper_[i];
3416                setStatus(i,atUpperBound);
3417              } else {
3418                setStatus(i,isFree);
3419              }
3420            }
3421          } else if (status==atUpperBound) {
3422            if (upper_[i]<1.0e20) {
3423              solution_[i]=upper_[i];
3424            } else {
3425              //printf("seq %d at upper of %g\n",i,upper_[i]);
3426              if (lower_[i]>-1.0e20) {
3427                solution_[i]=lower_[i];
3428                setStatus(i,atLowerBound);
3429              } else {
3430                setStatus(i,isFree);
3431              }
3432            }
3433          }
3434        }
3435      }
3436    }
3437  }
3438 
3439  if (what==63) {
3440    if (newArrays) {
3441      // get some arrays
3442      int iRow,iColumn;
3443      // these are "indexed" arrays so we always know where nonzeros are
3444      /**********************************************************
3445      rowArray_[3] is long enough for rows+columns
3446      *********************************************************/
3447      for (iRow=0;iRow<4;iRow++) {
3448        int length =numberRows2+factorization_->maximumPivots();
3449        if (iRow==3||objective_->type()>1)
3450          length += numberColumns_;
3451        if ((specialOptions_&65536)==0||!rowArray_[iRow]) {
3452          delete rowArray_[iRow];
3453          rowArray_[iRow]=new CoinIndexedVector();
3454        }
3455        rowArray_[iRow]->reserve(length);
3456      }
3457     
3458      for (iColumn=0;iColumn<2;iColumn++) {
3459        if ((specialOptions_&65536)==0||!columnArray_[iColumn]) {
3460          delete columnArray_[iColumn];
3461          columnArray_[iColumn]=new CoinIndexedVector();
3462        }
3463        if (!iColumn)
3464          columnArray_[iColumn]->reserve(numberColumns_);
3465        else
3466          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
3467      }
3468    } else {
3469      int iRow,iColumn;
3470      if (auxiliaryModel_) {
3471        for (iRow=0;iRow<4;iRow++) {
3472          assert (!rowArray_[iRow]);
3473          rowArray_[iRow]=auxiliaryModel_->rowArray_[iRow];
3474          // For safety
3475          memset(rowArray_[iRow]->denseVector(),0,rowArray_[iRow]->capacity()*sizeof(double));
3476        }
3477        for (iColumn=0;iColumn<2;iColumn++) {
3478          assert (!columnArray_[iColumn]);
3479          columnArray_[iColumn]=auxiliaryModel_->columnArray_[iColumn];
3480          // For safety
3481          memset(columnArray_[iColumn]->denseVector(),0,columnArray_[iColumn]->capacity()*sizeof(double));
3482        }
3483        // do matrices
3484        rowCopy_=auxiliaryModel_->rowCopy_;
3485        ClpMatrixBase * temp = matrix_;
3486        matrix_=auxiliaryModel_->matrix_;
3487        auxiliaryModel_->matrix_=temp;
3488      }
3489      for (iRow=0;iRow<4;iRow++) {
3490        rowArray_[iRow]->clear();
3491#ifndef NDEBUG
3492        int length =numberRows2+factorization_->maximumPivots();
3493        if (iRow==3||objective_->type()>1)
3494          length += numberColumns_;
3495        assert(rowArray_[iRow]->capacity()>=length);
3496        rowArray_[iRow]->checkClear();
3497#endif
3498      }
3499     
3500      for (iColumn=0;iColumn<2;iColumn++) {
3501        columnArray_[iColumn]->clear();
3502#ifndef NDEBUG
3503        int length =numberColumns_;
3504        if (iColumn)
3505          length=CoinMax(numberRows2,numberColumns_);
3506        assert(columnArray_[iColumn]->capacity()>=length);
3507        columnArray_[iColumn]->checkClear();
3508#endif
3509      }
3510    }   
3511  }
3512  if (problemStatus_==10) {
3513    problemStatus_=-1;
3514    handler_->setLogLevel(saveLevel); // switch back messages
3515  }
3516  if ((what&5)!=0) 
3517    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
3518  return goodMatrix;
3519}
3520void
3521ClpSimplex::deleteRim(int getRidOfFactorizationData)
3522{
3523  // Just possible empty problem
3524  int numberRows=numberRows_;
3525  int numberColumns=numberColumns_;
3526  if (!numberRows||!numberColumns) {
3527    numberRows=0;
3528    numberColumns=0;
3529  }
3530  if (!auxiliaryModel_) {
3531    int i;
3532    if (problemStatus_!=1&&problemStatus_!=2) {
3533      delete [] ray_;
3534      ray_=NULL;
3535    }
3536    // set upperOut_ to furthest away from bound so can use in dual for dualBound_
3537    upperOut_=1.0;
3538#if 0
3539    {
3540      int nBad=0;
3541      for (i=0;i<numberColumns;i++) {
3542        if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
3543          nBad++;
3544      }
3545      if (nBad)
3546        printf("yy %d basic fixed\n",nBad);
3547    }
3548#endif
3549    // ray may be null if in branch and bound
3550    if (rowScale_) {
3551      // Collect infeasibilities
3552      int numberPrimalScaled=0;
3553      int numberPrimalUnscaled=0;
3554      int numberDualScaled=0;
3555      int numberDualUnscaled=0;
3556      double scaleC = 1.0/objectiveScale_;
3557      double scaleR = 1.0/rhsScale_;
3558      if ((specialOptions_&131072)==0) {
3559        for (i=0;i<numberColumns;i++) {
3560          double scaleFactor = columnScale_[i];
3561          double valueScaled = columnActivityWork_[i];
3562          double lowerScaled = columnLowerWork_[i];
3563          double upperScaled = columnUpperWork_[i];
3564          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3565            if (valueScaled<lowerScaled-primalTolerance_||
3566                valueScaled>upperScaled+primalTolerance_)
3567              numberPrimalScaled++;
3568            else
3569              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3570          }
3571          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
3572          double value = columnActivity_[i];
3573          if (value<columnLower_[i]-primalTolerance_)
3574            numberPrimalUnscaled++;
3575          else if (value>columnUpper_[i]+primalTolerance_)
3576            numberPrimalUnscaled++;
3577          double valueScaledDual = reducedCostWork_[i];
3578          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3579            numberDualScaled++;
3580          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3581            numberDualScaled++;
3582          reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
3583          double valueDual = reducedCost_[i];
3584          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3585            numberDualUnscaled++;
3586          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3587            numberDualUnscaled++;
3588        }
3589        for (i=0;i<numberRows;i++) {
3590          double scaleFactor = rowScale_[i];
3591          double valueScaled = rowActivityWork_[i];
3592          double lowerScaled = rowLowerWork_[i];
3593          double upperScaled = rowUpperWork_[i];
3594          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3595            if (valueScaled<lowerScaled-primalTolerance_||
3596                valueScaled>upperScaled+primalTolerance_)
3597              numberPrimalScaled++;
3598            else
3599              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3600          }
3601          rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
3602          double value = rowActivity_[i];
3603          if (value<rowLower_[i]-primalTolerance_)
3604            numberPrimalUnscaled++;
3605          else if (value>rowUpper_[i]+primalTolerance_)
3606            numberPrimalUnscaled++;
3607          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
3608          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3609            numberDualScaled++;
3610          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3611            numberDualScaled++;
3612          dual_[i] *= scaleFactor*scaleC;
3613          double valueDual = dual_[i]; 
3614          if (rowObjective_)
3615            valueDual += rowObjective_[i];
3616          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3617            numberDualUnscaled++;
3618          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3619            numberDualUnscaled++;
3620        }
3621      } else {
3622        const double * inverseScale = columnScale_+numberColumns;
3623        for (i=0;i<numberColumns;i++) {
3624          double scaleFactor = columnScale_[i];
3625          double valueScaled = columnActivityWork_[i];
3626          double lowerScaled = columnLowerWork_[i];
3627          double upperScaled = columnUpperWork_[i];
3628          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3629            if (valueScaled<lowerScaled-primalTolerance_||
3630                valueScaled>upperScaled+primalTolerance_)
3631              numberPrimalScaled++;
3632            else
3633              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3634          }
3635          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
3636          double value = columnActivity_[i];
3637          if (value<columnLower_[i]-primalTolerance_)
3638            numberPrimalUnscaled++;
3639          else if (value>columnUpper_[i]+primalTolerance_)
3640            numberPrimalUnscaled++;
3641          double valueScaledDual = reducedCostWork_[i];
3642          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3643            numberDualScaled++;
3644          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3645            numberDualScaled++;
3646          reducedCost_[i] = (valueScaledDual*scaleC)*inverseScale[i];
3647          double valueDual = reducedCost_[i];
3648          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3649            numberDualUnscaled++;
3650          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3651            numberDualUnscaled++;
3652        }
3653        inverseScale = rowScale_+numberRows;
3654        for (i=0;i<numberRows;i++) {
3655          double scaleFactor = rowScale_[i];
3656          double valueScaled = rowActivityWork_[i];
3657          double lowerScaled = rowLowerWork_[i];
3658          double upperScaled = rowUpperWork_[i];
3659          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3660            if (valueScaled<lowerScaled-primalTolerance_||
3661                valueScaled>upperScaled+primalTolerance_)
3662              numberPrimalScaled++;
3663            else
3664              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3665          }
3666          rowActivity_[i] = (valueScaled*scaleR)*inverseScale[i];
3667          double value = rowActivity_[i];
3668          if (value<rowLower_[i]-primalTolerance_)
3669            numberPrimalUnscaled++;
3670          else if (value>rowUpper_[i]+primalTolerance_)
3671            numberPrimalUnscaled++;
3672          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
3673          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3674            numberDualScaled++;
3675          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3676            numberDualScaled++;
3677          dual_[i] *= scaleFactor*scaleC;
3678          double valueDual = dual_[i]; 
3679          if (rowObjective_)
3680            valueDual += rowObjective_[i];
3681          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3682            numberDualUnscaled++;
3683          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3684            numberDualUnscaled++;
3685        }
3686      }
3687      if (!problemStatus_&&!secondaryStatus_) {
3688        // See if we need to set secondary status
3689        if (numberPrimalUnscaled) {
3690          if (numberDualUnscaled) 
3691            secondaryStatus_=4;
3692          else
3693            secondaryStatus_=2;
3694        } else {
3695          if (numberDualUnscaled) 
3696            secondaryStatus_=3;
3697        }
3698      }
3699      if (problemStatus_==2) {
3700        for (i=0;i<numberColumns;i++) {
3701          ray_[i] *= columnScale_[i];
3702        }
3703      } else if (problemStatus_==1&&ray_) {
3704        for (i=0;i<numberRows;i++) {
3705          ray_[i] *= rowScale_[i];
3706        }
3707      }
3708    } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
3709      // Collect infeasibilities
3710      int numberPrimalScaled=0;
3711      int numberPrimalUnscaled=0;
3712      int numberDualScaled=0;
3713      int numberDualUnscaled=0;
3714      double scaleC = 1.0/objectiveScale_;
3715      double scaleR = 1.0/rhsScale_;
3716      for (i=0;i<numberColumns;i++) {
3717        double valueScaled = columnActivityWork_[i];
3718        double lowerScaled = columnLowerWork_[i];
3719        double upperScaled = columnUpperWork_[i];
3720        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3721          if (valueScaled<lowerScaled-primalTolerance_||
3722              valueScaled>upperScaled+primalTolerance_)
3723            numberPrimalScaled++;
3724          else
3725            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3726        }
3727        columnActivity_[i] = valueScaled*scaleR;
3728        double value = columnActivity_[i];
3729        if (value<columnLower_[i]-primalTolerance_)
3730          numberPrimalUnscaled++;
3731        else if (value>columnUpper_[i]+primalTolerance_)
3732          numberPrimalUnscaled++;
3733        double valueScaledDual = reducedCostWork_[i];
3734        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3735          numberDualScaled++;
3736        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3737          numberDualScaled++;
3738        reducedCost_[i] = valueScaledDual*scaleC;
3739        double valueDual = reducedCost_[i];
3740        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3741          numberDualUnscaled++;
3742        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3743          numberDualUnscaled++;
3744      }
3745      for (i=0;i<numberRows;i++) {
3746        double valueScaled = rowActivityWork_[i];
3747        double lowerScaled = rowLowerWork_[i];
3748        double upperScaled = rowUpperWork_[i];
3749        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3750          if (valueScaled<lowerScaled-primalTolerance_||
3751              valueScaled>upperScaled+primalTolerance_)
3752            numberPrimalScaled++;
3753          else
3754            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3755        }
3756        rowActivity_[i] = valueScaled*scaleR;
3757        double value = rowActivity_[i];
3758        if (value<rowLower_[i]-primalTolerance_)
3759          numberPrimalUnscaled++;
3760        else if (value>rowUpper_[i]+primalTolerance_)
3761          numberPrimalUnscaled++;
3762        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
3763        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3764          numberDualScaled++;
3765        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3766          numberDualScaled++;
3767        dual_[i] *= scaleC;
3768        double valueDual = dual_[i]; 
3769        if (rowObjective_)
3770          valueDual += rowObjective_[i];
3771        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3772          numberDualUnscaled++;
3773        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3774          numberDualUnscaled++;
3775      }
3776      if (!problemStatus_&&!secondaryStatus_) {
3777        // See if we need to set secondary status
3778        if (numberPrimalUnscaled) {
3779          if (numberDualUnscaled) 
3780            secondaryStatus_=4;
3781          else
3782            secondaryStatus_=2;
3783        } else {
3784          if (numberDualUnscaled) 
3785            secondaryStatus_=3;
3786        }
3787      }
3788    } else {
3789      if (columnActivityWork_) {
3790        for (i=0;i<numberColumns;i++) {
3791          double value = columnActivityWork_[i];
3792          double lower = columnLowerWork_[i];
3793          double upper = columnUpperWork_[i];
3794          if (lower>-1.0e20||upper<1.0e20) {
3795            if (value>lower&&value<upper)
3796              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3797          }
3798          columnActivity_[i] = columnActivityWork_[i];
3799          reducedCost_[i] = reducedCostWork_[i];
3800        }
3801        for (i=0;i<numberRows;i++) {
3802          double value = rowActivityWork_[i];
3803          double lower = rowLowerWork_[i];
3804          double upper = rowUpperWork_[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          rowActivity_[i] = rowActivityWork_[i];
3810        }
3811      }
3812    }
3813    // switch off scalefactor if auto
3814    if (automaticScale_) {
3815      rhsScale_=1.0;
3816      objectiveScale_=1.0;
3817    }
3818    if (optimizationDirection_!=1.0) {
3819      // and modify all dual signs
3820      for (i=0;i<numberColumns;i++) 
3821        reducedCost_[i] *= optimizationDirection_;
3822      for (i=0;i<numberRows;i++) 
3823        dual_[i] *= optimizationDirection_;
3824    }
3825  } else {
3826    // auxiliary model
3827    cost_=NULL;
3828    lower_=NULL;
3829    upper_=NULL;
3830    dj_=NULL;
3831    solution_=NULL;
3832    int iRow,iColumn;
3833    for (iRow=0;iRow<4;iRow++) {
3834      if (rowArray_[iRow]) {
3835        rowArray_[iRow]->clear();
3836        rowArray_[iRow]->checkClear();
3837      }
3838      rowArray_[iRow]=NULL;
3839    }
3840    for (iColumn=0;iColumn<2;iColumn++) {
3841      if (columnArray_[iColumn]) {
3842        columnArray_[iColumn]->clear();
3843        columnArray_[iColumn]->checkClear();
3844      }
3845      columnArray_[iColumn]=NULL;
3846    }
3847    rowCopy_=NULL;
3848    ClpMatrixBase * temp = matrix_;
3849    matrix_=auxiliaryModel_->matrix_;
3850    auxiliaryModel_->matrix_=temp;
3851    assert ((auxiliaryModel_->specialOptions_&16+32)==16+32);
3852    if (problemStatus_!=1&&problemStatus_!=10) {
3853      int i;
3854      if (auxiliaryModel_->rowScale_) {
3855        const double * scale = auxiliaryModel_->columnScale_;
3856        const double * inverseScale = scale + numberColumns;
3857        for (i=0;i<numberColumns;i++) {
3858          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i]*scale[i];
3859          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i]*inverseScale[i];
3860        }
3861        scale = auxiliaryModel_->rowScale_;
3862        inverseScale = scale + numberRows;
3863        for (i=0;i<numberRows;i++) {
3864          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
3865        }
3866      } else {
3867        for (i=0;i<numberColumns;i++) {
3868          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i];
3869          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i];
3870        }
3871        for (i=0;i<numberRows;i++) {
3872          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i];
3873        }
3874      }
3875      if (optimizationDirection_!=1.0) {
3876        // and modify reduced costs
3877        for (i=0;i<numberColumns;i++) 
3878          reducedCost_[i] *= optimizationDirection_;
3879      }
3880    } else if (problemStatus_==10) {
3881      int i;
3882      if (auxiliaryModel_->rowScale_) {
3883        const double * scale = auxiliaryModel_->columnScale_;
3884        const double * inverseScale = scale + numberColumns;
3885        for (i=0;i<numberColumns;i++) {
3886          double lower = auxiliaryModel_->columnLowerWork_[i];
3887          double upper = auxiliaryModel_->columnUpperWork_[i];
3888          double value = auxiliaryModel_->columnActivityWork_[i];
3889          if (value>lower&&value<upper) {
3890            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3891          }
3892          columnActivity_[i] = value*scale[i];
3893        }
3894        scale = auxiliaryModel_->rowScale_;
3895        inverseScale = scale + numberRows;
3896        for (i=0;i<numberRows;i++) {
3897          double lower = auxiliaryModel_->rowLowerWork_[i];
3898          double upper = auxiliaryModel_->rowUpperWork_[i];
3899          double value = auxiliaryModel_->rowActivityWork_[i];
3900          if (value>lower&&value<upper) {
3901            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3902          }
3903          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
3904        }
3905      } else {
3906        for (i=0;i<numberColumns;i++) {
3907          double lower = auxiliaryModel_->columnLowerWork_[i];
3908          double upper = auxiliaryModel_->columnUpperWork_[i];
3909          double value = auxiliaryModel_->columnActivityWork_[i];
3910          if (value>lower&&value<upper) {
3911            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3912          }
3913          columnActivity_[i] = value;
3914        }
3915        for (i=0;i<numberRows;i++) {
3916          double lower = auxiliaryModel_->rowLowerWork_[i];
3917          double upper = auxiliaryModel_->rowUpperWork_[i];
3918          double value = auxiliaryModel_->rowActivityWork_[i];
3919          if (value>lower&&value<upper) {
3920            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3921          }
3922          rowActivity_[i] = value;
3923        }
3924      }
3925    }
3926  }
3927  // scaling may have been turned off
3928  scalingFlag_ = abs(scalingFlag_);
3929  if(getRidOfFactorizationData>0) {
3930    gutsOfDelete(getRidOfFactorizationData+1);
3931  } else {
3932    // at least get rid of nonLinearCost_
3933    delete nonLinearCost_;
3934    nonLinearCost_=NULL;
3935  }
3936  // get rid of data
3937  matrix_->generalExpanded(this,13,scalingFlag_);
3938}
3939void 
3940ClpSimplex::setDualBound(double value)
3941{
3942  if (value>0.0)
3943    dualBound_=value;
3944}
3945void 
3946ClpSimplex::setInfeasibilityCost(double value)
3947{
3948  if (value>0.0)
3949    infeasibilityCost_=value;
3950}
3951void ClpSimplex::setNumberRefinements( int value) 
3952{
3953  if (value>=0&&value<10)
3954    numberRefinements_=value;
3955}
3956// Sets row pivot choice algorithm in dual
3957void 
3958ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
3959{
3960  delete dualRowPivot_;
3961  dualRowPivot_ = choice.clone(true);
3962}
3963// Sets row pivot choice algorithm in dual
3964void 
3965ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
3966{
3967  delete primalColumnPivot_;
3968  primalColumnPivot_ = choice.clone(true);
3969}
3970// Passes in factorization
3971void 
3972ClpSimplex::setFactorization( ClpFactorization & factorization)
3973{
3974  delete factorization_;
3975  factorization_= new ClpFactorization(factorization);
3976}
3977/* Perturbation:
3978   -50 to +50 - perturb by this power of ten (-6 sounds good)
3979   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
3980   101 - we are perturbed
3981   102 - don't try perturbing again
3982   default is 100
3983*/
3984void 
3985ClpSimplex::setPerturbation(int value)
3986{
3987  if(value<=100&&value >=-1000) {
3988    perturbation_=value;
3989  } 
3990}
3991// Sparsity on or off
3992bool 
3993ClpSimplex::sparseFactorization() const
3994{
3995  return factorization_->sparseThreshold()!=0;
3996}
3997void 
3998ClpSimplex::setSparseFactorization(bool value)
3999{
4000  if (value) {
4001    if (!factorization_->sparseThreshold())
4002      factorization_->goSparse();
4003  } else {
4004    factorization_->sparseThreshold(0);
4005  }
4006}
4007void checkCorrect(ClpSimplex * model,int iRow,
4008                  const double * element,const int * rowStart,const int * rowLength,
4009                  const int * column,
4010                  const double * columnLower_, const double * columnUpper_,
4011                  int infiniteUpperC,
4012                  int infiniteLowerC,
4013                  double &maximumUpC,
4014                  double &maximumDownC)
4015{
4016  int infiniteUpper = 0;
4017  int infiniteLower = 0;
4018  double maximumUp = 0.0;
4019  double maximumDown = 0.0;
4020  CoinBigIndex rStart = rowStart[iRow];
4021  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
4022  CoinBigIndex j;
4023  double large=1.0e15;
4024  int iColumn;
4025  // Compute possible lower and upper ranges
4026 
4027  for (j = rStart; j < rEnd; ++j) {
4028    double value=element[j];
4029    iColumn = column[j];
4030    if (value > 0.0) {
4031      if (columnUpper_[iColumn] >= large) {
4032        ++infiniteUpper;
4033      } else {
4034        maximumUp += columnUpper_[iColumn] * value;
4035      }
4036      if (columnLower_[iColumn] <= -large) {
4037        ++infiniteLower;
4038      } else {
4039        maximumDown += columnLower_[iColumn] * value;
4040      }
4041    } else if (value<0.0) {
4042      if (columnUpper_[iColumn] >= large) {
4043        ++infiniteLower;
4044      } else {
4045        maximumDown += columnUpper_[iColumn] * value;
4046      }
4047      if (columnLower_[iColumn] <= -large) {
4048        ++infiniteUpper;
4049      } else {
4050        maximumUp += columnLower_[iColumn] * value;
4051      }
4052    }
4053  }
4054  assert (infiniteLowerC==infiniteLower);
4055  assert (infiniteUpperC==infiniteUpper);
4056  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
4057    printf("row %d comp up %g, true up %g\n",iRow,
4058           maximumUpC,maximumUp);
4059  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
4060    printf("row %d comp down %g, true down %g\n",iRow,
4061           maximumDownC,maximumDown);
4062  maximumUpC=maximumUp;
4063  maximumDownC=maximumDown;
4064}
4065
4066/* Tightens primal bounds to make dual faster.  Unless
4067   fixed, bounds are slightly looser than they could be.
4068   This is to make dual go faster and is probably not needed
4069   with a presolve.  Returns non-zero if problem infeasible
4070
4071   Fudge for branch and bound - put bounds on columns of factor *
4072   largest value (at continuous) - should improve stability
4073   in branch and bound on infeasible branches (0.0 is off)
4074*/
4075int 
4076ClpSimplex::tightenPrimalBounds(double factor,int doTight)
4077{
4078 
4079  // Get a row copy in standard format
4080  CoinPackedMatrix copy;
4081  copy.reverseOrderedCopyOf(*matrix());
4082  // get matrix data pointers
4083  const int * column = copy.getIndices();
4084  const CoinBigIndex * rowStart = copy.getVectorStarts();
4085  const int * rowLength = copy.getVectorLengths(); 
4086  const double * element = copy.getElements();
4087  int numberChanged=1,iPass=0;
4088  double large = largeValue(); // treat bounds > this as infinite
4089#ifndef NDEBUG
4090  double large2= 1.0e10*large;
4091#endif
4092  int numberInfeasible=0;
4093  int totalTightened = 0;
4094
4095  double tolerance = primalTolerance();
4096
4097
4098  // Save column bounds
4099  double * saveLower = new double [numberColumns_];
4100  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
4101  double * saveUpper = new double [numberColumns_];
4102  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
4103
4104  int iRow, iColumn;
4105
4106  // If wanted - tighten column bounds using solution
4107  if (factor) {
4108    double largest=0.0;
4109    if (factor>0.0) {
4110      assert (factor>1.0);
4111      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4112        if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
4113          largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
4114        }
4115      }
4116      largest *= factor;
4117    } else {
4118      // absolute
4119       largest = - factor;
4120    }
4121    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4122      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
4123        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
4124        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
4125      }
4126    }
4127  }
4128#define MAXPASS 10
4129
4130  // Loop round seeing if we can tighten bounds
4131  // Would be faster to have a stack of possible rows
4132  // and we put altered rows back on stack
4133  int numberCheck=-1;
4134  while(numberChanged>numberCheck) {
4135
4136    numberChanged = 0; // Bounds tightened this pass
4137   
4138    if (iPass==MAXPASS) break;
4139    iPass++;
4140   
4141    for (iRow = 0; iRow < numberRows_; iRow++) {
4142
4143      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
4144
4145        // possible row
4146        int infiniteUpper = 0;
4147        int infiniteLower = 0;
4148        double maximumUp = 0.0;
4149        double maximumDown = 0.0;
4150        double newBound;
4151        CoinBigIndex rStart = rowStart[iRow];
4152        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
4153        CoinBigIndex j;
4154        // Compute possible lower and upper ranges
4155     
4156        for (j = rStart; j < rEnd; ++j) {
4157          double value=element[j];
4158          iColumn = column[j];
4159          if (value > 0.0) {
4160            if (columnUpper_[iColumn] >= large) {
4161              ++infiniteUpper;
4162            } else {
4163              maximumUp += columnUpper_[iColumn] * value;
4164            }
4165            if (columnLower_[iColumn] <= -large) {
4166              ++infiniteLower;
4167            } else {
4168              maximumDown += columnLower_[iColumn] * value;
4169            }
4170          } else if (value<0.0) {
4171            if (columnUpper_[iColumn] >= large) {
4172              ++infiniteLower;
4173            } else {
4174              maximumDown += columnUpper_[iColumn] * value;
4175            }
4176            if (columnLower_[iColumn] <= -large) {
4177              ++infiniteUpper;
4178            } else {
4179              maximumUp += columnLower_[iColumn] * value;
4180            }
4181          }
4182        }
4183        // Build in a margin of error
4184        maximumUp += 1.0e-8*fabs(maximumUp);
4185        maximumDown -= 1.0e-8*fabs(maximumDown);
4186        double maxUp = maximumUp+infiniteUpper*1.0e31;
4187        double maxDown = maximumDown-infiniteLower*1.0e31;
4188        if (maxUp <= rowUpper_[iRow] + tolerance && 
4189            maxDown >= rowLower_[iRow] - tolerance) {
4190         
4191          // Row is redundant - make totally free
4192          // NO - can't do this for postsolve
4193          // rowLower_[iRow]=-COIN_DBL_MAX;
4194          // rowUpper_[iRow]=COIN_DBL_MAX;
4195          //printf("Redundant row in presolveX %d\n",iRow);
4196
4197        } else {
4198          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
4199              maxDown > rowUpper_[iRow]+100.0*tolerance) {
4200            // problem is infeasible - exit at once
4201            numberInfeasible++;
4202            break;
4203          }
4204          double lower = rowLower_[iRow];
4205          double upper = rowUpper_[iRow];
4206          for (j = rStart; j < rEnd; ++j) {
4207            double value=element[j];
4208            iColumn = column[j];
4209            double nowLower = columnLower_[iColumn];
4210            double nowUpper = columnUpper_[iColumn];
4211            if (value > 0.0) {
4212              // positive value
4213              if (lower>-large) {
4214                if (!infiniteUpper) {
4215                  assert(nowUpper < large2);
4216                  newBound = nowUpper + 
4217                    (lower - maximumUp) / value;
4218                  // relax if original was large
4219                  if (fabs(maximumUp)>1.0e8)
4220                    newBound -= 1.0e-12*fabs(maximumUp);
4221                } else if (infiniteUpper==1&&nowUpper>large) {
4222                  newBound = (lower -maximumUp) / value;
4223                  // relax if original was large
4224                  if (fabs(maximumUp)>1.0e8)
4225                    newBound -= 1.0e-12*fabs(maximumUp);
4226                } else {
4227                  newBound = -COIN_DBL_MAX;
4228                }
4229                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
4230                  // Tighten the lower bound
4231                  numberChanged++;
4232                  // check infeasible (relaxed)
4233                  if (nowUpper < newBound) { 
4234                    if (nowUpper - newBound < 
4235                        -100.0*tolerance) 
4236                      numberInfeasible++;
4237                    else 
4238                      newBound=nowUpper;
4239                  }
4240                  columnLower_[iColumn] = newBound;
4241                  // adjust
4242                  double now;
4243                  if (nowLower<-large) {
4244                    now=0.0;
4245                    infiniteLower--;
4246                  } else {
4247                    now = nowLower;
4248                  }
4249                  maximumDown += (newBound-now) * value;
4250                  nowLower = newBound;
4251#ifdef DEBUG
4252                  checkCorrect(this,iRow,
4253                               element, rowStart, rowLength,
4254                               column,
4255                               columnLower_,  columnUpper_,
4256                               infiniteUpper,
4257                               infiniteLower,
4258                               maximumUp,
4259                               maximumDown);
4260#endif
4261                }
4262              } 
4263              if (upper <large) {
4264                if (!infiniteLower) {
4265                  assert(nowLower >- large2);
4266                  newBound = nowLower + 
4267                    (upper - maximumDown) / value;
4268                  // relax if original was large
4269                  if (fabs(maximumDown)>1.0e8)
4270                    newBound += 1.0e-12*fabs(maximumDown);
4271                } else if (infiniteLower==1&&nowLower<-large) {
4272                  newBound =   (upper - maximumDown) / value;
4273                  // relax if original was large
4274                  if (fabs(maximumDown)>1.0e8)
4275                    newBound += 1.0e-12*fabs(maximumDown);
4276                } else {
4277                  newBound = COIN_DBL_MAX;
4278                }
4279                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
4280                  // Tighten the upper bound
4281                  numberChanged++;
4282                  // check infeasible (relaxed)
4283                  if (nowLower > newBound) { 
4284                    if (newBound - nowLower < 
4285                        -100.0*tolerance) 
4286                      numberInfeasible++;
4287                    else 
4288                      newBound=nowLower;
4289                  }
4290                  columnUpper_[iColumn] = newBound;
4291                  // adjust
4292                  double now;
4293                  if (nowUpper>large) {
4294                    now=0.0;
4295                    infiniteUpper--;
4296                  } else {
4297                    now = nowUpper;
4298                  }
4299                  maximumUp += (newBound-now) * value;
4300                  nowUpper = newBound;
4301#ifdef DEBUG
4302                  checkCorrect(this,iRow,
4303                               element, rowStart, rowLength,
4304                               column,
4305                               columnLower_,  columnUpper_,
4306                               infiniteUpper,
4307                               infiniteLower,
4308                               maximumUp,
4309                               maximumDown);
4310#endif
4311                }
4312              }
4313            } else {
4314              // negative value
4315              if (lower>-large) {
4316                if (!infiniteUpper) {
4317                  assert(nowLower < large2);
4318                  newBound = nowLower + 
4319                    (lower - maximumUp) / value;
4320                  // relax if original was large
4321                  if (fabs(maximumUp)>1.0e8)
4322                    newBound += 1.0e-12*fabs(maximumUp);
4323                } else if (infiniteUpper==1&&nowLower<-large) {
4324                  newBound = (lower -maximumUp) / value;
4325                  // relax if original was large
4326                  if (fabs(maximumUp)>1.0e8)
4327                    newBound += 1.0e-12*fabs(maximumUp);
4328                } else {
4329                  newBound = COIN_DBL_MAX;
4330                }
4331                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
4332                  // Tighten the upper bound
4333                  numberChanged++;
4334                  // check infeasible (relaxed)
4335                  if (nowLower > newBound) { 
4336                    if (newBound - nowLower < 
4337                        -100.0*tolerance) 
4338                      numberInfeasible++;
4339                    else 
4340                      newBound=nowLower;
4341                  }
4342                  columnUpper_[iColumn] = newBound;
4343                  // adjust
4344                  double now;
4345                  if (nowUpper>large) {
4346                    now=0.0;
4347                    infiniteLower--;
4348                  } else {
4349                    now = nowUpper;
4350                  }
4351                  maximumDown += (newBound-now) * value;
4352                  nowUpper = newBound;
4353#ifdef DEBUG
4354                  checkCorrect(this,iRow,
4355                               element, rowStart, rowLength,
4356                               column,
4357                               columnLower_,  columnUpper_,
4358                               infiniteUpper,
4359                               infiniteLower,
4360                               maximumUp,
4361                               maximumDown);
4362#endif
4363                }
4364              }
4365              if (upper <large) {
4366                if (!infiniteLower) {
4367                  assert(nowUpper < large2);
4368                  newBound = nowUpper + 
4369                    (upper - maximumDown) / value;
4370                  // relax if original was large
4371                  if (fabs(maximumDown)>1.0e8)
4372                    newBound -= 1.0e-12*fabs(maximumDown);
4373                } else if (infiniteLower==1&&nowUpper>large) {
4374                  newBound =   (upper - maximumDown) / value;
4375                  // relax if original was large
4376                  if (fabs(maximumDown)>1.0e8)
4377                    newBound -= 1.0e-12*fabs(maximumDown);
4378                } else {
4379                  newBound = -COIN_DBL_MAX;
4380                }
4381                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
4382                  // Tighten the lower bound
4383                  numberChanged++;
4384                  // check infeasible (relaxed)
4385                  if (nowUpper < newBound) { 
4386                    if (nowUpper - newBound < 
4387                        -100.0*tolerance) 
4388                      numberInfeasible++;
4389                    else 
4390                      newBound=nowUpper;
4391                  }
4392                  columnLower_[iColumn] = newBound;
4393                  // adjust
4394                  double now;
4395                  if (nowLower<-large) {
4396                    now=0.0;
4397                    infiniteUpper--;
4398                  } else {
4399                    now = nowLower;
4400                  }
4401                  maximumUp += (newBound-now) * value;
4402                  nowLower = newBound;
4403#ifdef DEBUG
4404                  checkCorrect(this,iRow,
4405                               element, rowStart, rowLength,
4406                               column,
4407                               columnLower_,  columnUpper_,
4408                               infiniteUpper,
4409                               infiniteLower,
4410                               maximumUp,
4411                               maximumDown);
4412#endif
4413                }
4414              }
4415            }
4416          }
4417        }
4418      }
4419    }
4420    totalTightened += numberChanged;
4421    if (iPass==1)
4422      numberCheck=numberChanged>>4;
4423    if (numberInfeasible) break;
4424  }
4425  if (!numberInfeasible) {
4426    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
4427      <<totalTightened
4428      <<CoinMessageEol;
4429    // Set bounds slightly loose
4430    double useTolerance = 1.0e-3;
4431    if (doTight>0) {
4432      if (doTight>10) { 
4433        useTolerance=0.0;
4434      } else {
4435        while (doTight) {
4436          useTolerance *= 0.1;
4437          doTight--;
4438        }
4439      }
4440    }
4441    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4442      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
4443        // Make large bounds stay infinite
4444        if (saveUpper[iColumn]>1.0e30&&columnUpper_[iColumn]>1.0e10) {
4445          columnUpper_[iColumn]=COIN_DBL_MAX;
4446        }
4447        if (saveLower[iColumn]<-1.0e30&&columnLower_[iColumn]<-1.0e10) {
4448          columnLower_[iColumn]=-COIN_DBL_MAX;
4449        }
4450        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e-8) {
4451          // relax enough so will have correct dj
4452#if 1
4453          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
4454                                    columnLower_[iColumn]-100.0*useTolerance);
4455          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
4456                                    columnUpper_[iColumn]+100.0*useTolerance);
4457#else
4458          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
4459            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
4460              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
4461            } else {
4462              columnLower_[iColumn]=saveLower[iColumn];
4463              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
4464                                        saveLower[iColumn]+100.0*useTolerance);
4465            }
4466          } else {
4467            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
4468              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
4469            } else {
4470              columnUpper_[iColumn]=saveUpper[iColumn];
4471              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
4472                                        saveUpper[iColumn]-100.0*useTolerance);
4473            }
4474          }
4475#endif
4476        } else {
4477          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
4478            // relax a bit
4479            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
4480                                        saveUpper[iColumn]);
4481          }
4482          if (columnLower_[iColumn]>saveLower[iColumn]) {
4483            // relax a bit
4484            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
4485                                        saveLower[iColumn]);
4486          }
4487        }
4488      }
4489    }
4490  } else {
4491    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
4492      <<numberInfeasible
4493      <<CoinMessageEol;
4494    // restore column bounds
4495    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
4496    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
4497  }
4498  delete [] saveLower;
4499  delete [] saveUpper;
4500  return (numberInfeasible);
4501}
4502//#define SAVE_AND_RESTORE
4503// dual
4504#include "ClpSimplexDual.hpp"
4505#include "ClpSimplexPrimal.hpp"
4506#ifndef SAVE_AND_RESTORE
4507int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
4508#else
4509int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
4510{
4511  // May be empty problem
4512  if (numberRows_&&numberColumns_) {
4513    // Save on file for debug
4514    int returnCode;
4515    returnCode = saveModel("debug.sav");
4516    if (returnCode) {
4517      printf("** Unable to save model to debug.sav\n");
4518      abort();
4519    }
4520    ClpSimplex temp;
4521    returnCode=temp.restoreModel("debug.sav");
4522    if (returnCode) {
4523      printf("** Unable to restore model from debug.sav\n");
4524      abort();
4525    }
4526    temp.setLogLevel(handler_->logLevel());
4527    // Now do dual
4528    returnCode=temp.dualDebug(ifValuesPass,startFinishOptions);
4529    // Move status and solution back
4530    int numberTotal = numberRows_+numberColumns_;
4531    memcpy(status_,temp.statusArray(),numberTotal);
4532    memcpy(columnActivity_,temp.primalColumnSolution(),numberColumns_*sizeof(double));
4533    memcpy(rowActivity_,temp.primalRowSolution(),numberRows_*sizeof(double));
4534    memcpy(reducedCost_,temp.dualColumnSolution(),numberColumns_*sizeof(double));
4535    memcpy(dual_,temp.dualRowSolution(),numberRows_*sizeof(double));
4536    problemStatus_ = temp.problemStatus_;
4537    setObjectiveValue(temp.objectiveValue());
4538    setSumDualInfeasibilities(temp.sumDualInfeasibilities());
4539    setNumberDualInfeasibilities(temp.numberDualInfeasibilities());
4540    setSumPrimalInfeasibilities(temp.sumPrimalInfeasibilities());
4541    setNumberPrimalInfeasibilities(temp.numberPrimalInfeasibilities());
4542    setNumberIterations(temp.numberIterations());
4543    return returnCode;
4544  } else {
4545    // empty
4546    return dualDebug(ifValuesPass,startFinishOptions);
4547  }
4548}
4549int ClpSimplex::dualDebug (int ifValuesPass , int startFinishOptions)
4550#endif
4551{
4552  //double savedPivotTolerance = factorization_->pivotTolerance();
4553  int saveQuadraticActivated = objective_->activated();
4554  objective_->setActivated(0);
4555  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
4556  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
4557      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
4558      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
4559      additional data and have no destructor or (non-default) constructor.
4560
4561      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4562      in primal and dual.
4563
4564      As far as I can see this is perfectly safe.
4565  */
4566  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass, startFinishOptions);
4567  if ((specialOptions_&2048)!=0&&problemStatus_==10&&!numberPrimalInfeasibilities_
4568      &&sumDualInfeasibilities_<1000.0*dualTolerance_&&perturbation_>=100)
4569    problemStatus_=0; // ignore
4570  if (problemStatus_==10) {
4571    //printf("Cleaning up with primal\n");
4572    int savePerturbation = perturbation_;
4573    int saveLog = handler_->logLevel();
4574    //handler_->setLogLevel(63);
4575    perturbation_=100;
4576    bool denseFactorization = initialDenseFactorization();
4577    // It will be safe to allow dense
4578    setInitialDenseFactorization(true);
4579    // Allow for catastrophe
4580    int saveMax = intParam_[ClpMaxNumIteration];
4581    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
4582      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
4583    // check which algorithms allowed
4584    int dummy;
4585    if (problemStatus_==10)
4586      startFinishOptions |= 2;
4587    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
4588      returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
4589    else
4590      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
4591    if (problemStatus_==3&&numberIterations_<saveMax) {
4592      if (handler_->logLevel()==63)
4593        printf("looks like trouble - too many iterations in clean up - trying again\n");
4594      // flatten solution and try again
4595      int iRow,iColumn;
4596      for (iRow=0;iRow<numberRows_;iRow++) {
4597        if (getRowStatus(iRow)!=basic) {
4598          setRowStatus(iRow,superBasic);
4599          // but put to bound if close
4600          if (fabs(rowActivity_[iRow]-rowLower_[iRow])
4601              <=primalTolerance_) {
4602            rowActivity_[iRow]=rowLower_[iRow];
4603            setRowStatus(iRow,atLowerBound);
4604          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])
4605                     <=primalTolerance_) {
4606            rowActivity_[iRow]=rowUpper_[iRow];
4607            setRowStatus(iRow,atUpperBound);
4608          }
4609        }
4610      }
4611      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4612        if (getColumnStatus(iColumn)!=basic) {
4613          setColumnStatus(iColumn,superBasic);
4614          // but put to bound if close
4615          if (fabs(columnActivity_[iColumn]-columnLower_[iColumn])
4616              <=primalTolerance_) {
4617            columnActivity_[iColumn]=columnLower_[iColumn];
4618            setColumnStatus(iColumn,atLowerBound);
4619          } else if (fabs(columnActivity_[iColumn]
4620                          -columnUpper_[iColumn])
4621                     <=primalTolerance_) {
4622            columnActivity_[iColumn]=columnUpper_[iColumn];
4623            setColumnStatus(iColumn,atUpperBound);
4624          }
4625        }
4626      }
4627      problemStatus_=-1;
4628      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 
4629                                          2*numberRows_+numberColumns_,saveMax);
4630      perturbation_=savePerturbation;
4631      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
4632      if (problemStatus_==3&&numberIterations_<saveMax&& 
4633          handler_->logLevel()>0)
4634        printf("looks like real trouble - too many iterations in second clean up - giving up\n");
4635    }
4636    intParam_[ClpMaxNumIteration] = saveMax;
4637
4638    setInitialDenseFactorization(denseFactorization);
4639    perturbation_=savePerturbation;
4640    if (problemStatus_==10) { 
4641      if (!numberPrimalInfeasibilities_)
4642        problemStatus_=0;
4643      else
4644        problemStatus_=4;
4645    }
4646    handler_->setLogLevel(saveLog);
4647  }
4648  objective_->setActivated(saveQuadraticActivated);
4649  //factorization_->pivotTolerance(savedPivotTolerance);
4650  return returnCode;
4651}
4652// primal
4653int ClpSimplex::primal (int ifValuesPass , int startFinishOptions)
4654{
4655  //double savedPivotTolerance = factorization_->pivotTolerance();
4656#ifndef SLIM_CLP
4657  // See if nonlinear
4658  if (objective_->type()>1&&objective_->activated()) 
4659    return reducedGradient();
4660#endif 
4661  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
4662  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
4663      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
4664      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
4665      additional data and have no destructor or (non-default) constructor.
4666
4667      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4668      in primal and dual.
4669
4670      As far as I can see this is perfectly safe.
4671  */
4672  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass,startFinishOptions);
4673  if (problemStatus_==10) {
4674    //printf("Cleaning up with dual\n");
4675    int savePerturbation = perturbation_;
4676    perturbation_=100;
4677    bool denseFactorization = initialDenseFactorization();
4678    // It will be safe to allow dense
4679    setInitialDenseFactorization(true);
4680    // check which algorithms allowed
4681    int dummy;
4682    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0&&(specialOptions_&8192)==0) {
4683      double saveBound = dualBound_;
4684      // upperOut_ has largest away from bound
4685      dualBound_=CoinMin(CoinMax(2.0*upperOut_,1.0e8),dualBound_);
4686      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
4687      dualBound_=saveBound;
4688    } else {
4689      returnCode = ((ClpSimplexPrimal *) this)->primal(0,startFinishOptions);
4690    }
4691    setInitialDenseFactorization(denseFactorization);
4692    perturbation_=savePerturbation;
4693    if (problemStatus_==10) 
4694      problemStatus_=0;
4695  }
4696  //factorization_->pivotTolerance(savedPivotTolerance);
4697  return returnCode;
4698}
4699#ifndef SLIM_CLP
4700#include "ClpQuadraticObjective.hpp"
4701/* Dual ranging.
4702   This computes increase/decrease in cost for each given variable and corresponding
4703   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
4704   and numberColumns.. for artificials/slacks.
4705   For non-basic variables the sequence number will be that of the non-basic variables.
4706   
4707   Up to user to provide correct length arrays.
4708   
4709   Returns non-zero if infeasible unbounded etc
4710*/
4711#include "ClpSimplexOther.hpp"
4712int ClpSimplex::dualRanging(int numberCheck,const int * which,
4713                            double * costIncrease, int * sequenceIncrease,
4714                            double * costDecrease, int * sequenceDecrease,
4715                            double * valueIncrease, double * valueDecrease)
4716{
4717  int savePerturbation = perturbation_;
4718  perturbation_=100;
4719  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4720  if (problemStatus_==10) {
4721    //printf("Cleaning up with dual\n");
4722    bool denseFactorization = initialDenseFactorization();
4723    // It will be safe to allow dense
4724    setInitialDenseFactorization(true);
4725    // check which algorithms allowed
4726    int dummy;
4727    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
4728      // upperOut_ has largest away from bound
4729      double saveBound=dualBound_;
4730      if (upperOut_>0.0)
4731        dualBound_=2.0*upperOut_;
4732      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
4733      dualBound_=saveBound;
4734    } else {
4735      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4736    }
4737    setInitialDenseFactorization(denseFactorization);
4738    if (problemStatus_==10) 
4739      problemStatus_=0;
4740  }
4741  perturbation_=savePerturbation;
4742  if (problemStatus_||secondaryStatus_==6) {
4743    finish(); // get rid of arrays
4744    return 1; // odd status
4745  }
4746  ((ClpSimplexOther *) this)->dualRanging(numberCheck,which,
4747                                          costIncrease,sequenceIncrease,
4748                                          costDecrease,sequenceDecrease,
4749                                          valueIncrease, valueDecrease);
4750  finish(); // get rid of arrays
4751  return 0;
4752}
4753/* Primal ranging.
4754   This computes increase/decrease in value for each given variable and corresponding
4755   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
4756   and numberColumns.. for artificials/slacks.
4757   For basic variables the sequence number will be that of the basic variables.
4758   
4759   Up to user to providen correct length arrays.
4760   
4761   Returns non-zero if infeasible unbounded etc
4762*/
4763int ClpSimplex::primalRanging(int numberCheck,const int * which,
4764                  double * valueIncrease, int * sequenceIncrease,
4765                  double * valueDecrease, int * sequenceDecrease)
4766{
4767  int savePerturbation = perturbation_;
4768  perturbation_=100;
4769  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4770  if (problemStatus_==10) {
4771    //printf("Cleaning up with dual\n");
4772    bool denseFactorization = initialDenseFactorization();
4773    // It will be safe to allow dense
4774    setInitialDenseFactorization(true);
4775    // check which algorithms allowed
4776    int dummy;
4777    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
4778      // upperOut_ has largest away from bound
4779      double saveBound=dualBound_;
4780      if (upperOut_>0.0)
4781        dualBound_=2.0*upperOut_;
4782      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
4783      dualBound_=saveBound;
4784    } else {
4785      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4786    }
4787    setInitialDenseFactorization(denseFactorization);
4788    if (problemStatus_==10) 
4789      problemStatus_=0;
4790  }
4791  perturbation_=savePerturbation;
4792  if (problemStatus_||secondaryStatus_==6) {
4793    finish(); // get rid of arrays
4794    return 1; // odd status
4795  }
4796  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
4797                                          valueIncrease,sequenceIncrease,
4798                                          valueDecrease,sequenceDecrease);
4799  finish(); // get rid of arrays
4800  return 0;
4801}
4802/* Write the basis in MPS format to the specified file.
4803   If writeValues true writes values of structurals
4804   (and adds VALUES to end of NAME card)
4805   
4806   Row and column names may be null.
4807   formatType is
4808   <ul>
4809   <li> 0 - normal
4810   <li> 1 - extra accuracy
4811   <li> 2 - IEEE hex (later)
4812   </ul>
4813   
4814   Returns non-zero on I/O error
4815*/
4816int 
4817ClpSimplex::writeBasis(const char *filename,
4818                            bool writeValues,
4819                            int formatType) const
4820{
4821  return ((const ClpSimplexOther *) this)->writeBasis(filename,writeValues,
4822                                         formatType);
4823}
4824// Read a basis from the given filename
4825int 
4826ClpSimplex::readBasis(const char *filename)
4827{
4828  return ((ClpSimplexOther *) this)->readBasis(filename);
4829}
4830#include "ClpSimplexNonlinear.hpp"
4831/* Solves nonlinear problem using SLP - may be used as crash
4832   for other algorithms when number of iterations small
4833*/
4834int 
4835ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
4836{
4837  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
4838}
4839// Solves non-linear using reduced gradient
4840int ClpSimplex::reducedGradient(int phase)
4841{
4842  if (objective_->type()<2||!objective_->activated()) {
4843    // no quadratic part
4844    return primal(0);
4845  }
4846  // get feasible
4847  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
4848    objective_->setActivated(0);
4849    double saveDirection = optimizationDirection();
4850    setOptimizationDirection(0.0);
4851    primal(1);
4852    setOptimizationDirection(saveDirection);
4853    objective_->setActivated(1);
4854    // still infeasible
4855    if (numberPrimalInfeasibilities())
4856      return 0;
4857  }
4858  // Now enter method
4859  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
4860  return returnCode;
4861}
4862#include "ClpPredictorCorrector.hpp"
4863#include "ClpCholeskyBase.hpp"
4864// Preference is WSSMP, TAUCS, UFL (just ordering) then base
4865#if WSSMP_BARRIER
4866#include "ClpCholeskyWssmp.hpp"
4867#include "ClpCholeskyWssmpKKT.hpp"
4868#elif UFL_BARRIER
4869#include "ClpCholeskyUfl.hpp"
4870#elif TAUCS_BARRIER
4871#include "ClpCholeskyTaucs.hpp"
4872#endif
4873#include "ClpPresolve.hpp"
4874/* Solves using barrier (assumes you have good cholesky factor code).
4875   Does crossover to simplex if asked*/
4876int 
4877ClpSimplex::barrier(bool crossover)
4878{
4879  ClpSimplex * model2 = this;
4880  int savePerturbation=perturbation_;
4881  ClpInterior barrier;
4882  barrier.borrowModel(*model2);
4883  // See if quadratic objective
4884  ClpQuadraticObjective * quadraticObj = NULL;
4885  if (objective_->type()==2)
4886    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
4887  // If Quadratic we need KKT
4888  bool doKKT = (quadraticObj!=NULL);
4889  // Preference is WSSMP, UFL, TAUCS then base
4890#ifdef WSSMP_BARRIER
4891 if (!doKKT) {
4892   ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
4893   barrier.setCholesky(cholesky);
4894 } else {
4895  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
4896   ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
4897   barrier.setCholesky(cholesky);
4898 }
4899#elif UFL_BARRIER
4900 if (!doKKT) {
4901   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
4902   barrier.setCholesky(cholesky);
4903 } else {
4904   ClpCholeskyBase * cholesky = new ClpCholeskyBase();
4905   // not yetClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
4906   cholesky->setKKT(true);
4907   barrier.setCholesky(cholesky);
4908 }
4909#elif TAUCS_BARRIER
4910 assert (!doKKT);
4911 ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
4912 barrier.setCholesky(cholesky);
4913#else
4914 if (!doKKT) {
4915  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
4916  barrier.setCholesky(cholesky);
4917 } else {
4918   ClpCholeskyBase * cholesky = new ClpCholeskyBase();
4919   cholesky->setKKT(true);
4920   barrier.setCholesky(cholesky);
4921 }
4922#endif
4923  barrier.setDiagonalPerturbation(1.0e-14);
4924  int numberRows = model2->numberRows();
4925  int numberColumns = model2->numberColumns();
4926  int saveMaxIts = model2->maximumIterations();
4927  if (saveMaxIts<1000) {
4928    barrier.setMaximumBarrierIterations(saveMaxIts);
4929    model2->setMaximumIterations(1000000);
4930  }
4931  barrier.primalDual();
4932  int barrierStatus=barrier.status();
4933  double gap = barrier.complementarityGap();
4934  // get which variables are fixed
4935  double * saveLower=NULL;
4936  double * saveUpper=NULL;
4937  ClpPresolve pinfo2;
4938  ClpSimplex * saveModel2=NULL;
4939  int numberFixed = barrier.numberFixed();
4940  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
4941    // may as well do presolve
4942    int numberRows = barrier.numberRows();
4943    int numberColumns = barrier.numberColumns();
4944    int numberTotal = numberRows+numberColumns;
4945    saveLower = new double [numberTotal];
4946    saveUpper = new double [numberTotal];
4947    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
4948    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
4949    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
4950    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
4951    barrier.fixFixed();
4952    saveModel2=model2;
4953  }
4954  barrier.returnModel(*model2);
4955  double * rowPrimal = new double [numberRows];
4956  double * columnPrimal = new double [numberColumns];
4957  double * rowDual = new double [numberRows];
4958  double * columnDual = new double [numberColumns];
4959  // move solutions other way
4960  CoinMemcpyN(model2->primalRowSolution(),
4961              numberRows,rowPrimal);
4962  CoinMemcpyN(model2->dualRowSolution(),
4963              numberRows,rowDual);
4964  CoinMemcpyN(model2->primalColumnSolution(),
4965              numberColumns,columnPrimal);
4966  CoinMemcpyN(model2->dualColumnSolution(),
4967              numberColumns,columnDual);
4968  if (saveModel2) {
4969    // do presolve
4970    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
4971                                   false,5,true);
4972  }
4973  if (barrierStatus<4&&crossover) {
4974    // make sure no status left
4975    model2->createStatus();
4976    // solve
4977    model2->setPerturbation(100);
4978    // throw some into basis
4979    {
4980      int numberRows = model2->numberRows();
4981      int numberColumns = model2->numberColumns();
4982      double * dsort = new double[numberColumns];
4983      int * sort = new int[numberColumns];
4984      int n=0;
4985      const double * columnLower = model2->columnLower();
4986      const double * columnUpper = model2->columnUpper();
4987      const double * primalSolution = model2->primalColumnSolution();
4988      double tolerance = 10.0*primalTolerance_;
4989      int i;
4990      for ( i=0;i<numberRows;i++) 
4991        model2->setRowStatus(i,superBasic);
4992      for ( i=0;i<numberColumns;i++) {
4993        double distance = CoinMin(columnUpper[i]-primalSolution[i],
4994                              primalSolution[i]-columnLower[i]);
4995        if (distance>tolerance) {
4996          dsort[n]=-distance;
4997          sort[n++]=i;
4998          model2->setStatus(i,superBasic);
4999        } else if (distance>primalTolerance_) {
5000          model2->setStatus(i,superBasic);
5001        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
5002          model2->setStatus(i,atLowerBound);
5003        } else {
5004          model2->setStatus(i,atUpperBound);
5005        }
5006      }
5007      CoinSort_2(dsort,dsort+n,sort);
5008      n = CoinMin(numberRows,n);
5009      for ( i=0;i<n;i++) {
5010        int iColumn = sort[i];
5011        model2->setStatus(iColumn,basic);
5012      }
5013      delete [] sort;
5014      delete [] dsort;
5015    }
5016    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
5017      int numberRows = model2->numberRows();
5018      int numberColumns = model2->numberColumns();
5019      // just primal values pass
5020      double saveScale = model2->objectiveScale();
5021      model2->setObjectiveScale(1.0e-3);
5022      model2->primal(2);
5023      model2->setObjectiveScale(saveScale);
5024      // save primal solution and copy back dual
5025      CoinMemcpyN(model2->primalRowSolution(),
5026                  numberRows,rowPrimal);
5027      CoinMemcpyN(rowDual,
5028                  numberRows,model2->dualRowSolution());
5029      CoinMemcpyN(model2->primalColumnSolution(),
5030                  numberColumns,columnPrimal);
5031      CoinMemcpyN(columnDual,
5032                  numberColumns,model2->dualColumnSolution());
5033      //model2->primal(1);
5034      // clean up reduced costs and flag variables
5035      {
5036        double * dj = model2->dualColumnSolution();
5037        double * cost = model2->objective();
5038        double * saveCost = new double[numberColumns];
5039        memcpy(saveCost,cost,numberColumns*sizeof(double));
5040        double * saveLower = new double[numberColumns];
5041        double * lower = model2->columnLower();
5042        memcpy(saveLower,lower,numberColumns*sizeof(double));
5043        double * saveUpper = new double[numberColumns];
5044        double * upper = model2->columnUpper();
5045        memcpy(saveUpper,upper,numberColumns*sizeof(double));
5046        int i;
5047        double tolerance = 10.0*dualTolerance_;
5048        for ( i=0;i<numberColumns;i++) {
5049          if (model2->getStatus(i)==basic) {
5050            dj[i]=0.0;
5051          } else if (model2->getStatus(i)==atLowerBound) {
5052            if (optimizationDirection_*dj[i]<tolerance) {
5053              if (optimizationDirection_*dj[i]<0.0) {
5054                //if (dj[i]<-1.0e-3)
5055                //printf("bad dj at lb %d %g\n",i,dj[i]);
5056                cost[i] -= dj[i];
5057                dj[i]=0.0;
5058              }
5059            } else {
5060              upper[i]=lower[i];
5061            }
5062          } else if (model2->getStatus(i)==atUpperBound) {
5063            if (optimizationDirection_*dj[i]>tolerance) {
5064              if (optimizationDirection_*dj[i]>0.0) {
5065                //if (dj[i]>1.0e-3)
5066                //printf("bad dj at ub %d %g\n",i,dj[i]);
5067                cost[i] -= dj[i];
5068                dj[i]=0.0;
5069              }
5070            } else {
5071              lower[i]=upper[i];
5072            }
5073          }
5074        }
5075        // just dual values pass
5076        //model2->setLogLevel(63);
5077        //model2->setFactorizationFrequency(1);
5078        model2->dual(2);
5079        memcpy(cost,saveCost,numberColumns*sizeof(double));
5080        delete [] saveCost;
5081        memcpy(lower,saveLower,numberColumns*sizeof(double));
5082        delete [] saveLower;
5083        memcpy(upper,saveUpper,numberColumns*sizeof(double));
5084        delete [] saveUpper;
5085      }
5086      // and finish
5087      // move solutions
5088      CoinMemcpyN(rowPrimal,
5089                  numberRows,model2->primalRowSolution());
5090      CoinMemcpyN(columnPrimal,
5091                  numberColumns,model2->primalColumnSolution());
5092    }
5093//     double saveScale = model2->objectiveScale();
5094//     model2->setObjectiveScale(1.0e-3);
5095//     model2->primal(2);
5096//    model2->setObjectiveScale(saveScale);
5097    model2->primal(1);
5098  } else if (barrierStatus==4&&crossover) {
5099    // memory problems
5100    model2->setPerturbation(savePerturbation);
5101    model2->createStatus();
5102    model2->dual();
5103  }
5104  model2->setMaximumIterations(saveMaxIts);
5105  delete [] rowPrimal;
5106  delete [] columnPrimal;
5107  delete [] rowDual;
5108  delete [] columnDual;
5109  if (saveLower) {
5110    pinfo2.postsolve(true);
5111    delete model2;
5112    model2=saveModel2;
5113    int numberRows = model2->numberRows();
5114    int numberColumns = model2->numberColumns();
5115    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
5116    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
5117    delete [] saveLower;
5118    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
5119    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
5120    delete [] saveUpper;
5121    model2->primal(1);
5122  }
5123  model2->setPerturbation(savePerturbation);
5124  return model2->status();
5125}
5126/* For strong branching.  On input lower and upper are new bounds
5127   while on output they are objective function values (>1.0e50 infeasible).
5128   Return code is 0 if nothing interesting, -1 if infeasible both
5129   ways and +1 if infeasible one way (check values to see which one(s))
5130*/
5131int ClpSimplex::strongBranching(int numberVariables,const int * variables,
5132                                double * newLower, double * newUpper,
5133                                double ** outputSolution,
5134                                int * outputStatus, int * outputIterations,
5135                                bool stopOnFirstInfeasible,
5136                                bool alwaysFinish,
5137                                int startFinishOptions)
5138{
5139  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
5140                                                    newLower,  newUpper,outputSolution,
5141                                                    outputStatus, outputIterations,
5142                                                    stopOnFirstInfeasible,
5143                                                    alwaysFinish,startFinishOptions);
5144}
5145#endif
5146/* Borrow model.  This is so we dont have to copy large amounts
5147   of data around.  It assumes a derived class wants to overwrite
5148   an empty model with a real one - while it does an algorithm.
5149   This is same as ClpModel one, but sets scaling on etc. */
5150void 
5151ClpSimplex::borrowModel(ClpModel & otherModel) 
5152{
5153  ClpModel::borrowModel(otherModel);
5154  createStatus();
5155  //ClpDualRowSteepest steep1;
5156  //setDualRowPivotAlgorithm(steep1);
5157  //ClpPrimalColumnSteepest steep2;
5158  //setPrimalColumnPivotAlgorithm(steep2);
5159}
5160void 
5161ClpSimplex::borrowModel(ClpSimplex & otherModel) 
5162{
5163  ClpModel::borrowModel(otherModel);
5164  createStatus();
5165  dualBound_ = otherModel.dualBound_;
5166  dualTolerance_ = otherModel.dualTolerance_;
5167  primalTolerance_ = otherModel.primalTolerance_;
5168  delete dualRowPivot_;
5169  dualRowPivot_ = otherModel.dualRowPivot_->clone(true);
5170  delete primalColumnPivot_;
5171  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(true);
5172  perturbation_ = otherModel.perturbation_;
5173  specialOptions_ = otherModel.specialOptions_;
5174  automaticScale_ = otherModel.automaticScale_;
5175}
5176/// Saves scalars for ClpSimplex
5177typedef struct {
5178  double optimizationDirection;
5179  double dblParam[ClpLastDblParam];
5180  double objectiveValue;
5181  double dualBound;
5182  double dualTolerance;
5183  double primalTolerance;
5184  double sumDualInfeasibilities;
5185  double sumPrimalInfeasibilities;
5186  double infeasibilityCost;
5187  int numberRows;
5188  int numberColumns;
5189  int intParam[ClpLastIntParam];
5190  int numberIterations;
5191  int problemStatus;
5192  int maximumIterations;
5193  int lengthNames;
5194  int numberDualInfeasibilities;
5195  int numberDualInfeasibilitiesWithoutFree;
5196  int numberPrimalInfeasibilities;
5197  int numberRefinements;
5198  int scalingFlag;
5199  int algorithm;
5200  unsigned int specialOptions;
5201  int dualPivotChoice;
5202  int primalPivotChoice;
5203  int matrixStorageChoice;
5204} Clp_scalars;
5205#ifndef SLIM_NOIO
5206int outDoubleArray(double * array, int length, FILE * fp)
5207{
5208  int numberWritten;
5209  if (array&&length) {
5210    numberWritten = fwrite(&length,sizeof(int),1,fp);
5211    if (numberWritten!=1)
5212      return 1;
5213    numberWritten = fwrite(array,sizeof(double),length,fp);
5214    if (numberWritten!=length)
5215      return 1;
5216  } else {
5217    length = 0;
5218    numberWritten = fwrite(&length,sizeof(int),1,fp);
5219    if (numberWritten!=1)
5220      return 1;
5221  }
5222  return 0;
5223}
5224// Save model to file, returns 0 if success
5225int
5226ClpSimplex::saveModel(const char * fileName)
5227{
5228  FILE * fp = fopen(fileName,"wb");
5229  if (fp) {
5230    Clp_scalars scalars;
5231    CoinBigIndex numberWritten;
5232    // Fill in scalars
5233    scalars.optimizationDirection = optimizationDirection_;
5234    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
5235    scalars.objectiveValue = objectiveValue_;
5236    scalars.dualBound = dualBound_;
5237    scalars.dualTolerance = dualTolerance_;
5238    scalars.primalTolerance = primalTolerance_;
5239    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
5240    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
5241    scalars.infeasibilityCost = infeasibilityCost_;
5242    scalars.numberRows = numberRows_;
5243    scalars.numberColumns = numberColumns_;
5244    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
5245    scalars.numberIterations = numberIterations_;
5246    scalars.problemStatus = problemStatus_;
5247    scalars.maximumIterations = maximumIterations();
5248    scalars.lengthNames = lengthNames_;
5249    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
5250    scalars.numberDualInfeasibilitiesWithoutFree
5251      = numberDualInfeasibilitiesWithoutFree_;
5252    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
5253    scalars.numberRefinements = numberRefinements_;
5254    scalars.scalingFlag = scalingFlag_;
5255    scalars.algorithm = algorithm_;
5256    scalars.specialOptions = specialOptions_;
5257    scalars.dualPivotChoice = dualRowPivot_->type();
5258    scalars.primalPivotChoice = primalColumnPivot_->type();
5259    scalars.matrixStorageChoice = matrix_->type();
5260
5261    // put out scalars
5262    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
5263    if (numberWritten!=1)
5264      return 1;
5265    CoinBigIndex length;
5266#ifndef CLP_NO_STD
5267    int i;
5268    // strings
5269    for (i=0;i<ClpLastStrParam;i++) {
5270      length = strParam_[i].size();
5271      numberWritten = fwrite(&length,sizeof(int),1,fp);
5272      if (numberWritten!=1)
5273        return 1;
5274      if (length) {
5275        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
5276        if (numberWritten!=1)
5277          return 1;
5278      }
5279    }
5280#endif
5281    // arrays - in no particular order
5282    if (outDoubleArray(rowActivity_,numberRows_,fp))
5283        return 1;
5284    if (outDoubleArray(columnActivity_,numberColumns_,fp))
5285        return 1;
5286    if (outDoubleArray(dual_,numberRows_,fp))
5287        return 1;
5288    if (outDoubleArray(reducedCost_,numberColumns_,fp))
5289        return 1;
5290    if (outDoubleArray(rowLower_,numberRows_,fp))
5291        return 1;
5292    if (outDoubleArray(rowUpper_,numberRows_,fp))
5293        return 1;
5294    if (outDoubleArray(objective(),numberColumns_,fp))
5295        return 1;
5296    if (outDoubleArray(rowObjective_,numberRows_,fp))
5297        return 1;
5298    if (outDoubleArray(columnLower_,numberColumns_,fp))
5299        return 1;
5300    if (outDoubleArray(columnUpper_,numberColumns_,fp))
5301        return 1;
5302    if (ray_) {
5303      if (problemStatus_==1)
5304        if (outDoubleArray(ray_,numberRows_,fp))
5305          return 1;
5306      else if (problemStatus_==2)
5307        if (outDoubleArray(ray_,numberColumns_,fp))
5308          return 1;
5309      else
5310        if (outDoubleArray(NULL,0,fp))
5311          return 1;
5312    } else {
5313      if (outDoubleArray(NULL,0,fp))
5314        return 1;
5315    }
5316    if (status_&&(numberRows_+numberColumns_)>0) {
5317      length = numberRows_+numberColumns_;
5318      numberWritten = fwrite(&length,sizeof(int),1,fp);
5319      if (numberWritten!=1)
5320        return 1;
5321      numberWritten = fwrite(status_,sizeof(char),length, fp);
5322      if (numberWritten!=length)
5323        return 1;
5324    } else {
5325      length = 0;
5326      numberWritten = fwrite(&length,sizeof(int),1,fp);
5327      if (numberWritten!=1)
5328        return 1;
5329    }
5330#ifndef CLP_NO_STD
5331    if (lengthNames_) {
5332      char * array = 
5333        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
5334      char * put = array;
5335      CoinAssert (numberRows_ == (int) rowNames_.size());
5336      for (i=0;i<numberRows_;i++) {
5337        assert((int)rowNames_[i].size()<=lengthNames_);
5338        strcpy(put,rowNames_[i].c_str());
5339        put += lengthNames_+1;
5340      }
5341      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
5342      if (numberWritten!=numberRows_)
5343        return 1;
5344      put=array;
5345      CoinAssert (numberColumns_ == (int) columnNames_.size());
5346      for (i=0;i<numberColumns_;i++) {
5347        assert((int) columnNames_[i].size()<=lengthNames_);
5348        strcpy(put,columnNames_[i].c_str());
5349        put += lengthNames_+1;
5350      }
5351      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
5352      if (numberWritten!=numberColumns_)
5353        return 1;
5354      delete [] array;
5355    }
5356#endif
5357    // integers
5358    if (integerType_) {
5359      int marker=1;
5360      fwrite(&marker,sizeof(int),1,fp);
5361      numberWritten = fwrite(integerType_,1,numberColumns_,fp);
5362      if (numberWritten!=numberColumns_)
5363        return 1;
5364    } else {
5365      int marker=0;
5366      fwrite(&marker,sizeof(int),1,fp);
5367    }
5368    // just standard type at present
5369    assert (matrix_->type()==1);
5370    CoinAssert (matrix_->getNumCols() == numberColumns_);
5371    CoinAssert (matrix_->getNumRows() == numberRows_);
5372    // we are going to save with gaps
5373    length = matrix_->getVectorStarts()[numberColumns_-1]
5374      + matrix_->getVectorLengths()[numberColumns_-1];
5375    numberWritten = fwrite(&length,sizeof(int),1,fp);
5376    if (numberWritten!=1)
5377      return 1;
5378    numberWritten = fwrite(matrix_->getElements(),
5379                           sizeof(double),length,fp);
5380    if (numberWritten!=length)
5381      return 1;
5382    numberWritten = fwrite(matrix_->getIndices(),
5383                           sizeof(int),length,fp);
5384    if (numberWritten!=length)
5385      return 1;
5386    numberWritten = fwrite(matrix_->getVectorStarts(),
5387                           sizeof(int),numberColumns_+1,fp);
5388    if (numberWritten!=numberColumns_+1)
5389      return 1;
5390    numberWritten = fwrite(matrix_->getVectorLengths(),
5391                           sizeof(int),numberColumns_,fp);
5392    if (numberWritten!=numberColumns_)
5393      return 1;
5394    // finished
5395    fclose(fp);
5396    return 0;
5397  } else {
5398    return -1;
5399  }
5400}
5401
5402int inDoubleArray(double * &array, int length, FILE * fp)
5403{
5404  int numberRead;
5405  int length2;
5406  numberRead = fread(&length2,sizeof(int),1,fp);
5407  if (numberRead!=1)
5408    return 1;
5409  if (length2) {
5410    // lengths must match
5411    if (length!=length2)
5412      return 2;
5413    array = new double[length];
5414    numberRead = fread(array,sizeof(double),length,fp);
5415    if (numberRead!=length)
5416      return 1;
5417  } 
5418  return 0;
5419}
5420/* Restore model from file, returns 0 if success,
5421   deletes current model */
5422int 
5423ClpSimplex::restoreModel(const char * fileName)
5424{
5425  FILE * fp = fopen(fileName,"rb");
5426  if (fp) {
5427    // Get rid of current model
5428    // save event handler in case already set
5429    ClpEventHandler * handler = eventHandler_->clone();
5430    ClpModel::gutsOfDelete();
5431    eventHandler_ = handler;
5432    gutsOfDelete(0);
5433    int i;
5434    for (i=0;i<6;i++) {
5435      rowArray_[i]=NULL;
5436      columnArray_[i]=NULL;
5437    }
5438    // get an empty factorization so we can set tolerances etc
5439    getEmptyFactorization();
5440    // Say sparse
5441    factorization_->sparseThreshold(1);
5442    Clp_scalars scalars;
5443    CoinBigIndex numberRead;
5444
5445    // get scalars
5446    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
5447    if (numberRead!=1)
5448      return 1;
5449    // Fill in scalars
5450    optimizationDirection_ = scalars.optimizationDirection;
5451    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
5452    objectiveValue_ = scalars.objectiveValue;
5453    dualBound_ = scalars.dualBound;
5454    dualTolerance_ = scalars.dualTolerance;
5455    primalTolerance_ = scalars.primalTolerance;
5456    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
5457    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
5458    infeasibilityCost_ = scalars.infeasibilityCost;
5459    numberRows_ = scalars.numberRows;
5460    numberColumns_ = scalars.numberColumns;
5461    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
5462    numberIterations_ = scalars.numberIterations;
5463    problemStatus_ = scalars.problemStatus;
5464    setMaximumIterations(scalars.maximumIterations);
5465    lengthNames_ = scalars.lengthNames;
5466    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
5467    numberDualInfeasibilitiesWithoutFree_
5468      = scalars.numberDualInfeasibilitiesWithoutFree;
5469    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
5470    numberRefinements_ = scalars.numberRefinements;
5471    scalingFlag_ = scalars.scalingFlag;
5472    algorithm_ = scalars.algorithm;
5473    specialOptions_ = scalars.specialOptions;
5474    // strings
5475    CoinBigIndex length;
5476#ifndef CLP_NO_STD
5477    for (i=0;i<ClpLastStrParam;i++) {
5478      numberRead = fread(&length,sizeof(int),1,fp);
5479      if (numberRead!=1)
5480        return 1;
5481      if (length) {
5482        char * array = new char[length+1];
5483        numberRead =