source: branches/BSP/trunk/Clp/src/ClpSimplex.cpp @ 1087

Last change on this file since 1087 was 1087, checked in by ladanyi, 12 years ago

committing merging w/ trunk of Clp, Osi, Cbc

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