source: stable/1.8/Clp/src/ClpSimplex.cpp @ 1249

Last change on this file since 1249 was 1249, checked in by ladanyi, 12 years ago

Updated stable/1.8 to match trunk

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