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

Last change on this file since 1397 was 1397, checked in by forrest, 10 years ago

stuff for Marc and Stefan

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