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

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

stupid mistake in primal

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