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

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

add keyword for naive heuristic

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