source: trunk/ClpSimplex.cpp @ 702

Last change on this file since 702 was 702, checked in by forrest, 15 years ago

out assert

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