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

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

fix assert failure

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