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

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

for getting many solutions

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