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

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

sqrt to double

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