source: stable/1.11/Clp/src/ClpSimplex.cpp @ 1498

Last change on this file since 1498 was 1498, checked in by forrest, 10 years ago

fix serious - if rare - bug which says optimal when not

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