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

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

print only when feasible

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