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

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

fix bad write and optional parentheses

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