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

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

minor changes

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