source: stable/1.10/Clp/src/ClpSimplex.cpp @ 1359

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

modify asserts etc on difficult problem

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