source: trunk/ClpSimplex.cpp @ 701

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

messages and preprocessing

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 265.6 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      assert (scalingFlag_>0||!rowScale_);
2585      delete [] rowScale_;
2586      delete [] columnScale_;
2587      rowScale_=NULL;
2588      columnScale_=NULL;
2589    }
2590    if (scalingFlag_>0&&!rowScale_) {
2591      if (matrix_->scale(this))
2592        scalingFlag_=-scalingFlag_; // not scaled after all
2593      if (rowScale_&&automaticScale_) {
2594        // try automatic scaling
2595        double smallestObj=1.0e100;
2596        double largestObj=0.0;
2597        double largestRhs=0.0;
2598        const double * obj = objective();
2599        for (i=0;i<numberColumns_;i++) {
2600          double value = fabs(obj[i]);
2601          value *= columnScale_[i];
2602          if (value&&columnLower_[i]!=columnUpper_[i]) {
2603            smallestObj = CoinMin(smallestObj,value);
2604            largestObj = CoinMax(largestObj,value);
2605          }
2606          if (columnLower_[i]>0.0||columnUpper_[i]<0.0) {
2607            double scale = 1.0/columnScale_[i];
2608            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2609            if (columnLower_[i]>0)
2610              largestRhs=CoinMax(largestRhs,columnLower_[i]*scale);
2611            if (columnUpper_[i]<0.0)
2612              largestRhs=CoinMax(largestRhs,-columnUpper_[i]*scale);
2613          }
2614        }
2615        for (i=0;i<numberRows_;i++) {
2616          if (rowLower_[i]>0.0||rowUpper_[i]<0.0) {
2617            double scale = rowScale_[i];
2618            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2619            if (rowLower_[i]>0)
2620              largestRhs=CoinMax(largestRhs,rowLower_[i]*scale);
2621            if (rowUpper_[i]<0.0)
2622              largestRhs=CoinMax(largestRhs,-rowUpper_[i]*scale);
2623          }
2624        }
2625        printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2626        bool scalingDone=false;
2627        // look at element range
2628        double smallestNegative;
2629        double largestNegative;
2630        double smallestPositive;
2631        double largestPositive;
2632        matrix_->rangeOfElements(smallestNegative, largestNegative,
2633                                 smallestPositive, largestPositive);
2634        smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2635        largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2636        if (largestObj) {
2637          double ratio = largestObj/smallestObj;
2638          double scale=1.0;
2639          if (ratio<1.0e8) {
2640            // reasonable
2641            if (smallestObj<1.0e-4) {
2642              // may as well scale up
2643              scalingDone=true;
2644              scale=1.0e-3/smallestObj;
2645            } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2646              //done=true;
2647            } else {
2648              scalingDone=true;
2649              if (algorithm_<0) {
2650                scale = 1.0e6/largestObj;
2651              } else {
2652                scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2653              }
2654            }
2655          } else if (ratio<1.0e12) {
2656            // not so good
2657            if (smallestObj<1.0e-7) {
2658              // may as well scale up
2659              scalingDone=true;
2660              scale=1.0e-6/smallestObj;
2661            } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2662              //done=true;
2663            } else {
2664              scalingDone=true;
2665              if (algorithm_<0) {
2666                scale = 1.0e7/largestObj;
2667              } else {
2668                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2669              }
2670            }
2671          } else {
2672            // Really nasty problem
2673            if (smallestObj<1.0e-8) {
2674              // may as well scale up
2675              scalingDone=true;
2676              scale=1.0e-7/smallestObj;
2677              largestObj *= scale;
2678            } 
2679            if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2680              //done=true;
2681            } else {
2682              scalingDone=true;
2683              if (algorithm_<0) {
2684                scale = 1.0e7/largestObj;
2685              } else {
2686                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2687              }
2688            }
2689          }
2690          objectiveScale_=scale;
2691        }
2692        if (largestRhs>1.0e12) {
2693          scalingDone=true;
2694          rhsScale_=1.0e9/largestRhs;
2695        } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
2696          scalingDone=true;
2697          rhsScale_=1.0e6/largestRhs;
2698        } else {
2699          rhsScale_=1.0;
2700        }
2701        if (scalingDone) {
2702          handler_->message(CLP_RIM_SCALE,messages_)
2703            <<objectiveScale_<<rhsScale_
2704            <<CoinMessageEol;
2705        }
2706      }
2707    } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
2708      matrix_->scaleRowCopy(this);
2709    }
2710  }
2711  if (what==63) {
2712    if (newArrays) {
2713      delete [] cost_;
2714      // extra copy with original costs
2715      //cost_ = new double[2*numberTotal];
2716      cost_ = new double[numberTotal];
2717      delete [] lower_;
2718      delete [] upper_;
2719      lower_ = new double[numberColumns_+numberRows2];
2720      upper_ = new double[numberColumns_+numberRows2];
2721      delete [] dj_;
2722      dj_ = new double[numberRows2+numberColumns_];
2723      delete [] solution_;
2724      solution_ = new double[numberRows2+numberColumns_];
2725    } else if (auxiliaryModel_) {
2726      assert (!cost_);
2727      cost_=auxiliaryModel_->cost_;
2728      assert (!lower_);
2729      lower_=auxiliaryModel_->lower_;
2730      assert (!upper_);
2731      upper_=auxiliaryModel_->upper_;
2732      assert (!dj_);
2733      dj_=auxiliaryModel_->dj_;
2734      assert (!solution_);
2735      solution_=auxiliaryModel_->solution_;
2736      assert (!rowScale_);
2737      assert (!columnScale_);
2738    }
2739    reducedCostWork_ = dj_;
2740    rowReducedCost_ = dj_+numberColumns_;
2741    columnActivityWork_ = solution_;
2742    rowActivityWork_ = solution_+numberColumns_;
2743    objectiveWork_ = cost_;
2744    rowObjectiveWork_ = cost_+numberColumns_;
2745    rowLowerWork_ = lower_+numberColumns_;
2746    columnLowerWork_ = lower_;
2747    rowUpperWork_ = upper_+numberColumns_;
2748    columnUpperWork_ = upper_;
2749  }
2750  if ((what&4)!=0) {
2751    if (!auxiliaryModel_||(what==63&&(auxiliaryModel_->specialOptions_&4)==0)) {
2752      double direction = optimizationDirection_*objectiveScale_;
2753      const double * obj = objective();
2754      const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
2755      const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
2756      // and also scale by scale factors
2757      if (rowScale) {
2758        if (rowObjective_) {
2759          for (i=0;i<numberRows_;i++)
2760            rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
2761        } else {
2762          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2763        }
2764        // If scaled then do all columns later in one loop
2765        if (what!=63||auxiliaryModel_) {
2766          for (i=0;i<numberColumns_;i++) {
2767            CoinAssert(fabs(obj[i])<1.0e25);
2768            objectiveWork_[i] = obj[i]*direction*columnScale[i];
2769          }
2770        }
2771      } else {
2772        if (rowObjective_) {
2773          for (i=0;i<numberRows_;i++)
2774            rowObjectiveWork_[i] = rowObjective_[i]*direction;
2775        } else {
2776          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2777        }
2778        for (i=0;i<numberColumns_;i++) {
2779          CoinAssert(fabs(obj[i])<1.0e25);
2780          objectiveWork_[i] = obj[i]*direction;
2781        }
2782      }
2783      if (auxiliaryModel_) {
2784        // save costs
2785        memcpy(auxiliaryModel_->cost_+numberTotal,cost_,numberTotal*sizeof(double));
2786      }
2787    } else {
2788      // just copy
2789      memcpy(cost_,auxiliaryModel_->cost_+numberTotal,numberTotal*sizeof(double));
2790    }
2791  }
2792  if ((what&1)!=0) {
2793    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
2794    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
2795    // clean up any mismatches on infinity
2796    // and fix any variables with tiny gaps
2797    double primalTolerance=dblParam_[ClpPrimalTolerance];
2798    if (!auxiliaryModel_) {
2799      if(rowScale) {
2800        // If scaled then do all columns later in one loop
2801        if (what!=63) {
2802          for (i=0;i<numberColumns_;i++) {
2803            double multiplier = rhsScale_/columnScale[i];
2804            double lowerValue = columnLower_[i];
2805            double upperValue = columnUpper_[i];
2806            if (lowerValue>-1.0e20) {
2807              columnLowerWork_[i]=lowerValue*multiplier;
2808              if (upperValue>=1.0e20) {
2809                columnUpperWork_[i]=COIN_DBL_MAX;
2810              } else {
2811                columnUpperWork_[i]=upperValue*multiplier;
2812                if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2813                  if (columnLowerWork_[i]>=0.0) {
2814                    columnUpperWork_[i] = columnLowerWork_[i];
2815                  } else if (columnUpperWork_[i]<=0.0) {
2816                    columnLowerWork_[i] = columnUpperWork_[i];
2817                  } else {
2818                    columnUpperWork_[i] = 0.0;
2819                    columnLowerWork_[i] = 0.0;
2820                  }
2821                }
2822              }
2823            } else if (upperValue<1.0e20) {
2824              columnLowerWork_[i]=-COIN_DBL_MAX;
2825              columnUpperWork_[i]=upperValue*multiplier;
2826            } else {
2827              // free
2828              columnLowerWork_[i]=-COIN_DBL_MAX;
2829              columnUpperWork_[i]=COIN_DBL_MAX;
2830            }
2831          }
2832        }
2833        for (i=0;i<numberRows_;i++) {
2834          double multiplier = rhsScale_*rowScale[i];
2835          double lowerValue = rowLower_[i];
2836          double upperValue = rowUpper_[i];
2837          if (lowerValue>-1.0e20) {
2838            rowLowerWork_[i]=lowerValue*multiplier;
2839            if (upperValue>=1.0e20) {
2840              rowUpperWork_[i]=COIN_DBL_MAX;
2841            } else {
2842              rowUpperWork_[i]=upperValue*multiplier;
2843              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2844                if (rowLowerWork_[i]>=0.0) {
2845                  rowUpperWork_[i] = rowLowerWork_[i];
2846                } else if (rowUpperWork_[i]<=0.0) {
2847                  rowLowerWork_[i] = rowUpperWork_[i];
2848                } else {
2849                  rowUpperWork_[i] = 0.0;
2850                  rowLowerWork_[i] = 0.0;
2851                }
2852              }
2853            }
2854          } else if (upperValue<1.0e20) {
2855            rowLowerWork_[i]=-COIN_DBL_MAX;
2856            rowUpperWork_[i]=upperValue*multiplier;
2857          } else {
2858            // free
2859            rowLowerWork_[i]=-COIN_DBL_MAX;
2860            rowUpperWork_[i]=COIN_DBL_MAX;
2861          }
2862        }
2863      } else if (rhsScale_!=1.0) {
2864        for (i=0;i<numberColumns_;i++) {
2865          double lowerValue = columnLower_[i];
2866          double upperValue = columnUpper_[i];
2867          if (lowerValue>-1.0e20) {
2868            columnLowerWork_[i]=lowerValue*rhsScale_;
2869            if (upperValue>=1.0e20) {
2870              columnUpperWork_[i]=COIN_DBL_MAX;
2871            } else {
2872              columnUpperWork_[i]=upperValue*rhsScale_;
2873              if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2874                if (columnLowerWork_[i]>=0.0) {
2875                  columnUpperWork_[i] = columnLowerWork_[i];
2876                } else if (columnUpperWork_[i]<=0.0) {
2877                  columnLowerWork_[i] = columnUpperWork_[i];
2878                } else {
2879                  columnUpperWork_[i] = 0.0;
2880                  columnLowerWork_[i] = 0.0;
2881                }
2882              }
2883            }
2884          } else if (upperValue<1.0e20) {
2885            columnLowerWork_[i]=-COIN_DBL_MAX;
2886            columnUpperWork_[i]=upperValue*rhsScale_;
2887          } else {
2888            // free
2889            columnLowerWork_[i]=-COIN_DBL_MAX;
2890            columnUpperWork_[i]=COIN_DBL_MAX;
2891          }
2892        }
2893        for (i=0;i<numberRows_;i++) {
2894          double lowerValue = rowLower_[i];
2895          double upperValue = rowUpper_[i];
2896          if (lowerValue>-1.0e20) {
2897            rowLowerWork_[i]=lowerValue*rhsScale_;
2898            if (upperValue>=1.0e20) {
2899              rowUpperWork_[i]=COIN_DBL_MAX;
2900            } else {
2901              rowUpperWork_[i]=upperValue*rhsScale_;
2902              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2903                if (rowLowerWork_[i]>=0.0) {
2904                  rowUpperWork_[i] = rowLowerWork_[i];
2905                } else if (rowUpperWork_[i]<=0.0) {
2906                  rowLowerWork_[i] = rowUpperWork_[i];
2907                } else {
2908                  rowUpperWork_[i] = 0.0;
2909                  rowLowerWork_[i] = 0.0;
2910                }
2911              }
2912            }
2913          } else if (upperValue<1.0e20) {
2914            rowLowerWork_[i]=-COIN_DBL_MAX;
2915            rowUpperWork_[i]=upperValue*rhsScale_;
2916          } else {
2917            // free
2918            rowLowerWork_[i]=-COIN_DBL_MAX;
2919            rowUpperWork_[i]=COIN_DBL_MAX;
2920          }
2921        }
2922      } else {
2923        for (i=0;i<numberColumns_;i++) {
2924          double lowerValue = columnLower_[i];
2925          double upperValue = columnUpper_[i];
2926          if (lowerValue>-1.0e20) {
2927            columnLowerWork_[i]=lowerValue;
2928            if (upperValue>=1.0e20) {
2929              columnUpperWork_[i]=COIN_DBL_MAX;
2930            } else {
2931              columnUpperWork_[i]=upperValue;
2932              if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2933                if (columnLowerWork_[i]>=0.0) {
2934                  columnUpperWork_[i] = columnLowerWork_[i];
2935                } else if (columnUpperWork_[i]<=0.0) {
2936                  columnLowerWork_[i] = columnUpperWork_[i];
2937                } else {
2938                  columnUpperWork_[i] = 0.0;
2939                  columnLowerWork_[i] = 0.0;
2940                }
2941              }
2942            }
2943          } else if (upperValue<1.0e20) {
2944            columnLowerWork_[i]=-COIN_DBL_MAX;
2945            columnUpperWork_[i]=upperValue;
2946          } else {
2947            // free
2948            columnLowerWork_[i]=-COIN_DBL_MAX;
2949            columnUpperWork_[i]=COIN_DBL_MAX;
2950          }
2951        }
2952        for (i=0;i<numberRows_;i++) {
2953          double lowerValue = rowLower_[i];
2954          double upperValue = rowUpper_[i];
2955          if (lowerValue>-1.0e20) {
2956            rowLowerWork_[i]=lowerValue;
2957            if (upperValue>=1.0e20) {
2958              rowUpperWork_[i]=COIN_DBL_MAX;
2959            } else {
2960              rowUpperWork_[i]=upperValue;
2961              if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2962                if (rowLowerWork_[i]>=0.0) {
2963                  rowUpperWork_[i] = rowLowerWork_[i];
2964                } else if (rowUpperWork_[i]<=0.0) {
2965                  rowLowerWork_[i] = rowUpperWork_[i];
2966                } else {
2967                  rowUpperWork_[i] = 0.0;
2968                  rowLowerWork_[i] = 0.0;
2969                }
2970              }
2971            }
2972          } else if (upperValue<1.0e20) {
2973            rowLowerWork_[i]=-COIN_DBL_MAX;
2974            rowUpperWork_[i]=upperValue;
2975          } else {
2976            // free
2977            rowLowerWork_[i]=-COIN_DBL_MAX;
2978            rowUpperWork_[i]=COIN_DBL_MAX;
2979          }
2980        }
2981      }
2982    } else {
2983      // auxiliary model
2984      if (what!=63) {
2985        // just copy
2986        memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberTotal*sizeof(double));
2987        memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberTotal*sizeof(double));
2988      } else {
2989        assert (rhsScale_==1.0);
2990        assert (objectiveScale_==1.0);
2991        if ((auxiliaryModel_->specialOptions_&2)==0) {
2992          // need to do column bounds
2993          if (columnScale) {
2994            const double * inverseScale = columnScale+numberColumns_;
2995            // scaled
2996            for (i=0;i<numberColumns_;i++) {
2997              double value;
2998              value = columnLower_[i];
2999              if (value>-1.0e20) 
3000                value *= inverseScale[i];
3001              lower_[i] = value;
3002              value = columnUpper_[i];
3003              if (value<1.0e20) 
3004                value *= inverseScale[i];
3005              upper_[i] = value;
3006            }
3007          } else {
3008            for (i=0;i<numberColumns_;i++) {
3009              lower_[i] = columnLower_[i];
3010              upper_[i] = columnUpper_[i];
3011            }
3012          }
3013          // save
3014          memcpy(auxiliaryModel_->lower_+numberTotal,lower_,numberColumns_*sizeof(double));
3015          memcpy(auxiliaryModel_->upper_+numberTotal,upper_,numberColumns_*sizeof(double));
3016        } else {
3017          // just copy
3018          memcpy(lower_,auxiliaryModel_->lower_+numberTotal,numberColumns_*sizeof(double));
3019          memcpy(upper_,auxiliaryModel_->upper_+numberTotal,numberColumns_*sizeof(double));
3020        }
3021        if ((auxiliaryModel_->specialOptions_&1)==0) {
3022          // need to do row bounds
3023          if (rowScale) {
3024            // scaled
3025            for (i=0;i<numberRows_;i++) {
3026              double value;
3027              value = rowLower_[i];
3028              if (value>-1.0e20) 
3029                value *= rowScale[i];
3030              lower_[i+numberColumns_] = value;
3031              value = rowUpper_[i];
3032              if (value<1.0e20) 
3033                value *= rowScale[i];
3034              upper_[i+numberColumns_] = value;
3035            }
3036          } else {
3037            for (i=0;i<numberRows_;i++) {
3038              lower_[i+numberColumns_] = rowLower_[i];
3039              upper_[i+numberColumns_] = rowUpper_[i];
3040            }
3041          }
3042          // save
3043          memcpy(auxiliaryModel_->lower_+numberTotal+numberColumns_,
3044                 lower_+numberColumns_,numberRows_*sizeof(double));
3045          memcpy(auxiliaryModel_->upper_+numberTotal+numberColumns_,
3046                 upper_+numberColumns_,numberRows_*sizeof(double));
3047        } else {
3048          // just copy
3049          memcpy(lower_+numberColumns_,
3050                 auxiliaryModel_->lower_+numberTotal+numberColumns_,numberRows_*sizeof(double));
3051          memcpy(upper_+numberColumns_,
3052                 auxiliaryModel_->upper_+numberTotal+numberColumns_,numberRows_*sizeof(double));
3053        }
3054      }
3055      if (what==63) {
3056        // do basis
3057        assert ((auxiliaryModel_->specialOptions_&8)!=0);
3058        // clean solution trusting basis
3059        for ( i=0;i<numberTotal;i++) {
3060          Status status =getStatus(i);
3061          double value=solution_[i]; // may as well leave if basic (values pass)
3062          if (status!=basic) {
3063            if (upper_[i]==lower_[i]) {
3064              setStatus(i,isFixed);
3065              value=lower_[i];
3066            } else if (status==atLowerBound) {
3067              if (lower_[i]>-1.0e20) {
3068                value=lower_[i];
3069              } else {
3070                if (upper_[i]<1.0e20) {
3071                  value=upper_[i];
3072                  setStatus(i,atUpperBound);
3073                } else {
3074                  setStatus(i,isFree);
3075                }
3076              }
3077            } else if (status==atUpperBound) {
3078              if (upper_[i]<1.0e20) {
3079              value=upper_[i];
3080              } else {
3081                if (lower_[i]>-1.0e20) {
3082                  value=lower_[i];
3083                  setStatus(i,atLowerBound);
3084                } else {
3085                  setStatus(i,isFree);
3086                }
3087              }
3088            }
3089          }
3090          solution_[i]=value;
3091        }
3092      }
3093    }
3094  }
3095  if (what==63) {
3096    // move information to work arrays
3097    double direction = optimizationDirection_;
3098    // direction is actually scale out not scale in
3099    if (direction)
3100      direction = 1.0/direction;
3101    if (direction!=1.0&&(!auxiliaryModel_||(auxiliaryModel_->specialOptions_&8)==0)) {
3102      // reverse all dual signs
3103      for (i=0;i<numberColumns_;i++) 
3104        reducedCost_[i] *= direction;
3105      for (i=0;i<numberRows_;i++) 
3106        dual_[i] *= direction;
3107    }
3108    for (i=0;i<numberRows_+numberColumns_;i++) {
3109      setFakeBound(i,noFake);
3110    }
3111    if (!auxiliaryModel_) {
3112      if (rowScale_) {
3113        const double * obj = objective();
3114        double direction = optimizationDirection_*objectiveScale_;
3115        // clean up any mismatches on infinity
3116        // and fix any variables with tiny gaps
3117        double primalTolerance=dblParam_[ClpPrimalTolerance];
3118        // on entry
3119        for (i=0;i<numberColumns_;i++) {
3120          CoinAssert(fabs(obj[i])<1.0e25);
3121          double scaleFactor = columnScale_[i];
3122          double multiplier = rhsScale_/scaleFactor;
3123          scaleFactor *= direction;
3124          objectiveWork_[i] = obj[i]*scaleFactor;
3125          reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
3126          double lowerValue = columnLower_[i];
3127          double upperValue = columnUpper_[i];
3128          if (lowerValue>-1.0e20) {
3129            columnLowerWork_[i]=lowerValue*multiplier;
3130            if (upperValue>=1.0e20) {
3131              columnUpperWork_[i]=COIN_DBL_MAX;
3132            } else {
3133              columnUpperWork_[i]=upperValue*multiplier;
3134              if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
3135                if (columnLowerWork_[i]>=0.0) {
3136                  columnUpperWork_[i] = columnLowerWork_[i];
3137                } else if (columnUpperWork_[i]<=0.0) {
3138                  columnLowerWork_[i] = columnUpperWork_[i];
3139                } else {
3140                  columnUpperWork_[i] = 0.0;
3141                  columnLowerWork_[i] = 0.0;
3142                }
3143              }
3144            }
3145          } else if (upperValue<1.0e20) {
3146            columnLowerWork_[i]=-COIN_DBL_MAX;
3147            columnUpperWork_[i]=upperValue*multiplier;
3148          } else {
3149            // free
3150            columnLowerWork_[i]=-COIN_DBL_MAX;
3151            columnUpperWork_[i]=COIN_DBL_MAX;
3152          }
3153          double value = columnActivity_[i] * multiplier;
3154          if (fabs(value)>1.0e20) {
3155            //printf("bad value of %g for column %d\n",value,i);
3156            setColumnStatus(i,superBasic);
3157            if (columnUpperWork_[i]<0.0) {
3158              value=columnUpperWork_[i];
3159            } else if (columnLowerWork_[i]>0.0) {
3160              value=columnLowerWork_[i];
3161            } else {
3162              value=0.0;
3163            }
3164          }
3165          columnActivityWork_[i] = value;
3166        }
3167        for (i=0;i<numberRows_;i++) {
3168          dual_[i] /= rowScale_[i];
3169          dual_[i] *= objectiveScale_;
3170          rowReducedCost_[i] = dual_[i];
3171          double multiplier = rhsScale_*rowScale_[i];
3172          double value = rowActivity_[i] * multiplier;
3173          if (fabs(value)>1.0e20) {
3174            //printf("bad value of %g for row %d\n",value,i);
3175            setRowStatus(i,superBasic);
3176            if (rowUpperWork_[i]<0.0) {
3177              value=rowUpperWork_[i];
3178            } else if (rowLowerWork_[i]>0.0) {
3179              value=rowLowerWork_[i];
3180            } else {
3181              value=0.0;
3182            }
3183          }
3184          rowActivityWork_[i] = value;
3185        }
3186      } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
3187        // on entry
3188        for (i=0;i<numberColumns_;i++) {
3189          double value = columnActivity_[i];
3190          value *= rhsScale_;
3191          if (fabs(value)>1.0e20) {
3192            //printf("bad value of %g for column %d\n",value,i);
3193            setColumnStatus(i,superBasic);
3194            if (columnUpperWork_[i]<0.0) {
3195              value=columnUpperWork_[i];
3196            } else if (columnLowerWork_[i]>0.0) {
3197              value=columnLowerWork_[i];
3198            } else {
3199              value=0.0;
3200            }
3201          }
3202          columnActivityWork_[i] = value;
3203          reducedCostWork_[i] = reducedCost_[i]*objectiveScale_;
3204        }
3205        for (i=0;i<numberRows_;i++) {
3206          double value = rowActivity_[i];
3207          value *= rhsScale_;
3208          if (fabs(value)>1.0e20) {
3209            //printf("bad value of %g for row %d\n",value,i);
3210            setRowStatus(i,superBasic);
3211            if (rowUpperWork_[i]<0.0) {
3212              value=rowUpperWork_[i];
3213            } else if (rowLowerWork_[i]>0.0) {
3214              value=rowLowerWork_[i];
3215            } else {
3216              value=0.0;
3217            }
3218          }
3219          rowActivityWork_[i] = value;
3220          dual_[i] *= objectiveScale_;
3221          rowReducedCost_[i] = dual_[i];
3222        }
3223      } else {
3224        // on entry
3225        for (i=0;i<numberColumns_;i++) {
3226          double value = columnActivity_[i];
3227          if (fabs(value)>1.0e20) {
3228            //printf("bad value of %g for column %d\n",value,i);
3229            setColumnStatus(i,superBasic);
3230            if (columnUpperWork_[i]<0.0) {
3231              value=columnUpperWork_[i];
3232            } else if (columnLowerWork_[i]>0.0) {
3233              value=columnLowerWork_[i];
3234            } else {
3235              value=0.0;
3236            }
3237          }
3238          columnActivityWork_[i] = value;
3239          reducedCostWork_[i] = reducedCost_[i];
3240        }
3241        for (i=0;i<numberRows_;i++) {
3242          double value = rowActivity_[i];
3243          if (fabs(value)>1.0e20) {
3244            //printf("bad value of %g for row %d\n",value,i);
3245            setRowStatus(i,superBasic);
3246            if (rowUpperWork_[i]<0.0) {
3247              value=rowUpperWork_[i];
3248            } else if (rowLowerWork_[i]>0.0) {
3249              value=rowLowerWork_[i];
3250            } else {
3251              value=0.0;
3252            }
3253          }
3254          rowActivityWork_[i] = value;
3255          rowReducedCost_[i] = dual_[i];
3256        }
3257      }
3258    }
3259  }
3260 
3261  if (what==63&&doSanityCheck&&!auxiliaryModel_) {
3262    // check rim of problem okay
3263    if (!sanityCheck())
3264      goodMatrix=false;
3265  } 
3266  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
3267  // maybe we need to move scales to SimplexModel for factorization?
3268  if ((what==63&&!pivotVariable_)||(newArrays&&!keepPivots)) {
3269    delete [] pivotVariable_;
3270    pivotVariable_=new int[numberRows2];
3271    for (int i=0;i<numberRows2;i++)
3272      pivotVariable_[i]=-1;
3273  } else if (what==63&&!keepPivots) {
3274    // just reset
3275    for (int i=0;i<numberRows2;i++)
3276      pivotVariable_[i]=-1;
3277  } else if (what==63) {
3278    // check pivots
3279    for (int i=0;i<numberRows2;i++) {
3280      int iSequence = pivotVariable_[i];
3281      if (iSequence<numberRows_+numberColumns_&&
3282          getStatus(iSequence)!=basic) {
3283        keepPivots =false;
3284        break;
3285      }
3286    }
3287    if (!keepPivots) {
3288      // reset
3289      for (int i=0;i<numberRows2;i++)
3290        pivotVariable_[i]=-1;
3291    } else {
3292      // clean
3293      for (int i=0;i<numberColumns_+numberRows_;i++) {
3294        Status status =getStatus(i);
3295        if (status!=basic) {
3296          if (upper_[i]==lower_[i]) {
3297            setStatus(i,isFixed);
3298            solution_[i]=lower_[i];
3299          } else if (status==atLowerBound) {
3300            if (lower_[i]>-1.0e20) {
3301              solution_[i]=lower_[i];
3302            } else {
3303              //printf("seq %d at lower of %g\n",i,lower_[i]);
3304              if (upper_[i]<1.0e20) {
3305                solution_[i]=upper_[i];
3306                setStatus(i,atUpperBound);
3307              } else {
3308                setStatus(i,isFree);
3309              }
3310            }
3311          } else if (status==atUpperBound) {
3312            if (upper_[i]<1.0e20) {
3313              solution_[i]=upper_[i];
3314            } else {
3315              //printf("seq %d at upper of %g\n",i,upper_[i]);
3316              if (lower_[i]>-1.0e20) {
3317                solution_[i]=lower_[i];
3318                setStatus(i,atLowerBound);
3319              } else {
3320                setStatus(i,isFree);
3321              }
3322            }
3323          }
3324        }
3325      }
3326    }
3327  }
3328 
3329  if (what==63) {
3330    if (newArrays) {
3331      // get some arrays
3332      int iRow,iColumn;
3333      // these are "indexed" arrays so we always know where nonzeros are
3334      /**********************************************************
3335      rowArray_[3] is long enough for rows+columns
3336      *********************************************************/
3337      for (iRow=0;iRow<4;iRow++) {
3338        delete rowArray_[iRow];
3339        rowArray_[iRow]=new CoinIndexedVector();
3340        int length =numberRows2+factorization_->maximumPivots();
3341        if (iRow==3||objective_->type()>1)
3342          length += numberColumns_;
3343        rowArray_[iRow]->reserve(length);
3344      }
3345     
3346      for (iColumn=0;iColumn<2;iColumn++) {
3347        delete columnArray_[iColumn];
3348        columnArray_[iColumn]=new CoinIndexedVector();
3349        if (!iColumn)
3350          columnArray_[iColumn]->reserve(numberColumns_);
3351        else
3352          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
3353      }
3354    } else {
3355      int iRow,iColumn;
3356      if (auxiliaryModel_) {
3357        for (iRow=0;iRow<4;iRow++) {
3358          assert (!rowArray_[iRow]);
3359          rowArray_[iRow]=auxiliaryModel_->rowArray_[iRow];
3360          // For safety
3361          memset(rowArray_[iRow]->denseVector(),0,rowArray_[iRow]->capacity()*sizeof(double));
3362        }
3363        for (iColumn=0;iColumn<2;iColumn++) {
3364          assert (!columnArray_[iColumn]);
3365          columnArray_[iColumn]=auxiliaryModel_->columnArray_[iColumn];
3366          // For safety
3367          memset(columnArray_[iColumn]->denseVector(),0,columnArray_[iColumn]->capacity()*sizeof(double));
3368        }
3369        // do matrices
3370        rowCopy_=auxiliaryModel_->rowCopy_;
3371        ClpMatrixBase * temp = matrix_;
3372        matrix_=auxiliaryModel_->matrix_;
3373        auxiliaryModel_->matrix_=temp;
3374      }
3375      for (iRow=0;iRow<4;iRow++) {
3376        rowArray_[iRow]->clear();
3377#ifndef NDEBUG
3378        int length =numberRows2+factorization_->maximumPivots();
3379        if (iRow==3||objective_->type()>1)
3380          length += numberColumns_;
3381        assert(rowArray_[iRow]->capacity()==length);
3382        rowArray_[iRow]->checkClear();
3383#endif
3384      }
3385     
3386      for (iColumn=0;iColumn<2;iColumn++) {
3387        columnArray_[iColumn]->clear();
3388#ifndef NDEBUG
3389        int length =numberColumns_;
3390        if (iColumn)
3391          length=CoinMax(numberRows2,numberColumns_);
3392        assert(columnArray_[iColumn]->capacity()==length);
3393        columnArray_[iColumn]->checkClear();
3394#endif
3395      }
3396    }   
3397  }
3398  if (problemStatus_==10) {
3399    problemStatus_=-1;
3400    handler_->setLogLevel(saveLevel); // switch back messages
3401  }
3402  if ((what&5)!=0) 
3403    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
3404  return goodMatrix;
3405}
3406void
3407ClpSimplex::deleteRim(int getRidOfFactorizationData)
3408{
3409  if (!auxiliaryModel_) {
3410    int i;
3411    if (problemStatus_!=1&&problemStatus_!=2) {
3412      delete [] ray_;
3413      ray_=NULL;
3414    }
3415    // set upperOut_ to furthest away from bound so can use in dual for dualBound_
3416    upperOut_=1.0;
3417#if 0
3418    {
3419      int nBad=0;
3420      for (i=0;i<numberColumns_;i++) {
3421        if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
3422          nBad++;
3423      }
3424      if (nBad)
3425        printf("yy %d basic fixed\n",nBad);
3426    }
3427#endif
3428    // ray may be null if in branch and bound
3429    if (rowScale_) {
3430      // Collect infeasibilities
3431      int numberPrimalScaled=0;
3432      int numberPrimalUnscaled=0;
3433      int numberDualScaled=0;
3434      int numberDualUnscaled=0;
3435      double scaleC = 1.0/objectiveScale_;
3436      double scaleR = 1.0/rhsScale_;
3437      for (i=0;i<numberColumns_;i++) {
3438        double scaleFactor = columnScale_[i];
3439        double valueScaled = columnActivityWork_[i];
3440        double lowerScaled = columnLowerWork_[i];
3441        double upperScaled = columnUpperWork_[i];
3442        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3443          if (valueScaled<lowerScaled-primalTolerance_||
3444              valueScaled>upperScaled+primalTolerance_)
3445            numberPrimalScaled++;
3446          else
3447            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3448        }
3449        columnActivity_[i] = valueScaled*scaleFactor*scaleR;
3450        double value = columnActivity_[i];
3451        if (value<columnLower_[i]-primalTolerance_)
3452          numberPrimalUnscaled++;
3453        else if (value>columnUpper_[i]+primalTolerance_)
3454          numberPrimalUnscaled++;
3455        double valueScaledDual = reducedCostWork_[i];
3456        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3457          numberDualScaled++;
3458        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3459          numberDualScaled++;
3460        reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
3461        double valueDual = reducedCost_[i];
3462        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3463          numberDualUnscaled++;
3464        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3465          numberDualUnscaled++;
3466      }
3467      for (i=0;i<numberRows_;i++) {
3468        double scaleFactor = rowScale_[i];
3469        double valueScaled = rowActivityWork_[i];
3470        double lowerScaled = rowLowerWork_[i];
3471        double upperScaled = rowUpperWork_[i];
3472        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3473          if (valueScaled<lowerScaled-primalTolerance_||
3474              valueScaled>upperScaled+primalTolerance_)
3475            numberPrimalScaled++;
3476          else
3477            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3478        }
3479        rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
3480        double value = rowActivity_[i];
3481        if (value<rowLower_[i]-primalTolerance_)
3482          numberPrimalUnscaled++;
3483        else if (value>rowUpper_[i]+primalTolerance_)
3484          numberPrimalUnscaled++;
3485        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
3486        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3487          numberDualScaled++;
3488        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3489          numberDualScaled++;
3490        dual_[i] *= scaleFactor*scaleC;
3491        double valueDual = dual_[i]; 
3492        if (rowObjective_)
3493          valueDual += rowObjective_[i];
3494        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3495          numberDualUnscaled++;
3496        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3497          numberDualUnscaled++;
3498      }
3499      if (!problemStatus_&&!secondaryStatus_) {
3500        // See if we need to set secondary status
3501        if (numberPrimalUnscaled) {
3502          if (numberDualUnscaled) 
3503            secondaryStatus_=4;
3504          else
3505            secondaryStatus_=2;
3506        } else {
3507          if (numberDualUnscaled) 
3508            secondaryStatus_=3;
3509        }
3510      }
3511      if (problemStatus_==2) {
3512        for (i=0;i<numberColumns_;i++) {
3513          ray_[i] *= columnScale_[i];
3514        }
3515      } else if (problemStatus_==1&&ray_) {
3516        for (i=0;i<numberRows_;i++) {
3517          ray_[i] *= rowScale_[i];
3518        }
3519      }
3520    } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
3521      // Collect infeasibilities
3522      int numberPrimalScaled=0;
3523      int numberPrimalUnscaled=0;
3524      int numberDualScaled=0;
3525      int numberDualUnscaled=0;
3526      double scaleC = 1.0/objectiveScale_;
3527      double scaleR = 1.0/rhsScale_;
3528      for (i=0;i<numberColumns_;i++) {
3529        double valueScaled = columnActivityWork_[i];
3530        double lowerScaled = columnLowerWork_[i];
3531        double upperScaled = columnUpperWork_[i];
3532        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3533          if (valueScaled<lowerScaled-primalTolerance_||
3534              valueScaled>upperScaled+primalTolerance_)
3535            numberPrimalScaled++;
3536          else
3537            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3538        }
3539        columnActivity_[i] = valueScaled*scaleR;
3540        double value = columnActivity_[i];
3541        if (value<columnLower_[i]-primalTolerance_)
3542          numberPrimalUnscaled++;
3543        else if (value>columnUpper_[i]+primalTolerance_)
3544          numberPrimalUnscaled++;
3545        double valueScaledDual = reducedCostWork_[i];
3546        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3547          numberDualScaled++;
3548        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3549          numberDualScaled++;
3550        reducedCost_[i] = valueScaledDual*scaleC;
3551        double valueDual = reducedCost_[i];
3552        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3553          numberDualUnscaled++;
3554        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3555          numberDualUnscaled++;
3556      }
3557      for (i=0;i<numberRows_;i++) {
3558        double valueScaled = rowActivityWork_[i];
3559        double lowerScaled = rowLowerWork_[i];
3560        double upperScaled = rowUpperWork_[i];
3561        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
3562          if (valueScaled<lowerScaled-primalTolerance_||
3563              valueScaled>upperScaled+primalTolerance_)
3564            numberPrimalScaled++;
3565          else
3566            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
3567        }
3568        rowActivity_[i] = valueScaled*scaleR;
3569        double value = rowActivity_[i];
3570        if (value<rowLower_[i]-primalTolerance_)
3571          numberPrimalUnscaled++;
3572        else if (value>rowUpper_[i]+primalTolerance_)
3573          numberPrimalUnscaled++;
3574        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
3575        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
3576          numberDualScaled++;
3577        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
3578          numberDualScaled++;
3579        dual_[i] *= scaleC;
3580        double valueDual = dual_[i]; 
3581        if (rowObjective_)
3582          valueDual += rowObjective_[i];
3583        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
3584          numberDualUnscaled++;
3585        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
3586          numberDualUnscaled++;
3587      }
3588      if (!problemStatus_&&!secondaryStatus_) {
3589        // See if we need to set secondary status
3590        if (numberPrimalUnscaled) {
3591          if (numberDualUnscaled) 
3592            secondaryStatus_=4;
3593          else
3594            secondaryStatus_=2;
3595        } else {
3596          if (numberDualUnscaled) 
3597            secondaryStatus_=3;
3598        }
3599      }
3600    } else {
3601      if (columnActivityWork_) {
3602        for (i=0;i<numberColumns_;i++) {
3603          double value = columnActivityWork_[i];
3604          double lower = columnLowerWork_[i];
3605          double upper = columnUpperWork_[i];
3606          if (lower>-1.0e20||upper<1.0e20) {
3607            if (value>lower&&value<upper)
3608              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3609          }
3610          columnActivity_[i] = columnActivityWork_[i];
3611          reducedCost_[i] = reducedCostWork_[i];
3612        }
3613        for (i=0;i<numberRows_;i++) {
3614          double value = rowActivityWork_[i];
3615          double lower = rowLowerWork_[i];
3616          double upper = rowUpperWork_[i];
3617          if (lower>-1.0e20||upper<1.0e20) {
3618            if (value>lower&&value<upper)
3619              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3620          }
3621          rowActivity_[i] = rowActivityWork_[i];
3622        }
3623      }
3624    }
3625    // switch off scalefactor if auto
3626    if (automaticScale_) {
3627      rhsScale_=1.0;
3628      objectiveScale_=1.0;
3629    }
3630    if (optimizationDirection_!=1.0) {
3631      // and modify all dual signs
3632      for (i=0;i<numberColumns_;i++) 
3633        reducedCost_[i] *= optimizationDirection_;
3634      for (i=0;i<numberRows_;i++) 
3635        dual_[i] *= optimizationDirection_;
3636    }
3637  } else {
3638    // auxiliary model
3639    cost_=NULL;
3640    lower_=NULL;
3641    upper_=NULL;
3642    dj_=NULL;
3643    solution_=NULL;
3644    int iRow,iColumn;
3645    for (iRow=0;iRow<4;iRow++) {
3646      if (rowArray_[iRow]) {
3647        rowArray_[iRow]->clear();
3648        rowArray_[iRow]->checkClear();
3649      }
3650      rowArray_[iRow]=NULL;
3651    }
3652    for (iColumn=0;iColumn<2;iColumn++) {
3653      if (columnArray_[iColumn]) {
3654        columnArray_[iColumn]->clear();
3655        columnArray_[iColumn]->checkClear();
3656      }
3657      columnArray_[iColumn]=NULL;
3658    }
3659    rowCopy_=NULL;
3660    ClpMatrixBase * temp = matrix_;
3661    matrix_=auxiliaryModel_->matrix_;
3662    auxiliaryModel_->matrix_=temp;
3663    assert ((auxiliaryModel_->specialOptions_&16+32)==16+32);
3664    if (problemStatus_!=1&&problemStatus_!=10) {
3665      int i;
3666      if (auxiliaryModel_->rowScale_) {
3667        const double * scale = auxiliaryModel_->columnScale_;
3668        const double * inverseScale = scale + numberColumns_;
3669        for (i=0;i<numberColumns_;i++) {
3670          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i]*scale[i];
3671          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i]*inverseScale[i];
3672        }
3673        scale = auxiliaryModel_->rowScale_;
3674        inverseScale = scale + numberRows_;
3675        for (i=0;i<numberRows_;i++) {
3676          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
3677        }
3678      } else {
3679        for (i=0;i<numberColumns_;i++) {
3680          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i];
3681          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i];
3682        }
3683        for (i=0;i<numberRows_;i++) {
3684          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i];
3685        }
3686      }
3687      if (optimizationDirection_!=1.0) {
3688        // and modify reduced costs
3689        for (i=0;i<numberColumns_;i++) 
3690          reducedCost_[i] *= optimizationDirection_;
3691      }
3692    } else if (problemStatus_==10) {
3693      int i;
3694      if (auxiliaryModel_->rowScale_) {
3695        const double * scale = auxiliaryModel_->columnScale_;
3696        const double * inverseScale = scale + numberColumns_;
3697        for (i=0;i<numberColumns_;i++) {
3698          double lower = auxiliaryModel_->columnLowerWork_[i];
3699          double upper = auxiliaryModel_->columnUpperWork_[i];
3700          double value = auxiliaryModel_->columnActivityWork_[i];
3701          if (value>lower&&value<upper) {
3702            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3703          }
3704          columnActivity_[i] = value*scale[i];
3705        }
3706        scale = auxiliaryModel_->rowScale_;
3707        inverseScale = scale + numberRows_;
3708        for (i=0;i<numberRows_;i++) {
3709          double lower = auxiliaryModel_->rowLowerWork_[i];
3710          double upper = auxiliaryModel_->rowUpperWork_[i];
3711          double value = auxiliaryModel_->rowActivityWork_[i];
3712          if (value>lower&&value<upper) {
3713            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3714          }
3715          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
3716        }
3717      } else {
3718        for (i=0;i<numberColumns_;i++) {
3719          double lower = auxiliaryModel_->columnLowerWork_[i];
3720          double upper = auxiliaryModel_->columnUpperWork_[i];
3721          double value = auxiliaryModel_->columnActivityWork_[i];
3722          if (value>lower&&value<upper) {
3723            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3724          }
3725          columnActivity_[i] = value;
3726        }
3727        for (i=0;i<numberRows_;i++) {
3728          double lower = auxiliaryModel_->rowLowerWork_[i];
3729          double upper = auxiliaryModel_->rowUpperWork_[i];
3730          double value = auxiliaryModel_->rowActivityWork_[i];
3731          if (value>lower&&value<upper) {
3732            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
3733          }
3734          rowActivity_[i] = value;
3735        }
3736      }
3737    }
3738  }
3739  // scaling may have been turned off
3740  scalingFlag_ = abs(scalingFlag_);
3741  if(getRidOfFactorizationData>0) {
3742    gutsOfDelete(getRidOfFactorizationData+1);
3743  } else {
3744    // at least get rid of nonLinearCost_
3745    delete nonLinearCost_;
3746    nonLinearCost_=NULL;
3747  }
3748  // get rid of data
3749  matrix_->generalExpanded(this,13,scalingFlag_);
3750}
3751void 
3752ClpSimplex::setDualBound(double value)
3753{
3754  if (value>0.0)
3755    dualBound_=value;
3756}
3757void 
3758ClpSimplex::setInfeasibilityCost(double value)
3759{
3760  if (value>0.0)
3761    infeasibilityCost_=value;
3762}
3763void ClpSimplex::setNumberRefinements( int value) 
3764{
3765  if (value>=0&&value<10)
3766    numberRefinements_=value;
3767}
3768// Sets row pivot choice algorithm in dual
3769void 
3770ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
3771{
3772  delete dualRowPivot_;
3773  dualRowPivot_ = choice.clone(true);
3774}
3775// Sets row pivot choice algorithm in dual
3776void 
3777ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
3778{
3779  delete primalColumnPivot_;
3780  primalColumnPivot_ = choice.clone(true);
3781}
3782// Passes in factorization
3783void 
3784ClpSimplex::setFactorization( ClpFactorization & factorization)
3785{
3786  delete factorization_;
3787  factorization_= new ClpFactorization(factorization);
3788}
3789/* Perturbation:
3790   -50 to +50 - perturb by this power of ten (-6 sounds good)
3791   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
3792   101 - we are perturbed
3793   102 - don't try perturbing again
3794   default is 100
3795*/
3796void 
3797ClpSimplex::setPerturbation(int value)
3798{
3799  if(value<=100&&value >=-1000) {
3800    perturbation_=value;
3801  } 
3802}
3803// Sparsity on or off
3804bool 
3805ClpSimplex::sparseFactorization() const
3806{
3807  return factorization_->sparseThreshold()!=0;
3808}
3809void 
3810ClpSimplex::setSparseFactorization(bool value)
3811{
3812  if (value) {
3813    if (!factorization_->sparseThreshold())
3814      factorization_->goSparse();
3815  } else {
3816    factorization_->sparseThreshold(0);
3817  }
3818}
3819void checkCorrect(ClpSimplex * model,int iRow,
3820                  const double * element,const int * rowStart,const int * rowLength,
3821                  const int * column,
3822                  const double * columnLower_, const double * columnUpper_,
3823                  int infiniteUpperC,
3824                  int infiniteLowerC,
3825                  double &maximumUpC,
3826                  double &maximumDownC)
3827{
3828  int infiniteUpper = 0;
3829  int infiniteLower = 0;
3830  double maximumUp = 0.0;
3831  double maximumDown = 0.0;
3832  CoinBigIndex rStart = rowStart[iRow];
3833  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
3834  CoinBigIndex j;
3835  double large=1.0e15;
3836  int iColumn;
3837  // Compute possible lower and upper ranges
3838 
3839  for (j = rStart; j < rEnd; ++j) {
3840    double value=element[j];
3841    iColumn = column[j];
3842    if (value > 0.0) {
3843      if (columnUpper_[iColumn] >= large) {
3844        ++infiniteUpper;
3845      } else {
3846        maximumUp += columnUpper_[iColumn] * value;
3847      }
3848      if (columnLower_[iColumn] <= -large) {
3849        ++infiniteLower;
3850      } else {
3851        maximumDown += columnLower_[iColumn] * value;
3852      }
3853    } else if (value<0.0) {
3854      if (columnUpper_[iColumn] >= large) {
3855        ++infiniteLower;
3856      } else {
3857        maximumDown += columnUpper_[iColumn] * value;
3858      }
3859      if (columnLower_[iColumn] <= -large) {
3860        ++infiniteUpper;
3861      } else {
3862        maximumUp += columnLower_[iColumn] * value;
3863      }
3864    }
3865  }
3866  assert (infiniteLowerC==infiniteLower);
3867  assert (infiniteUpperC==infiniteUpper);
3868  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
3869    printf("row %d comp up %g, true up %g\n",iRow,
3870           maximumUpC,maximumUp);
3871  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
3872    printf("row %d comp down %g, true down %g\n",iRow,
3873           maximumDownC,maximumDown);
3874  maximumUpC=maximumUp;
3875  maximumDownC=maximumDown;
3876}
3877
3878/* Tightens primal bounds to make dual faster.  Unless
3879   fixed, bounds are slightly looser than they could be.
3880   This is to make dual go faster and is probably not needed
3881   with a presolve.  Returns non-zero if problem infeasible
3882
3883   Fudge for branch and bound - put bounds on columns of factor *
3884   largest value (at continuous) - should improve stability
3885   in branch and bound on infeasible branches (0.0 is off)
3886*/
3887int 
3888ClpSimplex::tightenPrimalBounds(double factor,int doTight)
3889{
3890 
3891  // Get a row copy in standard format
3892  CoinPackedMatrix copy;
3893  copy.reverseOrderedCopyOf(*matrix());
3894  // get matrix data pointers
3895  const int * column = copy.getIndices();
3896  const CoinBigIndex * rowStart = copy.getVectorStarts();
3897  const int * rowLength = copy.getVectorLengths(); 
3898  const double * element = copy.getElements();
3899  int numberChanged=1,iPass=0;
3900  double large = largeValue(); // treat bounds > this as infinite
3901#ifndef NDEBUG
3902  double large2= 1.0e10*large;
3903#endif
3904  int numberInfeasible=0;
3905  int totalTightened = 0;
3906
3907  double tolerance = primalTolerance();
3908
3909
3910  // Save column bounds
3911  double * saveLower = new double [numberColumns_];
3912  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
3913  double * saveUpper = new double [numberColumns_];
3914  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
3915
3916  int iRow, iColumn;
3917
3918  // If wanted - tighten column bounds using solution
3919  if (factor) {
3920    double largest=0.0;
3921    if (factor>0.0) {
3922      assert (factor>1.0);
3923      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3924        if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3925          largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
3926        }
3927      }
3928      largest *= factor;
3929    } else {
3930      // absolute
3931       largest = - factor;
3932    }
3933    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3934      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3935        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
3936        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
3937      }
3938    }
3939  }
3940#define MAXPASS 10
3941
3942  // Loop round seeing if we can tighten bounds
3943  // Would be faster to have a stack of possible rows
3944  // and we put altered rows back on stack
3945  int numberCheck=-1;
3946  while(numberChanged>numberCheck) {
3947
3948    numberChanged = 0; // Bounds tightened this pass
3949   
3950    if (iPass==MAXPASS) break;
3951    iPass++;
3952   
3953    for (iRow = 0; iRow < numberRows_; iRow++) {
3954
3955      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
3956
3957        // possible row
3958        int infiniteUpper = 0;
3959        int infiniteLower = 0;
3960        double maximumUp = 0.0;
3961        double maximumDown = 0.0;
3962        double newBound;
3963        CoinBigIndex rStart = rowStart[iRow];
3964        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
3965        CoinBigIndex j;
3966        // Compute possible lower and upper ranges
3967     
3968        for (j = rStart; j < rEnd; ++j) {
3969          double value=element[j];
3970          iColumn = column[j];
3971          if (value > 0.0) {
3972            if (columnUpper_[iColumn] >= large) {
3973              ++infiniteUpper;
3974            } else {
3975              maximumUp += columnUpper_[iColumn] * value;
3976            }
3977            if (columnLower_[iColumn] <= -large) {
3978              ++infiniteLower;
3979            } else {
3980              maximumDown += columnLower_[iColumn] * value;
3981            }
3982          } else if (value<0.0) {
3983            if (columnUpper_[iColumn] >= large) {
3984              ++infiniteLower;
3985            } else {
3986              maximumDown += columnUpper_[iColumn] * value;
3987            }
3988            if (columnLower_[iColumn] <= -large) {
3989              ++infiniteUpper;
3990            } else {
3991              maximumUp += columnLower_[iColumn] * value;
3992            }
3993          }
3994        }
3995        // Build in a margin of error
3996        maximumUp += 1.0e-8*fabs(maximumUp);
3997        maximumDown -= 1.0e-8*fabs(maximumDown);
3998        double maxUp = maximumUp+infiniteUpper*1.0e31;
3999        double maxDown = maximumDown-infiniteLower*1.0e31;
4000        if (maxUp <= rowUpper_[iRow] + tolerance && 
4001            maxDown >= rowLower_[iRow] - tolerance) {
4002         
4003          // Row is redundant - make totally free
4004          // NO - can't do this for postsolve
4005          // rowLower_[iRow]=-COIN_DBL_MAX;
4006          // rowUpper_[iRow]=COIN_DBL_MAX;
4007          //printf("Redundant row in presolveX %d\n",iRow);
4008
4009        } else {
4010          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
4011              maxDown > rowUpper_[iRow]+100.0*tolerance) {
4012            // problem is infeasible - exit at once
4013            numberInfeasible++;
4014            break;
4015          }
4016          double lower = rowLower_[iRow];
4017          double upper = rowUpper_[iRow];
4018          for (j = rStart; j < rEnd; ++j) {
4019            double value=element[j];
4020            iColumn = column[j];
4021            double nowLower = columnLower_[iColumn];
4022            double nowUpper = columnUpper_[iColumn];
4023            if (value > 0.0) {
4024              // positive value
4025              if (lower>-large) {
4026                if (!infiniteUpper) {
4027                  assert(nowUpper < large2);
4028                  newBound = nowUpper + 
4029                    (lower - maximumUp) / value;
4030                  // relax if original was large
4031                  if (fabs(maximumUp)>1.0e8)
4032                    newBound -= 1.0e-12*fabs(maximumUp);
4033                } else if (infiniteUpper==1&&nowUpper>large) {
4034                  newBound = (lower -maximumUp) / value;
4035                  // relax if original was large
4036                  if (fabs(maximumUp)>1.0e8)
4037                    newBound -= 1.0e-12*fabs(maximumUp);
4038                } else {
4039                  newBound = -COIN_DBL_MAX;
4040                }
4041                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
4042                  // Tighten the lower bound
4043                  columnLower_[iColumn] = newBound;
4044                  numberChanged++;
4045                  // check infeasible (relaxed)
4046                  if (nowUpper - newBound < 
4047                      -100.0*tolerance) {
4048                    numberInfeasible++;
4049                  }
4050                  // adjust
4051                  double now;
4052                  if (nowLower<-large) {
4053                    now=0.0;
4054                    infiniteLower--;
4055                  } else {
4056                    now = nowLower;
4057                  }
4058                  maximumDown += (newBound-now) * value;
4059                  nowLower = newBound;
4060#ifdef DEBUG
4061                  checkCorrect(this,iRow,
4062                               element, rowStart, rowLength,
4063                               column,
4064                               columnLower_,  columnUpper_,
4065                               infiniteUpper,
4066                               infiniteLower,
4067                               maximumUp,
4068                               maximumDown);
4069#endif
4070                }
4071              } 
4072              if (upper <large) {
4073                if (!infiniteLower) {
4074                  assert(nowLower >- large2);
4075                  newBound = nowLower + 
4076                    (upper - maximumDown) / value;
4077                  // relax if original was large
4078                  if (fabs(maximumDown)>1.0e8)
4079                    newBound += 1.0e-12*fabs(maximumDown);
4080                } else if (infiniteLower==1&&nowLower<-large) {
4081                  newBound =   (upper - maximumDown) / value;
4082                  // relax if original was large
4083                  if (fabs(maximumDown)>1.0e8)
4084                    newBound += 1.0e-12*fabs(maximumDown);
4085                } else {
4086                  newBound = COIN_DBL_MAX;
4087                }
4088                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
4089                  // Tighten the upper bound
4090                  columnUpper_[iColumn] = newBound;
4091                  numberChanged++;
4092                  // check infeasible (relaxed)
4093                  if (newBound - nowLower < 
4094                      -100.0*tolerance) {
4095                    numberInfeasible++;
4096                  }
4097                  // adjust
4098                  double now;
4099                  if (nowUpper>large) {
4100                    now=0.0;
4101                    infiniteUpper--;
4102                  } else {
4103                    now = nowUpper;
4104                  }
4105                  maximumUp += (newBound-now) * value;
4106                  nowUpper = newBound;
4107#ifdef DEBUG
4108                  checkCorrect(this,iRow,
4109                               element, rowStart, rowLength,
4110                               column,
4111                               columnLower_,  columnUpper_,
4112                               infiniteUpper,
4113                               infiniteLower,
4114                               maximumUp,
4115                               maximumDown);
4116#endif
4117                }
4118              }
4119            } else {
4120              // negative value
4121              if (lower>-large) {
4122                if (!infiniteUpper) {
4123                  assert(nowLower < large2);
4124                  newBound = nowLower + 
4125                    (lower - maximumUp) / value;
4126                  // relax if original was large
4127                  if (fabs(maximumUp)>1.0e8)
4128                    newBound += 1.0e-12*fabs(maximumUp);
4129                } else if (infiniteUpper==1&&nowLower<-large) {
4130                  newBound = (lower -maximumUp) / value;
4131                  // relax if original was large
4132                  if (fabs(maximumUp)>1.0e8)
4133                    newBound += 1.0e-12*fabs(maximumUp);
4134                } else {
4135                  newBound = COIN_DBL_MAX;
4136                }
4137                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
4138                  // Tighten the upper bound
4139                  columnUpper_[iColumn] = newBound;
4140                  numberChanged++;
4141                  // check infeasible (relaxed)
4142                  if (newBound - nowLower < 
4143                      -100.0*tolerance) {
4144                    numberInfeasible++;
4145                  }
4146                  // adjust
4147                  double now;
4148                  if (nowUpper>large) {
4149                    now=0.0;
4150                    infiniteLower--;
4151                  } else {
4152                    now = nowUpper;
4153                  }
4154                  maximumDown += (newBound-now) * value;
4155                  nowUpper = newBound;
4156#ifdef DEBUG
4157                  checkCorrect(this,iRow,
4158                               element, rowStart, rowLength,
4159                               column,
4160                               columnLower_,  columnUpper_,
4161                               infiniteUpper,
4162                               infiniteLower,
4163                               maximumUp,
4164                               maximumDown);
4165#endif
4166                }
4167              }
4168              if (upper <large) {
4169                if (!infiniteLower) {
4170                  assert(nowUpper < large2);
4171                  newBound = nowUpper + 
4172                    (upper - maximumDown) / value;
4173                  // relax if original was large
4174                  if (fabs(maximumDown)>1.0e8)
4175                    newBound -= 1.0e-12*fabs(maximumDown);
4176                } else if (infiniteLower==1&&nowUpper>large) {
4177                  newBound =   (upper - maximumDown) / value;
4178                  // relax if original was large
4179                  if (fabs(maximumDown)>1.0e8)
4180                    newBound -= 1.0e-12*fabs(maximumDown);
4181                } else {
4182                  newBound = -COIN_DBL_MAX;
4183                }
4184                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
4185                  // Tighten the lower bound
4186                  columnLower_[iColumn] = newBound;
4187                  numberChanged++;
4188                  // check infeasible (relaxed)
4189                  if (nowUpper - newBound < 
4190                      -100.0*tolerance) {
4191                    numberInfeasible++;
4192                  }
4193                  // adjust
4194                  double now;
4195                  if (nowLower<-large) {
4196                    now=0.0;
4197                    infiniteUpper--;
4198                  } else {
4199                    now = nowLower;
4200                  }
4201                  maximumUp += (newBound-now) * value;
4202                  nowLower = newBound;
4203#ifdef DEBUG
4204                  checkCorrect(this,iRow,
4205                               element, rowStart, rowLength,
4206                               column,
4207                               columnLower_,  columnUpper_,
4208                               infiniteUpper,
4209                               infiniteLower,
4210                               maximumUp,
4211                               maximumDown);
4212#endif
4213                }
4214              }
4215            }
4216          }
4217        }
4218      }
4219    }
4220    totalTightened += numberChanged;
4221    if (iPass==1)
4222      numberCheck=numberChanged>>4;
4223    if (numberInfeasible) break;
4224  }
4225  if (!numberInfeasible) {
4226    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
4227      <<totalTightened
4228      <<CoinMessageEol;
4229    // Set bounds slightly loose
4230    double useTolerance = 1.0e-3;
4231    if (doTight>0) {
4232      if (doTight>10) { 
4233        useTolerance=0.0;
4234      } else {
4235        while (doTight) {
4236          useTolerance *= 0.1;
4237          doTight--;
4238        }
4239      }
4240    }
4241    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4242      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
4243        // Make large bounds stay infinite
4244        if (saveUpper[iColumn]>1.0e30&&columnUpper_[iColumn]>1.0e10) {
4245          columnUpper_[iColumn]=COIN_DBL_MAX;
4246        }
4247        if (saveLower[iColumn]<-1.0e30&&columnLower_[iColumn]<-1.0e10) {
4248          columnLower_[iColumn]=-COIN_DBL_MAX;
4249        }
4250        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e-8) {
4251          // relax enough so will have correct dj
4252#if 1
4253          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
4254                                    columnLower_[iColumn]-100.0*useTolerance);
4255          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
4256                                    columnUpper_[iColumn]+100.0*useTolerance);
4257#else
4258          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
4259            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
4260              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
4261            } else {
4262              columnLower_[iColumn]=saveLower[iColumn];
4263              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
4264                                        saveLower[iColumn]+100.0*useTolerance);
4265            }
4266          } else {
4267            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
4268              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
4269            } else {
4270              columnUpper_[iColumn]=saveUpper[iColumn];
4271              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
4272                                        saveUpper[iColumn]-100.0*useTolerance);
4273            }
4274          }
4275#endif
4276        } else {
4277          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
4278            // relax a bit
4279            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
4280                                        saveUpper[iColumn]);
4281          }
4282          if (columnLower_[iColumn]>saveLower[iColumn]) {
4283            // relax a bit
4284            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
4285                                        saveLower[iColumn]);
4286          }
4287        }
4288      }
4289    }
4290  } else {
4291    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
4292      <<numberInfeasible
4293      <<CoinMessageEol;
4294    // restore column bounds
4295    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
4296    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
4297  }
4298  delete [] saveLower;
4299  delete [] saveUpper;
4300  return (numberInfeasible);
4301}
4302//#define SAVE_AND_RESTORE
4303// dual
4304#include "ClpSimplexDual.hpp"
4305#include "ClpSimplexPrimal.hpp"
4306#ifndef SAVE_AND_RESTORE
4307int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
4308#else
4309int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
4310{
4311  // May be empty problem
4312  if (numberRows_&&numberColumns_) {
4313    // Save on file for debug
4314    int returnCode;
4315    returnCode = saveModel("debug.sav");
4316    if (returnCode) {
4317      printf("** Unable to save model to debug.sav\n");
4318      abort();
4319    }
4320    ClpSimplex temp;
4321    returnCode=temp.restoreModel("debug.sav");
4322    if (returnCode) {
4323      printf("** Unable to restore model from debug.sav\n");
4324      abort();
4325    }
4326    temp.setLogLevel(handler_->logLevel());
4327    // Now do dual
4328    returnCode=temp.dualDebug(ifValuesPass,startFinishOptions);
4329    // Move status and solution back
4330    int numberTotal = numberRows_+numberColumns_;
4331    memcpy(status_,temp.statusArray(),numberTotal);
4332    memcpy(columnActivity_,temp.primalColumnSolution(),numberColumns_*sizeof(double));
4333    memcpy(rowActivity_,temp.primalRowSolution(),numberRows_*sizeof(double));
4334    memcpy(reducedCost_,temp.dualColumnSolution(),numberColumns_*sizeof(double));
4335    memcpy(dual_,temp.dualRowSolution(),numberRows_*sizeof(double));
4336    problemStatus_ = temp.problemStatus_;
4337    setObjectiveValue(temp.objectiveValue());
4338    setSumDualInfeasibilities(temp.sumDualInfeasibilities());
4339    setNumberDualInfeasibilities(temp.numberDualInfeasibilities());
4340    setSumPrimalInfeasibilities(temp.sumPrimalInfeasibilities());
4341    setNumberPrimalInfeasibilities(temp.numberPrimalInfeasibilities());
4342    setNumberIterations(temp.numberIterations());
4343    return returnCode;
4344  } else {
4345    // empty
4346    return dualDebug(ifValuesPass,startFinishOptions);
4347  }
4348}
4349int ClpSimplex::dualDebug (int ifValuesPass , int startFinishOptions)
4350#endif
4351{
4352  //double savedPivotTolerance = factorization_->pivotTolerance();
4353  int saveQuadraticActivated = objective_->activated();
4354  objective_->setActivated(0);
4355  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
4356  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
4357      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
4358      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
4359      additional data and have no destructor or (non-default) constructor.
4360
4361      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4362      in primal and dual.
4363
4364      As far as I can see this is perfectly safe.
4365  */
4366  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass, startFinishOptions);
4367  if ((specialOptions_&2048)!=0&&problemStatus_==10&&!numberPrimalInfeasibilities_
4368      &&sumDualInfeasibilities_<1000.0*dualTolerance_&&perturbation_>=100)
4369    problemStatus_=0; // ignore
4370  if (problemStatus_==10) {
4371    //printf("Cleaning up with primal\n");
4372    int savePerturbation = perturbation_;
4373    int saveLog = handler_->logLevel();
4374    //handler_->setLogLevel(63);
4375    perturbation_=100;
4376    bool denseFactorization = initialDenseFactorization();
4377    // It will be safe to allow dense
4378    setInitialDenseFactorization(true);
4379    // Allow for catastrophe
4380    int saveMax = intParam_[ClpMaxNumIteration];
4381    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
4382      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
4383    // check which algorithms allowed
4384    int dummy;
4385    if (problemStatus_==10)
4386      startFinishOptions |= 2;
4387    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
4388      returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
4389    else
4390      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
4391    if (problemStatus_==3&&numberIterations_<saveMax) {
4392      if (handler_->logLevel()==63)
4393        printf("looks like trouble - too many iterations in clean up - trying again\n");
4394      // flatten solution and try again
4395      int iRow,iColumn;
4396      for (iRow=0;iRow<numberRows_;iRow++) {
4397        if (getRowStatus(iRow)!=basic) {
4398          setRowStatus(iRow,superBasic);
4399          // but put to bound if close
4400          if (fabs(rowActivity_[iRow]-rowLower_[iRow])
4401              <=primalTolerance_) {
4402            rowActivity_[iRow]=rowLower_[iRow];
4403            setRowStatus(iRow,atLowerBound);
4404          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])
4405                     <=primalTolerance_) {
4406            rowActivity_[iRow]=rowUpper_[iRow];
4407            setRowStatus(iRow,atUpperBound);
4408          }
4409        }
4410      }
4411      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4412        if (getColumnStatus(iColumn)!=basic) {
4413          setColumnStatus(iColumn,superBasic);
4414          // but put to bound if close
4415          if (fabs(columnActivity_[iColumn]-columnLower_[iColumn])
4416              <=primalTolerance_) {
4417            columnActivity_[iColumn]=columnLower_[iColumn];
4418            setColumnStatus(iColumn,atLowerBound);
4419          } else if (fabs(columnActivity_[iColumn]
4420                          -columnUpper_[iColumn])
4421                     <=primalTolerance_) {
4422            columnActivity_[iColumn]=columnUpper_[iColumn];
4423            setColumnStatus(iColumn,atUpperBound);
4424          }
4425        }
4426      }
4427      problemStatus_=-1;
4428      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 
4429                                          2*numberRows_+numberColumns_,saveMax);
4430      perturbation_=savePerturbation;
4431      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
4432      if (problemStatus_==3&&numberIterations_<saveMax&& 
4433          handler_->logLevel()>0)
4434        printf("looks like real trouble - too many iterations in second clean up - giving up\n");
4435    }
4436    intParam_[ClpMaxNumIteration] = saveMax;
4437
4438    setInitialDenseFactorization(denseFactorization);
4439    perturbation_=savePerturbation;
4440    if (problemStatus_==10) 
4441      problemStatus_=0;
4442    handler_->setLogLevel(saveLog);
4443  }
4444  objective_->setActivated(saveQuadraticActivated);
4445  //factorization_->pivotTolerance(savedPivotTolerance);
4446  return returnCode;
4447}
4448// primal
4449int ClpSimplex::primal (int ifValuesPass , int startFinishOptions)
4450{
4451  //double savedPivotTolerance = factorization_->pivotTolerance();
4452#ifndef SLIM_CLP
4453  // See if nonlinear
4454  if (objective_->type()>1&&objective_->activated()) 
4455    return reducedGradient();
4456#endif 
4457  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
4458  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
4459      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
4460      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
4461      additional data and have no destructor or (non-default) constructor.
4462
4463      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4464      in primal and dual.
4465
4466      As far as I can see this is perfectly safe.
4467  */
4468  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass,startFinishOptions);
4469  if (problemStatus_==10) {
4470    //printf("Cleaning up with dual\n");
4471    int savePerturbation = perturbation_;
4472    perturbation_=100;
4473    bool denseFactorization = initialDenseFactorization();
4474    // It will be safe to allow dense
4475    setInitialDenseFactorization(true);
4476    // check which algorithms allowed
4477    int dummy;
4478    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0&&(specialOptions_&8192)==0) {
4479      double saveBound = dualBound_;
4480      // upperOut_ has largest away from bound
4481      dualBound_=CoinMin(CoinMax(2.0*upperOut_,1.0e8),dualBound_);
4482      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
4483      dualBound_=saveBound;
4484    } else {
4485      returnCode = ((ClpSimplexPrimal *) this)->primal(0,startFinishOptions);
4486    }
4487    setInitialDenseFactorization(denseFactorization);
4488    perturbation_=savePerturbation;
4489    if (problemStatus_==10) 
4490      problemStatus_=0;
4491  }
4492  //factorization_->pivotTolerance(savedPivotTolerance);
4493  return returnCode;
4494}
4495#ifndef SLIM_CLP
4496#include "ClpQuadraticObjective.hpp"
4497/* Dual ranging.
4498   This computes increase/decrease in cost for each given variable and corresponding
4499   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
4500   and numberColumns.. for artificials/slacks.
4501   For non-basic variables the sequence number will be that of the non-basic variables.
4502   
4503   Up to user to provide correct length arrays.
4504   
4505   Returns non-zero if infeasible unbounded etc
4506*/
4507#include "ClpSimplexOther.hpp"
4508int ClpSimplex::dualRanging(int numberCheck,const int * which,
4509                            double * costIncrease, int * sequenceIncrease,
4510                            double * costDecrease, int * sequenceDecrease)
4511{
4512  int savePerturbation = perturbation_;
4513  perturbation_=100;
4514  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4515  if (problemStatus_==10) {
4516    //printf("Cleaning up with dual\n");
4517    bool denseFactorization = initialDenseFactorization();
4518    // It will be safe to allow dense
4519    setInitialDenseFactorization(true);
4520    // check which algorithms allowed
4521    int dummy;
4522    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
4523      // upperOut_ has largest away from bound
4524      double saveBound=dualBound_;
4525      if (upperOut_>0.0)
4526        dualBound_=2.0*upperOut_;
4527      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
4528      dualBound_=saveBound;
4529    } else {
4530      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4531    }
4532    setInitialDenseFactorization(denseFactorization);
4533    if (problemStatus_==10) 
4534      problemStatus_=0;
4535  }
4536  perturbation_=savePerturbation;
4537  if (problemStatus_||secondaryStatus_==6) {
4538    finish(); // get rid of arrays
4539    return 1; // odd status
4540  }
4541  ((ClpSimplexOther *) this)->dualRanging(numberCheck,which,
4542                                          costIncrease,sequenceIncrease,
4543                                          costDecrease,sequenceDecrease);
4544  finish(); // get rid of arrays
4545  return 0;
4546}
4547/* Primal ranging.
4548   This computes increase/decrease in value for each given variable and corresponding
4549   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
4550   and numberColumns.. for artificials/slacks.
4551   For basic variables the sequence number will be that of the basic variables.
4552   
4553   Up to user to providen correct length arrays.
4554   
4555   Returns non-zero if infeasible unbounded etc
4556*/
4557int ClpSimplex::primalRanging(int numberCheck,const int * which,
4558                  double * valueIncrease, int * sequenceIncrease,
4559                  double * valueDecrease, int * sequenceDecrease)
4560{
4561  int savePerturbation = perturbation_;
4562  perturbation_=100;
4563  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4564  if (problemStatus_==10) {
4565    //printf("Cleaning up with dual\n");
4566    bool denseFactorization = initialDenseFactorization();
4567    // It will be safe to allow dense
4568    setInitialDenseFactorization(true);
4569    // check which algorithms allowed
4570    int dummy;
4571    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
4572      // upperOut_ has largest away from bound
4573      double saveBound=dualBound_;
4574      if (upperOut_>0.0)
4575        dualBound_=2.0*upperOut_;
4576      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
4577      dualBound_=saveBound;
4578    } else {
4579      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
4580    }
4581    setInitialDenseFactorization(denseFactorization);
4582    if (problemStatus_==10) 
4583      problemStatus_=0;
4584  }
4585  perturbation_=savePerturbation;
4586  if (problemStatus_||secondaryStatus_==6) {
4587    finish(); // get rid of arrays
4588    return 1; // odd status
4589  }
4590  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
4591                                          valueIncrease,sequenceIncrease,
4592                                          valueDecrease,sequenceDecrease);
4593  finish(); // get rid of arrays
4594  return 0;
4595}
4596/* Write the basis in MPS format to the specified file.
4597   If writeValues true writes values of structurals
4598   (and adds VALUES to end of NAME card)
4599   
4600   Row and column names may be null.
4601   formatType is
4602   <ul>
4603   <li> 0 - normal
4604   <li> 1 - extra accuracy
4605   <li> 2 - IEEE hex (later)
4606   </ul>
4607   
4608   Returns non-zero on I/O error
4609*/
4610int 
4611ClpSimplex::writeBasis(const char *filename,
4612                            bool writeValues,
4613                            int formatType) const
4614{
4615  return ((const ClpSimplexOther *) this)->writeBasis(filename,writeValues,
4616                                         formatType);
4617}
4618// Read a basis from the given filename
4619int 
4620ClpSimplex::readBasis(const char *filename)
4621{
4622  return ((ClpSimplexOther *) this)->readBasis(filename);
4623}
4624#include "ClpSimplexNonlinear.hpp"
4625/* Solves nonlinear problem using SLP - may be used as crash
4626   for other algorithms when number of iterations small
4627*/
4628int 
4629ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
4630{
4631  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
4632}
4633// Solves non-linear using reduced gradient
4634int ClpSimplex::reducedGradient(int phase)
4635{
4636  if (objective_->type()<2||!objective_->activated()) {
4637    // no quadratic part
4638    return primal(0);
4639  }
4640  // get feasible
4641  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
4642    objective_->setActivated(0);
4643    double saveDirection = optimizationDirection();
4644    setOptimizationDirection(0.0);
4645    primal(1);
4646    setOptimizationDirection(saveDirection);
4647    objective_->setActivated(1);
4648    // still infeasible
4649    if (numberPrimalInfeasibilities())
4650      return 0;
4651  }
4652  // Now enter method
4653  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
4654  return returnCode;
4655}
4656#include "ClpPredictorCorrector.hpp"
4657#include "ClpCholeskyBase.hpp"
4658// Preference is WSSMP, TAUCS, UFL (just ordering) then base
4659#if WSSMP_BARRIER
4660#include "ClpCholeskyWssmp.hpp"
4661#include "ClpCholeskyWssmpKKT.hpp"
4662#elif UFL_BARRIER
4663#include "ClpCholeskyUfl.hpp"
4664#elif TAUCS_BARRIER
4665#include "ClpCholeskyTaucs.hpp"
4666#endif
4667#include "ClpPresolve.hpp"
4668/* Solves using barrier (assumes you have good cholesky factor code).
4669   Does crossover to simplex if asked*/
4670int 
4671ClpSimplex::barrier(bool crossover)
4672{
4673  ClpSimplex * model2 = this;
4674  int savePerturbation=perturbation_;
4675  ClpInterior barrier;
4676  barrier.borrowModel(*model2);
4677  // See if quadratic objective
4678  ClpQuadraticObjective * quadraticObj = NULL;
4679  if (objective_->type()==2)
4680    quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
4681  // If Quadratic we need KKT
4682  bool doKKT = (quadraticObj!=NULL);
4683  // Preference is WSSMP, UFL, TAUCS then base
4684#ifdef WSSMP_BARRIER
4685 if (!doKKT) {
4686   ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
4687   barrier.setCholesky(cholesky);
4688 } else {
4689  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
4690   ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
4691   barrier.setCholesky(cholesky);
4692 }
4693#elif UFL_BARRIER
4694 if (!doKKT) {
4695   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
4696   barrier.setCholesky(cholesky);
4697 } else {
4698   ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
4699   cholesky->setKKT(true);
4700   barrier.setCholesky(cholesky);
4701 }
4702#elif TAUCS_BARRIER
4703 assert (!doKKT);
4704 ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
4705 barrier.setCholesky(cholesky);
4706#else
4707 if (!doKKT) {
4708  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
4709  barrier.setCholesky(cholesky);
4710 } else {
4711   ClpCholeskyBase * cholesky = new ClpCholeskyBase();
4712   cholesky->setKKT(true);
4713   barrier.setCholesky(cholesky);
4714 }
4715#endif
4716  barrier.setDiagonalPerturbation(1.0e-14);
4717  int numberRows = model2->numberRows();
4718  int numberColumns = model2->numberColumns();
4719  int saveMaxIts = model2->maximumIterations();
4720  if (saveMaxIts<1000) {
4721    barrier.setMaximumBarrierIterations(saveMaxIts);
4722    model2->setMaximumIterations(1000000);
4723  }
4724  barrier.primalDual();
4725  int barrierStatus=barrier.status();
4726  double gap = barrier.complementarityGap();
4727  // get which variables are fixed
4728  double * saveLower=NULL;
4729  double * saveUpper=NULL;
4730  ClpPresolve pinfo2;
4731  ClpSimplex * saveModel2=NULL;
4732  int numberFixed = barrier.numberFixed();
4733  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
4734    // may as well do presolve
4735    int numberRows = barrier.numberRows();
4736    int numberColumns = barrier.numberColumns();
4737    int numberTotal = numberRows+numberColumns;
4738    saveLower = new double [numberTotal];
4739    saveUpper = new double [numberTotal];
4740    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
4741    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
4742    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
4743    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
4744    barrier.fixFixed();
4745    saveModel2=model2;
4746  }
4747  barrier.returnModel(*model2);
4748  double * rowPrimal = new double [numberRows];
4749  double * columnPrimal = new double [numberColumns];
4750  double * rowDual = new double [numberRows];
4751  double * columnDual = new double [numberColumns];
4752  // move solutions other way
4753  CoinMemcpyN(model2->primalRowSolution(),
4754              numberRows,rowPrimal);
4755  CoinMemcpyN(model2->dualRowSolution(),
4756              numberRows,rowDual);
4757  CoinMemcpyN(model2->primalColumnSolution(),
4758              numberColumns,columnPrimal);
4759  CoinMemcpyN(model2->dualColumnSolution(),
4760              numberColumns,columnDual);
4761  if (saveModel2) {
4762    // do presolve
4763    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
4764                                   false,5,true);
4765  }
4766  if (barrierStatus<4&&crossover) {
4767    // make sure no status left
4768    model2->createStatus();
4769    // solve
4770    model2->setPerturbation(100);
4771    // throw some into basis
4772    {
4773      int numberRows = model2->numberRows();
4774      int numberColumns = model2->numberColumns();
4775      double * dsort = new double[numberColumns];
4776      int * sort = new int[numberColumns];
4777      int n=0;
4778      const double * columnLower = model2->columnLower();
4779      const double * columnUpper = model2->columnUpper();
4780      const double * primalSolution = model2->primalColumnSolution();
4781      double tolerance = 10.0*primalTolerance_;
4782      int i;
4783      for ( i=0;i<numberRows;i++) 
4784        model2->setRowStatus(i,superBasic);
4785      for ( i=0;i<numberColumns;i++) {
4786        double distance = CoinMin(columnUpper[i]-primalSolution[i],
4787                              primalSolution[i]-columnLower[i]);
4788        if (distance>tolerance) {
4789          dsort[n]=-distance;
4790          sort[n++]=i;
4791          model2->setStatus(i,superBasic);
4792        } else if (distance>primalTolerance_) {
4793          model2->setStatus(i,superBasic);
4794        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
4795          model2->setStatus(i,atLowerBound);
4796        } else {
4797          model2->setStatus(i,atUpperBound);
4798        }
4799      }
4800      CoinSort_2(dsort,dsort+n,sort);
4801      n = CoinMin(numberRows,n);
4802      for ( i=0;i<n;i++) {
4803        int iColumn = sort[i];
4804        model2->setStatus(iColumn,basic);
4805      }
4806      delete [] sort;
4807      delete [] dsort;
4808    }
4809    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
4810      int numberRows = model2->numberRows();
4811      int numberColumns = model2->numberColumns();
4812      // just primal values pass
4813      double saveScale = model2->objectiveScale();
4814      model2->setObjectiveScale(1.0e-3);
4815      model2->primal(2);
4816      model2->setObjectiveScale(saveScale);
4817      // save primal solution and copy back dual
4818      CoinMemcpyN(model2->primalRowSolution(),
4819                  numberRows,rowPrimal);
4820      CoinMemcpyN(rowDual,
4821                  numberRows,model2->dualRowSolution());
4822      CoinMemcpyN(model2->primalColumnSolution(),
4823                  numberColumns,columnPrimal);
4824      CoinMemcpyN(columnDual,
4825                  numberColumns,model2->dualColumnSolution());
4826      //model2->primal(1);
4827      // clean up reduced costs and flag variables
4828      {
4829        double * dj = model2->dualColumnSolution();
4830        double * cost = model2->objective();
4831        double * saveCost = new double[numberColumns];
4832        memcpy(saveCost,cost,numberColumns*sizeof(double));
4833        double * saveLower = new double[numberColumns];
4834        double * lower = model2->columnLower();
4835        memcpy(saveLower,lower,numberColumns*sizeof(double));
4836        double * saveUpper = new double[numberColumns];
4837        double * upper = model2->columnUpper();
4838        memcpy(saveUpper,upper,numberColumns*sizeof(double));
4839        int i;
4840        double tolerance = 10.0*dualTolerance_;
4841        for ( i=0;i<numberColumns;i++) {
4842          if (model2->getStatus(i)==basic) {
4843            dj[i]=0.0;
4844          } else if (model2->getStatus(i)==atLowerBound) {
4845            if (optimizationDirection_*dj[i]<tolerance) {
4846              if (optimizationDirection_*dj[i]<0.0) {
4847                //if (dj[i]<-1.0e-3)
4848                //printf("bad dj at lb %d %g\n",i,dj[i]);
4849                cost[i] -= dj[i];
4850                dj[i]=0.0;
4851              }
4852            } else {
4853              upper[i]=lower[i];
4854            }
4855          } else if (model2->getStatus(i)==atUpperBound) {
4856            if (optimizationDirection_*dj[i]>tolerance) {
4857              if (optimizationDirection_*dj[i]>0.0) {
4858                //if (dj[i]>1.0e-3)
4859                //printf("bad dj at ub %d %g\n",i,dj[i]);
4860                cost[i] -= dj[i];
4861                dj[i]=0.0;
4862              }
4863            } else {
4864              lower[i]=upper[i];
4865            }
4866          }
4867        }
4868        // just dual values pass
4869        //model2->setLogLevel(63);
4870        //model2->setFactorizationFrequency(1);
4871        model2->dual(2);
4872        memcpy(cost,saveCost,numberColumns*sizeof(double));
4873        delete [] saveCost;
4874        memcpy(lower,saveLower,numberColumns*sizeof(double));
4875        delete [] saveLower;
4876        memcpy(upper,saveUpper,numberColumns*sizeof(double));
4877        delete [] saveUpper;
4878      }
4879      // and finish
4880      // move solutions
4881      CoinMemcpyN(rowPrimal,
4882                  numberRows,model2->primalRowSolution());
4883      CoinMemcpyN(columnPrimal,
4884                  numberColumns,model2->primalColumnSolution());
4885    }
4886//     double saveScale = model2->objectiveScale();
4887//     model2->setObjectiveScale(1.0e-3);
4888//     model2->primal(2);
4889//    model2->setObjectiveScale(saveScale);
4890    model2->primal(1);
4891  } else if (barrierStatus==4&&crossover) {
4892    // memory problems
4893    model2->setPerturbation(savePerturbation);
4894    model2->createStatus();
4895    model2->dual();
4896  }
4897  model2->setMaximumIterations(saveMaxIts);
4898  delete [] rowPrimal;
4899  delete [] columnPrimal;
4900  delete [] rowDual;
4901  delete [] columnDual;
4902  if (saveLower) {
4903    pinfo2.postsolve(true);
4904    delete model2;
4905    model2=saveModel2;
4906    int numberRows = model2->numberRows();
4907    int numberColumns = model2->numberColumns();
4908    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
4909    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
4910    delete [] saveLower;
4911    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
4912    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
4913    delete [] saveUpper;
4914    model2->primal(1);
4915  }
4916  model2->setPerturbation(savePerturbation);
4917  return model2->status();
4918}
4919/* For strong branching.  On input lower and upper are new bounds
4920   while on output they are objective function values (>1.0e50 infeasible).
4921   Return code is 0 if nothing interesting, -1 if infeasible both
4922   ways and +1 if infeasible one way (check values to see which one(s))
4923*/
4924int ClpSimplex::strongBranching(int numberVariables,const int * variables,
4925                                double * newLower, double * newUpper,
4926                                double ** outputSolution,
4927                                int * outputStatus, int * outputIterations,
4928                                bool stopOnFirstInfeasible,
4929                                bool alwaysFinish,
4930                                int startFinishOptions)
4931{
4932  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
4933                                                    newLower,  newUpper,outputSolution,
4934                                                    outputStatus, outputIterations,
4935                                                    stopOnFirstInfeasible,
4936                                                    alwaysFinish,startFinishOptions);
4937}
4938#endif
4939/* Borrow model.  This is so we dont have to copy large amounts
4940   of data around.  It assumes a derived class wants to overwrite
4941   an empty model with a real one - while it does an algorithm.
4942   This is same as ClpModel one, but sets scaling on etc. */
4943void 
4944ClpSimplex::borrowModel(ClpModel & otherModel) 
4945{
4946  ClpModel::borrowModel(otherModel);
4947  createStatus();
4948  //ClpDualRowSteepest steep1;
4949  //setDualRowPivotAlgorithm(steep1);
4950  //ClpPrimalColumnSteepest steep2;
4951  //setPrimalColumnPivotAlgorithm(steep2);
4952}
4953void 
4954ClpSimplex::borrowModel(ClpSimplex & otherModel) 
4955{
4956  ClpModel::borrowModel(otherModel);
4957  createStatus();
4958  dualBound_ = otherModel.dualBound_;
4959  dualTolerance_ = otherModel.dualTolerance_;
4960  primalTolerance_ = otherModel.primalTolerance_;
4961  delete dualRowPivot_;
4962  dualRowPivot_ = otherModel.dualRowPivot_->clone(false);
4963  delete primalColumnPivot_;
4964  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(false);
4965  perturbation_ = otherModel.perturbation_;
4966  specialOptions_ = otherModel.specialOptions_;
4967  automaticScale_ = otherModel.automaticScale_;
4968}
4969/// Saves scalars for ClpSimplex
4970typedef struct {
4971  double optimizationDirection;
4972  double dblParam[ClpLastDblParam];
4973  double objectiveValue;
4974  double dualBound;
4975  double dualTolerance;
4976  double primalTolerance;
4977  double sumDualInfeasibilities;
4978  double sumPrimalInfeasibilities;
4979  double infeasibilityCost;
4980  int numberRows;
4981  int numberColumns;
4982  int intParam[ClpLastIntParam];
4983  int numberIterations;
4984  int problemStatus;
4985  int maximumIterations;
4986  int lengthNames;
4987  int numberDualInfeasibilities;
4988  int numberDualInfeasibilitiesWithoutFree;
4989  int numberPrimalInfeasibilities;
4990  int numberRefinements;
4991  int scalingFlag;
4992  int algorithm;
4993  unsigned int specialOptions;
4994  int dualPivotChoice;
4995  int primalPivotChoice;
4996  int matrixStorageChoice;
4997} Clp_scalars;
4998#ifndef SLIM_NOIO
4999int outDoubleArray(double * array, int length, FILE * fp)
5000{
5001  int numberWritten;
5002  if (array&&length) {
5003    numberWritten = fwrite(&length,sizeof(int),1,fp);
5004    if (numberWritten!=1)
5005      return 1;
5006    numberWritten = fwrite(array,sizeof(double),length,fp);
5007    if (numberWritten!=length)
5008      return 1;
5009  } else {
5010    length = 0;
5011    numberWritten = fwrite(&length,sizeof(int),1,fp);
5012    if (numberWritten!=1)
5013      return 1;
5014  }
5015  return 0;
5016}
5017// Save model to file, returns 0 if success
5018int
5019ClpSimplex::saveModel(const char * fileName)
5020{
5021  FILE * fp = fopen(fileName,"wb");
5022  if (fp) {
5023    Clp_scalars scalars;
5024    CoinBigIndex numberWritten;
5025    // Fill in scalars
5026    scalars.optimizationDirection = optimizationDirection_;
5027    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
5028    scalars.objectiveValue = objectiveValue_;
5029    scalars.dualBound = dualBound_;
5030    scalars.dualTolerance = dualTolerance_;
5031    scalars.primalTolerance = primalTolerance_;
5032    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
5033    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
5034    scalars.infeasibilityCost = infeasibilityCost_;
5035    scalars.numberRows = numberRows_;
5036    scalars.numberColumns = numberColumns_;
5037    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
5038    scalars.numberIterations = numberIterations_;
5039    scalars.problemStatus = problemStatus_;
5040    scalars.maximumIterations = maximumIterations();
5041    scalars.lengthNames = lengthNames_;
5042    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
5043    scalars.numberDualInfeasibilitiesWithoutFree
5044      = numberDualInfeasibilitiesWithoutFree_;
5045    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
5046    scalars.numberRefinements = numberRefinements_;
5047    scalars.scalingFlag = scalingFlag_;
5048    scalars.algorithm = algorithm_;
5049    scalars.specialOptions = specialOptions_;
5050    scalars.dualPivotChoice = dualRowPivot_->type();
5051    scalars.primalPivotChoice = primalColumnPivot_->type();
5052    scalars.matrixStorageChoice = matrix_->type();
5053
5054    // put out scalars
5055    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
5056    if (numberWritten!=1)
5057      return 1;
5058    CoinBigIndex length;
5059#ifndef CLP_NO_STD
5060    int i;
5061    // strings
5062    for (i=0;i<ClpLastStrParam;i++) {
5063      length = strParam_[i].size();
5064      numberWritten = fwrite(&length,sizeof(int),1,fp);
5065      if (numberWritten!=1)
5066        return 1;
5067      if (length) {
5068        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
5069        if (numberWritten!=1)
5070          return 1;
5071      }
5072    }
5073#endif
5074    // arrays - in no particular order
5075    if (outDoubleArray(rowActivity_,numberRows_,fp))
5076        return 1;
5077    if (outDoubleArray(columnActivity_,numberColumns_,fp))
5078        return 1;
5079    if (outDoubleArray(dual_,numberRows_,fp))
5080        return 1;
5081    if (outDoubleArray(reducedCost_,numberColumns_,fp))
5082        return 1;
5083    if (outDoubleArray(rowLower_,numberRows_,fp))
5084        return 1;
5085    if (outDoubleArray(rowUpper_,numberRows_,fp))
5086        return 1;
5087    if (outDoubleArray(objective(),numberColumns_,fp))
5088        return 1;
5089    if (outDoubleArray(rowObjective_,numberRows_,fp))
5090        return 1;
5091    if (outDoubleArray(columnLower_,numberColumns_,fp))
5092        return 1;
5093    if (outDoubleArray(columnUpper_,numberColumns_,fp))
5094        return 1;
5095    if (ray_) {
5096      if (problemStatus_==1)
5097        if (outDoubleArray(ray_,numberRows_,fp))
5098          return 1;
5099      else if (problemStatus_==2)
5100        if (outDoubleArray(ray_,numberColumns_,fp))
5101          return 1;
5102      else
5103        if (outDoubleArray(NULL,0,fp))
5104          return 1;
5105    } else {
5106      if (outDoubleArray(NULL,0,fp))
5107        return 1;
5108    }
5109    if (status_&&(numberRows_+numberColumns_)>0) {
5110      length = numberRows_+numberColumns_;
5111      numberWritten = fwrite(&length,sizeof(int),1,fp);
5112      if (numberWritten!=1)
5113        return 1;
5114      numberWritten = fwrite(status_,sizeof(char),length, fp);
5115      if (numberWritten!=length)
5116        return 1;
5117    } else {
5118      length = 0;
5119      numberWritten = fwrite(&length,sizeof(int),1,fp);
5120      if (numberWritten!=1)
5121        return 1;
5122    }
5123#ifndef CLP_NO_STD
5124    if (lengthNames_) {
5125      char * array = 
5126        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
5127      char * put = array;
5128      CoinAssert (numberRows_ == (int) rowNames_.size());
5129      for (i=0;i<numberRows_;i++) {
5130        assert((int)rowNames_[i].size()<=lengthNames_);
5131        strcpy(put,rowNames_[i].c_str());
5132        put += lengthNames_+1;
5133      }
5134      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
5135      if (numberWritten!=numberRows_)
5136        return 1;
5137      put=array;
5138      CoinAssert (numberColumns_ == (int) columnNames_.size());
5139      for (i=0;i<numberColumns_;i++) {
5140        assert((int) columnNames_[i].size()<=lengthNames_);
5141        strcpy(put,columnNames_[i].c_str());
5142        put += lengthNames_+1;
5143      }
5144      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
5145      if (numberWritten!=numberColumns_)
5146        return 1;
5147      delete [] array;
5148    }
5149#endif
5150    // just standard type at present
5151    assert (matrix_->type()==1);
5152    CoinAssert (matrix_->getNumCols() == numberColumns_);
5153    CoinAssert (matrix_->getNumRows() == numberRows_);
5154    // we are going to save with gaps
5155    length = matrix_->getVectorStarts()[numberColumns_-1]
5156      + matrix_->getVectorLengths()[numberColumns_-1];
5157    numberWritten = fwrite(&length,sizeof(int),1,fp);
5158    if (numberWritten!=1)
5159      return 1;
5160    numberWritten = fwrite(matrix_->getElements(),
5161                           sizeof(double),length,fp);
5162    if (numberWritten!=length)
5163      return 1;
5164    numberWritten = fwrite(matrix_->getIndices(),
5165                           sizeof(int),length,fp);
5166    if (numberWritten!=length)
5167      return 1;
5168    numberWritten = fwrite(matrix_->getVectorStarts(),
5169                           sizeof(int),numberColumns_+1,fp);
5170    if (numberWritten!=numberColumns_+1)
5171      return 1;
5172    numberWritten = fwrite(matrix_->getVectorLengths(),
5173                           sizeof(int),numberColumns_,fp);
5174    if (numberWritten!=numberColumns_)
5175      return 1;
5176    // finished
5177    fclose(fp);
5178    return 0;
5179  } else {
5180    return -1;
5181  }
5182}
5183
5184int inDoubleArray(double * &array, int length, FILE * fp)
5185{
5186  int numberRead;
5187  int length2;
5188  numberRead = fread(&length2,sizeof(int),1,fp);
5189  if (numberRead!=1)
5190    return 1;
5191  if (length2) {
5192    // lengths must match
5193    if (length!=length2)
5194      return 2;
5195    array = new double[length];
5196    numberRead = fread(array,sizeof(double),length,fp);
5197    if (numberRead!=length)
5198      return 1;
5199  } 
5200  return 0;
5201}
5202/* Restore model from file, returns 0 if success,
5203   deletes current model */
5204int 
5205ClpSimplex::restoreModel(const char * fileName)
5206{
5207  FILE * fp = fopen(fileName,"rb");
5208  if (fp) {
5209    // Get rid of current model
5210    // save event handler in case already set
5211    ClpEventHandler * handler = eventHandler_->clone();
5212    ClpModel::gutsOfDelete();
5213    eventHandler_ = handler;
5214    gutsOfDelete(0);
5215    int i;
5216    for (i=0;i<6;i++) {
5217      rowArray_[i]=NULL;
5218      columnArray_[i]=NULL;
5219    }
5220    // get an empty factorization so we can set tolerances etc
5221    factorization_ = new ClpFactorization();
5222    // Say sparse
5223    factorization_->sparseThreshold(1);
5224    Clp_scalars scalars;
5225    CoinBigIndex numberRead;
5226
5227    // get scalars
5228    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
5229    if (numberRead!=1)
5230      return 1;
5231    // Fill in scalars
5232    optimizationDirection_ = scalars.optimizationDirection;
5233    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
5234    objectiveValue_ = scalars.objectiveValue;
5235    dualBound_ = scalars.dualBound;
5236    dualTolerance_ = scalars.dualTolerance;
5237    primalTolerance_ = scalars.primalTolerance;
5238    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
5239    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
5240    infeasibilityCost_ = scalars.infeasibilityCost;
5241    numberRows_ = scalars.numberRows;
5242    numberColumns_ = scalars.numberColumns;
5243    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
5244    numberIterations_ = scalars.numberIterations;
5245    problemStatus_ = scalars.problemStatus;
5246    setMaximumIterations(scalars.maximumIterations);
5247    lengthNames_ = scalars.lengthNames;
5248    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
5249    numberDualInfeasibilitiesWithoutFree_
5250      = scalars.numberDualInfeasibilitiesWithoutFree;
5251    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
5252    numberRefinements_ = scalars.numberRefinements;
5253    scalingFlag_ = scalars.scalingFlag;
5254    algorithm_ = scalars.algorithm;
5255    specialOptions_ = scalars.specialOptions;
5256    // strings
5257    CoinBigIndex length;
5258#ifndef CLP_NO_STD
5259    for (i=0;i<ClpLastStrParam;i++) {
5260      numberRead = fread(&length,sizeof(int),1,fp);
5261      if (numberRead!=1)
5262        return 1;
5263      if (length) {
5264        char * array = new char[length+1];
5265        numberRead = fread(array,length,1,fp);
5266        if (numberRead!=1)
5267          return 1;
5268        array[length]='\0';
5269        strParam_[i]=array;
5270        delete [] array;
5271      }
5272    }
5273#endif
5274    // arrays - in no particular order
5275    if (inDoubleArray(rowActivity_,numberRows_,fp))
5276        return 1;
5277    if (inDoubleArray(columnActivity_,numberColumns_,fp))
5278        return 1;
5279    if (inDoubleArray(dual_,numberRows_,fp))
5280        return 1;
5281    if (inDoubleArray(reducedCost_,numberColumns_,fp))
5282        return 1;
5283    if (inDoubleArray(rowLower_,numberRows_,fp))
5284        return 1;
5285    if (inDoubleArray(rowUpper_,numberRows_,fp))
5286        return 1;
5287    double * objective;
5288    if (inDoubleArray(objective,numberColumns_,fp))
5289        return 1;
5290    delete objective_;
5291    objective_ = new ClpLinearObjective(objective,numberColumns_);
5292    delete [] objective;
5293    if (inDoubleArray(rowObjective_,numberRows_,fp))
5294        return 1;
5295    if (inDoubleArray(columnLower_,numberColumns_,fp))
5296        return 1;
5297    if (inDoubleArray(columnUpper_,numberColumns_,fp))
5298        return 1;
5299    if (problemStatus_==1) {
5300      if (inDoubleArray(ray_,numberRows_,fp))
5301        return 1;
5302    } else if (problemStatus_==2) {
5303      if (inDoubleArray(ray_,numberColumns_,fp))
5304        return 1;
5305    } else {
5306      // ray should be null
5307      numberRead = fread(&length,sizeof(int),1,fp);
5308      if (numberRead!=1)
5309        return 1;
5310      if (length)
5311        return 2;
5312    }
5313    delete [] status_;
5314    status_=NULL;
5315    // status region
5316    numberRead = fread(&length,sizeof(int),1,fp);
5317    if (numberRead!=1)
5318        return 1;
5319    if (length) {
5320      if (length!=numberRows_+numberColumns_)
5321        return 1;
5322      status_ = new char unsigned[length];
5323      numberRead = fread(status_,sizeof(char),length, fp);
5324      if (numberRead!=length)
5325        return 1;
5326    }
5327#ifndef CLP_NO_STD
5328    if (lengthNames_) {
5329      char * array = 
5330        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
5331      char * get = array;
5332      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
5333      if (numberRead!=numberRows_)
5334        return 1;
5335      rowNames_ = std::vector<std::string> ();
5336      rowNames_.resize(numberRows_);
5337      for (i=0;i<numberRows_;i++) {
5338        rowNames_.push_back(get);
5339        get += lengthNames_+1;
5340      }
5341      get = array;
5342      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
5343      if (numberRead!=numberColumns_)
5344        return 1;
5345      columnNames_ = std::vector<std::string> ();
5346      columnNames_.resize(numberColumns_);
5347      for (i=0;i<numberColumns_;i++) {
5348        columnNames_.push_back(get);
5349        get += lengthNames_+1;
5350      }
5351      delete [] array;
5352    }
5353#endif
5354    // Pivot choices
5355    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
5356    delete dualRowPivot_;
5357    switch ((scalars.dualPivotChoice&63)) {
5358    default:
5359      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
5360    case 1:
5361      // Dantzig
5362      dualRowPivot_ = new ClpDualRowDantzig();
5363      break;
5364    case 2:
5365      // Steepest - use mode
5366      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
5367      break;
5368    }
5369    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
5370    delete primalColumnPivot_;
5371    switch ((scalars.primalPivotChoice&63)) {
5372    default:
5373      printf("Need another primalPivot case %d\n",
5374             scalars.primalPivotChoice&63);
5375    case 1:
5376      // Dantzig
5377      primalColumnPivot_ = new ClpPrimalColumnDantzig();
5378      break;
5379    case 2:
5380      // Steepest - use mode
5381      primalColumnPivot_
5382        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
5383      break;
5384    }
5385    assert(scalars.matrixStorageChoice==1);
5386    delete matrix_;
5387    // get arrays
5388    numberRead = fread(&length,sizeof(int),1,fp);
5389    if (numberRead!=1)
5390      return 1;
5391    double * elements = new double[length];
5392    int * indices = new int[length];
5393    CoinBigIndex * starts = new CoinBigIndex[numberColumns_+1];
5394    int * lengths = new int[numberColumns_];
5395    numberRead = fread(elements, sizeof(double),length,fp);
5396    if (numberRead!=length)
5397      return 1;
5398    numberRead = fread(indices, sizeof(int),length,fp);
5399    if (numberRead!=length)
5400      return 1;
5401    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
5402    if (numberRead!=numberColumns_+1)
5403      return 1;
5404    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
5405    if (numberRead!=numberColumns_)
5406      return 1;
5407    // assign matrix
5408    CoinPackedMatrix * matrix = new CoinPackedMatrix();
5409    matrix->setExtraGap(0.0);
5410    matrix->setExtraMajor(0.0);
5411    // Pack down
5412    length=0;
5413    for (i=0;i<numberColumns_;i++) {
5414      int start = starts[i];
5415      starts[i]=length;
5416      for (CoinBigIndex j=start;j<start+lengths[i];j++) {
5417        elements[length]=elements[j];
5418        indices[length++]=indices[j];
5419      }
5420      lengths[i]=length-starts[i];
5421    }
5422    starts[numberColumns_]=length;
5423    matrix->assignMatrix(true, numberRows_, numberColumns_,
5424                         length, elements, indices, starts, lengths);
5425    // and transfer to Clp
5426    matrix_ = new ClpPackedMatrix(matrix);
5427    // finished
5428    fclose(fp);
5429    return 0;
5430  } else {
5431    return -1;
5432  }
5433  return 0;
5434}
5435#endif
5436// value of incoming variable (in Dual)
5437double 
5438ClpSimplex::valueIncomingDual() const
5439{
5440  // Need value of incoming for list of infeasibilities as may be infeasible
5441  double valueIncoming = (dualOut_/alpha_)*directionOut_;
5442  if (directionIn_==-1)
5443    valueIncoming = upperIn_-valueIncoming;
5444  else
5445    valueIncoming = lowerIn_-valueIncoming;
5446  return valueIncoming;
5447}
5448// Sanity check on input data - returns true if okay
5449bool 
5450ClpSimplex::sanityCheck()
5451{
5452  // bad if empty
5453  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
5454    int infeasNumber[2];
5455    double infeasSum[2];
5456    problemStatus_=emptyProblem(infeasNumber,infeasSum,false);
5457    numberDualInfeasibilities_=infeasNumber[0];
5458    sumDualInfeasibilities_=infeasSum[0];
5459    numberPrimalInfeasibilities_=infeasNumber[1];
5460    sumPrimalInfeasibilities_=infeasSum[1];
5461    return false;
5462  }
5463  int numberBad ;
5464  double largestBound, smallestBound, minimumGap;
5465  double smallestObj, largestObj;
5466  int firstBad;
5467  int modifiedBounds=0;
5468  int i;
5469  numberBad=0;
5470  firstBad=-1;
5471  minimumGap=1.0e100;
5472  smallestBound=1.0e100;
5473  largestBound=0.0;
5474  smallestObj=1.0e100;
5475  largestObj=0.0;
5476  // If bounds are too close - fix
5477  double fixTolerance = primalTolerance_;
5478  if (fixTolerance<2.0e-8)
5479    fixTolerance *= 1.1;
5480  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
5481    double value;
5482    value = fabs(cost_[i]);
5483    if (value>1.0e50) {
5484      numberBad++;
5485      if (firstBad<0)
5486        firstBad=i;
5487    } else if (value) {
5488      if (value>largestObj)
5489        largestObj=value;
5490      if (value<smallestObj)
5491        smallestObj=value;
5492    }
5493    value=upper_[i]-lower_[i];
5494    if (value<-primalTolerance_) {
5495      numberBad++;
5496      if (firstBad<0)
5497        firstBad=i;
5498    } else if (value<=fixTolerance) {
5499      if (value) {
5500        // modify
5501        upper_[i] = lower_[i];
5502        modifiedBounds++;
5503      }
5504    } else {
5505      if (value<minimumGap)
5506        minimumGap=value;
5507    }
5508    if (lower_[i]>-1.0e100&&lower_[i]) {
5509      value = fabs(lower_[i]);
5510      if (value>largestBound)
5511        largestBound=value;
5512      if (value<smallestBound)
5513        smallestBound=value;
5514    }
5515    if (upper_[i]<1.0e100&&upper_[i]) {
5516      value = fabs(upper_[i]);
5517      if (value>largestBound)
5518        largestBound=value;
5519      if (value<smallestBound)
5520        smallestBound=value;
5521    }
5522  }
5523  if (largestBound)
5524    handler_->message(CLP_RIMSTATISTICS3,messages_)
5525      <<smallestBound
5526      <<largestBound
5527      <<minimumGap
5528      <<CoinMessageEol;
5529  minimumGap=1.0e100;
5530  smallestBound=1.0e100;
5531  largestBound=0.0;
5532  for (i