source: trunk/ClpSimplex.cpp @ 748

Last change on this file since 748 was 748, checked in by forrest, 14 years ago

for writeLp mainly

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