source: stable/1.7/Clp/src/ClpSimplex.cpp @ 1213

Last change on this file since 1213 was 1213, checked in by forrest, 13 years ago

for Lou Hafer to update Osi stable

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