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

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

correct solution investigation

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 350.2 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  //#define COMPUTE_INT_INFEAS
1595#ifdef COMPUTE_INT_INFEAS
1596  if (algorithm_>0&&integerType_&&!nonLinearCost_->numberInfeasibilities()) {
1597    if (fabs(theta_)>1.0e-6||!numberIterations_) {
1598      int numberFixed=0;
1599      int numberUnsat=0;
1600      int numberSat=0;
1601      double sumUnsat=0.0;
1602      double tolerance = 10.0*primalTolerance_;
1603      double mostAway=0.0;
1604      int iAway=-1;
1605      for (int i=0;i<numberColumns_;i++) {
1606        if (upper_[i]>lower_[i]) {
1607          double value = solution_[i];
1608          if (value>lower_[i]+tolerance&&
1609              value<upper_[i]-tolerance) {
1610            if (columnScale_)
1611              value *= columnScale_[i];
1612            double closest = floor(value+0.5);
1613            // problem may be perturbed so relax test
1614            if (fabs(value-closest)>1.0e-4) {
1615              numberUnsat++;
1616              sumUnsat += fabs(value-closest);
1617              if (mostAway<fabs(value-closest)) {
1618                mostAway=fabs(value-closest);
1619                iAway=i;
1620              }
1621            } else {
1622              numberSat++;
1623            }
1624          } else {
1625            numberSat++;
1626          }
1627        } else {
1628          numberFixed++;
1629        }
1630      }
1631      printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
1632             numberIterations_,numberUnsat,sumUnsat,mostAway,numberFixed,numberSat);
1633    }
1634  }
1635#endif
1636  // save value of incoming and outgoing
1637  double oldIn = solution_[sequenceIn_];
1638  double oldOut = solution_[sequenceOut_];
1639  numberIterations_++;
1640  changeMade_++; // something has happened
1641  // incoming variable
1642  if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
1643    handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
1644      <<directionOut_
1645      <<directionIn_<<theta_
1646      <<dualOut_<<dualIn_<<alpha_
1647      <<CoinMessageEol;
1648    if (getStatus(sequenceIn_)==isFree) {
1649      handler_->message(CLP_SIMPLEX_FREEIN,messages_)
1650        <<sequenceIn_
1651        <<CoinMessageEol;
1652    }
1653  }
1654  // change of incoming
1655  char rowcol[]={'R','C'};
1656  if (pivotRow_>=0)
1657    pivotVariable_[pivotRow_]=sequenceIn();
1658  if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20)
1659    progressFlag_ |= 2; // making real progress
1660  solution_[sequenceIn_]=valueIn_;
1661  if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
1662    progressFlag_ |= 1; // making real progress
1663  if (sequenceIn_!=sequenceOut_) {
1664    if (alphaAccuracy_>0.0) {
1665      double value = fabs(alpha_);
1666      if (value>1.0)
1667        alphaAccuracy_ *= value;
1668      else
1669        alphaAccuracy_ /= value;
1670    }
1671    //assert( getStatus(sequenceOut_)== basic);
1672    setStatus(sequenceIn_,basic);
1673    if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
1674      // As Nonlinear costs may have moved bounds (to more feasible)
1675      // Redo using value
1676      if (fabs(valueOut_-lower_[sequenceOut_])<fabs(valueOut_-upper_[sequenceOut_])) {
1677        // going to lower
1678        setStatus(sequenceOut_,atLowerBound);
1679        oldOut = lower_[sequenceOut_];
1680      } else {
1681        // going to upper
1682        setStatus(sequenceOut_,atUpperBound);
1683        oldOut = upper_[sequenceOut_];
1684      }
1685    } else {
1686      // fixed
1687      setStatus(sequenceOut_,isFixed);
1688    }
1689    solution_[sequenceOut_]=valueOut_;
1690  } else {
1691    //if (objective_->type()<2)
1692    //assert (fabs(theta_)>1.0e-13);
1693    // flip from bound to bound
1694    // As Nonlinear costs may have moved bounds (to more feasible)
1695    // Redo using value
1696    if (fabs(valueIn_-lower_[sequenceIn_])<fabs(valueIn_-upper_[sequenceIn_])) {
1697      // as if from upper bound
1698      setStatus(sequenceIn_, atLowerBound);
1699    } else {
1700      // as if from lower bound
1701      setStatus(sequenceIn_, atUpperBound);
1702    }
1703  }
1704 
1705  // Update hidden stuff e.g. effective RHS and gub
1706  matrix_->updatePivot(this,oldIn,oldOut);
1707  objectiveValue_ += objectiveChange/(objectiveScale_*rhsScale_);
1708  if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
1709    handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
1710      <<numberIterations_<<objectiveValue()
1711      <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
1712      <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
1713    handler_->printing(algorithm_<0)<<dualOut_<<theta_;
1714    handler_->printing(algorithm_>0)<<dualIn_<<theta_;
1715    handler_->message()<<CoinMessageEol;
1716  }
1717  if (hitMaximumIterations())
1718    return 2;
1719#if 1
1720  //if (numberIterations_>14000)
1721  //handler_->setLogLevel(63);
1722  //if (numberIterations_>24000)
1723  //exit(77);
1724  // check for small cycles
1725  int in = sequenceIn_;
1726  int out = sequenceOut_;
1727  matrix_->correctSequence(this,in,out);
1728  int cycle=progress_.cycle(in,out,
1729                            directionIn_,directionOut_);
1730  if (cycle>0&&objective_->type()<2) {
1731    //if (cycle>0) {
1732    if (handler_->logLevel()>=63)
1733      printf("Cycle of %d\n",cycle);
1734    // reset
1735    progress_.startCheck();
1736    double random = randomNumberGenerator_.randomDouble();
1737    int extra = (int) (9.999*random);
1738    int off[]={1,1,1,1,2,2,2,3,3,4};
1739    if (factorization_->pivots()>cycle) {
1740      forceFactorization_=CoinMax(1,cycle-off[extra]);
1741    } else {
1742      // need to reject something
1743      int iSequence;
1744      if (algorithm_>0)
1745        iSequence = sequenceIn_;
1746      else
1747        iSequence = sequenceOut_;
1748      char x = isColumn(iSequence) ? 'C' :'R';
1749      if (handler_->logLevel()>=63)
1750        handler_->message(CLP_SIMPLEX_FLAG,messages_)
1751          <<x<<sequenceWithin(iSequence)
1752          <<CoinMessageEol;
1753      setFlagged(iSequence);
1754      //printf("flagging %d\n",iSequence);
1755    }
1756    return 1;
1757  }
1758#endif
1759  // only time to re-factorize if one before real time
1760  // this is so user won't be surprised that maximumPivots has exact meaning
1761  int numberPivots=factorization_->pivots();
1762  int maximumPivots = factorization_->maximumPivots();
1763  int numberDense = factorization_->numberDense();
1764  if (numberPivots==maximumPivots||
1765      maximumPivots<2) {
1766    // If dense then increase
1767    if (maximumPivots>100&&numberDense>1.5*maximumPivots) {
1768      factorization_->maximumPivots(numberDense);
1769      dualRowPivot_->maximumPivotsChanged();
1770      primalColumnPivot_->maximumPivotsChanged();
1771      // and redo arrays
1772      for (int iRow=0;iRow<4;iRow++) {
1773        int length =rowArray_[iRow]->capacity()+numberDense-maximumPivots;
1774        rowArray_[iRow]->reserve(length);
1775      }
1776    }
1777    return 1;
1778  } else if (factorization_->timeToRefactorize()) {
1779    //printf("ret after %d pivots\n",factorization_->pivots());
1780    return 1;
1781  } else if (forceFactorization_>0&&
1782             factorization_->pivots()==forceFactorization_) {
1783    // relax
1784    forceFactorization_ = (3+5*forceFactorization_)/4;
1785    if (forceFactorization_>factorization_->maximumPivots())
1786      forceFactorization_ = -1; //off
1787    return 1;
1788  } else if (numberIterations_>1000+10*(numberRows_+(numberColumns_>>2))) {
1789    double random = randomNumberGenerator_.randomDouble();
1790    int maxNumber = (forceFactorization_<0) ? maximumPivots : CoinMin(forceFactorization_,maximumPivots);
1791    if (factorization_->pivots()>=random*maxNumber) {
1792      return 1;
1793    } else if (numberIterations_>1000000+10*(numberRows_+(numberColumns_>>2))&&
1794               numberIterations_<1001000+10*(numberRows_+(numberColumns_>>2))) {
1795      return 1;
1796    } else {
1797      // carry on iterating
1798      return 0;
1799    }
1800  } else {
1801    // carry on iterating
1802    return 0;
1803  }
1804}
1805// Copy constructor.
1806ClpSimplex::ClpSimplex(const ClpSimplex &rhs,int scalingMode) :
1807  ClpModel(rhs,scalingMode),
1808  columnPrimalInfeasibility_(0.0),
1809  rowPrimalInfeasibility_(0.0),
1810  columnPrimalSequence_(-2),
1811  rowPrimalSequence_(-2), 
1812  columnDualInfeasibility_(0.0),
1813  rowDualInfeasibility_(0.0),
1814  moreSpecialOptions_(2),
1815  baseIteration_(0),
1816  primalToleranceToGetOptimal_(-1.0),
1817  remainingDualInfeasibility_(0.0),
1818  largeValue_(1.0e15),
1819  largestPrimalError_(0.0),
1820  largestDualError_(0.0),
1821  alphaAccuracy_(-1.0),
1822  dualBound_(1.0e10),
1823  alpha_(0.0),
1824  theta_(0.0),
1825  lowerIn_(0.0),
1826  valueIn_(0.0),
1827  upperIn_(-COIN_DBL_MAX),
1828  dualIn_(0.0),
1829  lowerOut_(-1),
1830  valueOut_(-1),
1831  upperOut_(-1),
1832  dualOut_(-1),
1833  dualTolerance_(0.0),
1834  primalTolerance_(0.0),
1835  sumDualInfeasibilities_(0.0),
1836  sumPrimalInfeasibilities_(0.0),
1837  infeasibilityCost_(1.0e10),
1838  sumOfRelaxedDualInfeasibilities_(0.0),
1839  sumOfRelaxedPrimalInfeasibilities_(0.0),
1840  acceptablePivot_(1.0e-8),
1841  lower_(NULL),
1842  rowLowerWork_(NULL),
1843  columnLowerWork_(NULL),
1844  upper_(NULL),
1845  rowUpperWork_(NULL),
1846  columnUpperWork_(NULL),
1847  cost_(NULL),
1848  rowObjectiveWork_(NULL),
1849  objectiveWork_(NULL),
1850  sequenceIn_(-1),
1851  directionIn_(-1),
1852  sequenceOut_(-1),
1853  directionOut_(-1),
1854  pivotRow_(-1),
1855  lastGoodIteration_(-100),
1856  dj_(NULL),
1857  rowReducedCost_(NULL),
1858  reducedCostWork_(NULL),
1859  solution_(NULL),
1860  rowActivityWork_(NULL),
1861  columnActivityWork_(NULL),
1862  auxiliaryModel_(NULL),
1863  numberDualInfeasibilities_(0),
1864  numberDualInfeasibilitiesWithoutFree_(0),
1865  numberPrimalInfeasibilities_(100),
1866  numberRefinements_(0),
1867  pivotVariable_(NULL),
1868  factorization_(NULL),
1869  savedSolution_(NULL),
1870  numberTimesOptimal_(0),
1871  disasterArea_(NULL),
1872  changeMade_(1),
1873  algorithm_(0),
1874  forceFactorization_(-1),
1875  perturbation_(100),
1876  nonLinearCost_(NULL),
1877  lastBadIteration_(-999999),
1878  lastFlaggedIteration_(-999999),
1879  numberFake_(0),
1880  numberChanged_(0),
1881  progressFlag_(0),
1882  firstFree_(-1),
1883  numberExtraRows_(0),
1884  maximumBasic_(0),
1885  dontFactorizePivots_(0),
1886  incomingInfeasibility_(1.0),
1887  allowedInfeasibility_(10.0),
1888  automaticScale_(0),
1889  baseModel_(NULL)
1890{
1891  int i;
1892  for (i=0;i<6;i++) {
1893    rowArray_[i]=NULL;
1894    columnArray_[i]=NULL;
1895  }
1896  for (i=0;i<4;i++) {
1897    spareIntArray_[i]=0;
1898    spareDoubleArray_[i]=0.0;
1899  }
1900  saveStatus_=NULL;
1901  factorization_ = NULL;
1902  dualRowPivot_ = NULL;
1903  primalColumnPivot_ = NULL;
1904  gutsOfDelete(0);
1905  delete nonLinearCost_;
1906  nonLinearCost_ = NULL;
1907  gutsOfCopy(rhs);
1908  solveType_=1; // say simplex based life form
1909}
1910// Copy constructor from model
1911ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
1912  ClpModel(rhs,scalingMode),
1913  columnPrimalInfeasibility_(0.0),
1914  rowPrimalInfeasibility_(0.0),
1915  columnPrimalSequence_(-2),
1916  rowPrimalSequence_(-2), 
1917  columnDualInfeasibility_(0.0),
1918  rowDualInfeasibility_(0.0),
1919  moreSpecialOptions_(2),
1920  baseIteration_(0),
1921  primalToleranceToGetOptimal_(-1.0),
1922  remainingDualInfeasibility_(0.0),
1923  largeValue_(1.0e15),
1924  largestPrimalError_(0.0),
1925  largestDualError_(0.0),
1926  alphaAccuracy_(-1.0),
1927  dualBound_(1.0e10),
1928  alpha_(0.0),
1929  theta_(0.0),
1930  lowerIn_(0.0),
1931  valueIn_(0.0),
1932  upperIn_(-COIN_DBL_MAX),
1933  dualIn_(0.0),
1934  lowerOut_(-1),
1935  valueOut_(-1),
1936  upperOut_(-1),
1937  dualOut_(-1),
1938  dualTolerance_(0.0),
1939  primalTolerance_(0.0),
1940  sumDualInfeasibilities_(0.0),
1941  sumPrimalInfeasibilities_(0.0),
1942  infeasibilityCost_(1.0e10),
1943  sumOfRelaxedDualInfeasibilities_(0.0),
1944  sumOfRelaxedPrimalInfeasibilities_(0.0),
1945  acceptablePivot_(1.0e-8),
1946  lower_(NULL),
1947  rowLowerWork_(NULL),
1948  columnLowerWork_(NULL),
1949  upper_(NULL),
1950  rowUpperWork_(NULL),
1951  columnUpperWork_(NULL),
1952  cost_(NULL),
1953  rowObjectiveWork_(NULL),
1954  objectiveWork_(NULL),
1955  sequenceIn_(-1),
1956  directionIn_(-1),
1957  sequenceOut_(-1),
1958  directionOut_(-1),
1959  pivotRow_(-1),
1960  lastGoodIteration_(-100),
1961  dj_(NULL),
1962  rowReducedCost_(NULL),
1963  reducedCostWork_(NULL),
1964  solution_(NULL),
1965  rowActivityWork_(NULL),
1966  columnActivityWork_(NULL),
1967  auxiliaryModel_(NULL),
1968  numberDualInfeasibilities_(0),
1969  numberDualInfeasibilitiesWithoutFree_(0),
1970  numberPrimalInfeasibilities_(100),
1971  numberRefinements_(0),
1972  pivotVariable_(NULL),
1973  factorization_(NULL),
1974  savedSolution_(NULL),
1975  numberTimesOptimal_(0),
1976  disasterArea_(NULL),
1977  changeMade_(1),
1978  algorithm_(0),
1979  forceFactorization_(-1),
1980  perturbation_(100),
1981  nonLinearCost_(NULL),
1982  lastBadIteration_(-999999),
1983  lastFlaggedIteration_(-999999),
1984  numberFake_(0),
1985  numberChanged_(0),
1986  progressFlag_(0),
1987  firstFree_(-1),
1988  numberExtraRows_(0),
1989  maximumBasic_(0),
1990  dontFactorizePivots_(0),
1991  incomingInfeasibility_(1.0),
1992  allowedInfeasibility_(10.0),
1993  automaticScale_(0),
1994  baseModel_(NULL)
1995{
1996  int i;
1997  for (i=0;i<6;i++) {
1998    rowArray_[i]=NULL;
1999    columnArray_[i]=NULL;
2000  }
2001  for (i=0;i<4;i++) {
2002    spareIntArray_[i]=0;
2003    spareDoubleArray_[i]=0.0;
2004  }
2005  saveStatus_=NULL;
2006  // get an empty factorization so we can set tolerances etc
2007  getEmptyFactorization();
2008  // say Steepest pricing
2009  dualRowPivot_ = new ClpDualRowSteepest();
2010  // say Steepest pricing
2011  primalColumnPivot_ = new ClpPrimalColumnSteepest();
2012  solveType_=1; // say simplex based life form
2013 
2014}
2015// Assignment operator. This copies the data
2016ClpSimplex & 
2017ClpSimplex::operator=(const ClpSimplex & rhs)
2018{
2019  if (this != &rhs) {
2020    gutsOfDelete(0);
2021    delete nonLinearCost_;
2022    nonLinearCost_ = NULL;
2023    ClpModel::operator=(rhs);
2024    gutsOfCopy(rhs);
2025  }
2026  return *this;
2027}
2028void 
2029ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
2030{
2031  assert (numberRows_==rhs.numberRows_);
2032  assert (numberColumns_==rhs.numberColumns_);
2033  numberExtraRows_ = rhs.numberExtraRows_;
2034  maximumBasic_ = rhs.maximumBasic_;
2035  dontFactorizePivots_ = rhs.dontFactorizePivots_;
2036  int numberRows2 = numberRows_+numberExtraRows_;
2037  auxiliaryModel_ = NULL;
2038  moreSpecialOptions_ = rhs.moreSpecialOptions_;
2039  if ((whatsChanged_&1)!=0) {
2040    int numberTotal = numberColumns_+numberRows2;
2041    if ((specialOptions_&65536)!=0&&maximumRows_>=0) {
2042      assert (maximumInternalRows_>=numberRows2);
2043      assert (maximumInternalColumns_>=numberColumns_);
2044      numberTotal = 2*(maximumInternalColumns_+ maximumInternalRows_);
2045    }
2046    lower_ = ClpCopyOfArray(rhs.lower_,numberTotal);
2047    rowLowerWork_ = lower_+numberColumns_;
2048    columnLowerWork_ = lower_;
2049    upper_ = ClpCopyOfArray(rhs.upper_,numberTotal);
2050    rowUpperWork_ = upper_+numberColumns_;
2051    columnUpperWork_ = upper_;
2052    cost_ = ClpCopyOfArray(rhs.cost_,numberTotal);
2053    objectiveWork_ = cost_;
2054    rowObjectiveWork_ = cost_+numberColumns_;
2055    dj_ = ClpCopyOfArray(rhs.dj_,numberTotal);
2056    if (dj_) {
2057      reducedCostWork_ = dj_;
2058      rowReducedCost_ = dj_+numberColumns_;
2059    }
2060    solution_ = ClpCopyOfArray(rhs.solution_,numberTotal);
2061    if (solution_) {
2062      columnActivityWork_ = solution_;
2063      rowActivityWork_ = solution_+numberColumns_;
2064    }
2065    if (rhs.pivotVariable_) {
2066      pivotVariable_ = new int[numberRows2];
2067      CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
2068    } else {
2069      pivotVariable_=NULL;
2070    }
2071    savedSolution_ = ClpCopyOfArray(rhs.savedSolution_,numberTotal);
2072    int i;
2073    for (i=0;i<6;i++) {
2074      rowArray_[i]=NULL;
2075      if (rhs.rowArray_[i]) 
2076        rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
2077      columnArray_[i]=NULL;
2078      if (rhs.columnArray_[i]) 
2079        columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
2080    }
2081    if (rhs.saveStatus_) {
2082      saveStatus_ = ClpCopyOfArray( rhs.saveStatus_,numberTotal);
2083    }
2084  } else {
2085    lower_ = NULL;
2086    rowLowerWork_ = NULL;
2087    columnLowerWork_ = NULL; 
2088    upper_ = NULL; 
2089    rowUpperWork_ = NULL; 
2090    columnUpperWork_ = NULL;
2091    cost_ = NULL; 
2092    objectiveWork_ = NULL; 
2093    rowObjectiveWork_ = NULL; 
2094    dj_ = NULL; 
2095    reducedCostWork_ = NULL;
2096    rowReducedCost_ = NULL;
2097    solution_ = NULL; 
2098    columnActivityWork_ = NULL; 
2099    rowActivityWork_ = NULL; 
2100    pivotVariable_=NULL;
2101    savedSolution_ = NULL;
2102    int i;
2103    for (i=0;i<6;i++) {
2104      rowArray_[i]=NULL;
2105      columnArray_[i]=NULL;
2106    }
2107    saveStatus_ = NULL;
2108  }
2109  if (rhs.factorization_) {
2110    delete factorization_;
2111    factorization_ = new ClpFactorization(*rhs.factorization_);
2112  } else {
2113    factorization_=NULL;
2114  }
2115  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
2116  columnPrimalSequence_ = rhs.columnPrimalSequence_;
2117  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
2118  rowPrimalSequence_ = rhs.rowPrimalSequence_;
2119  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
2120  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
2121  baseIteration_ = rhs.baseIteration_;
2122  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
2123  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
2124  largeValue_ = rhs.largeValue_;
2125  largestPrimalError_ = rhs.largestPrimalError_;
2126  largestDualError_ = rhs.largestDualError_;
2127  alphaAccuracy_ = rhs.alphaAccuracy_;
2128  dualBound_ = rhs.dualBound_;
2129  alpha_ = rhs.alpha_;
2130  theta_ = rhs.theta_;
2131  lowerIn_ = rhs.lowerIn_;
2132  valueIn_ = rhs.valueIn_;
2133  upperIn_ = rhs.upperIn_;
2134  dualIn_ = rhs.dualIn_;
2135  sequenceIn_ = rhs.sequenceIn_;
2136  directionIn_ = rhs.directionIn_;
2137  lowerOut_ = rhs.lowerOut_;
2138  valueOut_ = rhs.valueOut_;
2139  upperOut_ = rhs.upperOut_;
2140  dualOut_ = rhs.dualOut_;
2141  sequenceOut_ = rhs.sequenceOut_;
2142  directionOut_ = rhs.directionOut_;
2143  pivotRow_ = rhs.pivotRow_;
2144  lastGoodIteration_ = rhs.lastGoodIteration_;
2145  numberRefinements_ = rhs.numberRefinements_;
2146  dualTolerance_ = rhs.dualTolerance_;
2147  primalTolerance_ = rhs.primalTolerance_;
2148  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
2149  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
2150  numberDualInfeasibilitiesWithoutFree_ = 
2151    rhs.numberDualInfeasibilitiesWithoutFree_;
2152  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
2153  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
2154  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
2155  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
2156  numberTimesOptimal_ = rhs.numberTimesOptimal_;
2157  disasterArea_ = NULL;
2158  changeMade_ = rhs.changeMade_;
2159  algorithm_ = rhs.algorithm_;
2160  forceFactorization_ = rhs.forceFactorization_;
2161  perturbation_ = rhs.perturbation_;
2162  infeasibilityCost_ = rhs.infeasibilityCost_;
2163  lastBadIteration_ = rhs.lastBadIteration_;
2164  lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
2165  numberFake_ = rhs.numberFake_;
2166  numberChanged_ = rhs.numberChanged_;
2167  progressFlag_ = rhs.progressFlag_;
2168  firstFree_ = rhs.firstFree_;
2169  incomingInfeasibility_ = rhs.incomingInfeasibility_;
2170  allowedInfeasibility_ = rhs.allowedInfeasibility_;
2171  automaticScale_ = rhs.automaticScale_;
2172  if (rhs.baseModel_) {
2173    baseModel_ = new ClpSimplex(*rhs.baseModel_);
2174  } else {
2175    baseModel_ = NULL;
2176  }
2177  progress_=rhs.progress_;
2178  for (int i=0;i<4;i++) {
2179    spareIntArray_[i]=rhs.spareIntArray_[i];
2180    spareDoubleArray_[i]=rhs.spareDoubleArray_[i];
2181  }
2182  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
2183  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
2184  acceptablePivot_ = rhs.acceptablePivot_;
2185  if (rhs.nonLinearCost_!=NULL)
2186    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
2187  else
2188    nonLinearCost_=NULL;
2189  solveType_=rhs.solveType_;
2190}
2191// type == 0 do everything, most + pivot data, 2 factorization data as well
2192void 
2193ClpSimplex::gutsOfDelete(int type)
2194{
2195  if (!type||(specialOptions_&65536)==0) {
2196    maximumInternalColumns_ = -1;
2197    maximumInternalRows_ = -1;
2198    delete [] lower_;
2199    lower_=NULL;
2200    rowLowerWork_=NULL;
2201    columnLowerWork_=NULL;
2202    delete [] upper_;
2203    upper_=NULL;
2204    rowUpperWork_=NULL;
2205    columnUpperWork_=NULL;
2206    delete [] cost_;
2207    cost_=NULL;
2208    objectiveWork_=NULL;
2209    rowObjectiveWork_=NULL;
2210    delete [] dj_;
2211    dj_=NULL;
2212    reducedCostWork_=NULL;
2213    rowReducedCost_=NULL;
2214    delete [] solution_;
2215    solution_=NULL;
2216    rowActivityWork_=NULL;
2217    columnActivityWork_=NULL;
2218    delete [] savedSolution_;
2219    savedSolution_ = NULL;
2220  }
2221  if ((specialOptions_&2)==0) {
2222    delete nonLinearCost_;
2223    nonLinearCost_ = NULL;
2224  }
2225  int i;
2226  if ((specialOptions_&65536)==0) {
2227    for (i=0;i<6;i++) {
2228      delete rowArray_[i];
2229      rowArray_[i]=NULL;
2230      delete columnArray_[i];
2231      columnArray_[i]=NULL;
2232    }
2233  }
2234  delete rowCopy_;
2235  rowCopy_=NULL;
2236  delete [] saveStatus_;
2237  saveStatus_=NULL;
2238  if (!type) {
2239    // delete everything
2240    delete auxiliaryModel_;
2241    auxiliaryModel_ = NULL;
2242    setEmptyFactorization();
2243    delete [] pivotVariable_;
2244    pivotVariable_=NULL;
2245    delete dualRowPivot_;
2246    dualRowPivot_ = NULL;
2247    delete primalColumnPivot_;
2248    primalColumnPivot_ = NULL;
2249    delete baseModel_;
2250    baseModel_=NULL;
2251  } else {
2252    // delete any size information in methods
2253    if (type>1) {
2254      factorization_->clearArrays();
2255      delete [] pivotVariable_;
2256      pivotVariable_=NULL;
2257    }
2258    dualRowPivot_->clearArrays();
2259    primalColumnPivot_->clearArrays();
2260  }
2261}
2262// This sets largest infeasibility and most infeasible
2263void 
2264ClpSimplex::checkPrimalSolution(const double * rowActivities,
2265                                        const double * columnActivities)
2266{
2267  double * solution;
2268  int iRow,iColumn;
2269
2270  objectiveValue_ = 0.0;
2271  // now look at primal solution
2272  solution = rowActivityWork_;
2273  sumPrimalInfeasibilities_=0.0;
2274  numberPrimalInfeasibilities_=0;
2275  double primalTolerance = primalTolerance_;
2276  double relaxedTolerance=primalTolerance_;
2277  // we can't really trust infeasibilities if there is primal error
2278  double error = CoinMin(1.0e-2,largestPrimalError_);
2279  // allow tolerance at least slightly bigger than standard
2280  relaxedTolerance = relaxedTolerance +  error;
2281  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2282  for (iRow=0;iRow<numberRows_;iRow++) {
2283    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
2284    double infeasibility=0.0;
2285    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
2286    if (solution[iRow]>rowUpperWork_[iRow]) {
2287      infeasibility=solution[iRow]-rowUpperWork_[iRow];
2288    } else if (solution[iRow]<rowLowerWork_[iRow]) {
2289      infeasibility=rowLowerWork_[iRow]-solution[iRow];
2290    }
2291    if (infeasibility>primalTolerance) {
2292      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2293      if (infeasibility>relaxedTolerance) 
2294        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2295      numberPrimalInfeasibilities_ ++;
2296    }
2297    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
2298  }
2299  // Check any infeasibilities from dynamic rows
2300  matrix_->primalExpanded(this,2);
2301  solution = columnActivityWork_;
2302  if (!matrix_->rhsOffset(this)) {
2303    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2304      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2305      double infeasibility=0.0;
2306      objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
2307      if (solution[iColumn]>columnUpperWork_[iColumn]) {
2308        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
2309      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
2310        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
2311      }
2312      if (infeasibility>primalTolerance) {
2313        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2314        if (infeasibility>relaxedTolerance) 
2315          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2316        numberPrimalInfeasibilities_ ++;
2317      }
2318      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
2319    }
2320  } else {
2321    // as we are using effective rhs we only check basics
2322    // But we do need to get objective
2323    objectiveValue_ += innerProduct(objectiveWork_,numberColumns_,solution);
2324    for (int j=0;j<numberRows_;j++) {
2325      int iColumn = pivotVariable_[j];
2326      //assert (fabs(solution[iColumn])<1.0e15||getColumnStatus(iColumn) == basic);
2327      double infeasibility=0.0;
2328      if (solution[iColumn]>columnUpperWork_[iColumn]) {
2329        infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
2330      } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
2331        infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
2332      }
2333      if (infeasibility>primalTolerance) {
2334        sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2335        if (infeasibility>relaxedTolerance) 
2336          sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
2337        numberPrimalInfeasibilities_ ++;
2338      }
2339      infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
2340    }
2341  }
2342  objectiveValue_ += objective_->nonlinearOffset();
2343  objectiveValue_ /= (objectiveScale_*rhsScale_);
2344}
2345void 
2346ClpSimplex::checkDualSolution()
2347{
2348
2349  int iRow,iColumn;
2350  sumDualInfeasibilities_=0.0;
2351  numberDualInfeasibilities_=0;
2352  numberDualInfeasibilitiesWithoutFree_=0;
2353  if (matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2) {
2354    // pretend we found dual infeasibilities
2355    sumOfRelaxedDualInfeasibilities_ = 1.0;
2356    sumDualInfeasibilities_=1.0;
2357    numberDualInfeasibilities_=1;
2358    return;
2359  }
2360  int firstFreePrimal = -1;
2361  int firstFreeDual = -1;
2362  int numberSuperBasicWithDj=0;
2363  double relaxedTolerance=dualTolerance_;
2364  // we can't really trust infeasibilities if there is dual error
2365  double error = CoinMin(1.0e-2,largestDualError_);
2366  // allow tolerance at least slightly bigger than standard
2367  relaxedTolerance = relaxedTolerance +  error;
2368  sumOfRelaxedDualInfeasibilities_ = 0.0;
2369
2370  // Check any djs from dynamic rows
2371  matrix_->dualExpanded(this,NULL,NULL,3);
2372  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
2373  objectiveValue_ = 0.0;
2374  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2375    objectiveValue_ += objectiveWork_[iColumn]*columnActivityWork_[iColumn];
2376    if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
2377      // not basic
2378      double distanceUp = columnUpperWork_[iColumn]-
2379        columnActivityWork_[iColumn];
2380      double distanceDown = columnActivityWork_[iColumn] -
2381        columnLowerWork_[iColumn];
2382      if (distanceUp>primalTolerance_) {
2383        double value = reducedCostWork_[iColumn];
2384        // Check if "free"
2385        if (distanceDown>primalTolerance_) {
2386          if (fabs(value)>1.0e2*relaxedTolerance) {
2387            numberSuperBasicWithDj++;
2388            if (firstFreeDual<0)
2389              firstFreeDual = iColumn;
2390          }
2391          if (firstFreePrimal<0)
2392            firstFreePrimal = iColumn;
2393        }
2394        // should not be negative
2395        if (value<0.0) {
2396          value = - value;
2397          if (value>dualTolerance_) {
2398            if (getColumnStatus(iColumn) != isFree) {
2399              numberDualInfeasibilitiesWithoutFree_ ++;
2400              sumDualInfeasibilities_ += value-dualTolerance_;
2401              if (value>relaxedTolerance) 
2402                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2403              numberDualInfeasibilities_ ++;
2404            } else {
2405              // free so relax a lot
2406              value *= 0.01;
2407              if (value>dualTolerance_) {
2408                sumDualInfeasibilities_ += value-dualTolerance_;
2409                if (value>relaxedTolerance) 
2410                  sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2411                numberDualInfeasibilities_ ++;
2412              }
2413            }
2414          }
2415        }
2416      }
2417      if (distanceDown>primalTolerance_) {
2418        double value = reducedCostWork_[iColumn];
2419        // should not be positive
2420        if (value>0.0) {
2421          if (value>dualTolerance_) {
2422            sumDualInfeasibilities_ += value-dualTolerance_;
2423            if (value>relaxedTolerance) 
2424              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2425            numberDualInfeasibilities_ ++;
2426            if (getColumnStatus(iColumn) != isFree) 
2427              numberDualInfeasibilitiesWithoutFree_ ++;
2428            // maybe we can make feasible by increasing tolerance
2429          }
2430        }
2431      }
2432    }
2433  }
2434  for (iRow=0;iRow<numberRows_;iRow++) {
2435    objectiveValue_ += rowActivityWork_[iRow]*rowObjectiveWork_[iRow];
2436    if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
2437      // not basic
2438      double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
2439      double distanceDown = rowActivityWork_[iRow] -rowLowerWork_[iRow];
2440      if (distanceUp>primalTolerance_) {
2441        double value = rowReducedCost_[iRow];
2442        // Check if "free"
2443        if (distanceDown>primalTolerance_) {
2444          if (fabs(value)>1.0e2*relaxedTolerance) {
2445            numberSuperBasicWithDj++;
2446            if (firstFreeDual<0)
2447              firstFreeDual = iRow+numberColumns_;
2448          }
2449          if (firstFreePrimal<0)
2450            firstFreePrimal = iRow+numberColumns_;
2451        }
2452        // should not be negative
2453        if (value<0.0) {
2454          value = - value;
2455          if (value>dualTolerance_) {
2456            sumDualInfeasibilities_ += value-dualTolerance_;
2457            if (value>relaxedTolerance) 
2458              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2459            numberDualInfeasibilities_ ++;
2460            if (getRowStatus(iRow) != isFree) 
2461              numberDualInfeasibilitiesWithoutFree_ ++;
2462          }
2463        }
2464      }
2465      if (distanceDown>primalTolerance_) {
2466        double value = rowReducedCost_[iRow];
2467        // should not be positive
2468        if (value>0.0) {
2469          if (value>dualTolerance_) {
2470            sumDualInfeasibilities_ += value-dualTolerance_;
2471            if (value>relaxedTolerance) 
2472              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
2473            numberDualInfeasibilities_ ++;
2474            if (getRowStatus(iRow) != isFree) 
2475              numberDualInfeasibilitiesWithoutFree_ ++;
2476            // maybe we can make feasible by increasing tolerance
2477          }
2478        }
2479      }
2480    }
2481  }
2482  if (algorithm_<0&&firstFreeDual>=0) {
2483    // dual
2484    firstFree_ = firstFreeDual;
2485  } else if (numberSuperBasicWithDj||
2486             (progress_.lastIterationNumber(0)<=0)) {
2487    firstFree_=firstFreePrimal;
2488  }
2489  objectiveValue_ += objective_->nonlinearOffset();
2490  objectiveValue_ /= (objectiveScale_*rhsScale_);
2491}
2492/* This sets sum and number of infeasibilities (Dual and Primal) */
2493void 
2494ClpSimplex::checkBothSolutions()
2495{
2496  if ((matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2)||
2497      matrix_->rhsOffset(this)) {
2498    // old way
2499    checkPrimalSolution(rowActivityWork_,columnActivityWork_);
2500    checkDualSolution();
2501    return;
2502  }
2503  int iSequence;
2504
2505  objectiveValue_ = 0.0;
2506  sumPrimalInfeasibilities_=0.0;
2507  numberPrimalInfeasibilities_=0;
2508  double primalTolerance = primalTolerance_;
2509  double relaxedToleranceP=primalTolerance_;
2510  // we can't really trust infeasibilities if there is primal error
2511  double error = CoinMin(1.0e-2,largestPrimalError_);
2512  // allow tolerance at least slightly bigger than standard
2513  relaxedToleranceP = relaxedToleranceP +  error;
2514  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2515  sumDualInfeasibilities_=0.0;
2516  numberDualInfeasibilities_=0;
2517  double dualTolerance=dualTolerance_;
2518  double relaxedToleranceD=dualTolerance;
2519  // we can't really trust infeasibilities if there is dual error
2520  error = CoinMin(1.0e-2,largestDualError_);
2521  // allow tolerance at least slightly bigger than standard
2522  relaxedToleranceD = relaxedToleranceD +  error;
2523  sumOfRelaxedDualInfeasibilities_ = 0.0;
2524
2525  // Check any infeasibilities from dynamic rows
2526  matrix_->primalExpanded(this,2);
2527  // Check any djs from dynamic rows
2528  matrix_->dualExpanded(this,NULL,NULL,3);
2529  int numberDualInfeasibilitiesFree= 0;
2530  int firstFreePrimal = -1;
2531  int firstFreeDual = -1;
2532  int numberSuperBasicWithDj=0;
2533
2534  int numberTotal = numberRows_ + numberColumns_;
2535  for (iSequence=0;iSequence<numberTotal;iSequence++) {
2536    double value = solution_[iSequence];
2537#ifdef COIN_DEBUG
2538    if (fabs(value)>1.0e20)
2539      printf("%d values %g %g %g - status %d\n",iSequence,lower_[iSequence],
2540             solution_[iSequence],upper_[iSequence],status_[iSequence]);
2541#endif
2542    objectiveValue_ += value*cost_[iSequence];
2543    double distanceUp =upper_[iSequence]-value;
2544    double distanceDown =value-lower_[iSequence];
2545    if (distanceUp<-primalTolerance) {
2546      double infeasibility=-distanceUp;
2547      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2548      if (infeasibility>relaxedToleranceP) 
2549        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedToleranceP;
2550      numberPrimalInfeasibilities_ ++;
2551    } else if (distanceDown<-primalTolerance) {
2552      double infeasibility=-distanceDown;
2553      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
2554      if (infeasibility>relaxedToleranceP) 
2555        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedToleranceP;
2556      numberPrimalInfeasibilities_ ++;
2557    } else {
2558      // feasible (so could be free)
2559      if (getStatus(iSequence) != basic&&!flagged(iSequence)) {
2560        // not basic
2561        double djValue = dj_[iSequence];
2562        if (distanceDown<primalTolerance) {
2563          if (distanceUp>primalTolerance&&djValue<-dualTolerance) {
2564            sumDualInfeasibilities_ -= djValue+dualTolerance;
2565            if (djValue<-relaxedToleranceD) 
2566              sumOfRelaxedDualInfeasibilities_ -= djValue+relaxedToleranceD;
2567            numberDualInfeasibilities_ ++;
2568          } 
2569        } else if (distanceUp<primalTolerance) {
2570          if (djValue>dualTolerance) {
2571            sumDualInfeasibilities_ += djValue-dualTolerance;
2572            if (djValue>relaxedToleranceD) 
2573              sumOfRelaxedDualInfeasibilities_ += djValue-relaxedToleranceD;
2574            numberDualInfeasibilities_ ++;
2575          }
2576        } else {
2577          // may be free
2578          djValue *= 0.01;
2579          if (fabs(djValue)>dualTolerance) {
2580            if (getStatus(iSequence) == isFree) 
2581              numberDualInfeasibilitiesFree++;
2582            sumDualInfeasibilities_ += fabs(djValue)-dualTolerance;
2583            numberDualInfeasibilities_ ++;
2584            if (fabs(djValue)>relaxedToleranceD) {
2585              sumOfRelaxedDualInfeasibilities_ += value-relaxedToleranceD;
2586              numberSuperBasicWithDj++;
2587              if (firstFreeDual<0)
2588                firstFreeDual = iSequence;
2589            }
2590          }
2591          if (firstFreePrimal<0)
2592            firstFreePrimal = iSequence;
2593        }
2594      }
2595    }
2596  }
2597  objectiveValue_ += objective_->nonlinearOffset();
2598  objectiveValue_ /= (objectiveScale_*rhsScale_);
2599  numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_-
2600    numberDualInfeasibilitiesFree;
2601  if (algorithm_<0&&firstFreeDual>=0) {
2602    // dual
2603    firstFree_ = firstFreeDual;
2604  } else if (numberSuperBasicWithDj||
2605             (progress_.lastIterationNumber(0)<=0)) {
2606    firstFree_=firstFreePrimal;
2607  }
2608}
2609/* Adds multiple of a column into an array */
2610void 
2611ClpSimplex::add(double * array,
2612                int sequence, double multiplier) const
2613{
2614  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2615    //slack
2616    array [sequence-numberColumns_] -= multiplier;
2617  } else {
2618    // column
2619    matrix_->add(this,array,sequence,multiplier);
2620  }
2621}
2622/*
2623  Unpacks one column of the matrix into indexed array
2624*/
2625void 
2626ClpSimplex::unpack(CoinIndexedVector * rowArray) const
2627{
2628  rowArray->clear();
2629  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2630    //slack
2631    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
2632  } else {
2633    // column
2634    matrix_->unpack(this,rowArray,sequenceIn_);
2635  }
2636}
2637void 
2638ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence) const
2639{
2640  rowArray->clear();
2641  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2642    //slack
2643    rowArray->insert(sequence-numberColumns_,-1.0);
2644  } else {
2645    // column
2646    matrix_->unpack(this,rowArray,sequence);
2647  }
2648}
2649/*
2650  Unpacks one column of the matrix into indexed array
2651*/
2652void 
2653ClpSimplex::unpackPacked(CoinIndexedVector * rowArray) 
2654{
2655  rowArray->clear();
2656  if (sequenceIn_>=numberColumns_&&sequenceIn_<numberColumns_+numberRows_) {
2657    //slack
2658    int * index = rowArray->getIndices();
2659    double * array = rowArray->denseVector();
2660    array[0]=-1.0;
2661    index[0]=sequenceIn_-numberColumns_;
2662    rowArray->setNumElements(1);
2663    rowArray->setPackedMode(true);
2664  } else {
2665    // column
2666    matrix_->unpackPacked(this,rowArray,sequenceIn_);
2667  }
2668}
2669void 
2670ClpSimplex::unpackPacked(CoinIndexedVector * rowArray,int sequence)
2671{
2672  rowArray->clear();
2673  if (sequence>=numberColumns_&&sequence<numberColumns_+numberRows_) {
2674    //slack
2675    int * index = rowArray->getIndices();
2676    double * array = rowArray->denseVector();
2677    array[0]=-1.0;
2678    index[0]=sequence-numberColumns_;
2679    rowArray->setNumElements(1);
2680    rowArray->setPackedMode(true);
2681  } else {
2682    // column
2683    matrix_->unpackPacked(this,rowArray,sequence);
2684  }
2685}
2686//static int scale_times[]={0,0,0,0};
2687bool
2688ClpSimplex::createRim(int what,bool makeRowCopy, int startFinishOptions)
2689{
2690  bool goodMatrix=true;
2691  spareIntArray_[0]=0;
2692  if (!matrix_->canGetRowCopy())
2693    makeRowCopy=false; // switch off row copy if can't produce
2694  int saveLevel=handler_->logLevel();
2695  // Arrays will be there and correct size unless what is 63
2696  bool newArrays = (what==63);
2697  // We may be restarting with same size
2698  bool keepPivots=false;
2699  if (startFinishOptions==-1) {
2700    startFinishOptions=0;
2701    keepPivots=true;
2702  }
2703  bool oldMatrix = ((startFinishOptions&4)!=0&&(whatsChanged_&1)!=0);
2704  if (auxiliaryModel_) {
2705    if (auxiliaryModel_->numberRows_==numberRows_&&
2706        auxiliaryModel_->numberColumns_==numberColumns_&&
2707        (whatsChanged_&511)==511) {
2708      oldMatrix=true; 
2709    } else {
2710      deleteAuxiliaryModel();
2711      oldMatrix=false;
2712    }
2713  }
2714  if (what==63) {
2715    pivotRow_=-1; 
2716    if (!status_)
2717      createStatus();
2718    if (oldMatrix)
2719      newArrays=false;
2720    if (problemStatus_==10) {
2721      handler_->setLogLevel(0); // switch off messages
2722      if (rowArray_[0]) {
2723        // stuff is still there
2724        oldMatrix=true;
2725        newArrays=false;
2726        keepPivots=true;
2727        for (int iRow=0;iRow<4;iRow++) {
2728          rowArray_[iRow]->clear();
2729        }
2730        for (int iColumn=0;iColumn<2;iColumn++) {
2731          columnArray_[iColumn]->clear();
2732        }
2733      }
2734    } else if (factorization_) {
2735      // match up factorization messages
2736      if (handler_->logLevel()<3)
2737        factorization_->messageLevel(0);
2738      else
2739        factorization_->messageLevel(CoinMax(3,factorization_->messageLevel()));
2740      /* Faster to keep pivots rather than re-scan matrix.  Matrix may have changed
2741         i.e. oldMatrix false but okay as long as same number rows and status array exists
2742      */
2743      if ((startFinishOptions&2)!=0&&factorization_->numberRows()==numberRows_&&status_)
2744        keepPivots=true;
2745    }
2746    numberExtraRows_ = matrix_->generalExpanded(this,2,maximumBasic_);
2747    if (numberExtraRows_&&newArrays) {
2748      // make sure status array large enough
2749      assert (status_);
2750      int numberOld = numberRows_+numberColumns_;
2751      int numberNew = numberRows_+numberColumns_+numberExtraRows_;
2752      unsigned char * newStatus = new unsigned char [numberNew];
2753      memset(newStatus+numberOld,0,numberExtraRows_);
2754      CoinMemcpyN(status_,numberOld,newStatus);
2755      delete [] status_;
2756      status_=newStatus;
2757    }
2758  }
2759  int numberRows2 = numberRows_+numberExtraRows_;
2760  int numberTotal = numberRows2+numberColumns_;
2761  if ((specialOptions_&65536)!=0&&!auxiliaryModel_) {
2762    assert (!numberExtraRows_);
2763    if (!cost_||numberRows2>maximumInternalRows_||
2764        numberColumns_>maximumInternalColumns_) {
2765      newArrays=true;
2766      keepPivots=false;
2767      printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
2768             numberRows_,maximumRows_,maximumInternalRows_);
2769      int oldMaximumRows=maximumInternalRows_;
2770      int oldMaximumColumns=maximumInternalColumns_;
2771      if (cost_) {
2772        if (numberRows2>maximumInternalRows_) 
2773          maximumInternalRows_ = numberRows2;
2774        if (numberColumns_>maximumInternalColumns_) 
2775          maximumInternalColumns_ = numberColumns_;
2776      } else {
2777        maximumInternalRows_ = numberRows2;
2778        maximumInternalColumns_ = numberColumns_;
2779      }
2780      assert(maximumInternalRows_ == maximumRows_);
2781      assert(maximumInternalColumns_ == maximumColumns_);
2782      printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
2783             numberRows_,maximumRows_,maximumInternalRows_);
2784      int numberTotal2 = (maximumInternalRows_+maximumInternalColumns_)*2;
2785      delete [] cost_;
2786      cost_ = new double[numberTotal2];
2787      delete [] lower_;
2788      delete [] upper_;
2789      lower_ = new double[numberTotal2];
2790      upper_ = new double[numberTotal2];
2791      delete [] dj_;
2792      dj_ = new double[numberTotal2];
2793      delete [] solution_;
2794      solution_ = new double[numberTotal2];
2795      // ***** should be non NULL but seems to be too much
2796      //printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
2797      if (savedRowScale_) {
2798        assert (oldMaximumRows>0);
2799        double * temp;
2800        temp = new double [4*maximumRows_];
2801        CoinFillN(temp,4*maximumRows_,1.0);
2802        CoinMemcpyN(savedRowScale_,numberRows_,temp);
2803        CoinMemcpyN(savedRowScale_+oldMaximumRows,numberRows_,temp+maximumRows_);
2804        CoinMemcpyN(savedRowScale_+2*oldMaximumRows,numberRows_,temp+2*maximumRows_);
2805        CoinMemcpyN(savedRowScale_+3*oldMaximumRows,numberRows_,temp+3*maximumRows_);
2806        delete [] savedRowScale_;
2807        savedRowScale_ = temp;
2808        temp = new double [4*maximumColumns_];
2809        CoinFillN(temp,4*maximumColumns_,1.0);
2810        CoinMemcpyN(savedColumnScale_,numberColumns_,temp);
2811        CoinMemcpyN(savedColumnScale_+oldMaximumColumns,numberColumns_,temp+maximumColumns_);
2812        CoinMemcpyN(savedColumnScale_+2*oldMaximumColumns,numberColumns_,temp+2*maximumColumns_);
2813        CoinMemcpyN(savedColumnScale_+3*oldMaximumColumns,numberColumns_,temp+3*maximumColumns_);
2814        delete [] savedColumnScale_;
2815        savedColumnScale_ = temp;
2816      }
2817    }
2818  }
2819  int i;
2820  bool doSanityCheck=true;
2821  if (what==63&&!auxiliaryModel_) {
2822    // We may want to switch stuff off for speed
2823    if ((specialOptions_&256)!=0)
2824      makeRowCopy=false; // no row copy
2825    if ((specialOptions_&128)!=0)
2826      doSanityCheck=false; // no sanity check
2827    //check matrix
2828    if (!matrix_)
2829      matrix_=new ClpPackedMatrix();
2830    int checkType=(doSanityCheck) ? 15 : 14;
2831    if (oldMatrix)
2832      checkType = 14;
2833    if ((specialOptions_&COIN_CBC_USING_CLP)!=0)
2834      checkType -= 4; // don't check for duplicates
2835    if (!matrix_->allElementsInRange(this,smallElement_,1.0e20,checkType)) {
2836      problemStatus_=4;
2837      secondaryStatus_=8;
2838      //goodMatrix= false;
2839      return false;
2840    }
2841    if (makeRowCopy&&!oldMatrix) {
2842      delete rowCopy_;
2843      // may return NULL if can't give row copy
2844      rowCopy_ = matrix_->reverseOrderedCopy();
2845    }
2846#if 0
2847    if (what==63&&(specialOptions_&131072)!=0) {
2848      int k=rowScale_ ? 1 : 0;
2849      if (oldMatrix)
2850        k+=2;
2851      scale_times[k]++;
2852      if ((scale_times[0]+scale_times[1]+scale_times[2]+scale_times[3])%1000==0)
2853        printf("scale counts %d %d %d %d\n",
2854               scale_times[0],scale_times[1],scale_times[2],scale_times[3]);
2855    }
2856#endif
2857    // do scaling if needed
2858    if (!oldMatrix&&scalingFlag_<0) {
2859      if (scalingFlag_<0&&rowScale_) {
2860        //if (handler_->logLevel()>0)
2861          printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
2862                 scalingFlag_);
2863        scalingFlag_=0;
2864      }
2865      delete [] rowScale_;
2866      delete [] columnScale_;
2867      rowScale_=NULL;
2868      columnScale_=NULL;
2869    }
2870    inverseRowScale_ = NULL;
2871    inverseColumnScale_ = NULL;
2872    if (scalingFlag_>0&&!rowScale_) {
2873      if ((specialOptions_&131072)!=0) {
2874        inverseRowScale_ = rowScale_+numberRows2;
2875        inverseColumnScale_ = columnScale_+numberColumns_;
2876      }
2877      if ((specialOptions_&65536)!=0&&!auxiliaryModel_) {
2878        assert (!rowScale_);
2879        rowScale_ = savedRowScale_;
2880        columnScale_ = savedColumnScale_;
2881        // put back original
2882        if (savedRowScale_) {
2883          inverseRowScale_ = savedRowScale_+maximumInternalRows_;
2884          inverseColumnScale_ = savedColumnScale_+maximumInternalColumns_;
2885          CoinMemcpyN(savedRowScale_+2*maximumInternalRows_,
2886                      numberRows2,savedRowScale_);
2887          CoinMemcpyN(savedRowScale_+3*maximumInternalRows_,
2888                      numberRows2,inverseRowScale_);
2889          CoinMemcpyN(savedColumnScale_+2*maximumColumns_,
2890                      numberColumns_,savedColumnScale_);
2891          CoinMemcpyN(savedColumnScale_+3*maximumColumns_,
2892                      numberColumns_,inverseColumnScale_);
2893        }
2894      }
2895      if (matrix_->scale(this))
2896        scalingFlag_=-scalingFlag_; // not scaled after all
2897      if ((specialOptions_&131072)!=0&&rowScale_&&!savedRowScale_) {
2898        inverseRowScale_ = rowScale_+numberRows2;
2899        inverseColumnScale_ = columnScale_+numberColumns_;
2900      }
2901      if (rowScale_&&automaticScale_) {
2902        // try automatic scaling
2903        double smallestObj=1.0e100;
2904        double largestObj=0.0;
2905        double largestRhs=0.0;
2906        const double * obj = objective();
2907        for (i=0;i<numberColumns_;i++) {
2908          double value = fabs(obj[i]);
2909          value *= columnScale_[i];
2910          if (value&&columnLower_[i]!=columnUpper_[i]) {
2911            smallestObj = CoinMin(smallestObj,value);
2912            largestObj = CoinMax(largestObj,value);
2913          }
2914          if (columnLower_[i]>0.0||columnUpper_[i]<0.0) {
2915            double scale = 1.0/columnScale_[i];
2916            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2917            if (columnLower_[i]>0)
2918              largestRhs=CoinMax(largestRhs,columnLower_[i]*scale);
2919            if (columnUpper_[i]<0.0)
2920              largestRhs=CoinMax(largestRhs,-columnUpper_[i]*scale);
2921          }
2922        }
2923        for (i=0;i<numberRows_;i++) {
2924          if (rowLower_[i]>0.0||rowUpper_[i]<0.0) {
2925            double scale = rowScale_[i];
2926            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2927            if (rowLower_[i]>0)
2928              largestRhs=CoinMax(largestRhs,rowLower_[i]*scale);
2929            if (rowUpper_[i]<0.0)
2930              largestRhs=CoinMax(largestRhs,-rowUpper_[i]*scale);
2931          }
2932        }
2933        printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2934        bool scalingDone=false;
2935        // look at element range
2936        double smallestNegative;
2937        double largestNegative;
2938        double smallestPositive;
2939        double largestPositive;
2940        matrix_->rangeOfElements(smallestNegative, largestNegative,
2941                                 smallestPositive, largestPositive);
2942        smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2943        largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2944        if (largestObj) {
2945          double ratio = largestObj/smallestObj;
2946          double scale=1.0;
2947          if (ratio<1.0e8) {
2948            // reasonable
2949            if (smallestObj<1.0e-4) {
2950              // may as well scale up
2951              scalingDone=true;
2952              scale=1.0e-3/smallestObj;
2953            } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2954              //done=true;
2955            } else {
2956              scalingDone=true;
2957              if (algorithm_<0) {
2958                scale = 1.0e6/largestObj;
2959              } else {
2960                scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2961              }
2962            }
2963          } else if (ratio<1.0e12) {
2964            // not so good
2965            if (smallestObj<1.0e-7) {
2966              // may as well scale up
2967              scalingDone=true;
2968              scale=1.0e-6/smallestObj;
2969            } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2970              //done=true;
2971            } else {
2972              scalingDone=true;
2973              if (algorithm_<0) {
2974                scale = 1.0e7/largestObj;
2975              } else {
2976                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2977              }
2978            }
2979          } else {
2980            // Really nasty problem
2981            if (smallestObj<1.0e-8) {
2982              // may as well scale up
2983              scalingDone=true;
2984              scale=1.0e-7/smallestObj;
2985              largestObj *= scale;
2986            } 
2987            if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2988              //done=true;
2989            } else {
2990              scalingDone=true;
2991              if (algorithm_<0) {
2992                scale = 1.0e7/largestObj;
2993              } else {
2994                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2995              }
2996            }
2997          }
2998          objectiveScale_=scale;
2999        }
3000        if (largestRhs>1.0e12) {
3001          scalingDone=true;
3002          rhsScale_=1.0e9/largestRhs;
3003        } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
3004          scalingDone=true;
3005          rhsScale_=1.0e6/largestRhs;
3006        } else {
3007          rhsScale_=1.0;
3008        }
3009        if (scalingDone) {
3010          handler_->message(CLP_RIM_SCALE,messages_)
3011            <<objectiveScale_<<rhsScale_
3012            <<CoinMessageEol;
3013        }
3014      }
3015    } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
3016      matrix_->scaleRowCopy(this);
3017    }
3018    // See if we can try for faster row copy
3019    if (makeRowCopy&&!oldMatrix) {
3020      ClpPackedMatrix* clpMatrix =
3021        dynamic_cast< ClpPackedMatrix*>(matrix_);
3022      if (clpMatrix&&numberThreads_) 
3023        clpMatrix->specialRowCopy(this,rowCopy_);
3024      if (clpMatrix)
3025        clpMatrix->specialColumnCopy(this);
3026    }
3027  }
3028  if (what==63) {
3029    if (newArrays&&(specialOptions_&65536)==0) {
3030      delete [] cost_;
3031      cost_ = new double[numberTotal];
3032      delete [] lower_;
3033      delete [] upper_;
3034      lower_ = new double[numberTotal];
3035      upper_ = new double[numberTotal];
3036      delete [] dj_;
3037      dj_ = new double[numberTotal];
3038      delete [] solution_;
3039      solution_ = new double[numberTotal];
3040    } else if (auxiliaryModel_) {
3041      assert (!cost_);
3042      cost_=auxiliaryModel_->cost_;
3043      assert (!lower_);
3044      lower_=auxiliaryModel_->lower_;
3045      assert (!upper_);
3046      upper_=auxiliaryModel_->upper_;
3047      assert (!dj_);
3048      dj_=auxiliaryModel_->dj_;
3049      assert (!solution_);
3050      solution_=auxiliaryModel_->solution_;
3051      assert (!rowScale_);
3052      assert (!columnScale_);
3053    }
3054    reducedCostWork_ = dj_;
3055    rowReducedCost_ = dj_+numberColumns_;
3056    columnActivityWork_ = solution_;
3057    rowActivityWork_ = solution_+numberColumns_;
3058    objectiveWork_ = cost_;
3059    rowObjectiveWork_ = cost_+numberColumns_;
3060    rowLowerWork_ = lower_+numberColumns_;
3061    columnLowerWork_ = lower_;
3062    rowUpperWork_ = upper_+numberColumns_;
3063    columnUpperWork_ = upper_;
3064  }
3065  if ((what&4)!=0) {
3066    if (!auxiliaryModel_||(what==63&&(auxiliaryModel_->specialOptions_&4)==0)) {
3067      double direction = optimizationDirection_*objectiveScale_;
3068      const double * obj = objective();
3069      const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
3070      const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
3071      // and also scale by scale factors
3072      if (rowScale) {
3073        if (rowObjective_) {
3074          for (i=0;i<numberRows_;i++)
3075            rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
3076        } else {
3077          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
3078        }
3079        // If scaled then do all columns later in one loop
3080        if (what!=63||auxiliaryModel_) {
3081          for (i=0;i<numberColumns_;i++) {
3082            CoinAssert(fabs(obj[i])<1.0e25);
3083            objectiveWork_[i] = obj[i]*direction*columnScale[i];
3084          }
3085        }
3086      } else {
3087        if (rowObjective_) {
3088          for (i=0;i<numberRows_;i++)
3089            rowObjectiveWork_[i] = rowObjective_[i]*direction;
3090        } else {
3091          memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
3092        }
3093        for (i=0;i<numberColumns_;i++) {
3094          CoinAssert(fabs(obj[i])<1.0e25);
3095          objectiveWork_[i] = obj[i]*direction;
3096        }
3097      }
3098      if (auxiliaryModel_) {
3099        // save costs
3100        CoinMemcpyN(cost_,numberTotal,auxiliaryModel_->cost_+numberTotal);
3101      }
3102    } else {
3103      // just copy
3104      CoinMemcpyN(auxiliaryModel_->cost_+numberTotal,numberTotal,cost_);
3105    }
3106  }
3107  if ((what&1)!=0) {
3108    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
3109    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
3110    // clean up any mismatches on infinity
3111    // and fix any variables with tiny gaps
3112    double primalTolerance=dblParam_[ClpPrimalTolerance];
3113    if (!auxiliaryModel_) {
3114      if(rowScale) {
3115        // If scaled then do all columns later in one loop
3116        if (what!=63) {
3117          if (!inverseColumnScale_) {
3118            for (i=0;i<numberColumns_;i++) {
3119              double multiplier = rhsScale_/columnScale[i];
3120              double lowerValue = columnLower_[i];
3121              double upperValue = columnUpper_[i];
3122              if (lowerValue>-1.0e20) {
3123                columnLowerWork_[i]=lowerValue*multiplier;
3124                if (upperValue>=1.0e20) {
3125                  columnUpperWork_[i]=COIN_DBL_MAX;
3126                } else {
3127                  columnUpperWork_[i]=upperValue*multiplier;
3128                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3129                    if (columnLowerWork_[i]>=0.0) {
3130                      columnUpperWork_[i] = columnLowerWork_[i];
3131                    } else if (columnUpperWork_[i]<=0.0) {
3132                      columnLowerWork_[i] = columnUpperWork_[i];
3133                    } else {
3134                      columnUpperWork_[i] = 0.0;
3135                      columnLowerWork_[i] = 0.0;
3136                    }
3137                  }
3138                }
3139              } else if (upperValue<1.0e20) {
3140                columnLowerWork_[i]=-COIN_DBL_MAX;
3141                columnUpperWork_[i]=upperValue*multiplier;
3142              } else {
3143                // free
3144                columnLowerWork_[i]=-COIN_DBL_MAX;
3145                columnUpperWork_[i]=COIN_DBL_MAX;
3146              }
3147            }
3148          } else {
3149            assert (inverseColumnScale_);
3150            const double * inverseScale = inverseColumnScale_;
3151            for (i=0;i<numberColumns_;i++) {
3152              double multiplier = rhsScale_*inverseScale[i];
3153              double lowerValue = columnLower_[i];
3154              double upperValue = columnUpper_[i];
3155              if (lowerValue>-1.0e20) {
3156                columnLowerWork_[i]=lowerValue*multiplier;
3157                if (upperValue>=1.0e20) {
3158                  columnUpperWork_[i]=COIN_DBL_MAX;
3159                } else {
3160                  columnUpperWork_[i]=upperValue*multiplier;
3161                  if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3162                    if (columnLowerWork_[i]>=0.0) {
3163                      columnUpperWork_[i] = columnLowerWork_[i];
3164                    } else if (columnUpperWork_[i]<=0.0) {
3165                      columnLowerWork_[i] = columnUpperWork_[i];
3166                    } else {
3167                      columnUpperWork_[i] = 0.0;
3168                      columnLowerWork_[i] = 0.0;
3169                    }
3170                  }
3171                }
3172              } else if (upperValue<1.0e20) {
3173                columnLowerWork_[i]=-COIN_DBL_MAX;
3174                columnUpperWork_[i]=upperValue*multiplier;
3175              } else {
3176                // free
3177                columnLowerWork_[i]=-COIN_DBL_MAX;
3178                columnUpperWork_[i]=COIN_DBL_MAX;
3179              }
3180            }
3181          }
3182        }
3183        for (i=0;i<numberRows_;i++) {
3184          double multiplier = rhsScale_*rowScale[i];
3185          double lowerValue = rowLower_[i];
3186          double upperValue = rowUpper_[i];
3187          if (lowerValue>-1.0e20) {
3188            rowLowerWork_[i]=lowerValue*multiplier;
3189            if (upperValue>=1.0e20) {
3190              rowUpperWork_[i]=COIN_DBL_MAX;
3191            } else {
3192              rowUpperWork_[i]=upperValue*multiplier;
3193              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
3194                if (rowLowerWork_[i]>=0.0) {
3195                  rowUpperWork_[i] = rowLowerWork_[i];
3196                } else if (rowUpperWork_[i]<=0.0) {
3197                  rowLowerWork_[i] = rowUpperWork_[i];
3198                } else {
3199                  rowUpperWork_[i] = 0.0;
3200                  rowLowerWork_[i] = 0.0;
3201                }
3202              }
3203            }
3204          } else if (upperValue<1.0e20) {
3205            rowLowerWork_[i]=-COIN_DBL_MAX;
3206            rowUpperWork_[i]=upperValue*multiplier;
3207          } else {
3208            // free
3209            rowLowerWork_[i]=-COIN_DBL_MAX;
3210            rowUpperWork_[i]=COIN_DBL_MAX;
3211          }
3212        }
3213      } else if (rhsScale_!=1.0) {
3214        for (i=0;i<numberColumns_;i++) {
3215          double lowerValue = columnLower_[i];
3216          double upperValue = columnUpper_[i];
3217          if (lowerValue>-1.0e20) {
3218            columnLowerWork_[i]=lowerValue*rhsScale_;
3219            if (upperValue>=1.0e20) {
3220              columnUpperWork_[i]=COIN_DBL_MAX;
3221            } else {
3222              columnUpperWork_[i]=upperValue*rhsScale_;
3223              if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3224                if (columnLowerWork_[i]>=0.0) {
3225                  columnUpperWork_[i] = columnLowerWork_[i];
3226                } else if (columnUpperWork_[i]<=0.0) {
3227                  columnLowerWork_[i] = columnUpperWork_[i];
3228                } else {
3229                  columnUpperWork_[i] = 0.0;
3230                  columnLowerWork_[i] = 0.0;
3231                }
3232              }
3233            }
3234          } else if (upperValue<1.0e20) {
3235            columnLowerWork_[i]=-COIN_DBL_MAX;
3236            columnUpperWork_[i]=upperValue*rhsScale_;
3237          } else {
3238            // free
3239            columnLowerWork_[i]=-COIN_DBL_MAX;
3240            columnUpperWork_[i]=COIN_DBL_MAX;
3241          }
3242        }
3243        for (i=0;i<numberRows_;i++) {
3244          double lowerValue = rowLower_[i];
3245          double upperValue = rowUpper_[i];
3246          if (lowerValue>-1.0e20) {
3247            rowLowerWork_[i]=lowerValue*rhsScale_;
3248            if (upperValue>=1.0e20) {
3249              rowUpperWork_[i]=COIN_DBL_MAX;
3250            } else {
3251              rowUpperWork_[i]=upperValue*rhsScale_;
3252              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
3253                if (rowLowerWork_[i]>=0.0) {
3254                  rowUpperWork_[i] = rowLowerWork_[i];
3255                } else if (rowUpperWork_[i]<=0.0) {
3256                  rowLowerWork_[i] = rowUpperWork_[i];
3257                } else {
3258                  rowUpperWork_[i] = 0.0;
3259                  rowLowerWork_[i] = 0.0;
3260                }
3261              }
3262            }
3263          } else if (upperValue<1.0e20) {
3264            rowLowerWork_[i]=-COIN_DBL_MAX;
3265            rowUpperWork_[i]=upperValue*rhsScale_;
3266          } else {
3267            // free
3268            rowLowerWork_[i]=-COIN_DBL_MAX;
3269            rowUpperWork_[i]=COIN_DBL_MAX;
3270          }
3271        }
3272      } else {
3273        for (i=0;i<numberColumns_;i++) {
3274          double lowerValue = columnLower_[i];
3275          double upperValue = columnUpper_[i];
3276          if (lowerValue>-1.0e20) {
3277            columnLowerWork_[i]=lowerValue;
3278            if (upperValue>=1.0e20) {
3279              columnUpperWork_[i]=COIN_DBL_MAX;
3280            } else {
3281              columnUpperWork_[i]=upperValue;
3282              if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3283                if (columnLowerWork_[i]>=0.0) {
3284                  columnUpperWork_[i] = columnLowerWork_[i];
3285                } else if (columnUpperWork_[i]<=0.0) {
3286                  columnLowerWork_[i] = columnUpperWork_[i];
3287                } else {
3288                  columnUpperWork_[i] = 0.0;
3289                  columnLowerWork_[i] = 0.0;
3290                }
3291              }
3292            }
3293          } else if (upperValue<1.0e20) {
3294            columnLowerWork_[i]=-COIN_DBL_MAX;
3295            columnUpperWork_[i]=upperValue;
3296          } else {
3297            // free
3298            columnLowerWork_[i]=-COIN_DBL_MAX;
3299            columnUpperWork_[i]=COIN_DBL_MAX;
3300          }
3301        }
3302        for (i=0;i<numberRows_;i++) {
3303          double lowerValue = rowLower_[i];
3304          double upperValue = rowUpper_[i];
3305          if (lowerValue>-1.0e20) {
3306            rowLowerWork_[i]=lowerValue;
3307            if (upperValue>=1.0e20) {
3308              rowUpperWork_[i]=COIN_DBL_MAX;
3309            } else {
3310              rowUpperWork_[i]=upperValue;
3311              if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
3312                if (rowLowerWork_[i]>=0.0) {
3313                  rowUpperWork_[i] = rowLowerWork_[i];
3314                } else if (rowUpperWork_[i]<=0.0) {
3315                  rowLowerWork_[i] = rowUpperWork_[i];
3316                } else {
3317                  rowUpperWork_[i] = 0.0;
3318                  rowLowerWork_[i] = 0.0;
3319                }
3320              }
3321            }
3322          } else if (upperValue<1.0e20) {
3323            rowLowerWork_[i]=-COIN_DBL_MAX;
3324            rowUpperWork_[i]=upperValue;
3325          } else {
3326            // free
3327            rowLowerWork_[i]=-COIN_DBL_MAX;
3328            rowUpperWork_[i]=COIN_DBL_MAX;
3329          }
3330        }
3331      }
3332    } else {
3333      // auxiliary model
3334      if (what!=63) {
3335        // just copy
3336        CoinMemcpyN(auxiliaryModel_->lower_+numberTotal,numberTotal,lower_);
3337        CoinMemcpyN(auxiliaryModel_->upper_+numberTotal,numberTotal,upper_);
3338      } else {
3339        assert (rhsScale_==1.0);
3340        assert (objectiveScale_==1.0);
3341        if ((auxiliaryModel_->specialOptions_&2)==0) {
3342          // need to do column bounds
3343          if (columnScale) {
3344            const double * inverseScale = columnScale+numberColumns_;
3345            // scaled
3346            for (i=0;i<numberColumns_;i++) {
3347              double value;
3348              value = columnLower_[i];
3349              if (value>-1.0e20) 
3350                value *= inverseScale[i];
3351              lower_[i] = value;
3352              value = columnUpper_[i];
3353              if (value<1.0e20) 
3354                value *= inverseScale[i];
3355              upper_[i] = value;
3356            }
3357          } else {
3358            for (i=0;i<numberColumns_;i++) {
3359              lower_[i] = columnLower_[i];
3360              upper_[i] = columnUpper_[i];
3361            }
3362          }
3363          // save
3364          CoinMemcpyN(lower_,numberColumns_,auxiliaryModel_->lower_+numberTotal);
3365          CoinMemcpyN(upper_,numberColumns_,auxiliaryModel_->upper_+numberTotal);
3366        } else {
3367          // just copy
3368          CoinMemcpyN(auxiliaryModel_->lower_+numberTotal,numberColumns_,lower_);
3369          CoinMemcpyN(auxiliaryModel_->upper_+numberTotal,numberColumns_,upper_);
3370        }
3371        if ((auxiliaryModel_->specialOptions_&1)==0) {
3372          // need to do row bounds
3373          if (rowScale) {
3374            // scaled
3375            for (i=0;i<numberRows_;i++) {
3376              double value;
3377              value = rowLower_[i];
3378              if (value>-1.0e20) 
3379                value *= rowScale[i];
3380              lower_[i+numberColumns_] = value;
3381              value = rowUpper_[i];
3382              if (value<1.0e20) 
3383                value *= rowScale[i];
3384              upper_[i+numberColumns_] = value;
3385            }
3386          } else {
3387            for (i=0;i<numberRows_;i++) {
3388              lower_[i+numberColumns_] = rowLower_[i];
3389              upper_[i+numberColumns_] = rowUpper_[i];
3390            }
3391          }
3392          // save
3393          CoinMemcpyN(                 lower_+numberColumns_,numberRows_,auxiliaryModel_->lower_+numberTotal+numberColumns_);
3394          CoinMemcpyN(                 upper_+numberColumns_,numberRows_,auxiliaryModel_->upper_+numberTotal+numberColumns_);
3395        } else {
3396          // just copy
3397          CoinMemcpyN(                 auxiliaryModel_->lower_+numberTotal+numberColumns_,
3398                       numberRows_,lower_+numberColumns_);
3399          CoinMemcpyN(                 auxiliaryModel_->upper_+numberTotal+numberColumns_,
3400                       numberRows_,upper_+numberColumns_);
3401        }
3402      }
3403      if (what==63) {
3404        // do basis
3405        assert ((auxiliaryModel_->specialOptions_&8)!=0);
3406        // clean solution trusting basis
3407        for ( i=0;i<numberTotal;i++) {
3408          Status status =getStatus(i);
3409          double value=solution_[i]; // may as well leave if basic (values pass)
3410          if (status!=basic) {
3411            if (upper_[i]==lower_[i]) {
3412              setStatus(i,isFixed);
3413              value=lower_[i];
3414            } else if (status==atLowerBound) {
3415              if (lower_[i]>-1.0e20) {
3416                value=lower_[i];
3417              } else {
3418                if (upper_[i]<1.0e20) {
3419                  value=upper_[i];
3420                  setStatus(i,atUpperBound);
3421                } else {
3422                  setStatus(i,isFree);
3423                }
3424              }
3425            } else if (status==atUpperBound) {
3426              if (upper_[i]<1.0e20) {
3427              value=upper_[i];
3428              } else {
3429                if (lower_[i]>-1.0e20) {
3430                  value=lower_[i];
3431                  setStatus(i,atLowerBound);
3432                } else {
3433                  setStatus(i,isFree);
3434                }
3435              }
3436            }
3437          }
3438          solution_[i]=value;
3439        }
3440      }
3441    }
3442  }
3443  if (what==63) {
3444    // move information to work arrays
3445    double direction = optimizationDirection_;
3446    // direction is actually scale out not scale in
3447    if (direction)
3448      direction = 1.0/direction;
3449    if (direction!=1.0&&(!auxiliaryModel_||(auxiliaryModel_->specialOptions_&8)==0)) {
3450      // reverse all dual signs
3451      for (i=0;i<numberColumns_;i++) 
3452        reducedCost_[i] *= direction;
3453      for (i=0;i<numberRows_;i++) 
3454        dual_[i] *= direction;
3455    }
3456    for (i=0;i<numberRows_+numberColumns_;i++) {
3457      setFakeBound(i,noFake);
3458    }
3459    if (!auxiliaryModel_) {
3460      if (rowScale_) {
3461        const double * obj = objective();
3462        double direction = optimizationDirection_*objectiveScale_;
3463        // clean up any mismatches on infinity
3464        // and fix any variables with tiny gaps
3465        double primalTolerance=dblParam_[ClpPrimalTolerance];
3466        // on entry
3467        if (!inverseColumnScale_) {
3468          for (i=0;i<numberColumns_;i++) {
3469            CoinAssert(fabs(obj[i])<1.0e25);
3470            double scaleFactor = columnScale_[i];
3471            double multiplier = rhsScale_/scaleFactor;
3472            scaleFactor *= direction;
3473            objectiveWork_[i] = obj[i]*scaleFactor;
3474            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
3475            double lowerValue = columnLower_[i];
3476            double upperValue = columnUpper_[i];
3477            if (lowerValue>-1.0e20) {
3478              columnLowerWork_[i]=lowerValue*multiplier;
3479              if (upperValue>=1.0e20) {
3480                columnUpperWork_[i]=COIN_DBL_MAX;
3481              } else {
3482                columnUpperWork_[i]=upperValue*multiplier;
3483                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3484                  if (columnLowerWork_[i]>=0.0) {
3485                    columnUpperWork_[i] = columnLowerWork_[i];
3486                  } else if (columnUpperWork_[i]<=0.0) {
3487                    columnLowerWork_[i] = columnUpperWork_[i];
3488                  } else {
3489                    columnUpperWork_[i] = 0.0;
3490                    columnLowerWork_[i] = 0.0;
3491                  }
3492                }
3493              }
3494            } else if (upperValue<1.0e20) {
3495              columnLowerWork_[i]=-COIN_DBL_MAX;
3496              columnUpperWork_[i]=upperValue*multiplier;
3497            } else {
3498              // free
3499              columnLowerWork_[i]=-COIN_DBL_MAX;
3500              columnUpperWork_[i]=COIN_DBL_MAX;
3501            }
3502            double value = columnActivity_[i] * multiplier;
3503            if (fabs(value)>1.0e20) {
3504              //printf("bad value of %g for column %d\n",value,i);
3505              setColumnStatus(i,superBasic);
3506              if (columnUpperWork_[i]<0.0) {
3507                value=columnUpperWork_[i];
3508              } else if (columnLowerWork_[i]>0.0) {
3509                value=columnLowerWork_[i];
3510              } else {
3511                value=0.0;
3512              }
3513            }
3514            columnActivityWork_[i] = value;
3515          }
3516          for (i=0;i<numberRows_;i++) {
3517            dual_[i] /= rowScale_[i];
3518            dual_[i] *= objectiveScale_;
3519            rowReducedCost_[i] = dual_[i];
3520            double multiplier = rhsScale_*rowScale_[i];
3521            double value = rowActivity_[i] * multiplier;
3522            if (fabs(value)>1.0e20) {
3523              //printf("bad value of %g for row %d\n",value,i);
3524              setRowStatus(i,superBasic);
3525              if (rowUpperWork_[i]<0.0) {
3526                value=rowUpperWork_[i];
3527              } else if (rowLowerWork_[i]>0.0) {
3528                value=rowLowerWork_[i];
3529              } else {
3530                value=0.0;
3531              }
3532            }
3533            rowActivityWork_[i] = value;
3534          }
3535        } else {
3536          assert (inverseColumnScale_);
3537          const double * inverseScale = inverseColumnScale_;
3538          for (i=0;i<numberColumns_;i++) {
3539            CoinAssert(fabs(obj[i])<1.0e25);
3540            double scaleFactor = columnScale_[i];
3541            double multiplier = rhsScale_*inverseScale[i];
3542            scaleFactor *= direction;
3543            objectiveWork_[i] = obj[i]*scaleFactor;
3544            reducedCostWork_[i] = reducedCost_[i]*scaleFactor;
3545            double lowerValue = columnLower_[i];
3546            double upperValue = columnUpper_[i];
3547            if (lowerValue>-1.0e20) {
3548              columnLowerWork_[i]=lowerValue*multiplier;
3549              if (upperValue>=1.0e20) {
3550                columnUpperWork_[i]=COIN_DBL_MAX;
3551              } else {
3552                columnUpperWork_[i]=upperValue*multiplier;
3553                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3554                  if (columnLowerWork_[i]>=0.0) {
3555                    columnUpperWork_[i] = columnLowerWork_[i];
3556                  } else if (columnUpperWork_[i]<=0.0) {
3557                    columnLowerWork_[i] = columnUpperWork_[i];
3558                  } else {
3559                    columnUpperWork_[i] = 0.0;
3560                    columnLowerWork_[i] = 0.0;
3561                  }
3562                }
3563              }
3564            } else if (upperValue<1.0e20) {
3565              columnLowerWork_[i]=-COIN_DBL_MAX;
3566              columnUpperWork_[i]=upperValue*multiplier;
3567            } else {
3568              // free
3569              columnLowerWork_[i]=-COIN_DBL_MAX;
3570              columnUpperWork_[i]=COIN_DBL_MAX;
3571            }
3572            double value = columnActivity_[i] * multiplier;
3573            if (fabs(value)>1.0e20) {
3574              //printf("bad value of %g for column %d\n",value,i);
3575              setColumnStatus(i,superBasic);
3576              if (columnUpperWork_[i]<0.0) {
3577                value=columnUpperWork_[i];
3578              } else if (columnLowerWork_[i]>0.0) {
3579                value=columnLowerWork_[i];
3580              } else {
3581                value=0.0;
3582              }
3583            }
3584            columnActivityWork_[i] = value;
3585          }
3586          inverseScale = inverseRowScale_;
3587          for (i=0;i<numberRows_;i++) {
3588            dual_[i] *= inverseScale[i];
3589            dual_[i] *= objectiveScale_;
3590            rowReducedCost_[i] = dual_[i];
3591            double multiplier = rhsScale_*rowScale_[i];
3592            double value = rowActivity_[i] * multiplier;
3593            if (fabs(value)>1.0e20) {
3594              //printf("bad value of %g for row %d\n",value,i);
3595              setRowStatus(i,superBasic);
3596              if (rowUpperWork_[i]<0.0) {
3597                value=rowUpperWork_[i];
3598              } else if (rowLowerWork_[i]>0.0) {
3599                value=rowLowerWork_[i];
3600              } else {
3601                value=0.0;
3602              }
3603            }
3604            rowActivityWork_[i] = value;
3605          }
3606        }
3607      } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
3608        // on entry
3609        for (i=0;i<numberColumns_;i++) {
3610          double value = columnActivity_[i];
3611          value *= rhsScale_;
3612          if (fabs(value)>1.0e20) {
3613            //printf("bad value of %g for column %d\n",value,i);
3614            setColumnStatus(i,superBasic);
3615            if (columnUpperWork_[i]<0.0) {
3616              value=columnUpperWork_[i];
3617            } else if (columnLowerWork_[i]>0.0) {
3618              value=columnLowerWork_[i];
3619            } else {
3620              value=0.0;
3621            }
3622          }
3623          columnActivityWork_[i] = value;
3624          reducedCostWork_[i] = reducedCost_[i]*objectiveScale_;
3625        }
3626        for (i=0;i<numberRows_;i++) {
3627          double value = rowActivity_[i];
3628          value *= rhsScale_;
3629          if (fabs(value)>1.0e20) {
3630            //printf("bad value of %g for row %d\n",value,i);
3631            setRowStatus(i,superBasic);
3632            if (rowUpperWork_[i]<0.0) {
3633              value=rowUpperWork_[i];
3634            } else if (rowLowerWork_[i]>0.0) {
3635              value=rowLowerWork_[i];
3636            } else {
3637              value=0.0;
3638            }
3639          }
3640          rowActivityWork_[i] = value;
3641          dual_[i] *= objectiveScale_;
3642          rowReducedCost_[i] = dual_[i];
3643        }
3644      } else {
3645        // on entry
3646        for (i=0;i<numberColumns_;i++) {
3647          double value = columnActivity_[i];
3648          if (fabs(value)>1.0e20) {
3649            //printf("bad value of %g for column %d\n",value,i);
3650            setColumnStatus(i,superBasic);
3651            if (columnUpperWork_[i]<0.0) {
3652              value=columnUpperWork_[i];
3653            } else if (columnLowerWork_[i]>0.0) {
3654              value=columnLowerWork_[i];
3655            } else {
3656              value=0.0;
3657            }
3658          }
3659          columnActivityWork_[i] = value;
3660          reducedCostWork_[i] = reducedCost_[i];
3661        }
3662        for (i=0;i<numberRows_;i++) {
3663          double value = rowActivity_[i];
3664          if (fabs(value)>1.0e20) {
3665            //printf("bad value of %g for row %d\n",value,i);
3666            setRowStatus(i,superBasic);
3667            if (rowUpperWork_[i]<0.0) {
3668              value=rowUpperWork_[i];
3669            } else if (rowLowerWork_[i]>0.0) {
3670              value=rowLowerWork_[i];
3671            } else {
3672              value=0.0;
3673            }
3674          }
3675          rowActivityWork_[i] = value;
3676          rowReducedCost_[i] = dual_[i];
3677        }
3678      }
3679    }
3680  }
3681 
3682  if (what==63&&doSanityCheck&&!auxiliaryModel_) {
3683    // check rim of problem okay
3684    if (!sanityCheck())
3685      goodMatrix=false;
3686  } 
3687  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
3688  // maybe we need to move scales to SimplexModel for factorization?
3689  if ((what==63&&!pivotVariable_)||(newArrays&&!keepPivots)) {
3690    delete [] pivotVariable_;
3691    pivotVariable_=new int[numberRows2];
3692    for (int i=0;i<numberRows2;i++)
3693      pivotVariable_[i]=-1;
3694  } else if (what==63&&!keepPivots) {
3695    // just reset
3696    for (int i=0;i<numberRows2;i++)
3697      pivotVariable_[i]=-1;
3698  } else if (what==63) {
3699    // check pivots
3700    for (int i=0;i<numberRows2;i++) {
3701      int iSequence = pivotVariable_[i];
3702      if (iSequence<numberRows_+numberColumns_&&
3703          getStatus(iSequence)!=basic) {
3704        keepPivots =false;
3705        break;
3706      }
3707    }
3708    if (!keepPivots) {
3709      // reset
3710      for (int i=0;i<numberRows2;i++)
3711        pivotVariable_[i]=-1;
3712    } else {
3713      // clean
3714      for (int i=0;i<numberColumns_+numberRows_;i++) {
3715        Status status =getStatus(i);
3716        if (status!=basic) {
3717          if (upper_[i]==lower_[i]) {
3718            setStatus(i,isFixed);
3719            solution_[i]=lower_[i];
3720          } else if (status==atLowerBound) {
3721            if (lower_[i]>-1.0e20) {
3722              solution_[i]=lower_[i];
3723            } else {
3724              //printf("seq %d at lower of %g\n",i,lower_[i]);
3725              if (upper_[i]<1.0e20) {
3726                solution_[i]=upper_[i];
3727                setStatus(i,atUpperBound);
3728              } else {
3729                setStatus(i,isFree);
3730              }
3731            }
3732          } else if (status==atUpperBound) {
3733            if (upper_[i]<1.0e20) {
3734              solution_[i]=upper_[i];
3735            } else {
3736              //printf("seq %d at upper of %g\n",i,upper_[i]);
3737              if (lower_[i]>-1.0e20) {
3738                solution_[i]=lower_[i];
3739                setStatus(i,atLowerBound);
3740              } else {
3741                setStatus(i,isFree);
3742              }
3743            }
3744          } else if (status==isFixed&&upper_[i]>lower_[i]) {
3745            // was fixed - not now
3746            if (solution_[i]<=lower_[i]) {
3747              setStatus(i,atLowerBound);
3748            } else if (solution_[i]>=upper_[i]) {
3749              setStatus(i,atUpperBound);
3750            } else {
3751              setStatus(i,superBasic);
3752            }
3753          }
3754        }
3755      }
3756    }
3757  }
3758 
3759  if (what==63) {
3760    if (newArrays) {
3761      // get some arrays
3762      int iRow,iColumn;
3763      // these are "indexed" arrays so we always know where nonzeros are
3764      /**********************************************************
3765      rowArray_[3] is long enough for rows+columns
3766      *********************************************************/
3767      for (iRow=0;iRow<4;iRow++) {
3768        int length =numberRows2+factorization_->maximumPivots();
3769        if (iRow==3||objective_->type()>1)
3770          length += numberColumns_;
3771        if ((specialOptions_&65536)==0||!rowArray_[iRow]) {
3772          delete rowArray_[iRow];
3773          rowArray_[iRow]=new CoinIndexedVector();
3774        }
3775        rowArray_[iRow]->reserve(length);
3776      }
3777     
3778      for (iColumn=0;iColumn<2;iColumn++) {
3779        if ((specialOptions_&65536)==0||!columnArray_[iColumn]) {
3780          delete columnArray_[iColumn];
3781          columnArray_[iColumn]=new CoinIndexedVector();
3782        }
3783        if (!iColumn)
3784          columnArray_[iColumn]->reserve(numberColumns_);
3785        else
3786          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
3787      }
3788    } else {
3789      int iRow,iColumn;
3790      if (auxiliaryModel_) {
3791        for (iRow=0;iRow<4;iRow++) {
3792          assert (!rowArray_[iRow]);
3793          rowArray_[iRow]=auxiliaryModel_->rowArray_[iRow];
3794          // For safety
3795          memset(rowArray_[iRow]->denseVector(),0,rowArray_[iRow]->capacity()*sizeof(double));
3796        }
3797        for (iColumn=0;iColumn<2;iColumn++) {
3798          assert (!columnArray_[iColumn]);
3799          columnArray_[iColumn]=auxiliaryModel_->columnArray_[iColumn];
3800          // For safety
3801          memset(columnArray_[iColumn]->denseVector(),0,columnArray_[iColumn]->capacity()*sizeof(double));
3802        }
3803        // do matrices
3804        rowCopy_=auxiliaryModel_->rowCopy_;
3805        ClpMatrixBase * temp = matrix_;
3806        matrix_=auxiliaryModel_->matrix_;
3807        auxiliaryModel_->matrix_=temp;
3808      }
3809      for (iRow=0;iRow<4;iRow++) {
3810        rowArray_[iRow]->clear();
3811#ifndef NDEBUG
3812        int length =numberRows2+factorization_->maximumPivots();
3813        if (iRow==3||objective_->type()>1)
3814          length += numberColumns_;
3815        assert(rowArray_[iRow]->capacity()>=length);
3816        rowArray_[iRow]->checkClear();
3817#endif
3818      }
3819     
3820      for (iColumn=0;iColumn<2;iColumn++) {
3821        columnArray_[iColumn]->clear();
3822#ifndef NDEBUG
3823        int length =numberColumns_;
3824        if (iColumn)
3825          length=CoinMax(numberRows2,numberColumns_);
3826        assert(columnArray_[iColumn]->capacity()>=length);
3827        columnArray_[iColumn]->checkClear();
3828#endif
3829      }
3830    }   
3831  }
3832  if (problemStatus_==10) {
3833    problemStatus_=-1;
3834    handler_->setLogLevel(saveLevel); // switch back messages
3835  }
3836  if ((what&5)!=0) 
3837    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
3838  if (goodMatrix&&(specialOptions_&65536)!=0&&!auxiliaryModel_) {
3839    int save = maximumColumns_+maximumRows_;
3840    CoinMemcpyN(cost_,numberTotal,cost_+save);
3841    CoinMemcpyN(lower_,numberTotal,lower_+save);
3842    CoinMemcpyN(upper_,numberTotal,upper_+save);
3843    CoinMemcpyN(dj_,numberTotal,dj_+save);
3844    CoinMemcpyN(solution_,numberTotal,solution_+save);
3845    if (rowScale_&&!savedRowScale_) {
3846      double * temp;
3847      temp = new double [4*maximumRows_];
3848      CoinFillN(temp,4*maximumRows_,1.0);
3849      CoinMemcpyN(rowScale_,numberRows2,temp);
3850      CoinMemcpyN(rowScale_+numberRows2,numberRows2,temp+maximumRows_);
3851      CoinMemcpyN(rowScale_,numberRows2,temp+2*maximumRows_);
3852      CoinMemcpyN(rowScale_+numberRows2,numberRows2,temp+3*maximumRows_);
3853      delete [] rowScale_;
3854      savedRowScale_ = temp;
3855      rowScale_ = savedRowScale_;
3856      inverseRowScale_ = savedRowScale_+maximumInternalRows_;
3857      temp = new double [4*maximumColumns_];
3858      CoinFillN(temp,4*maximumColumns_,1.0);
3859      CoinMemcpyN(columnScale_,numberColumns_,temp);
3860      CoinMemcpyN(columnScale_+numberColumns_,numberColumns_,temp+maximumColumns_);
3861      CoinMemcpyN(columnScale_,numberColumns_,temp+2*maximumColumns_);
3862      CoinMemcpyN(columnScale_+numberColumns_,numberColumns_,temp+3*maximumColumns_);
3863      delete [] columnScale_;
3864      savedColumnScale_ = temp;
3865      columnScale_ = savedColumnScale_;
3866      inverseColumnScale_ = savedColumnScale_+maximumInternalColumns_;
3867    }
3868  }
3869  return goodMatrix;
3870}
3871// Does rows and columns
3872void
3873ClpSimplex::createRim1(bool initial)
3874{
3875  int i;
3876  int numberRows2 = numberRows_+numberExtraRows_;
3877  int numberTotal = numberRows2+numberColumns_;
3878  //assert (!auxiliaryModel_);
3879  if (!auxiliaryModel_) {
3880    if ((specialOptions_&65536)!=0&&true) {
3881      assert (!initial);
3882      int save = maximumColumns_+maximumRows_;
3883      CoinMemcpyN(lower_+save,numberTotal,lower_);
3884      CoinMemcpyN(upper_+save,numberTotal,upper_);
3885      return;
3886    }
3887    const double * rowScale = rowScale_;
3888    const double * columnScale = columnScale_;
3889    // clean up any mismatches on infinity
3890    // and fix any variables with tiny gaps
3891    double primalTolerance=dblParam_[ClpPrimalTolerance];
3892    if(rowScale) {
3893      // If scaled then do all columns later in one loop
3894      if (!initial) {
3895        if (!inverseColumnScale_) {
3896          for (i=0;i<numberColumns_;i++) {
3897            double multiplier = rhsScale_/columnScale[i];
3898            double lowerValue = columnLower_[i];
3899            double upperValue = columnUpper_[i];
3900            if (lowerValue>-1.0e20) {
3901              columnLowerWork_[i]=lowerValue*multiplier;
3902              if (upperValue>=1.0e20) {
3903                columnUpperWork_[i]=COIN_DBL_MAX;
3904              } else {
3905                columnUpperWork_[i]=upperValue*multiplier;
3906                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3907                  if (columnLowerWork_[i]>=0.0) {
3908                    columnUpperWork_[i] = columnLowerWork_[i];
3909                  } else if (columnUpperWork_[i]<=0.0) {
3910                    columnLowerWork_[i] = columnUpperWork_[i];
3911                  } else {
3912                    columnUpperWork_[i] = 0.0;
3913                    columnLowerWork_[i] = 0.0;
3914                  }
3915                }
3916              }
3917            } else if (upperValue<1.0e20) {
3918              columnLowerWork_[i]=-COIN_DBL_MAX;
3919              columnUpperWork_[i]=upperValue*multiplier;
3920            } else {
3921              // free
3922              columnLowerWork_[i]=-COIN_DBL_MAX;
3923              columnUpperWork_[i]=COIN_DBL_MAX;
3924            }
3925          }
3926        } else {
3927          assert (inverseColumnScale_);
3928          const double * inverseScale = inverseColumnScale_;
3929          for (i=0;i<numberColumns_;i++) {
3930            double multiplier = rhsScale_*inverseScale[i];
3931            double lowerValue = columnLower_[i];
3932            double upperValue = columnUpper_[i];
3933            if (lowerValue>-1.0e20) {
3934              columnLowerWork_[i]=lowerValue*multiplier;
3935              if (upperValue>=1.0e20) {
3936                columnUpperWork_[i]=COIN_DBL_MAX;
3937              } else {
3938                columnUpperWork_[i]=upperValue*multiplier;
3939                if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
3940                  if (columnLowerWork_[i]>=0.0) {
3941                    columnUpperWork_[i] = columnLowerWork_[i];
3942                  } else if (columnUpperWork_[i]<=0.0) {
3943                    columnLowerWork_[i] = columnUpperWork_[i];
3944                  } else {
3945                    columnUpperWork_[i] = 0.0;
3946                    columnLowerWork_[i] = 0.0;
3947                  }
3948                }
3949              }
3950            } else if (upperValue<1.0e20) {
3951              columnLowerWork_[i]=-COIN_DBL_MAX;
3952              columnUpperWork_[i]=upperValue*multiplier;
3953            } else {
3954              // free
3955              columnLowerWork_[i]=-COIN_DBL_MAX;
3956              columnUpperWork_[i]=COIN_DBL_MAX;
3957            }
3958          }
3959        }
3960      }
3961      for (i=0;i<numberRows_;i++) {
3962        double multiplier = rhsScale_*rowScale[i];
3963        double lowerValue = rowLower_[i];
3964        double upperValue = rowUpper_[i];
3965        if (lowerValue>-1.0e20) {
3966          rowLowerWork_[i]=lowerValue*multiplier;
3967          if (upperValue>=1.0e20) {
3968            rowUpperWork_[i]=COIN_DBL_MAX;
3969          } else {
3970            rowUpperWork_[i]=upperValue*multiplier;
3971            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
3972              if (rowLowerWork_[i]>=0.0) {
3973                rowUpperWork_[i] = rowLowerWork_[i];
3974              } else if (rowUpperWork_[i]<=0.0) {
3975                rowLowerWork_[i] = rowUpperWork_[i];
3976              } else {
3977                rowUpperWork_[i] = 0.0;
3978                rowLowerWork_[i] = 0.0;
3979              }
3980            }
3981          }
3982        } else if (upperValue<1.0e20) {
3983          rowLowerWork_[i]=-COIN_DBL_MAX;
3984          rowUpperWork_[i]=upperValue*multiplier;
3985        } else {
3986          // free
3987          rowLowerWork_[i]=-COIN_DBL_MAX;
3988          rowUpperWork_[i]=COIN_DBL_MAX;
3989        }
3990      }
3991    } else if (rhsScale_!=1.0) {
3992      for (i=0;i<numberColumns_;i++) {
3993        double lowerValue = columnLower_[i];
3994        double upperValue = columnUpper_[i];
3995        if (lowerValue>-1.0e20) {
3996          columnLowerWork_[i]=lowerValue*rhsScale_;
3997          if (upperValue>=1.0e20) {
3998            columnUpperWork_[i]=COIN_DBL_MAX;
3999          } else {
4000            columnUpperWork_[i]=upperValue*rhsScale_;
4001            if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
4002              if (columnLowerWork_[i]>=0.0) {
4003                columnUpperWork_[i] = columnLowerWork_[i];
4004              } else if (columnUpperWork_[i]<=0.0) {
4005                columnLowerWork_[i] = columnUpperWork_[i];
4006              } else {
4007                columnUpperWork_[i] = 0.0;
4008                columnLowerWork_[i] = 0.0;
4009              }
4010            }
4011          }
4012        } else if (upperValue<1.0e20) {
4013          columnLowerWork_[i]=-COIN_DBL_MAX;
4014          columnUpperWork_[i]=upperValue*rhsScale_;
4015        } else {
4016          // free
4017          columnLowerWork_[i]=-COIN_DBL_MAX;
4018          columnUpperWork_[i]=COIN_DBL_MAX;
4019        }
4020      }
4021      for (i=0;i<numberRows_;i++) {
4022        double lowerValue = rowLower_[i];
4023        double upperValue = rowUpper_[i];
4024        if (lowerValue>-1.0e20) {
4025          rowLowerWork_[i]=lowerValue*rhsScale_;
4026          if (upperValue>=1.0e20) {
4027            rowUpperWork_[i]=COIN_DBL_MAX;
4028          } else {
4029            rowUpperWork_[i]=upperValue*rhsScale_;
4030            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
4031              if (rowLowerWork_[i]>=0.0) {
4032                rowUpperWork_[i] = rowLowerWork_[i];
4033              } else if (rowUpperWork_[i]<=0.0) {
4034                rowLowerWork_[i] = rowUpperWork_[i];
4035              } else {
4036                rowUpperWork_[i] = 0.0;
4037                rowLowerWork_[i] = 0.0;
4038              }
4039            }
4040          }
4041        } else if (upperValue<1.0e20) {
4042          rowLowerWork_[i]=-COIN_DBL_MAX;
4043          rowUpperWork_[i]=upperValue*rhsScale_;
4044        } else {
4045          // free
4046          rowLowerWork_[i]=-COIN_DBL_MAX;
4047          rowUpperWork_[i]=COIN_DBL_MAX;
4048        }
4049      }
4050    } else {
4051      for (i=0;i<numberColumns_;i++) {
4052        double lowerValue = columnLower_[i];
4053        double upperValue = columnUpper_[i];
4054        if (lowerValue>-1.0e20) {
4055          columnLowerWork_[i]=lowerValue;
4056          if (upperValue>=1.0e20) {
4057            columnUpperWork_[i]=COIN_DBL_MAX;
4058          } else {
4059            columnUpperWork_[i]=upperValue;
4060            if (fabs(columnUpperWork_[i]-columnLowerWork_[i])<=primalTolerance) {
4061              if (columnLowerWork_[i]>=0.0) {
4062                columnUpperWork_[i] = columnLowerWork_[i];
4063              } else if (columnUpperWork_[i]<=0.0) {
4064                columnLowerWork_[i] = columnUpperWork_[i];
4065              } else {
4066                columnUpperWork_[i] = 0.0;
4067                columnLowerWork_[i] = 0.0;
4068              }
4069            }
4070          }
4071        } else if (upperValue<1.0e20) {
4072          columnLowerWork_[i]=-COIN_DBL_MAX;
4073          columnUpperWork_[i]=upperValue;
4074        } else {
4075          // free
4076          columnLowerWork_[i]=-COIN_DBL_MAX;
4077          columnUpperWork_[i]=COIN_DBL_MAX;
4078        }
4079      }
4080      for (i=0;i<numberRows_;i++) {
4081        double lowerValue = rowLower_[i];
4082        double upperValue = rowUpper_[i];
4083        if (lowerValue>-1.0e20) {
4084          rowLowerWork_[i]=lowerValue;
4085          if (upperValue>=1.0e20) {
4086            rowUpperWork_[i]=COIN_DBL_MAX;
4087          } else {
4088            rowUpperWork_[i]=upperValue;
4089            if (fabs(rowUpperWork_[i]-rowLowerWork_[i])<=primalTolerance) {
4090              if (rowLowerWork_[i]>=0.0) {
4091                rowUpperWork_[i] = rowLowerWork_[i];
4092              } else if (rowUpperWork_[i]<=0.0) {
4093                rowLowerWork_[i] = rowUpperWork_[i];
4094              } else {
4095                rowUpperWork_[i] = 0.0;
4096                rowLowerWork_[i] = 0.0;
4097              }
4098            }
4099          }
4100        } else if (upperValue<1.0e20) {
4101          rowLowerWork_[i]=-COIN_DBL_MAX;
4102          rowUpperWork_[i]=upperValue;
4103        } else {
4104          // free
4105          rowLowerWork_[i]=-COIN_DBL_MAX;
4106          rowUpperWork_[i]=COIN_DBL_MAX;
4107        }
4108      }
4109    }
4110  } else {
4111    // auxiliary model
4112    if (!initial) {
4113      // just copy
4114      CoinMemcpyN(auxiliaryModel_->lower_+numberTotal,numberTotal,lower_);
4115      CoinMemcpyN(auxiliaryModel_->upper_+numberTotal,numberTotal,upper_);
4116    } else {
4117      // clean up any mismatches on infinity
4118      // and fix any variables with tiny gaps
4119      const double * rowScale = auxiliaryModel_->rowScale_ ;
4120      const double * columnScale = auxiliaryModel_->columnScale_ ;
4121      assert (rhsScale_==1.0);
4122      assert (objectiveScale_==1.0);
4123      if ((auxiliaryModel_->specialOptions_&2)==0) {
4124        // need to do column bounds
4125        if (columnScale) {
4126          const double * inverseScale = columnScale+numberColumns_;
4127          // scaled
4128          for (i=0;i<numberColumns_;i++) {
4129            double value;
4130            value = columnLower_[i];
4131            if (value>-1.0e20) 
4132              value *= inverseScale[i];
4133            lower_[i] = value;
4134            value = columnUpper_[i];
4135            if (value<1.0e20) 
4136              value *= inverseScale[i];
4137            upper_[i] = value;
4138          }
4139        } else {
4140          for (i=0;i<numberColumns_;i++) {
4141            lower_[i] = columnLower_[i];
4142            upper_[i] = columnUpper_[i];
4143          }
4144        }
4145        // save
4146 CoinMemcpyN(lower_,numberColumns_,auxiliaryModel_->lower_+numberTotal);
4147 CoinMemcpyN(upper_,numberColumns_,auxiliaryModel_->upper_+numberTotal);
4148      } else {
4149        // just copy
4150 CoinMemcpyN(auxiliaryModel_->lower_+numberTotal,numberColumns_,lower_);
4151 CoinMemcpyN(auxiliaryModel_->upper_+numberTotal,numberColumns_,upper_);
4152      }
4153      if ((auxiliaryModel_->specialOptions_&1)==0) {
4154        // need to do row bounds
4155        if (rowScale) {
4156          // scaled
4157          for (i=0;i<numberRows_;i++) {
4158            double value;
4159            value = rowLower_[i];
4160            if (value>-1.0e20) 
4161              value *= rowScale[i];
4162            lower_[i+numberColumns_] = value;
4163            value = rowUpper_[i];
4164            if (value<1.0e20) 
4165              value *= rowScale[i];
4166            upper_[i+numberColumns_] = value;
4167          }
4168        } else {
4169          for (i=0;i<numberRows_;i++) {
4170            lower_[i+numberColumns_] = rowLower_[i];
4171            upper_[i+numberColumns_] = rowUpper_[i];
4172          }
4173        }
4174        // save
4175 CoinMemcpyN(          lower_+numberColumns_,numberRows_,auxiliaryModel_->lower_+numberTotal+numberColumns_);
4176 CoinMemcpyN(          upper_+numberColumns_,numberRows_,auxiliaryModel_->upper_+numberTotal+numberColumns_);
4177      } else {
4178        // just copy
4179 CoinMemcpyN(          auxiliaryModel_->lower_+numberTotal+numberColumns_,numberRows_,
4180              lower_+numberColumns_);
4181 CoinMemcpyN(          auxiliaryModel_->upper_+numberTotal+numberColumns_,numberRows_,
4182              upper_+numberColumns_);
4183      }
4184    }
4185    if (initial) {
4186      // do basis
4187      assert ((auxiliaryModel_->specialOptions_&8)!=0);
4188      // clean solution trusting basis
4189      for ( i=0;i<numberTotal;i++) {
4190        Status status =getStatus(i);
4191        double value=solution_[i]; // may as well leave if basic (values pass)
4192        if (status!=basic) {
4193          if (upper_[i]==lower_[i]) {
4194            setStatus(i,isFixed);
4195            value=lower_[i];
4196          } else if (status==atLowerBound) {
4197            if (lower_[i]>-1.0e20) {
4198              value=lower_[i];
4199            } else {
4200              if (upper_[i]<1.0e20) {
4201                value=upper_[i];
4202                setStatus(i,atUpperBound);
4203              } else {
4204                setStatus(i,isFree);
4205              }
4206            }
4207          } else if (status==atUpperBound) {
4208            if (upper_[i]<1.0e20) {
4209              value=upper_[i];
4210            } else {
4211              if (lower_[i]>-1.0e20) {
4212                value=lower_[i];
4213                setStatus(i,atLowerBound);
4214              } else {
4215                setStatus(i,isFree);
4216              }
4217            }
4218          }
4219        }
4220        solution_[i]=value;
4221      }
4222    }
4223  }
4224#ifndef NDEBUG
4225  if ((specialOptions_&65536)!=0&&false) {
4226    assert (!initial);
4227    int save = maximumColumns_+maximumRows_;
4228    for (int i=0;i<numberTotal;i++) {
4229      assert (fabs(lower_[i]-lower_[i+save])<1.0e-5);
4230      assert (fabs(upper_[i]-upper_[i+save])<1.0e-5);
4231    }
4232  }
4233#endif
4234}
4235// Does objective
4236void
4237ClpSimplex::createRim4(bool initial)
4238{
4239  int i;
4240  int numberRows2 = numberRows_+numberExtraRows_;
4241  int numberTotal = numberRows2+numberColumns_;
4242  //assert (!auxiliaryModel_);
4243  if (!auxiliaryModel_||(initial&&(auxiliaryModel_->specialOptions_&4)==0)) {
4244    if ((specialOptions_&65536)!=0&&true) {
4245      assert (!initial);
4246      int save = maximumColumns_+maximumRows_;
4247      CoinMemcpyN(cost_+save,numberTotal,cost_);
4248      return;
4249    }
4250    double direction = optimizationDirection_*objectiveScale_;
4251    const double * obj = objective();
4252    const double * rowScale = auxiliaryModel_ ? auxiliaryModel_->rowScale_ : rowScale_;
4253    const double * columnScale = auxiliaryModel_ ? auxiliaryModel_->columnScale_ : columnScale_;
4254    // and also scale by scale factors
4255    if (rowScale) {
4256      if (rowObjective_) {
4257        for (i=0;i<numberRows_;i++)
4258          rowObjectiveWork_[i] = rowObjective_[i]*direction/rowScale[i];
4259      } else {
4260        memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
4261      }
4262      // If scaled then do all columns later in one loop
4263      if (!initial||auxiliaryModel_) {
4264        for (i=0;i<numberColumns_;i++) {
4265          CoinAssert(fabs(obj[i])<1.0e25);
4266          objectiveWork_[i] = obj[i]*direction*columnScale[i];
4267        }
4268      }
4269    } else {
4270      if (rowObjective_) {
4271        for (i=0;i<numberRows_;i++)
4272          rowObjectiveWork_[i] = rowObjective_[i]*direction;
4273      } else {
4274        memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
4275      }
4276      for (i=0;i<numberColumns_;i++) {
4277        CoinAssert(fabs(obj[i])<1.0e25);
4278        objectiveWork_[i] = obj[i]*direction;
4279      }
4280    }
4281    if (auxiliaryModel_) {
4282      // save costs
4283      CoinMemcpyN(cost_,numberTotal,auxiliaryModel_->cost_+numberTotal);
4284    }
4285  } else {
4286    // just copy
4287    CoinMemcpyN(auxiliaryModel_->cost_+numberTotal,numberTotal,cost_);
4288  }
4289}
4290// Does rows and columns and objective
4291void
4292ClpSimplex::createRim5(bool initial)
4293{
4294  createRim4(initial);
4295  createRim1(initial);
4296}
4297void
4298ClpSimplex::deleteRim(int getRidOfFactorizationData)
4299{
4300  // Just possible empty problem
4301  int numberRows=numberRows_;
4302  int numberColumns=numberColumns_;
4303  if (!numberRows||!numberColumns) {
4304    numberRows=0;
4305    if (objective_->type()<2)
4306      numberColumns=0;
4307  }
4308  if (!auxiliaryModel_) {
4309    int i;
4310    if (problemStatus_!=1&&problemStatus_!=2) {
4311      delete [] ray_;
4312      ray_=NULL;
4313    }
4314    // set upperOut_ to furthest away from bound so can use in dual for dualBound_
4315    upperOut_=1.0;
4316#if 0
4317    {
4318      int nBad=0;
4319      for (i=0;i<numberColumns;i++) {
4320        if (lower_[i]==upper_[i]&&getColumnStatus(i)==basic)
4321          nBad++;
4322      }
4323      if (nBad)
4324        printf("yy %d basic fixed\n",nBad);
4325    }
4326#endif
4327    // ray may be null if in branch and bound
4328    if (rowScale_) {
4329      // Collect infeasibilities
4330      int numberPrimalScaled=0;
4331      int numberPrimalUnscaled=0;
4332      int numberDualScaled=0;
4333      int numberDualUnscaled=0;
4334      double scaleC = 1.0/objectiveScale_;
4335      double scaleR = 1.0/rhsScale_;
4336      if (!inverseColumnScale_) {
4337        for (i=0;i<numberColumns;i++) {
4338          double scaleFactor = columnScale_[i];
4339          double valueScaled = columnActivityWork_[i];
4340          double lowerScaled = columnLowerWork_[i];
4341          double upperScaled = columnUpperWork_[i];
4342          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4343            if (valueScaled<lowerScaled-primalTolerance_||
4344                valueScaled>upperScaled+primalTolerance_)
4345              numberPrimalScaled++;
4346            else
4347              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4348          }
4349          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
4350          double value = columnActivity_[i];
4351          if (value<columnLower_[i]-primalTolerance_)
4352            numberPrimalUnscaled++;
4353          else if (value>columnUpper_[i]+primalTolerance_)
4354            numberPrimalUnscaled++;
4355          double valueScaledDual = reducedCostWork_[i];
4356          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4357            numberDualScaled++;
4358          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4359            numberDualScaled++;
4360          reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
4361          double valueDual = reducedCost_[i];
4362          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4363            numberDualUnscaled++;
4364          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4365            numberDualUnscaled++;
4366        }
4367        for (i=0;i<numberRows;i++) {
4368          double scaleFactor = rowScale_[i];
4369          double valueScaled = rowActivityWork_[i];
4370          double lowerScaled = rowLowerWork_[i];
4371          double upperScaled = rowUpperWork_[i];
4372          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4373            if (valueScaled<lowerScaled-primalTolerance_||
4374                valueScaled>upperScaled+primalTolerance_)
4375              numberPrimalScaled++;
4376            else
4377              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4378          }
4379          rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
4380          double value = rowActivity_[i];
4381          if (value<rowLower_[i]-primalTolerance_)
4382            numberPrimalUnscaled++;
4383          else if (value>rowUpper_[i]+primalTolerance_)
4384            numberPrimalUnscaled++;
4385          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
4386          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4387            numberDualScaled++;
4388          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4389            numberDualScaled++;
4390          dual_[i] *= scaleFactor*scaleC;
4391          double valueDual = dual_[i]; 
4392          if (rowObjective_)
4393            valueDual += rowObjective_[i];
4394          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4395            numberDualUnscaled++;
4396          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4397            numberDualUnscaled++;
4398        }
4399      } else {
4400        assert (inverseColumnScale_);
4401        const double * inverseScale = inverseColumnScale_;
4402        for (i=0;i<numberColumns;i++) {
4403          double scaleFactor = columnScale_[i];
4404          double valueScaled = columnActivityWork_[i];
4405          double lowerScaled = columnLowerWork_[i];
4406          double upperScaled = columnUpperWork_[i];
4407          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4408            if (valueScaled<lowerScaled-primalTolerance_||
4409                valueScaled>upperScaled+primalTolerance_)
4410              numberPrimalScaled++;
4411            else
4412              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4413          }
4414          columnActivity_[i] = valueScaled*scaleFactor*scaleR;
4415          double value = columnActivity_[i];
4416          if (value<columnLower_[i]-primalTolerance_)
4417            numberPrimalUnscaled++;
4418          else if (value>columnUpper_[i]+primalTolerance_)
4419            numberPrimalUnscaled++;
4420          double valueScaledDual = reducedCostWork_[i];
4421          if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4422            numberDualScaled++;
4423          if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4424            numberDualScaled++;
4425          reducedCost_[i] = (valueScaledDual*scaleC)*inverseScale[i];
4426          double valueDual = reducedCost_[i];
4427          if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4428            numberDualUnscaled++;
4429          if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4430            numberDualUnscaled++;
4431        }
4432        inverseScale = inverseRowScale_;
4433        for (i=0;i<numberRows;i++) {
4434          double scaleFactor = rowScale_[i];
4435          double valueScaled = rowActivityWork_[i];
4436          double lowerScaled = rowLowerWork_[i];
4437          double upperScaled = rowUpperWork_[i];
4438          if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4439            if (valueScaled<lowerScaled-primalTolerance_||
4440                valueScaled>upperScaled+primalTolerance_)
4441              numberPrimalScaled++;
4442            else
4443              upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4444          }
4445          rowActivity_[i] = (valueScaled*scaleR)*inverseScale[i];
4446          double value = rowActivity_[i];
4447          if (value<rowLower_[i]-primalTolerance_)
4448            numberPrimalUnscaled++;
4449          else if (value>rowUpper_[i]+primalTolerance_)
4450            numberPrimalUnscaled++;
4451          double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
4452          if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4453            numberDualScaled++;
4454          if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4455            numberDualScaled++;
4456          dual_[i] *= scaleFactor*scaleC;
4457          double valueDual = dual_[i]; 
4458          if (rowObjective_)
4459            valueDual += rowObjective_[i];
4460          if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4461            numberDualUnscaled++;
4462          if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4463            numberDualUnscaled++;
4464        }
4465      }
4466      if (!problemStatus_&&!secondaryStatus_) {
4467        // See if we need to set secondary status
4468        if (numberPrimalUnscaled) {
4469          if (numberDualUnscaled) 
4470            secondaryStatus_=4;
4471          else
4472            secondaryStatus_=2;
4473        } else {
4474          if (numberDualUnscaled) 
4475            secondaryStatus_=3;
4476        }
4477      }
4478      if (problemStatus_==2) {
4479        for (i=0;i<numberColumns;i++) {
4480          ray_[i] *= columnScale_[i];
4481        }
4482      } else if (problemStatus_==1&&ray_) {
4483        for (i=0;i<numberRows;i++) {
4484          ray_[i] *= rowScale_[i];
4485        }
4486      }
4487    } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
4488      // Collect infeasibilities
4489      int numberPrimalScaled=0;
4490      int numberPrimalUnscaled=0;
4491      int numberDualScaled=0;
4492      int numberDualUnscaled=0;
4493      double scaleC = 1.0/objectiveScale_;
4494      double scaleR = 1.0/rhsScale_;
4495      for (i=0;i<numberColumns;i++) {
4496        double valueScaled = columnActivityWork_[i];
4497        double lowerScaled = columnLowerWork_[i];
4498        double upperScaled = columnUpperWork_[i];
4499        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4500          if (valueScaled<lowerScaled-primalTolerance_||
4501              valueScaled>upperScaled+primalTolerance_)
4502            numberPrimalScaled++;
4503          else
4504            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4505        }
4506        columnActivity_[i] = valueScaled*scaleR;
4507        double value = columnActivity_[i];
4508        if (value<columnLower_[i]-primalTolerance_)
4509          numberPrimalUnscaled++;
4510        else if (value>columnUpper_[i]+primalTolerance_)
4511          numberPrimalUnscaled++;
4512        double valueScaledDual = reducedCostWork_[i];
4513        if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4514          numberDualScaled++;
4515        if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4516          numberDualScaled++;
4517        reducedCost_[i] = valueScaledDual*scaleC;
4518        double valueDual = reducedCost_[i];
4519        if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4520          numberDualUnscaled++;
4521        if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4522          numberDualUnscaled++;
4523      }
4524      for (i=0;i<numberRows;i++) {
4525        double valueScaled = rowActivityWork_[i];
4526        double lowerScaled = rowLowerWork_[i];
4527        double upperScaled = rowUpperWork_[i];
4528        if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
4529          if (valueScaled<lowerScaled-primalTolerance_||
4530              valueScaled>upperScaled+primalTolerance_)
4531            numberPrimalScaled++;
4532          else
4533            upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
4534        }
4535        rowActivity_[i] = valueScaled*scaleR;
4536        double value = rowActivity_[i];
4537        if (value<rowLower_[i]-primalTolerance_)
4538          numberPrimalUnscaled++;
4539        else if (value>rowUpper_[i]+primalTolerance_)
4540          numberPrimalUnscaled++;
4541        double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
4542        if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
4543          numberDualScaled++;
4544        if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
4545          numberDualScaled++;
4546        dual_[i] *= scaleC;
4547        double valueDual = dual_[i]; 
4548        if (rowObjective_)
4549          valueDual += rowObjective_[i];
4550        if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
4551          numberDualUnscaled++;
4552        if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
4553          numberDualUnscaled++;
4554      }
4555      if (!problemStatus_&&!secondaryStatus_) {
4556        // See if we need to set secondary status
4557        if (numberPrimalUnscaled) {
4558          if (numberDualUnscaled) 
4559            secondaryStatus_=4;
4560          else
4561            secondaryStatus_=2;
4562        } else {
4563          if (numberDualUnscaled) 
4564            secondaryStatus_=3;
4565        }
4566      }
4567    } else {
4568      if (columnActivityWork_) {
4569        for (i=0;i<numberColumns;i++) {
4570          double value = columnActivityWork_[i];
4571          double lower = columnLowerWork_[i];
4572          double upper = columnUpperWork_[i];
4573          if (lower>-1.0e20||upper<1.0e20) {
4574            if (value>lower&&value<upper)
4575              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4576          }
4577          columnActivity_[i] = columnActivityWork_[i];
4578          reducedCost_[i] = reducedCostWork_[i];
4579        }
4580        for (i=0;i<numberRows;i++) {
4581          double value = rowActivityWork_[i];
4582          double lower = rowLowerWork_[i];
4583          double upper = rowUpperWork_[i];
4584          if (lower>-1.0e20||upper<1.0e20) {
4585            if (value>lower&&value<upper)
4586              upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4587          }
4588          rowActivity_[i] = rowActivityWork_[i];
4589        }
4590      }
4591    }
4592    // switch off scalefactor if auto
4593    if (automaticScale_) {
4594      rhsScale_=1.0;
4595      objectiveScale_=1.0;
4596    }
4597    if (optimizationDirection_!=1.0) {
4598      // and modify all dual signs
4599      for (i=0;i<numberColumns;i++) 
4600        reducedCost_[i] *= optimizationDirection_;
4601      for (i=0;i<numberRows;i++) 
4602        dual_[i] *= optimizationDirection_;
4603    }
4604  } else {
4605    // auxiliary model
4606    cost_=NULL;
4607    lower_=NULL;
4608    upper_=NULL;
4609    dj_=NULL;
4610    solution_=NULL;
4611    int iRow,iColumn;
4612    for (iRow=0;iRow<4;iRow++) {
4613      if (rowArray_[iRow]) {
4614        rowArray_[iRow]->clear();
4615        //rowArray_[iRow]->checkClear();
4616      }
4617      rowArray_[iRow]=NULL;
4618    }
4619    for (iColumn=0;iColumn<2;iColumn++) {
4620      if (columnArray_[iColumn]) {
4621        columnArray_[iColumn]->clear();
4622        //columnArray_[iColumn]->checkClear();
4623      }
4624      columnArray_[iColumn]=NULL;
4625    }
4626    rowCopy_=NULL;
4627    ClpMatrixBase * temp = matrix_;
4628    matrix_=auxiliaryModel_->matrix_;
4629    auxiliaryModel_->matrix_=temp;
4630    assert ((auxiliaryModel_->specialOptions_&(16+32))==16+32);
4631    if (problemStatus_!=1&&problemStatus_!=10) {
4632      int i;
4633      if (auxiliaryModel_->rowScale_) {
4634        const double * scale = auxiliaryModel_->columnScale_;
4635        const double * inverseScale = scale + numberColumns;
4636        for (i=0;i<numberColumns;i++) {
4637          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i]*scale[i];
4638          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i]*inverseScale[i];
4639        }
4640        scale = auxiliaryModel_->rowScale_;
4641        inverseScale = scale + numberRows;
4642        for (i=0;i<numberRows;i++) {
4643          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
4644        }
4645      } else {
4646        for (i=0;i<numberColumns;i++) {
4647          columnActivity_[i] = auxiliaryModel_->columnActivityWork_[i];
4648          reducedCost_[i] = auxiliaryModel_->reducedCostWork_[i];
4649        }
4650        for (i=0;i<numberRows;i++) {
4651          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i];
4652        }
4653      }
4654      if (optimizationDirection_!=1.0) {
4655        // and modify reduced costs
4656        for (i=0;i<numberColumns;i++) 
4657          reducedCost_[i] *= optimizationDirection_;
4658      }
4659    } else if (problemStatus_==10) {
4660      int i;
4661      if (auxiliaryModel_->rowScale_) {
4662        const double * scale = auxiliaryModel_->columnScale_;
4663        const double * inverseScale = scale + numberColumns;
4664        for (i=0;i<numberColumns;i++) {
4665          double lower = auxiliaryModel_->columnLowerWork_[i];
4666          double upper = auxiliaryModel_->columnUpperWork_[i];
4667          double value = auxiliaryModel_->columnActivityWork_[i];
4668          if (value>lower&&value<upper) {
4669            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4670          }
4671          columnActivity_[i] = value*scale[i];
4672        }
4673        scale = auxiliaryModel_->rowScale_;
4674        inverseScale = scale + numberRows;
4675        for (i=0;i<numberRows;i++) {
4676          double lower = auxiliaryModel_->rowLowerWork_[i];
4677          double upper = auxiliaryModel_->rowUpperWork_[i];
4678          double value = auxiliaryModel_->rowActivityWork_[i];
4679          if (value>lower&&value<upper) {
4680            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4681          }
4682          rowActivity_[i] = auxiliaryModel_->rowActivityWork_[i]*inverseScale[i];
4683        }
4684      } else {
4685        for (i=0;i<numberColumns;i++) {
4686          double lower = auxiliaryModel_->columnLowerWork_[i];
4687          double upper = auxiliaryModel_->columnUpperWork_[i];
4688          double value = auxiliaryModel_->columnActivityWork_[i];
4689          if (value>lower&&value<upper) {
4690            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4691          }
4692          columnActivity_[i] = value;
4693        }
4694        for (i=0;i<numberRows;i++) {
4695          double lower = auxiliaryModel_->rowLowerWork_[i];
4696          double upper = auxiliaryModel_->rowUpperWork_[i];
4697          double value = auxiliaryModel_->rowActivityWork_[i];
4698          if (value>lower&&value<upper) {
4699            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
4700          }
4701          rowActivity_[i] = value;
4702        }
4703      }
4704    }
4705  }
4706  // scaling may have been turned off
4707  scalingFlag_ = abs(scalingFlag_);
4708  if(getRidOfFactorizationData>0) {
4709    gutsOfDelete(getRidOfFactorizationData+1);
4710  } else {
4711    // at least get rid of nonLinearCost_
4712    delete nonLinearCost_;
4713    nonLinearCost_=NULL;
4714  }
4715  if (!rowObjective_&&problemStatus_==0&&objective_->type()==1&&
4716      numberRows&&numberColumns) {
4717  // Redo objective value
4718    double objectiveValue =0.0;
4719    const double * cost = objective();
4720    for (int i=0;i<numberColumns;i++) {
4721      double value = columnActivity_[i];
4722      objectiveValue += value*cost[i];
4723    }
4724    //if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4725    //printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4726    objectiveValue_ = objectiveValue*optimizationDirection();
4727  }
4728  // get rid of data
4729  matrix_->generalExpanded(this,13,scalingFlag_);
4730}
4731void 
4732ClpSimplex::setDualBound(double value)
4733{
4734  if (value>0.0)
4735    dualBound_=value;
4736}
4737void 
4738ClpSimplex::setInfeasibilityCost(double value)
4739{
4740  if (value>0.0)
4741    infeasibilityCost_=value;
4742}
4743void ClpSimplex::setNumberRefinements( int value) 
4744{
4745  if (value>=0&&value<10)
4746    numberRefinements_=value;
4747}
4748// Sets row pivot choice algorithm in dual
4749void 
4750ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
4751{
4752  delete dualRowPivot_;
4753  dualRowPivot_ = choice.clone(true);
4754}
4755// Sets row pivot choice algorithm in dual
4756void 
4757ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
4758{
4759  delete primalColumnPivot_;
4760  primalColumnPivot_ = choice.clone(true);
4761}
4762// Passes in factorization
4763void 
4764ClpSimplex::setFactorization( ClpFactorization & factorization)
4765{
4766  delete factorization_;
4767  factorization_= new ClpFactorization(factorization);
4768}
4769// Copies in factorization to existing one
4770void 
4771ClpSimplex::copyFactorization( ClpFactorization & factorization)
4772{
4773  *factorization_= factorization;
4774}
4775/* Perturbation:
4776   -50 to +50 - perturb by this power of ten (-6 sounds good)
4777   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4778   101 - we are perturbed
4779   102 - don't try perturbing again
4780   default is 100
4781*/
4782void 
4783ClpSimplex::setPerturbation(int value)
4784{
4785  if(value<=100&&value >=-1000) {
4786    perturbation_=value;
4787  } 
4788}
4789// Sparsity on or off
4790bool 
4791ClpSimplex::sparseFactorization() const
4792{
4793  return factorization_->sparseThreshold()!=0;
4794}
4795void 
4796ClpSimplex::setSparseFactorization(bool value)
4797{
4798  if (value) {
4799    if (!factorization_->sparseThreshold())
4800      factorization_->goSparse();
4801  } else {
4802    factorization_->sparseThreshold(0);
4803  }
4804}
4805void checkCorrect(ClpSimplex * model,int iRow,
4806                  const double * element,const int * rowStart,const int * rowLength,
4807                  const int * column,
4808                  const double * columnLower_, const double * columnUpper_,
4809                  int infiniteUpperC,
4810                  int infiniteLowerC,
4811                  double &maximumUpC,
4812                  double &maximumDownC)
4813{
4814  int infiniteUpper = 0;
4815  int infiniteLower = 0;
4816  double maximumUp = 0.0;
4817  double maximumDown = 0.0;
4818  CoinBigIndex rStart = rowStart[iRow];
4819  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
4820  CoinBigIndex j;
4821  double large=1.0e15;
4822  int iColumn;
4823  // Compute possible lower and upper ranges
4824 
4825  for (j = rStart; j < rEnd; ++j) {
4826    double value=element[j];
4827    iColumn = column[j];
4828    if (value > 0.0) {
4829      if (columnUpper_[iColumn] >= large) {
4830        ++infiniteUpper;
4831      } else {
4832        maximumUp += columnUpper_[iColumn] * value;
4833      }
4834      if (columnLower_[iColumn] <= -large) {
4835        ++infiniteLower;
4836      } else {
4837        maximumDown += columnLower_[iColumn] * value;
4838      }
4839    } else if (value<0.0) {
4840      if (columnUpper_[iColumn] >= large) {
4841        ++infiniteLower;
4842      } else {
4843        maximumDown += columnUpper_[iColumn] * value;
4844      }
4845      if (columnLower_[iColumn] <= -large) {
4846        ++infiniteUpper;
4847      } else {
4848        maximumUp += columnLower_[iColumn] * value;
4849      }
4850    }
4851  }
4852  assert (infiniteLowerC==infiniteLower);
4853  assert (infiniteUpperC==infiniteUpper);
4854  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
4855    printf("row %d comp up %g, true up %g\n",iRow,
4856           maximumUpC,maximumUp);
4857  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
4858    printf("row %d comp down %g, true down %g\n",iRow,
4859           maximumDownC,maximumDown);
4860  maximumUpC=maximumUp;
4861  maximumDownC=maximumDown;
4862}
4863
4864/* Tightens primal bounds to make dual faster.  Unless
4865   fixed, bounds are slightly looser than they could be.
4866   This is to make dual go faster and is probably not needed
4867   with a presolve.  Returns non-zero if problem infeasible
4868
4869   Fudge for branch and bound - put bounds on columns of factor *
4870   largest value (at continuous) - should improve stability
4871   in branch and bound on infeasible branches (0.0 is off)
4872*/
4873int 
4874ClpSimplex::tightenPrimalBounds(double factor,int doTight,bool tightIntegers)
4875{
4876 
4877  // Get a row copy in standard format
4878  CoinPackedMatrix copy;
4879  copy.reverseOrderedCopyOf(*matrix());
4880  // Matrix may have been created so get rid of it
4881  matrix_->releasePackedMatrix();
4882  // get matrix data pointers
4883  const int * column = copy.getIndices();
4884  const CoinBigIndex * rowStart = copy.getVectorStarts();
4885  const int * rowLength = copy.getVectorLengths(); 
4886  const double * element = copy.getElements();
4887  int numberChanged=1,iPass=0;
4888  double large = largeValue(); // treat bounds > this as infinite
4889#ifndef NDEBUG
4890  double large2= 1.0e10*large;
4891#endif
4892  int numberInfeasible=0;
4893  int totalTightened = 0;
4894
4895  double tolerance = primalTolerance();
4896
4897
4898  // Save column bounds
4899  double * saveLower = new double [numberColumns_];
4900  CoinMemcpyN(columnLower_,numberColumns_,saveLower);
4901  double * saveUpper = new double [numberColumns_];
4902  CoinMemcpyN(columnUpper_,numberColumns_,saveUpper);
4903
4904  int iRow, iColumn;
4905
4906  // If wanted - tighten column bounds using solution
4907  if (factor) {
4908    double largest=0.0;
4909    if (factor>0.0) {
4910      assert (factor>1.0);
4911      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4912        if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
4913          largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
4914        }
4915      }
4916      largest *= factor;
4917    } else {
4918      // absolute
4919       largest = - factor;
4920    }
4921    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4922      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
4923        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
4924        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
4925      }
4926    }
4927  }
4928#define MAXPASS 10
4929
4930  // Loop round seeing if we can tighten bounds
4931  // Would be faster to have a stack of possible rows
4932  // and we put altered rows back on stack
4933  int numberCheck=-1;
4934  while(numberChanged>numberCheck) {
4935
4936    numberChanged = 0; // Bounds tightened this pass
4937   
4938    if (iPass==MAXPASS) break;
4939    iPass++;
4940   
4941    for (iRow = 0; iRow < numberRows_; iRow++) {
4942
4943      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
4944
4945        // possible row
4946        int infiniteUpper = 0;
4947        int infiniteLower = 0;
4948        double maximumUp = 0.0;
4949        double maximumDown = 0.0;
4950        double newBound;
4951        CoinBigIndex rStart = rowStart[iRow];
4952        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
4953        CoinBigIndex j;
4954        // Compute possible lower and upper ranges
4955     
4956        for (j = rStart; j < rEnd; ++j) {
4957          double value=element[j];
4958          iColumn = column[j];
4959          if (value > 0.0) {
4960            if (columnUpper_[iColumn] >= large) {
4961              ++infiniteUpper;
4962            } else {
4963              maximumUp += columnUpper_[iColumn] * value;
4964            }
4965            if (columnLower_[iColumn] <= -large) {
4966              ++infiniteLower;
4967            } else {
4968              maximumDown += columnLower_[iColumn] * value;
4969            }
4970          } else if (value<0.0) {
4971            if (columnUpper_[iColumn] >= large) {
4972              ++infiniteLower;
4973            } else {
4974              maximumDown += columnUpper_[iColumn] * value;
4975            }
4976            if (columnLower_[iColumn] <= -large) {
4977              ++infiniteUpper;
4978            } else {
4979              maximumUp += columnLower_[iColumn] * value;
4980            }
4981          }
4982        }
4983        // Build in a margin of error
4984        maximumUp += 1.0e-8*fabs(maximumUp);
4985        maximumDown -= 1.0e-8*fabs(maximumDown);
4986        double maxUp = maximumUp+infiniteUpper*1.0e31;
4987        double maxDown = maximumDown-infiniteLower*1.0e31;
4988        if (maxUp <= rowUpper_[iRow] + tolerance && 
4989            maxDown >= rowLower_[iRow] - tolerance) {
4990         
4991          // Row is redundant - make totally free
4992          // NO - can't do this for postsolve
4993          // rowLower_[iRow]=-COIN_DBL_MAX;
4994          // rowUpper_[iRow]=COIN_DBL_MAX;
4995          //printf("Redundant row in presolveX %d\n",iRow);
4996
4997        } else {
4998          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
4999              maxDown > rowUpper_[iRow]+100.0*tolerance) {
5000            // problem is infeasible - exit at once
5001            numberInfeasible++;
5002            break;
5003          }
5004          double lower = rowLower_[iRow];
5005          double upper = rowUpper_[iRow];
5006          for (j = rStart; j < rEnd; ++j) {
5007            double value=element[j];
5008            iColumn = column[j];
5009            double nowLower = columnLower_[iColumn];
5010            double nowUpper = columnUpper_[iColumn];
5011            if (value > 0.0) {
5012              // positive value
5013              if (lower>-large) {
5014                if (!infiniteUpper) {
5015                  assert(nowUpper < large2);
5016                  newBound = nowUpper + 
5017                    (lower - maximumUp) / value;
5018                  // relax if original was large
5019                  if (fabs(maximumUp)>1.0e8)
5020                    newBound -= 1.0e-12*fabs(maximumUp);
5021                } else if (infiniteUpper==1&&nowUpper>large) {
5022                  newBound = (lower -maximumUp) / value;
5023                  // relax if original was large
5024                  if (fabs(maximumUp)>1.0e8)
5025                    newBound -= 1.0e-12*fabs(maximumUp);
5026                } else {
5027                  newBound = -COIN_DBL_MAX;
5028                }
5029                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
5030                  // Tighten the lower bound
5031                  numberChanged++;
5032                  // check infeasible (relaxed)
5033                  if (nowUpper < newBound) { 
5034                    if (nowUpper - newBound < 
5035                        -100.0*tolerance) 
5036                      numberInfeasible++;
5037                    else 
5038                      newBound=nowUpper;
5039                  }
5040                  columnLower_[iColumn] = newBound;
5041                  // adjust
5042                  double now;
5043                  if (nowLower<-large) {
5044                    now=0.0;
5045                    infiniteLower--;
5046                  } else {
5047                    now = nowLower;
5048                  }
5049                  maximumDown += (newBound-now) * value;
5050                  nowLower = newBound;
5051#ifdef DEBUG
5052                  checkCorrect(this,iRow,
5053                               element, rowStart, rowLength,
5054                               column,
5055                               columnLower_,  columnUpper_,
5056                               infiniteUpper,
5057                               infiniteLower,
5058                               maximumUp,
5059                               maximumDown);
5060#endif
5061                }
5062              } 
5063              if (upper <large) {
5064                if (!infiniteLower) {
5065                  assert(nowLower >- large2);
5066                  newBound = nowLower + 
5067                    (upper - maximumDown) / value;
5068                  // relax if original was large
5069                  if (fabs(maximumDown)>1.0e8)
5070                    newBound += 1.0e-12*fabs(maximumDown);
5071                } else if (infiniteLower==1&&nowLower<-large) {
5072                  newBound =   (upper - maximumDown) / value;
5073                  // relax if original was large
5074                  if (fabs(maximumDown)>1.0e8)
5075                    newBound += 1.0e-12*fabs(maximumDown);
5076                } else {
5077                  newBound = COIN_DBL_MAX;
5078                }
5079                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
5080                  // Tighten the upper bound
5081                  numberChanged++;
5082                  // check infeasible (relaxed)
5083                  if (nowLower > newBound) { 
5084                    if (newBound - nowLower < 
5085                        -100.0*tolerance) 
5086                      numberInfeasible++;
5087                    else 
5088                      newBound=nowLower;
5089                  }
5090                  columnUpper_[iColumn] = newBound;
5091                  // adjust
5092                  double now;
5093                  if (nowUpper>large) {
5094                    now=0.0;
5095                    infiniteUpper--;
5096                  } else {
5097                    now = nowUpper;
5098                  }
5099                  maximumUp += (newBound-now) * value;
5100                  nowUpper = newBound;
5101#ifdef DEBUG
5102                  checkCorrect(this,iRow,
5103                               element, rowStart, rowLength,
5104                               column,
5105                               columnLower_,  columnUpper_,
5106                               infiniteUpper,
5107                               infiniteLower,
5108                               maximumUp,
5109                               maximumDown);
5110#endif
5111                }
5112              }
5113            } else {
5114              // negative value
5115              if (lower>-large) {
5116                if (!infiniteUpper) {
5117                  assert(nowLower < large2);
5118                  newBound = nowLower + 
5119                    (lower - maximumUp) / value;
5120                  // relax if original was large
5121                  if (fabs(maximumUp)>1.0e8)
5122                    newBound += 1.0e-12*fabs(maximumUp);
5123                } else if (infiniteUpper==1&&nowLower<-large) {
5124                  newBound = (lower -maximumUp) / value;
5125                  // relax if original was large
5126                  if (fabs(maximumUp)>1.0e8)
5127                    newBound += 1.0e-12*fabs(maximumUp);
5128                } else {
5129                  newBound = COIN_DBL_MAX;
5130                }
5131                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
5132                  // Tighten the upper bound
5133                  numberChanged++;
5134                  // check infeasible (relaxed)
5135                  if (nowLower > newBound) { 
5136                    if (newBound - nowLower < 
5137                        -100.0*tolerance) 
5138                      numberInfeasible++;
5139                    else 
5140                      newBound=nowLower;
5141                  }
5142                  columnUpper_[iColumn] = newBound;
5143                  // adjust
5144                  double now;
5145                  if (nowUpper>large) {
5146                    now=0.0;
5147                    infiniteLower--;
5148                  } else {
5149                    now = nowUpper;
5150                  }
5151                  maximumDown += (newBound-now) * value;
5152                  nowUpper = newBound;
5153#ifdef DEBUG
5154                  checkCorrect(this,iRow,
5155                               element, rowStart, rowLength,
5156                               column,
5157                               columnLower_,  columnUpper_,
5158                               infiniteUpper,
5159                               infiniteLower,
5160                               maximumUp,
5161                               maximumDown);
5162#endif
5163                }
5164              }
5165              if (upper <large) {
5166                if (!infiniteLower) {
5167                  assert(nowUpper < large2);
5168                  newBound = nowUpper + 
5169                    (upper - maximumDown) / value;
5170                  // relax if original was large
5171                  if (fabs(maximumDown)>1.0e8)
5172                    newBound -= 1.0e-12*fabs(maximumDown);
5173                } else if (infiniteLower==1&&nowUpper>large) {
5174                  newBound =   (upper - maximumDown) / value;
5175                  // relax if original was large
5176                  if (fabs(maximumDown)>1.0e8)
5177                    newBound -= 1.0e-12*fabs(maximumDown);
5178                } else {
5179                  newBound = -COIN_DBL_MAX;
5180                }
5181                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
5182                  // Tighten the lower bound
5183                  numberChanged++;
5184                  // check infeasible (relaxed)
5185                  if (nowUpper < newBound) { 
5186                    if (nowUpper - newBound < 
5187                        -100.0*tolerance) 
5188                      numberInfeasible++;
5189                    else 
5190                      newBound=nowUpper;
5191                  }
5192                  columnLower_[iColumn] = newBound;
5193                  // adjust
5194                  double now;
5195                  if (nowLower<-large) {
5196                    now=0.0;
5197                    infiniteUpper--;
5198                  } else {
5199                    now = nowLower;
5200                  }
5201                  maximumUp += (newBound-now) * value;
5202                  nowLower = newBound;
5203#ifdef DEBUG
5204                  checkCorrect(this,iRow,
5205                               element, rowStart, rowLength,
5206                               column,
5207                               columnLower_,  columnUpper_,
5208                               infiniteUpper,
5209                               infiniteLower,
5210                               maximumUp,
5211                               maximumDown);
5212#endif
5213                }
5214              }
5215            }
5216          }
5217        }
5218      }
5219    }
5220    totalTightened += numberChanged;
5221    if (iPass==1)
5222      numberCheck=numberChanged>>4;
5223    if (numberInfeasible) break;
5224  }
5225  if (!numberInfeasible) {
5226    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
5227      <<totalTightened
5228      <<CoinMessageEol;
5229    // Set bounds slightly loose
5230    double useTolerance = 1.0e-3;
5231    if (doTight>0) {
5232      if (doTight>10) { 
5233        useTolerance=0.0;
5234      } else {
5235        while (doTight) {
5236          useTolerance *= 0.1;
5237          doTight--;
5238        }
5239      }
5240    }
5241    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5242      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
5243        // Make large bounds stay infinite
5244        if (saveUpper[iColumn]>1.0e30&&columnUpper_[iColumn]>1.0e10) {
5245          columnUpper_[iColumn]=COIN_DBL_MAX;
5246        }
5247        if (saveLower[iColumn]<-1.0e30&&columnLower_[iColumn]<-1.0e10) {
5248          columnLower_[iColumn]=-COIN_DBL_MAX;
5249        }
5250        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e-8) {
5251          // relax enough so will have correct dj
5252#if 1
5253          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
5254                                    columnLower_[iColumn]-100.0*useTolerance);
5255          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
5256                                    columnUpper_[iColumn]+100.0*useTolerance);
5257#else
5258          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
5259            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
5260              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
5261            } else {
5262              columnLower_[iColumn]=saveLower[iColumn];
5263              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
5264                                        saveLower[iColumn]+100.0*useTolerance);
5265            }
5266          } else {
5267            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
5268              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
5269            } else {
5270              columnUpper_[iColumn]=saveUpper[iColumn];
5271              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
5272                                        saveUpper[iColumn]-100.0*useTolerance);
5273            }
5274          }
5275#endif
5276        } else {
5277          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
5278            // relax a bit
5279            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
5280                                        saveUpper[iColumn]);
5281          }
5282          if (columnLower_[iColumn]>saveLower[iColumn]) {
5283            // relax a bit
5284            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
5285                                        saveLower[iColumn]);
5286          }
5287        }
5288      }
5289    }
5290    if (tightIntegers&&integerType_) {
5291      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5292        if (integerType_[iColumn]) {
5293          double value;
5294          value = floor(columnLower_[iColumn]+0.5);
5295          if (fabs(value-columnLower_[iColumn])>primalTolerance_)
5296            value = ceil(columnLower_[iColumn]);
5297          columnLower_[iColumn]=value;
5298          value = floor(columnUpper_[iColumn]+0.5);
5299          if (fabs(value-columnUpper_[iColumn])>primalTolerance_)
5300            value = floor(columnUpper_[iColumn]);
5301          columnUpper_[iColumn]=value;
5302          if (columnLower_[iColumn]>columnUpper_[iColumn])
5303            numberInfeasible++;
5304        }
5305      }
5306      if (numberInfeasible) {
5307        handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
5308          <<numberInfeasible
5309          <<CoinMessageEol;
5310        // restore column bounds
5311 CoinMemcpyN(saveLower,numberColumns_,columnLower_);
5312 CoinMemcpyN(saveUpper,numberColumns_,columnUpper_);
5313      }
5314    }
5315  } else {
5316    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
5317      <<numberInfeasible
5318      <<CoinMessageEol;
5319    // restore column bounds
5320    CoinMemcpyN(saveLower,numberColumns_,columnLower_);
5321    CoinMemcpyN(saveUpper,numberColumns_,columnUpper_);
5322  }
5323  delete [] saveLower;
5324  delete [] saveUpper;
5325  return (numberInfeasible);
5326}
5327//#define SAVE_AND_RESTORE
5328// dual
5329#include "ClpSimplexDual.hpp"
5330#include "ClpSimplexPrimal.hpp"
5331#ifndef SAVE_AND_RESTORE
5332int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
5333#else
5334int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
5335{
5336  // May be empty problem
5337  if (numberRows_&&numberColumns_) {
5338    // Save on file for debug
5339    int returnCode;
5340    returnCode = saveModel("debug.sav");
5341    if (returnCode) {
5342      printf("** Unable to save model to debug.sav\n");
5343      abort();
5344    }
5345    ClpSimplex temp;
5346    returnCode=temp.restoreModel("debug.sav");
5347    if (returnCode) {
5348      printf("** Unable to restore model from debug.sav\n");
5349      abort();
5350    }
5351    temp.setLogLevel(handler_->logLevel());
5352    // Now do dual
5353    returnCode=temp.dualDebug(ifValuesPass,startFinishOptions);
5354    // Move status and solution back
5355    int numberTotal = numberRows_+numberColumns_;
5356    CoinMemcpyN(temp.statusArray(),numberTotal,status_);
5357    CoinMemcpyN(temp.primalColumnSolution(),numberColumns_,columnActivity_);
5358    CoinMemcpyN(temp.primalRowSolution(),numberRows_,rowActivity_);
5359    CoinMemcpyN(temp.dualColumnSolution(),numberColumns_,reducedCost_);
5360    CoinMemcpyN(temp.dualRowSolution(),numberRows_,dual_);
5361    problemStatus_ = temp.problemStatus_;
5362    setObjectiveValue(temp.objectiveValue());
5363    setSumDualInfeasibilities(temp.sumDualInfeasibilities());
5364    setNumberDualInfeasibilities(temp.numberDualInfeasibilities());
5365    setSumPrimalInfeasibilities(temp.sumPrimalInfeasibilities());
5366    setNumberPrimalInfeasibilities(temp.numberPrimalInfeasibilities());
5367    setNumberIterations(temp.numberIterations());
5368    onStopped(); // set secondary status if stopped
5369    return returnCode;
5370  } else {
5371    // empty
5372    return dualDebug(ifValuesPass,startFinishOptions);
5373  }
5374}
5375int ClpSimplex::dualDebug (int ifValuesPass , int startFinishOptions)
5376#endif
5377{
5378  //double savedPivotTolerance = factorization_->pivotTolerance();
5379  int saveQuadraticActivated = objective_->activated();
5380  objective_->setActivated(0);
5381  ClpObjective * saveObjective = objective_;
5382  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
5383  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
5384      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
5385      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
5386      additional data and have no destructor or (non-default) constructor.
5387
5388      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
5389      in primal and dual.
5390
5391      As far as I can see this is perfectly safe.
5392  */
5393  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass, startFinishOptions);
5394  //int lastAlgorithm = -1;
5395  if ((specialOptions_&2048)!=0&&problemStatus_==10&&!numberPrimalInfeasibilities_
5396      &&sumDualInfeasibilities_<1000.0*dualTolerance_&&perturbation_>=100)
5397    problemStatus_=0; // ignore
5398  if (problemStatus_==10) {
5399    //printf("Cleaning up with primal\n");
5400#ifdef COIN_DEVELOP
5401    int saveNumberIterations=numberIterations_;
5402#endif
5403    //lastAlgorithm=1;
5404    int savePerturbation = perturbation_;
5405    int saveLog = handler_->logLevel();
5406    //handler_->setLogLevel(63);
5407    perturbation_=100;
5408    bool denseFactorization = initialDenseFactorization();
5409    // It will be safe to allow dense
5410    setInitialDenseFactorization(true);
5411    // Allow for catastrophe
5412    int saveMax = intParam_[ClpMaxNumIteration];
5413    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
5414      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
5415    // check which algorithms allowed
5416    int dummy;
5417    if (problemStatus_==10&&saveObjective==objective_)
5418      startFinishOptions |= 2;
5419    baseIteration_=numberIterations_;
5420    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
5421      returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
5422    else
5423      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
5424    baseIteration_=0;
5425    if (saveObjective != objective_) {
5426      // We changed objective to see if infeasible
5427      delete objective_;
5428      objective_=saveObjective;
5429      if (!problemStatus_) {
5430        // carry on
5431        returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
5432      }
5433    }
5434    if (problemStatus_==3&&numberIterations_<saveMax) {
5435#ifndef COIN_DEVELOP
5436      if (handler_->logLevel()==63)
5437#endif
5438        printf("looks like trouble - too many iterations in clean up - trying again\n");
5439      // flatten solution and try again
5440      int iRow,iColumn;
5441      for (iRow=0;iRow<numberRows_;iRow++) {
5442        if (getRowStatus(iRow)!=basic) {
5443          setRowStatus(iRow,superBasic);
5444          // but put to bound if close
5445          if (fabs(rowActivity_[iRow]-rowLower_[iRow])