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

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

changes to try and make faster

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