source: trunk/ClpSimplex.cpp @ 554

Last change on this file since 554 was 554, checked in by forrest, 15 years ago

add cleanup

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