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

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

out abort in crash

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