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

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

this may be a mistake

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