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

Last change on this file since 1212 was 1212, checked in by forrest, 12 years ago

for Francisco and change internal solution check

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