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

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

try and improve stability

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