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

Last change on this file since 1034 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

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