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

Last change on this file since 798 was 798, checked in by ladanyi, 14 years ago

finished switch to subversion

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