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

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

add ids

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