source: trunk/Clp/src/ClpSimplex.cpp @ 799

Last change on this file since 799 was 799, checked in by andreasw, 14 years ago

undid last commit (patches incorrectly applied)

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