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

Last change on this file since 1342 was 1342, checked in by forrest, 11 years ago

better to clear rowarrays on singular pivot

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