source: trunk/ClpSimplex.cpp @ 752

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

check bounds

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