source: trunk/ClpSimplex.cpp @ 560

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

CoinAssert? and check column numbers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 219.6 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    const double * obj = objective();
2422    for (i=0;i<numberColumns_;i++) {
2423      cost_[i]=obj[i];
2424      CoinAssert(fabs(cost_[i])<1.0e25);
2425    }
2426    if (rowObjective_)
2427      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
2428    else
2429      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
2430    // and initialize changes to zero
2431    //memset(cost_+nTotal,0,nTotal*sizeof(double));
2432  }
2433  // do scaling if needed
2434  if ((what&16)!=0) {
2435    if (!oldMatrix) {
2436      assert (scalingFlag_>0||!rowScale_);
2437      delete [] rowScale_;
2438      delete [] columnScale_;
2439      rowScale_=NULL;
2440      columnScale_=NULL;
2441    }
2442    if (scalingFlag_>0&&!rowScale_) {
2443      if (matrix_->scale(this))
2444        scalingFlag_=-scalingFlag_; // not scaled after all
2445      if (rowScale_&&automaticScale_) {
2446        // try automatic scaling
2447        double smallestObj=1.0e100;
2448        double largestObj=0.0;
2449        double largestRhs=0.0;
2450        for (i=0;i<numberColumns_+numberRows_;i++) {
2451          double value = fabs(cost_[i]);
2452          if (i<numberColumns_)
2453            value *= columnScale_[i];
2454          if (value&&lower_[i]!=upper_[i]) {
2455            smallestObj = CoinMin(smallestObj,value);
2456            largestObj = CoinMax(largestObj,value);
2457          }
2458          if (lower_[i]>0.0||upper_[i]<0.0) {
2459            double scale;
2460            if (i<numberColumns_)
2461              scale = 1.0/columnScale_[i];
2462            else
2463              scale = rowScale_[i-numberColumns_];
2464            //printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
2465            if (lower_[i]>0)
2466              largestRhs=CoinMax(largestRhs,lower_[i]*scale);
2467            if (upper_[i]<0.0)
2468              largestRhs=CoinMax(largestRhs,-upper_[i]*scale);
2469          }
2470        }
2471        printf("small obj %g, large %g - rhs %g\n",smallestObj,largestObj,largestRhs);
2472        bool scalingDone=false;
2473        // look at element range
2474        double smallestNegative;
2475        double largestNegative;
2476        double smallestPositive;
2477        double largestPositive;
2478        matrix_->rangeOfElements(smallestNegative, largestNegative,
2479                                 smallestPositive, largestPositive);
2480        smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2481        largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2482        if (largestObj) {
2483          double ratio = largestObj/smallestObj;
2484          double scale=1.0;
2485          if (ratio<1.0e8) {
2486            // reasonable
2487            if (smallestObj<1.0e-4) {
2488              // may as well scale up
2489              scalingDone=true;
2490              scale=1.0e-3/smallestObj;
2491            } else if (largestObj<1.0e6||(algorithm_>0&&largestObj<1.0e-4*infeasibilityCost_)) {
2492              //done=true;
2493            } else {
2494              scalingDone=true;
2495              if (algorithm_<0) {
2496                scale = 1.0e6/largestObj;
2497              } else {
2498                scale = CoinMax(1.0e6,1.0e-4*infeasibilityCost_)/largestObj;
2499              }
2500            }
2501          } else if (ratio<1.0e12) {
2502            // not so good
2503            if (smallestObj<1.0e-7) {
2504              // may as well scale up
2505              scalingDone=true;
2506              scale=1.0e-6/smallestObj;
2507            } else if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2508              //done=true;
2509            } else {
2510              scalingDone=true;
2511              if (algorithm_<0) {
2512                scale = 1.0e7/largestObj;
2513              } else {
2514                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2515              }
2516            }
2517          } else {
2518            // Really nasty problem
2519            if (smallestObj<1.0e-8) {
2520              // may as well scale up
2521              scalingDone=true;
2522              scale=1.0e-7/smallestObj;
2523              largestObj *= scale;
2524            } 
2525            if (largestObj<1.0e7||(algorithm_>0&&largestObj<1.0e-3*infeasibilityCost_)) {
2526              //done=true;
2527            } else {
2528              scalingDone=true;
2529              if (algorithm_<0) {
2530                scale = 1.0e7/largestObj;
2531              } else {
2532                scale = CoinMax(1.0e7,1.0e-3*infeasibilityCost_)/largestObj;
2533              }
2534            }
2535          }
2536          objectiveScale_=scale;
2537        }
2538        if (largestRhs>1.0e12) {
2539          scalingDone=true;
2540          rhsScale_=1.0e9/largestRhs;
2541        } else if (largestPositive>1.0e-14*smallestPositive&&largestRhs>1.0e6) {
2542          scalingDone=true;
2543          rhsScale_=1.0e6/largestRhs;
2544        } else {
2545          rhsScale_=1.0;
2546        }
2547        if (scalingDone) {
2548          handler_->message(CLP_RIM_SCALE,messages_)
2549            <<objectiveScale_<<rhsScale_
2550            <<CoinMessageEol;
2551        }
2552      }
2553    } else if (makeRowCopy&&scalingFlag_>0&&!oldMatrix) {
2554      matrix_->scaleRowCopy(this);
2555    }
2556  }
2557  if ((what&4)!=0) {
2558    double direction = optimizationDirection_*objectiveScale_;
2559    // but also scale by scale factors
2560    // not really sure about row scaling
2561    if (!rowScale_) {
2562      if (direction!=1.0) {
2563        for (i=0;i<numberRows_;i++)
2564          rowObjectiveWork_[i] *= direction;
2565        for (i=0;i<numberColumns_;i++)
2566          objectiveWork_[i] *= direction;
2567      }
2568    } else {
2569      for (i=0;i<numberRows_;i++)
2570        rowObjectiveWork_[i] *= direction/rowScale_[i];
2571      for (i=0;i<numberColumns_;i++)
2572        objectiveWork_[i] *= direction*columnScale_[i];
2573    }
2574  }
2575  if ((what&(1+32))!=0) {
2576    if(rowScale_) {
2577      for (i=0;i<numberColumns_;i++) {
2578        double multiplier = rhsScale_/columnScale_[i];
2579        if (columnLowerWork_[i]>-1.0e50)
2580          columnLowerWork_[i] *= multiplier;
2581        if (columnUpperWork_[i]<1.0e50)
2582          columnUpperWork_[i] *= multiplier;
2583       
2584      }
2585      for (i=0;i<numberRows_;i++) {
2586        double multiplier = rhsScale_*rowScale_[i];
2587        if (rowLowerWork_[i]>-1.0e50)
2588          rowLowerWork_[i] *= multiplier;
2589        if (rowUpperWork_[i]<1.0e50)
2590          rowUpperWork_[i] *= multiplier;
2591      }
2592    } else if (rhsScale_!=1.0) {
2593      for (i=0;i<numberColumns_+numberRows_;i++) {
2594        if (lower_[i]>-1.0e50)
2595          lower_[i] *= rhsScale_;
2596        if (upper_[i]<1.0e50)
2597          upper_[i] *= rhsScale_;
2598       
2599      }
2600    }
2601  }
2602  if ((what&(8+32))!=0) {
2603    if (rowScale_) {
2604      // on entry
2605      for (i=0;i<numberColumns_;i++) {
2606        columnActivityWork_[i] /= columnScale_[i];
2607        columnActivityWork_[i] *= rhsScale_;
2608        reducedCostWork_[i] *= objectiveScale_*columnScale_[i];
2609      }
2610      for (i=0;i<numberRows_;i++) {
2611        rowActivityWork_[i] *= rowScale_[i]*rhsScale_;
2612        dual_[i] /= rowScale_[i];
2613        dual_[i] *= objectiveScale_;
2614        rowReducedCost_[i] = dual_[i];
2615      }
2616    } else if (objectiveScale_!=1.0||rhsScale_!=1.0) {
2617      // on entry
2618      for (i=0;i<numberColumns_;i++) {
2619        columnActivityWork_[i] *= rhsScale_;
2620        reducedCostWork_[i] *= objectiveScale_;
2621      }
2622      for (i=0;i<numberRows_;i++) {
2623        rowActivityWork_[i] *= rhsScale_;
2624        dual_[i] *= objectiveScale_;
2625        rowReducedCost_[i] = dual_[i];
2626      }
2627    }
2628  }
2629 
2630  if ((what&16)!=0&&doSanityCheck) {
2631    // check rim of problem okay
2632    if (!sanityCheck())
2633      goodMatrix=false;
2634  } 
2635  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
2636  // maybe we need to move scales to SimplexModel for factorization?
2637  if (((what&8)!=0&&!pivotVariable_)||(newArrays&&!keepPivots)) {
2638    delete [] pivotVariable_;
2639    pivotVariable_=new int[numberRows2];
2640    for (int i=0;i<numberRows2;i++)
2641      pivotVariable_[i]=-1;
2642  } else if (what==63&&!keepPivots) {
2643    // just reset
2644    for (int i=0;i<numberRows2;i++)
2645      pivotVariable_[i]=-1;
2646  } else if (what==63) {
2647    // check pivots
2648    for (int i=0;i<numberRows2;i++) {
2649      int iSequence = pivotVariable_[i];
2650      if (iSequence<numberRows_+numberColumns_&&
2651          getStatus(iSequence)!=basic) {
2652        keepPivots =false;
2653        break;
2654      }
2655    }
2656    if (!keepPivots) {
2657      // reset
2658      for (int i=0;i<numberRows2;i++)
2659        pivotVariable_[i]=-1;
2660    } else {
2661      // clean
2662      for (int i=0;i<numberColumns_+numberRows_;i++) {
2663        Status status =getStatus(i);
2664        if (status!=basic) {
2665          if (upper_[i]==lower_[i]) {
2666            setStatus(i,isFixed);
2667            solution_[i]=lower_[i];
2668          } else if (status==atLowerBound) {
2669            if (lower_[i]>-1.0e20) {
2670              solution_[i]=lower_[i];
2671            } else {
2672              //printf("seq %d at lower of %g\n",i,lower_[i]);
2673              if (upper_[i]<1.0e20) {
2674                solution_[i]=upper_[i];
2675                setStatus(i,atUpperBound);
2676              } else {
2677                setStatus(i,isFree);
2678              }
2679            }
2680          } else if (status==atUpperBound) {
2681            if (upper_[i]<1.0e20) {
2682              solution_[i]=upper_[i];
2683            } else {
2684              //printf("seq %d at upper of %g\n",i,upper_[i]);
2685              if (lower_[i]>-1.0e20) {
2686                solution_[i]=lower_[i];
2687                setStatus(i,atLowerBound);
2688              } else {
2689                setStatus(i,isFree);
2690              }
2691            }
2692          }
2693        }
2694      }
2695    }
2696  }
2697
2698  if (what==63) {
2699    if (newArrays) {
2700      // get some arrays
2701      int iRow,iColumn;
2702      // these are "indexed" arrays so we always know where nonzeros are
2703      /**********************************************************
2704      rowArray_[3] is long enough for rows+columns
2705      *********************************************************/
2706      for (iRow=0;iRow<4;iRow++) {
2707        delete rowArray_[iRow];
2708        rowArray_[iRow]=new CoinIndexedVector();
2709        int length =numberRows2+factorization_->maximumPivots();
2710        if (iRow==3||objective_->type()>1)
2711          length += numberColumns_;
2712        rowArray_[iRow]->reserve(length);
2713      }
2714     
2715      for (iColumn=0;iColumn<2;iColumn++) {
2716        delete columnArray_[iColumn];
2717        columnArray_[iColumn]=new CoinIndexedVector();
2718        if (!iColumn)
2719          columnArray_[iColumn]->reserve(numberColumns_);
2720        else
2721          columnArray_[iColumn]->reserve(CoinMax(numberRows2,numberColumns_));
2722      }
2723    }   
2724  }
2725  double primalTolerance=dblParam_[ClpPrimalTolerance];
2726  if ((what&1)!=0&&doSanityCheck) {
2727    // fix any variables with tiny gaps
2728    for (i=0;i<numberColumns_;i++) {
2729      if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
2730        if (columnLowerWork_[i]>=0.0) {
2731          columnUpperWork_[i] = columnLowerWork_[i];
2732        } else if (columnUpperWork_[i]<=0.0) {
2733          columnLowerWork_[i] = columnUpperWork_[i];
2734        } else {
2735          columnUpperWork_[i] = 0.0;
2736          columnLowerWork_[i] = 0.0;
2737        }
2738      }
2739    }
2740    for (i=0;i<numberRows_;i++) {
2741      if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
2742        if (rowLowerWork_[i]>=0.0) {
2743          rowUpperWork_[i] = rowLowerWork_[i];
2744        } else if (rowUpperWork_[i]<=0.0) {
2745          rowLowerWork_[i] = rowUpperWork_[i];
2746        } else {
2747          rowUpperWork_[i] = 0.0;
2748          rowLowerWork_[i] = 0.0;
2749        }
2750      }
2751    }
2752  }
2753  if (problemStatus_==10) {
2754    problemStatus_=-1;
2755    handler_->setLogLevel(saveLevel); // switch back messages
2756  }
2757  if ((what&5)!=0) 
2758    matrix_->generalExpanded(this,9,what); // update costs and bounds if necessary
2759  return goodMatrix;
2760}
2761void
2762ClpSimplex::deleteRim(int getRidOfFactorizationData)
2763{
2764  int i;
2765  if (problemStatus_!=1&&problemStatus_!=2) {
2766    delete [] ray_;
2767    ray_=NULL;
2768  }
2769  // set upperOut_ to furthest away from bound so can use in dual for dualBound_
2770  upperOut_=1.0;
2771  // ray may be null if in branch and bound
2772  if (rowScale_) {
2773    // Collect infeasibilities
2774    int numberPrimalScaled=0;
2775    int numberPrimalUnscaled=0;
2776    int numberDualScaled=0;
2777    int numberDualUnscaled=0;
2778    double scaleC = 1.0/objectiveScale_;
2779    double scaleR = 1.0/rhsScale_;
2780    for (i=0;i<numberColumns_;i++) {
2781      double scaleFactor = columnScale_[i];
2782      double valueScaled = columnActivityWork_[i];
2783      double lowerScaled = columnLowerWork_[i];
2784      double upperScaled = columnUpperWork_[i];
2785      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2786        if (valueScaled<lowerScaled-primalTolerance_||
2787            valueScaled>upperScaled+primalTolerance_)
2788          numberPrimalScaled++;
2789        else
2790          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2791      }
2792      columnActivity_[i] = valueScaled*scaleFactor*scaleR;
2793      double value = columnActivity_[i];
2794      if (value<columnLower_[i]-primalTolerance_)
2795        numberPrimalUnscaled++;
2796      else if (value>columnUpper_[i]+primalTolerance_)
2797        numberPrimalUnscaled++;
2798      double valueScaledDual = reducedCostWork_[i];
2799      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2800        numberDualScaled++;
2801      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2802        numberDualScaled++;
2803      reducedCost_[i] = (valueScaledDual*scaleC)/scaleFactor;
2804      double valueDual = reducedCost_[i];
2805      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2806        numberDualUnscaled++;
2807      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2808        numberDualUnscaled++;
2809    }
2810    for (i=0;i<numberRows_;i++) {
2811      double scaleFactor = rowScale_[i];
2812      double valueScaled = rowActivityWork_[i];
2813      double lowerScaled = rowLowerWork_[i];
2814      double upperScaled = rowUpperWork_[i];
2815      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2816        if (valueScaled<lowerScaled-primalTolerance_||
2817            valueScaled>upperScaled+primalTolerance_)
2818          numberPrimalScaled++;
2819        else
2820          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2821      }
2822      rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
2823      double value = rowActivity_[i];
2824      if (value<rowLower_[i]-primalTolerance_)
2825        numberPrimalUnscaled++;
2826      else if (value>rowUpper_[i]+primalTolerance_)
2827        numberPrimalUnscaled++;
2828      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2829      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2830        numberDualScaled++;
2831      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2832        numberDualScaled++;
2833      dual_[i] = valueScaledDual*scaleFactor*scaleC;
2834      double valueDual = dual_[i]; 
2835      if (rowObjective_)
2836        valueDual += rowObjective_[i];
2837      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2838        numberDualUnscaled++;
2839      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2840        numberDualUnscaled++;
2841    }
2842    if (!problemStatus_&&!secondaryStatus_) {
2843      // See if we need to set secondary status
2844      if (numberPrimalUnscaled) {
2845        if (numberDualUnscaled) 
2846          secondaryStatus_=4;
2847        else
2848          secondaryStatus_=2;
2849      } else {
2850        if (numberDualUnscaled) 
2851          secondaryStatus_=3;
2852      }
2853    }
2854    if (problemStatus_==2) {
2855      for (i=0;i<numberColumns_;i++) {
2856        ray_[i] *= columnScale_[i];
2857      }
2858    } else if (problemStatus_==1&&ray_) {
2859      for (i=0;i<numberRows_;i++) {
2860        ray_[i] *= rowScale_[i];
2861      }
2862    }
2863  } else if (rhsScale_!=1.0||objectiveScale_!=1.0) {
2864    // Collect infeasibilities
2865    int numberPrimalScaled=0;
2866    int numberPrimalUnscaled=0;
2867    int numberDualScaled=0;
2868    int numberDualUnscaled=0;
2869    double scaleC = 1.0/objectiveScale_;
2870    double scaleR = 1.0/rhsScale_;
2871    for (i=0;i<numberColumns_;i++) {
2872      double valueScaled = columnActivityWork_[i];
2873      double lowerScaled = columnLowerWork_[i];
2874      double upperScaled = columnUpperWork_[i];
2875      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2876        if (valueScaled<lowerScaled-primalTolerance_||
2877            valueScaled>upperScaled+primalTolerance_)
2878          numberPrimalScaled++;
2879        else
2880          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2881      }
2882      columnActivity_[i] = valueScaled*scaleR;
2883      double value = columnActivity_[i];
2884      if (value<columnLower_[i]-primalTolerance_)
2885        numberPrimalUnscaled++;
2886      else if (value>columnUpper_[i]+primalTolerance_)
2887        numberPrimalUnscaled++;
2888      double valueScaledDual = reducedCostWork_[i];
2889      if (valueScaled>columnLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2890        numberDualScaled++;
2891      if (valueScaled<columnUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2892        numberDualScaled++;
2893      reducedCost_[i] = valueScaledDual*scaleC;
2894      double valueDual = reducedCost_[i];
2895      if (value>columnLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2896        numberDualUnscaled++;
2897      if (value<columnUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2898        numberDualUnscaled++;
2899    }
2900    for (i=0;i<numberRows_;i++) {
2901      double valueScaled = rowActivityWork_[i];
2902      double lowerScaled = rowLowerWork_[i];
2903      double upperScaled = rowUpperWork_[i];
2904      if (lowerScaled>-1.0e20||upperScaled<1.0e20) {
2905        if (valueScaled<lowerScaled-primalTolerance_||
2906            valueScaled>upperScaled+primalTolerance_)
2907          numberPrimalScaled++;
2908        else
2909          upperOut_ = CoinMax(upperOut_,CoinMin(valueScaled-lowerScaled,upperScaled-valueScaled));
2910      }
2911      rowActivity_[i] = valueScaled*scaleR;
2912      double value = rowActivity_[i];
2913      if (value<rowLower_[i]-primalTolerance_)
2914        numberPrimalUnscaled++;
2915      else if (value>rowUpper_[i]+primalTolerance_)
2916        numberPrimalUnscaled++;
2917      double valueScaledDual = dual_[i]+rowObjectiveWork_[i];;
2918      if (valueScaled>rowLowerWork_[i]+primalTolerance_&&valueScaledDual>dualTolerance_)
2919        numberDualScaled++;
2920      if (valueScaled<rowUpperWork_[i]-primalTolerance_&&valueScaledDual<-dualTolerance_)
2921        numberDualScaled++;
2922      dual_[i] = valueScaledDual*scaleC;
2923      double valueDual = dual_[i]; 
2924      if (rowObjective_)
2925        valueDual += rowObjective_[i];
2926      if (value>rowLower_[i]+primalTolerance_&&valueDual>dualTolerance_)
2927        numberDualUnscaled++;
2928      if (value<rowUpper_[i]-primalTolerance_&&valueDual<-dualTolerance_)
2929        numberDualUnscaled++;
2930    }
2931    if (!problemStatus_&&!secondaryStatus_) {
2932      // See if we need to set secondary status
2933      if (numberPrimalUnscaled) {
2934        if (numberDualUnscaled) 
2935          secondaryStatus_=4;
2936        else
2937          secondaryStatus_=2;
2938      } else {
2939        if (numberDualUnscaled) 
2940          secondaryStatus_=3;
2941      }
2942    }
2943  } else {
2944    if (columnActivityWork_) {
2945      for (i=0;i<numberColumns_;i++) {
2946        double value = columnActivityWork_[i];
2947        double lower = columnLowerWork_[i];
2948        double upper = columnUpperWork_[i];
2949        if (lower>-1.0e20||upper<1.0e20) {
2950          if (value>lower&&value<upper)
2951            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2952        }
2953        columnActivity_[i] = columnActivityWork_[i];
2954        reducedCost_[i] = reducedCostWork_[i];
2955      }
2956      for (i=0;i<numberRows_;i++) {
2957        double value = rowActivityWork_[i];
2958        double lower = rowLowerWork_[i];
2959        double upper = rowUpperWork_[i];
2960        if (lower>-1.0e20||upper<1.0e20) {
2961          if (value>lower&&value<upper)
2962            upperOut_ = CoinMax(upperOut_,CoinMin(value-lower,upper-value));
2963        }
2964        rowActivity_[i] = rowActivityWork_[i];
2965      }
2966    }
2967  }
2968  // switch off scalefactor if auto
2969  if (automaticScale_) {
2970    rhsScale_=1.0;
2971    objectiveScale_=1.0;
2972  }
2973  if (optimizationDirection_!=1.0) {
2974    // and modify all dual signs
2975    for (i=0;i<numberColumns_;i++) 
2976      reducedCost_[i] *= optimizationDirection_;
2977    for (i=0;i<numberRows_;i++) 
2978      dual_[i] *= optimizationDirection_;
2979  }
2980  // scaling may have been turned off
2981  scalingFlag_ = abs(scalingFlag_);
2982  if(getRidOfFactorizationData>0) {
2983    gutsOfDelete(getRidOfFactorizationData+1);
2984  } else {
2985    // at least get rid of nonLinearCost_
2986    delete nonLinearCost_;
2987    nonLinearCost_=NULL;
2988  }
2989  // get rid of data
2990  matrix_->generalExpanded(this,13,scalingFlag_);
2991}
2992void 
2993ClpSimplex::setDualBound(double value)
2994{
2995  if (value>0.0)
2996    dualBound_=value;
2997}
2998void 
2999ClpSimplex::setInfeasibilityCost(double value)
3000{
3001  if (value>0.0)
3002    infeasibilityCost_=value;
3003}
3004void ClpSimplex::setNumberRefinements( int value) 
3005{
3006  if (value>=0&&value<10)
3007    numberRefinements_=value;
3008}
3009// Sets row pivot choice algorithm in dual
3010void 
3011ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
3012{
3013  delete dualRowPivot_;
3014  dualRowPivot_ = choice.clone(true);
3015}
3016// Sets row pivot choice algorithm in dual
3017void 
3018ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
3019{
3020  delete primalColumnPivot_;
3021  primalColumnPivot_ = choice.clone(true);
3022}
3023// Passes in factorization
3024void 
3025ClpSimplex::setFactorization( ClpFactorization & factorization)
3026{
3027  delete factorization_;
3028  factorization_= new ClpFactorization(factorization);
3029}
3030/* Perturbation:
3031   -50 to +50 - perturb by this power of ten (-6 sounds good)
3032   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
3033   101 - we are perturbed
3034   102 - don't try perturbing again
3035   default is 100
3036*/
3037void 
3038ClpSimplex::setPerturbation(int value)
3039{
3040  if(value<=100&&value >=-1000) {
3041    perturbation_=value;
3042  } 
3043}
3044// Sparsity on or off
3045bool 
3046ClpSimplex::sparseFactorization() const
3047{
3048  return factorization_->sparseThreshold()!=0;
3049}
3050void 
3051ClpSimplex::setSparseFactorization(bool value)
3052{
3053  if (value) {
3054    if (!factorization_->sparseThreshold())
3055      factorization_->goSparse();
3056  } else {
3057    factorization_->sparseThreshold(0);
3058  }
3059}
3060void checkCorrect(ClpSimplex * model,int iRow,
3061                  const double * element,const int * rowStart,const int * rowLength,
3062                  const int * column,
3063                  const double * columnLower_, const double * columnUpper_,
3064                  int infiniteUpperC,
3065                  int infiniteLowerC,
3066                  double &maximumUpC,
3067                  double &maximumDownC)
3068{
3069  int infiniteUpper = 0;
3070  int infiniteLower = 0;
3071  double maximumUp = 0.0;
3072  double maximumDown = 0.0;
3073  CoinBigIndex rStart = rowStart[iRow];
3074  CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
3075  CoinBigIndex j;
3076  double large=1.0e15;
3077  int iColumn;
3078  // Compute possible lower and upper ranges
3079 
3080  for (j = rStart; j < rEnd; ++j) {
3081    double value=element[j];
3082    iColumn = column[j];
3083    if (value > 0.0) {
3084      if (columnUpper_[iColumn] >= large) {
3085        ++infiniteUpper;
3086      } else {
3087        maximumUp += columnUpper_[iColumn] * value;
3088      }
3089      if (columnLower_[iColumn] <= -large) {
3090        ++infiniteLower;
3091      } else {
3092        maximumDown += columnLower_[iColumn] * value;
3093      }
3094    } else if (value<0.0) {
3095      if (columnUpper_[iColumn] >= large) {
3096        ++infiniteLower;
3097      } else {
3098        maximumDown += columnUpper_[iColumn] * value;
3099      }
3100      if (columnLower_[iColumn] <= -large) {
3101        ++infiniteUpper;
3102      } else {
3103        maximumUp += columnLower_[iColumn] * value;
3104      }
3105    }
3106  }
3107  assert (infiniteLowerC==infiniteLower);
3108  assert (infiniteUpperC==infiniteUpper);
3109  if (fabs(maximumUp-maximumUpC)>1.0e-12*CoinMax(fabs(maximumUp),fabs(maximumUpC)))
3110    printf("row %d comp up %g, true up %g\n",iRow,
3111           maximumUpC,maximumUp);
3112  if (fabs(maximumDown-maximumDownC)>1.0e-12*CoinMax(fabs(maximumDown),fabs(maximumDownC)))
3113    printf("row %d comp down %g, true down %g\n",iRow,
3114           maximumDownC,maximumDown);
3115  maximumUpC=maximumUp;
3116  maximumDownC=maximumDown;
3117}
3118
3119/* Tightens primal bounds to make dual faster.  Unless
3120   fixed, bounds are slightly looser than they could be.
3121   This is to make dual go faster and is probably not needed
3122   with a presolve.  Returns non-zero if problem infeasible
3123
3124   Fudge for branch and bound - put bounds on columns of factor *
3125   largest value (at continuous) - should improve stability
3126   in branch and bound on infeasible branches (0.0 is off)
3127*/
3128int 
3129ClpSimplex::tightenPrimalBounds(double factor)
3130{
3131 
3132  // Get a row copy in standard format
3133  CoinPackedMatrix copy;
3134  copy.reverseOrderedCopyOf(*matrix());
3135  // get matrix data pointers
3136  const int * column = copy.getIndices();
3137  const CoinBigIndex * rowStart = copy.getVectorStarts();
3138  const int * rowLength = copy.getVectorLengths(); 
3139  const double * element = copy.getElements();
3140  int numberChanged=1,iPass=0;
3141  double large = largeValue(); // treat bounds > this as infinite
3142#ifndef NDEBUG
3143  double large2= 1.0e10*large;
3144#endif
3145  int numberInfeasible=0;
3146  int totalTightened = 0;
3147
3148  double tolerance = primalTolerance();
3149
3150
3151  // Save column bounds
3152  double * saveLower = new double [numberColumns_];
3153  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
3154  double * saveUpper = new double [numberColumns_];
3155  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
3156
3157  int iRow, iColumn;
3158
3159  // If wanted - tighten column bounds using solution
3160  if (factor) {
3161    assert (factor>1.0);
3162    double largest=0.0;
3163    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3164      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3165        largest = CoinMax(largest,fabs(columnActivity_[iColumn]));
3166      }
3167    }
3168    largest *= factor;
3169    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3170      if (columnUpper_[iColumn]-columnLower_[iColumn]>tolerance) {
3171        columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn],largest);
3172        columnLower_[iColumn] = CoinMax(columnLower_[iColumn],-largest);
3173      }
3174    }
3175  }
3176#define MAXPASS 10
3177
3178  // Loop round seeing if we can tighten bounds
3179  // Would be faster to have a stack of possible rows
3180  // and we put altered rows back on stack
3181  int numberCheck=-1;
3182  while(numberChanged>numberCheck) {
3183
3184    numberChanged = 0; // Bounds tightened this pass
3185   
3186    if (iPass==MAXPASS) break;
3187    iPass++;
3188   
3189    for (iRow = 0; iRow < numberRows_; iRow++) {
3190
3191      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
3192
3193        // possible row
3194        int infiniteUpper = 0;
3195        int infiniteLower = 0;
3196        double maximumUp = 0.0;
3197        double maximumDown = 0.0;
3198        double newBound;
3199        CoinBigIndex rStart = rowStart[iRow];
3200        CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
3201        CoinBigIndex j;
3202        // Compute possible lower and upper ranges
3203     
3204        for (j = rStart; j < rEnd; ++j) {
3205          double value=element[j];
3206          iColumn = column[j];
3207          if (value > 0.0) {
3208            if (columnUpper_[iColumn] >= large) {
3209              ++infiniteUpper;
3210            } else {
3211              maximumUp += columnUpper_[iColumn] * value;
3212            }
3213            if (columnLower_[iColumn] <= -large) {
3214              ++infiniteLower;
3215            } else {
3216              maximumDown += columnLower_[iColumn] * value;
3217            }
3218          } else if (value<0.0) {
3219            if (columnUpper_[iColumn] >= large) {
3220              ++infiniteLower;
3221            } else {
3222              maximumDown += columnUpper_[iColumn] * value;
3223            }
3224            if (columnLower_[iColumn] <= -large) {
3225              ++infiniteUpper;
3226            } else {
3227              maximumUp += columnLower_[iColumn] * value;
3228            }
3229          }
3230        }
3231        // Build in a margin of error
3232        maximumUp += 1.0e-8*fabs(maximumUp);
3233        maximumDown -= 1.0e-8*fabs(maximumDown);
3234        double maxUp = maximumUp+infiniteUpper*1.0e31;
3235        double maxDown = maximumDown-infiniteLower*1.0e31;
3236        if (maxUp <= rowUpper_[iRow] + tolerance && 
3237            maxDown >= rowLower_[iRow] - tolerance) {
3238         
3239          // Row is redundant - make totally free
3240          // NO - can't do this for postsolve
3241          // rowLower_[iRow]=-COIN_DBL_MAX;
3242          // rowUpper_[iRow]=COIN_DBL_MAX;
3243          //printf("Redundant row in presolveX %d\n",iRow);
3244
3245        } else {
3246          if (maxUp < rowLower_[iRow] -100.0*tolerance ||
3247              maxDown > rowUpper_[iRow]+100.0*tolerance) {
3248            // problem is infeasible - exit at once
3249            numberInfeasible++;
3250            break;
3251          }
3252          double lower = rowLower_[iRow];
3253          double upper = rowUpper_[iRow];
3254          for (j = rStart; j < rEnd; ++j) {
3255            double value=element[j];
3256            iColumn = column[j];
3257            double nowLower = columnLower_[iColumn];
3258            double nowUpper = columnUpper_[iColumn];
3259            if (value > 0.0) {
3260              // positive value
3261              if (lower>-large) {
3262                if (!infiniteUpper) {
3263                  assert(nowUpper < large2);
3264                  newBound = nowUpper + 
3265                    (lower - maximumUp) / value;
3266                  // relax if original was large
3267                  if (fabs(maximumUp)>1.0e8)
3268                    newBound -= 1.0e-12*fabs(maximumUp);
3269                } else if (infiniteUpper==1&&nowUpper>large) {
3270                  newBound = (lower -maximumUp) / value;
3271                  // relax if original was large
3272                  if (fabs(maximumUp)>1.0e8)
3273                    newBound -= 1.0e-12*fabs(maximumUp);
3274                } else {
3275                  newBound = -COIN_DBL_MAX;
3276                }
3277                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3278                  // Tighten the lower bound
3279                  columnLower_[iColumn] = newBound;
3280                  numberChanged++;
3281                  // check infeasible (relaxed)
3282                  if (nowUpper - newBound < 
3283                      -100.0*tolerance) {
3284                    numberInfeasible++;
3285                  }
3286                  // adjust
3287                  double now;
3288                  if (nowLower<-large) {
3289                    now=0.0;
3290                    infiniteLower--;
3291                  } else {
3292                    now = nowLower;
3293                  }
3294                  maximumDown += (newBound-now) * value;
3295                  nowLower = newBound;
3296#ifdef DEBUG
3297                  checkCorrect(this,iRow,
3298                               element, rowStart, rowLength,
3299                               column,
3300                               columnLower_,  columnUpper_,
3301                               infiniteUpper,
3302                               infiniteLower,
3303                               maximumUp,
3304                               maximumDown);
3305#endif
3306                }
3307              } 
3308              if (upper <large) {
3309                if (!infiniteLower) {
3310                  assert(nowLower >- large2);
3311                  newBound = nowLower + 
3312                    (upper - maximumDown) / value;
3313                  // relax if original was large
3314                  if (fabs(maximumDown)>1.0e8)
3315                    newBound += 1.0e-12*fabs(maximumDown);
3316                } else if (infiniteLower==1&&nowLower<-large) {
3317                  newBound =   (upper - maximumDown) / value;
3318                  // relax if original was large
3319                  if (fabs(maximumDown)>1.0e8)
3320                    newBound += 1.0e-12*fabs(maximumDown);
3321                } else {
3322                  newBound = COIN_DBL_MAX;
3323                }
3324                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3325                  // Tighten the upper bound
3326                  columnUpper_[iColumn] = newBound;
3327                  numberChanged++;
3328                  // check infeasible (relaxed)
3329                  if (newBound - nowLower < 
3330                      -100.0*tolerance) {
3331                    numberInfeasible++;
3332                  }
3333                  // adjust
3334                  double now;
3335                  if (nowUpper>large) {
3336                    now=0.0;
3337                    infiniteUpper--;
3338                  } else {
3339                    now = nowUpper;
3340                  }
3341                  maximumUp += (newBound-now) * value;
3342                  nowUpper = newBound;
3343#ifdef DEBUG
3344                  checkCorrect(this,iRow,
3345                               element, rowStart, rowLength,
3346                               column,
3347                               columnLower_,  columnUpper_,
3348                               infiniteUpper,
3349                               infiniteLower,
3350                               maximumUp,
3351                               maximumDown);
3352#endif
3353                }
3354              }
3355            } else {
3356              // negative value
3357              if (lower>-large) {
3358                if (!infiniteUpper) {
3359                  assert(nowLower < large2);
3360                  newBound = nowLower + 
3361                    (lower - maximumUp) / value;
3362                  // relax if original was large
3363                  if (fabs(maximumUp)>1.0e8)
3364                    newBound += 1.0e-12*fabs(maximumUp);
3365                } else if (infiniteUpper==1&&nowLower<-large) {
3366                  newBound = (lower -maximumUp) / value;
3367                  // relax if original was large
3368                  if (fabs(maximumUp)>1.0e8)
3369                    newBound += 1.0e-12*fabs(maximumUp);
3370                } else {
3371                  newBound = COIN_DBL_MAX;
3372                }
3373                if (newBound < nowUpper - 1.0e-12&&newBound<large) {
3374                  // Tighten the upper bound
3375                  columnUpper_[iColumn] = newBound;
3376                  numberChanged++;
3377                  // check infeasible (relaxed)
3378                  if (newBound - nowLower < 
3379                      -100.0*tolerance) {
3380                    numberInfeasible++;
3381                  }
3382                  // adjust
3383                  double now;
3384                  if (nowUpper>large) {
3385                    now=0.0;
3386                    infiniteLower--;
3387                  } else {
3388                    now = nowUpper;
3389                  }
3390                  maximumDown += (newBound-now) * value;
3391                  nowUpper = newBound;
3392#ifdef DEBUG
3393                  checkCorrect(this,iRow,
3394                               element, rowStart, rowLength,
3395                               column,
3396                               columnLower_,  columnUpper_,
3397                               infiniteUpper,
3398                               infiniteLower,
3399                               maximumUp,
3400                               maximumDown);
3401#endif
3402                }
3403              }
3404              if (upper <large) {
3405                if (!infiniteLower) {
3406                  assert(nowUpper < large2);
3407                  newBound = nowUpper + 
3408                    (upper - maximumDown) / value;
3409                  // relax if original was large
3410                  if (fabs(maximumDown)>1.0e8)
3411                    newBound -= 1.0e-12*fabs(maximumDown);
3412                } else if (infiniteLower==1&&nowUpper>large) {
3413                  newBound =   (upper - maximumDown) / value;
3414                  // relax if original was large
3415                  if (fabs(maximumDown)>1.0e8)
3416                    newBound -= 1.0e-12*fabs(maximumDown);
3417                } else {
3418                  newBound = -COIN_DBL_MAX;
3419                }
3420                if (newBound > nowLower + 1.0e-12&&newBound>-large) {
3421                  // Tighten the lower bound
3422                  columnLower_[iColumn] = newBound;
3423                  numberChanged++;
3424                  // check infeasible (relaxed)
3425                  if (nowUpper - newBound < 
3426                      -100.0*tolerance) {
3427                    numberInfeasible++;
3428                  }
3429                  // adjust
3430                  double now;
3431                  if (nowLower<-large) {
3432                    now=0.0;
3433                    infiniteUpper--;
3434                  } else {
3435                    now = nowLower;
3436                  }
3437                  maximumUp += (newBound-now) * value;
3438                  nowLower = newBound;
3439#ifdef DEBUG
3440                  checkCorrect(this,iRow,
3441                               element, rowStart, rowLength,
3442                               column,
3443                               columnLower_,  columnUpper_,
3444                               infiniteUpper,
3445                               infiniteLower,
3446                               maximumUp,
3447                               maximumDown);
3448#endif
3449                }
3450              }
3451            }
3452          }
3453        }
3454      }
3455    }
3456    totalTightened += numberChanged;
3457    if (iPass==1)
3458      numberCheck=numberChanged>>4;
3459    if (numberInfeasible) break;
3460  }
3461  if (!numberInfeasible) {
3462    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
3463      <<totalTightened
3464      <<CoinMessageEol;
3465    // Set bounds slightly loose
3466    double useTolerance = 1.0e-3;
3467    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3468      if (saveUpper[iColumn]>saveLower[iColumn]+useTolerance) {
3469        if (columnUpper_[iColumn]-columnLower_[iColumn]<useTolerance+1.0e8) {
3470          // relax enough so will have correct dj
3471#if 1
3472          columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3473                                    columnLower_[iColumn]-100.0*useTolerance);
3474          columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3475                                    columnUpper_[iColumn]+100.0*useTolerance);
3476#else
3477          if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) {
3478            if (columnUpper_[iColumn]- 100.0*useTolerance>saveLower[iColumn]) {
3479              columnLower_[iColumn]=columnUpper_[iColumn]-100.0*useTolerance;
3480            } else {
3481              columnLower_[iColumn]=saveLower[iColumn];
3482              columnUpper_[iColumn]=CoinMin(saveUpper[iColumn],
3483                                        saveLower[iColumn]+100.0*useTolerance);
3484            }
3485          } else {
3486            if (columnLower_[iColumn]+100.0*useTolerance<saveUpper[iColumn]) {
3487              columnUpper_[iColumn]=columnLower_[iColumn]+100.0*useTolerance;
3488            } else {
3489              columnUpper_[iColumn]=saveUpper[iColumn];
3490              columnLower_[iColumn]=CoinMax(saveLower[iColumn],
3491                                        saveUpper[iColumn]-100.0*useTolerance);
3492            }
3493          }
3494#endif
3495        } else {
3496          if (columnUpper_[iColumn]<saveUpper[iColumn]) {
3497            // relax a bit
3498            columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn]+100.0*useTolerance,
3499                                        saveUpper[iColumn]);
3500          }
3501          if (columnLower_[iColumn]>saveLower[iColumn]) {
3502            // relax a bit
3503            columnLower_[iColumn] = CoinMax(columnLower_[iColumn]-100.0*useTolerance,
3504                                        saveLower[iColumn]);
3505          }
3506        }
3507      }
3508    }
3509  } else {
3510    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
3511      <<numberInfeasible
3512      <<CoinMessageEol;
3513    // restore column bounds
3514    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
3515    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
3516  }
3517  delete [] saveLower;
3518  delete [] saveUpper;
3519  return (numberInfeasible);
3520}
3521// dual
3522#include "ClpSimplexDual.hpp"
3523#include "ClpSimplexPrimal.hpp"
3524int ClpSimplex::dual (int ifValuesPass , int startFinishOptions)
3525{
3526  int saveQuadraticActivated = objective_->activated();
3527  objective_->setActivated(0);
3528  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
3529  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3530      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3531      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3532      additional data and have no destructor or (non-default) constructor.
3533
3534      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3535      in primal and dual.
3536
3537      As far as I can see this is perfectly safe.
3538  */
3539  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass, startFinishOptions);
3540  if ((specialOptions_&2048)!=0&&problemStatus_==10&&!numberPrimalInfeasibilities_
3541      &&sumDualInfeasibilities_<1000.0*dualTolerance_&&perturbation_>=100)
3542    problemStatus_=0; // ignore
3543  if (problemStatus_==10) {
3544    //printf("Cleaning up with primal\n");
3545    int savePerturbation = perturbation_;
3546    int saveLog = handler_->logLevel();
3547    //handler_->setLogLevel(63);
3548    perturbation_=100;
3549    bool denseFactorization = initialDenseFactorization();
3550    // It will be safe to allow dense
3551    setInitialDenseFactorization(true);
3552    // Allow for catastrophe
3553    int saveMax = intParam_[ClpMaxNumIteration];
3554    if (intParam_[ClpMaxNumIteration]>100000+numberIterations_)
3555      intParam_[ClpMaxNumIteration] = numberIterations_ + 1000 + 2*numberRows_+numberColumns_;
3556    // check which algorithms allowed
3557    int dummy;
3558    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
3559      returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
3560    else
3561      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
3562    if (problemStatus_==3&&numberIterations_<saveMax) {
3563      if (handler_->logLevel()==63)
3564        printf("looks like trouble - too many iterations in clean up - trying again\n");
3565      // flatten solution and try again
3566      int iRow,iColumn;
3567      for (iRow=0;iRow<numberRows_;iRow++) {
3568        if (getRowStatus(iRow)!=basic) {
3569          setRowStatus(iRow,superBasic);
3570          // but put to bound if close
3571          if (fabs(rowActivity_[iRow]-rowLower_[iRow])
3572              <=primalTolerance_) {
3573            rowActivity_[iRow]=rowLower_[iRow];
3574            setRowStatus(iRow,atLowerBound);
3575          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])
3576                     <=primalTolerance_) {
3577            rowActivity_[iRow]=rowUpper_[iRow];
3578            setRowStatus(iRow,atUpperBound);
3579          }
3580        }
3581      }
3582      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3583        if (getColumnStatus(iColumn)!=basic) {
3584          setColumnStatus(iColumn,superBasic);
3585          // but put to bound if close
3586          if (fabs(columnActivity_[iColumn]-columnLower_[iColumn])
3587              <=primalTolerance_) {
3588            columnActivity_[iColumn]=columnLower_[iColumn];
3589            setColumnStatus(iColumn,atLowerBound);
3590          } else if (fabs(columnActivity_[iColumn]
3591                          -columnUpper_[iColumn])
3592                     <=primalTolerance_) {
3593            columnActivity_[iColumn]=columnUpper_[iColumn];
3594            setColumnStatus(iColumn,atUpperBound);
3595          }
3596        }
3597      }
3598      problemStatus_=-1;
3599      intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 
3600                                          2*numberRows_+numberColumns_,saveMax);
3601      perturbation_=savePerturbation;
3602      returnCode = ((ClpSimplexPrimal *) this)->primal(0);
3603      if (problemStatus_==3&&numberIterations_<saveMax&& 
3604          handler_->logLevel()>0)
3605        printf("looks like real trouble - too many iterations in second clean up - giving up\n");
3606    }
3607    intParam_[ClpMaxNumIteration] = saveMax;
3608
3609    setInitialDenseFactorization(denseFactorization);
3610    perturbation_=savePerturbation;
3611    if (problemStatus_==10) 
3612      problemStatus_=0;
3613    handler_->setLogLevel(saveLog);
3614  }
3615  objective_->setActivated(saveQuadraticActivated);
3616  return returnCode;
3617}
3618#include "ClpQuadraticObjective.hpp"
3619// primal
3620int ClpSimplex::primal (int ifValuesPass , int startFinishOptions)
3621{
3622  // See if nonlinear
3623  if (objective_->type()>1&&objective_->activated()) 
3624    return reducedGradient();
3625 
3626  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
3627  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
3628      Classes ClpSimplexDual, ClpSimplexPrimal, (ClpSimplexNonlinear)
3629      and ClpSimplexOther all exist and inherit from ClpSimplex but have no
3630      additional data and have no destructor or (non-default) constructor.
3631
3632      This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
3633      in primal and dual.
3634
3635      As far as I can see this is perfectly safe.
3636  */
3637  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass,startFinishOptions);
3638  if (problemStatus_==10) {
3639    //printf("Cleaning up with dual\n");
3640    int savePerturbation = perturbation_;
3641    perturbation_=100;
3642    bool denseFactorization = initialDenseFactorization();
3643    // It will be safe to allow dense
3644    setInitialDenseFactorization(true);
3645    // check which algorithms allowed
3646    int dummy;
3647    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0&&(specialOptions_&8192)==0) {
3648      double saveBound = dualBound_;
3649      // upperOut_ has largest away from bound
3650      dualBound_=CoinMin(2.0*upperOut_,dualBound_);
3651      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
3652      dualBound_=saveBound;
3653    } else {
3654      returnCode = ((ClpSimplexPrimal *) this)->primal(0,startFinishOptions);
3655    }
3656    setInitialDenseFactorization(denseFactorization);
3657    perturbation_=savePerturbation;
3658    if (problemStatus_==10) 
3659      problemStatus_=0;
3660  }
3661  return returnCode;
3662}
3663/* Dual ranging.
3664   This computes increase/decrease in cost for each given variable and corresponding
3665   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3666   and numberColumns.. for artificials/slacks.
3667   For non-basic variables the sequence number will be that of the non-basic variables.
3668   
3669   Up to user to provide correct length arrays.
3670   
3671   Returns non-zero if infeasible unbounded etc
3672*/
3673#include "ClpSimplexOther.hpp"
3674int ClpSimplex::dualRanging(int numberCheck,const int * which,
3675                            double * costIncrease, int * sequenceIncrease,
3676                            double * costDecrease, int * sequenceDecrease)
3677{
3678  int savePerturbation = perturbation_;
3679  perturbation_=100;
3680  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3681  if (problemStatus_==10) {
3682    //printf("Cleaning up with dual\n");
3683    bool denseFactorization = initialDenseFactorization();
3684    // It will be safe to allow dense
3685    setInitialDenseFactorization(true);
3686    // check which algorithms allowed
3687    int dummy;
3688    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3689      // upperOut_ has largest away from bound
3690      double saveBound=dualBound_;
3691      if (upperOut_>0.0)
3692        dualBound_=2.0*upperOut_;
3693      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3694      dualBound_=saveBound;
3695    } else {
3696      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3697    }
3698    setInitialDenseFactorization(denseFactorization);
3699    if (problemStatus_==10) 
3700      problemStatus_=0;
3701  }
3702  perturbation_=savePerturbation;
3703  if (problemStatus_||secondaryStatus_==6) {
3704    finish(); // get rid of arrays
3705    return 1; // odd status
3706  }
3707  ((ClpSimplexOther *) this)->dualRanging(numberCheck,which,
3708                                          costIncrease,sequenceIncrease,
3709                                          costDecrease,sequenceDecrease);
3710  finish(); // get rid of arrays
3711  return 0;
3712}
3713/* Primal ranging.
3714   This computes increase/decrease in value for each given variable and corresponding
3715   sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
3716   and numberColumns.. for artificials/slacks.
3717   For basic variables the sequence number will be that of the basic variables.
3718   
3719   Up to user to providen correct length arrays.
3720   
3721   Returns non-zero if infeasible unbounded etc
3722*/
3723int ClpSimplex::primalRanging(int numberCheck,const int * which,
3724                  double * valueIncrease, int * sequenceIncrease,
3725                  double * valueDecrease, int * sequenceDecrease)
3726{
3727  int savePerturbation = perturbation_;
3728  perturbation_=100;
3729  int returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3730  if (problemStatus_==10) {
3731    //printf("Cleaning up with dual\n");
3732    bool denseFactorization = initialDenseFactorization();
3733    // It will be safe to allow dense
3734    setInitialDenseFactorization(true);
3735    // check which algorithms allowed
3736    int dummy;
3737    if ((matrix_->generalExpanded(this,4,dummy)&2)!=0) {
3738      // upperOut_ has largest away from bound
3739      double saveBound=dualBound_;
3740      if (upperOut_>0.0)
3741        dualBound_=2.0*upperOut_;
3742      returnCode = ((ClpSimplexDual *) this)->dual(0,1);
3743      dualBound_=saveBound;
3744    } else {
3745      returnCode = ((ClpSimplexPrimal *) this)->primal(0,1);
3746    }
3747    setInitialDenseFactorization(denseFactorization);
3748    if (problemStatus_==10) 
3749      problemStatus_=0;
3750  }
3751  perturbation_=savePerturbation;
3752  if (problemStatus_||secondaryStatus_==6) {
3753    finish(); // get rid of arrays
3754    return 1; // odd status
3755  }
3756  ((ClpSimplexOther *) this)->primalRanging(numberCheck,which,
3757                                          valueIncrease,sequenceIncrease,
3758                                          valueDecrease,sequenceDecrease);
3759  finish(); // get rid of arrays
3760  return 0;
3761}
3762/* Write the basis in MPS format to the specified file.
3763   If writeValues true writes values of structurals
3764   (and adds VALUES to end of NAME card)
3765   
3766   Row and column names may be null.
3767   formatType is
3768   <ul>
3769   <li> 0 - normal
3770   <li> 1 - extra accuracy
3771   <li> 2 - IEEE hex (later)
3772   </ul>
3773   
3774   Returns non-zero on I/O error
3775*/
3776int 
3777ClpSimplex::writeBasis(const char *filename,
3778                            bool writeValues,
3779                            int formatType) const
3780{
3781  return ((const ClpSimplexOther *) this)->writeBasis(filename,writeValues,
3782                                         formatType);
3783}
3784// Read a basis from the given filename
3785int 
3786ClpSimplex::readBasis(const char *filename)
3787{
3788  return ((ClpSimplexOther *) this)->readBasis(filename);
3789}
3790#include "ClpSimplexNonlinear.hpp"
3791/* Solves nonlinear problem using SLP - may be used as crash
3792   for other algorithms when number of iterations small
3793*/
3794int 
3795ClpSimplex::nonlinearSLP(int numberPasses, double deltaTolerance)
3796{
3797  return ((ClpSimplexNonlinear *) this)->primalSLP(numberPasses,deltaTolerance);
3798}
3799// Solves non-linear using reduced gradient
3800int ClpSimplex::reducedGradient(int phase)
3801{
3802  if (objective_->type()<2||!objective_->activated()) {
3803    // no quadratic part
3804    return primal(0);
3805  }
3806  // get feasible
3807  if ((this->status()<0||numberPrimalInfeasibilities())&&phase==0) {
3808    objective_->setActivated(0);
3809    double saveDirection = optimizationDirection();
3810    setOptimizationDirection(0.0);
3811    primal(1);
3812    setOptimizationDirection(saveDirection);
3813    objective_->setActivated(1);
3814    // still infeasible
3815    if (numberPrimalInfeasibilities())
3816      return 0;
3817  }
3818  // Now enter method
3819  int returnCode = ((ClpSimplexNonlinear *) this)->primal();
3820  return returnCode;
3821}
3822#include "ClpPredictorCorrector.hpp"
3823#include "ClpCholeskyBase.hpp"
3824// Preference is WSSMP, TAUCS, UFL (just ordering) then base
3825#if WSSMP_BARRIER
3826#include "ClpCholeskyWssmp.hpp"
3827#include "ClpCholeskyWssmpKKT.hpp"
3828#elif UFL_BARRIER
3829#include "ClpCholeskyUfl.hpp"
3830#elif TAUCS_BARRIER
3831#include "ClpCholeskyTaucs.hpp"
3832#endif
3833#include "ClpPresolve.hpp"
3834/* Solves using barrier (assumes you have good cholesky factor code).
3835   Does crossover to simplex if asked*/
3836int 
3837ClpSimplex::barrier(bool crossover)
3838{
3839  ClpSimplex * model2 = this;
3840  int savePerturbation=perturbation_;
3841  ClpInterior barrier;
3842  barrier.borrowModel(*model2);
3843// Preference is WSSMP, UFL, TAUCS then base
3844#ifdef WSSMP_BARRIER
3845  ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100,model2->numberRows()/10));
3846  //ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp();
3847  //ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100,model2->numberRows()/10));
3848#elif UFL_BARRIER
3849  ClpCholeskyUfl * cholesky = new ClpCholeskyUfl();
3850#elif TAUCS_BARRIER
3851  ClpCholeskyTaucs * cholesky = new ClpCholeskyTaucs();
3852#else
3853  ClpCholeskyBase * cholesky = new ClpCholeskyBase();
3854#endif
3855  barrier.setCholesky(cholesky);
3856  int numberRows = model2->numberRows();
3857  int numberColumns = model2->numberColumns();
3858  int saveMaxIts = model2->maximumIterations();
3859  if (saveMaxIts<1000) {
3860    barrier.setMaximumBarrierIterations(saveMaxIts);
3861    model2->setMaximumIterations(1000000);
3862  }
3863  barrier.primalDual();
3864  int barrierStatus=barrier.status();
3865  double gap = barrier.complementarityGap();
3866  // get which variables are fixed
3867  double * saveLower=NULL;
3868  double * saveUpper=NULL;
3869  ClpPresolve pinfo2;
3870  ClpSimplex * saveModel2=NULL;
3871  int numberFixed = barrier.numberFixed();
3872  if (numberFixed*20>barrier.numberRows()&&numberFixed>5000&&crossover&&0) {
3873    // may as well do presolve
3874    int numberRows = barrier.numberRows();
3875    int numberColumns = barrier.numberColumns();
3876    int numberTotal = numberRows+numberColumns;
3877    saveLower = new double [numberTotal];
3878    saveUpper = new double [numberTotal];
3879    memcpy(saveLower,barrier.columnLower(),numberColumns*sizeof(double));
3880    memcpy(saveLower+numberColumns,barrier.rowLower(),numberRows*sizeof(double));
3881    memcpy(saveUpper,barrier.columnUpper(),numberColumns*sizeof(double));
3882    memcpy(saveUpper+numberColumns,barrier.rowUpper(),numberRows*sizeof(double));
3883    barrier.fixFixed();
3884    saveModel2=model2;
3885  }
3886  barrier.returnModel(*model2);
3887  double * rowPrimal = new double [numberRows];
3888  double * columnPrimal = new double [numberColumns];
3889  double * rowDual = new double [numberRows];
3890  double * columnDual = new double [numberColumns];
3891  // move solutions other way
3892  CoinMemcpyN(model2->primalRowSolution(),
3893              numberRows,rowPrimal);
3894  CoinMemcpyN(model2->dualRowSolution(),
3895              numberRows,rowDual);
3896  CoinMemcpyN(model2->primalColumnSolution(),
3897              numberColumns,columnPrimal);
3898  CoinMemcpyN(model2->dualColumnSolution(),
3899              numberColumns,columnDual);
3900  if (saveModel2) {
3901    // do presolve
3902    model2 = pinfo2.presolvedModel(*model2,1.0e-8,
3903                                   false,5,true);
3904  }
3905  if (barrierStatus<4&&crossover) {
3906    // make sure no status left
3907    model2->createStatus();
3908    // solve
3909    model2->setPerturbation(100);
3910    // throw some into basis
3911    {
3912      int numberRows = model2->numberRows();
3913      int numberColumns = model2->numberColumns();
3914      double * dsort = new double[numberColumns];
3915      int * sort = new int[numberColumns];
3916      int n=0;
3917      const double * columnLower = model2->columnLower();
3918      const double * columnUpper = model2->columnUpper();
3919      const double * primalSolution = model2->primalColumnSolution();
3920      double tolerance = 10.0*primalTolerance_;
3921      int i;
3922      for ( i=0;i<numberRows;i++) 
3923        model2->setRowStatus(i,superBasic);
3924      for ( i=0;i<numberColumns;i++) {
3925        double distance = CoinMin(columnUpper[i]-primalSolution[i],
3926                              primalSolution[i]-columnLower[i]);
3927        if (distance>tolerance) {
3928          dsort[n]=-distance;
3929          sort[n++]=i;
3930          model2->setStatus(i,superBasic);
3931        } else if (distance>primalTolerance_) {
3932          model2->setStatus(i,superBasic);
3933        } else if (primalSolution[i]<=columnLower[i]+primalTolerance_) {
3934          model2->setStatus(i,atLowerBound);
3935        } else {
3936          model2->setStatus(i,atUpperBound);
3937        }
3938      }
3939      CoinSort_2(dsort,dsort+n,sort);
3940      n = CoinMin(numberRows,n);
3941      for ( i=0;i<n;i++) {
3942        int iColumn = sort[i];
3943        model2->setStatus(iColumn,basic);
3944      }
3945      delete [] sort;
3946      delete [] dsort;
3947    }
3948    if (gap<1.0e-3*((double) (numberRows+numberColumns))) {
3949      int numberRows = model2->numberRows();
3950      int numberColumns = model2->numberColumns();
3951      // just primal values pass
3952      double saveScale = model2->objectiveScale();
3953      model2->setObjectiveScale(1.0e-3);
3954      model2->primal(2);
3955      model2->setObjectiveScale(saveScale);
3956      // save primal solution and copy back dual
3957      CoinMemcpyN(model2->primalRowSolution(),
3958                  numberRows,rowPrimal);
3959      CoinMemcpyN(rowDual,
3960                  numberRows,model2->dualRowSolution());
3961      CoinMemcpyN(model2->primalColumnSolution(),
3962                  numberColumns,columnPrimal);
3963      CoinMemcpyN(columnDual,
3964                  numberColumns,model2->dualColumnSolution());
3965      //model2->primal(1);
3966      // clean up reduced costs and flag variables
3967      {
3968        double * dj = model2->dualColumnSolution();
3969        double * cost = model2->objective();
3970        double * saveCost = new double[numberColumns];
3971        memcpy(saveCost,cost,numberColumns*sizeof(double));
3972        double * saveLower = new double[numberColumns];
3973        double * lower = model2->columnLower();
3974        memcpy(saveLower,lower,numberColumns*sizeof(double));
3975        double * saveUpper = new double[numberColumns];
3976        double * upper = model2->columnUpper();
3977        memcpy(saveUpper,upper,numberColumns*sizeof(double));
3978        int i;
3979        double tolerance = 10.0*dualTolerance_;
3980        for ( i=0;i<numberColumns;i++) {
3981          if (model2->getStatus(i)==basic) {
3982            dj[i]=0.0;
3983          } else if (model2->getStatus(i)==atLowerBound) {
3984            if (optimizationDirection_*dj[i]<tolerance) {
3985              if (optimizationDirection_*dj[i]<0.0) {
3986                //if (dj[i]<-1.0e-3)
3987                //printf("bad dj at lb %d %g\n",i,dj[i]);
3988                cost[i] -= dj[i];
3989                dj[i]=0.0;
3990              }
3991            } else {
3992              upper[i]=lower[i];
3993            }
3994          } else if (model2->getStatus(i)==atUpperBound) {
3995            if (optimizationDirection_*dj[i]>tolerance) {
3996              if (optimizationDirection_*dj[i]>0.0) {
3997                //if (dj[i]>1.0e-3)
3998                //printf("bad dj at ub %d %g\n",i,dj[i]);
3999                cost[i] -= dj[i];
4000                dj[i]=0.0;
4001              }
4002            } else {
4003              lower[i]=upper[i];
4004            }
4005          }
4006        }
4007        // just dual values pass
4008        //model2->setLogLevel(63);
4009        //model2->setFactorizationFrequency(1);
4010        model2->dual(2);
4011        memcpy(cost,saveCost,numberColumns*sizeof(double));
4012        delete [] saveCost;
4013        memcpy(lower,saveLower,numberColumns*sizeof(double));
4014        delete [] saveLower;
4015        memcpy(upper,saveUpper,numberColumns*sizeof(double));
4016        delete [] saveUpper;
4017      }
4018      // and finish
4019      // move solutions
4020      CoinMemcpyN(rowPrimal,
4021                  numberRows,model2->primalRowSolution());
4022      CoinMemcpyN(columnPrimal,
4023                  numberColumns,model2->primalColumnSolution());
4024    }
4025//     double saveScale = model2->objectiveScale();
4026//     model2->setObjectiveScale(1.0e-3);
4027//     model2->primal(2);
4028//    model2->setObjectiveScale(saveScale);
4029    model2->primal(1);
4030  } else if (barrierStatus==4&&crossover) {
4031    // memory problems
4032    model2->setPerturbation(savePerturbation);
4033    model2->createStatus();
4034    model2->dual();
4035  }
4036  model2->setMaximumIterations(saveMaxIts);
4037  delete [] rowPrimal;
4038  delete [] columnPrimal;
4039  delete [] rowDual;
4040  delete [] columnDual;
4041  if (saveLower) {
4042    pinfo2.postsolve(true);
4043    delete model2;
4044    model2=saveModel2;
4045    int numberRows = model2->numberRows();
4046    int numberColumns = model2->numberColumns();
4047    memcpy(model2->columnLower(),saveLower,numberColumns*sizeof(double));
4048    memcpy(model2->rowLower(),saveLower+numberColumns,numberRows*sizeof(double));
4049    delete [] saveLower;
4050    memcpy(model2->columnUpper(),saveUpper,numberColumns*sizeof(double));
4051    memcpy(model2->rowUpper(),saveUpper+numberColumns,numberRows*sizeof(double));
4052    delete [] saveUpper;
4053    model2->primal(1);
4054  }
4055  model2->setPerturbation(savePerturbation);
4056  return model2->status();
4057}
4058/* For strong branching.  On input lower and upper are new bounds
4059   while on output they are objective function values (>1.0e50 infeasible).
4060   Return code is 0 if nothing interesting, -1 if infeasible both
4061   ways and +1 if infeasible one way (check values to see which one(s))
4062*/
4063int ClpSimplex::strongBranching(int numberVariables,const int * variables,
4064                                double * newLower, double * newUpper,
4065                                double ** outputSolution,
4066                                int * outputStatus, int * outputIterations,
4067                                bool stopOnFirstInfeasible,
4068                                bool alwaysFinish,
4069                                int startFinishOptions)
4070{
4071  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
4072                                                    newLower,  newUpper,outputSolution,
4073                                                    outputStatus, outputIterations,
4074                                                    stopOnFirstInfeasible,
4075                                                    alwaysFinish,startFinishOptions);
4076}
4077/* Borrow model.  This is so we dont have to copy large amounts
4078   of data around.  It assumes a derived class wants to overwrite
4079   an empty model with a real one - while it does an algorithm.
4080   This is same as ClpModel one, but sets scaling on etc. */
4081void 
4082ClpSimplex::borrowModel(ClpModel & otherModel) 
4083{
4084  ClpModel::borrowModel(otherModel);
4085  createStatus();
4086  //ClpDualRowSteepest steep1;
4087  //setDualRowPivotAlgorithm(steep1);
4088  //ClpPrimalColumnSteepest steep2;
4089  //setPrimalColumnPivotAlgorithm(steep2);
4090}
4091void 
4092ClpSimplex::borrowModel(ClpSimplex & otherModel) 
4093{
4094  ClpModel::borrowModel(otherModel);
4095  createStatus();
4096  dualBound_ = otherModel.dualBound_;
4097  dualTolerance_ = otherModel.dualTolerance_;
4098  primalTolerance_ = otherModel.primalTolerance_;
4099  delete dualRowPivot_;
4100  dualRowPivot_ = otherModel.dualRowPivot_->clone(false);
4101  delete primalColumnPivot_;
4102  primalColumnPivot_ = otherModel.primalColumnPivot_->clone(false);
4103  perturbation_ = otherModel.perturbation_;
4104  specialOptions_ = otherModel.specialOptions_;
4105  automaticScale_ = otherModel.automaticScale_;
4106}
4107/// Saves scalars for ClpSimplex
4108typedef struct {
4109  double optimizationDirection;
4110  double dblParam[ClpLastDblParam];
4111  double objectiveValue;
4112  double dualBound;
4113  double dualTolerance;
4114  double primalTolerance;
4115  double sumDualInfeasibilities;
4116  double sumPrimalInfeasibilities;
4117  double infeasibilityCost;
4118  int numberRows;
4119  int numberColumns;
4120  int intParam[ClpLastIntParam];
4121  int numberIterations;
4122  int problemStatus;
4123  int maximumIterations;
4124  int lengthNames;
4125  int numberDualInfeasibilities;
4126  int numberDualInfeasibilitiesWithoutFree;
4127  int numberPrimalInfeasibilities;
4128  int numberRefinements;
4129  int scalingFlag;
4130  int algorithm;
4131  unsigned int specialOptions;
4132  int dualPivotChoice;
4133  int primalPivotChoice;
4134  int matrixStorageChoice;
4135} Clp_scalars;
4136
4137int outDoubleArray(double * array, int length, FILE * fp)
4138{
4139  int numberWritten;
4140  if (array&&length) {
4141    numberWritten = fwrite(&length,sizeof(int),1,fp);
4142    if (numberWritten!=1)
4143      return 1;
4144    numberWritten = fwrite(array,sizeof(double),length,fp);
4145    if (numberWritten!=length)
4146      return 1;
4147  } else {
4148    length = 0;
4149    numberWritten = fwrite(&length,sizeof(int),1,fp);
4150    if (numberWritten!=1)
4151      return 1;
4152  }
4153  return 0;
4154}
4155// Save model to file, returns 0 if success
4156int
4157ClpSimplex::saveModel(const char * fileName)
4158{
4159  FILE * fp = fopen(fileName,"wb");
4160  if (fp) {
4161    Clp_scalars scalars;
4162    int i;
4163    CoinBigIndex numberWritten;
4164    // Fill in scalars
4165    scalars.optimizationDirection = optimizationDirection_;
4166    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
4167    scalars.objectiveValue = objectiveValue_;
4168    scalars.dualBound = dualBound_;
4169    scalars.dualTolerance = dualTolerance_;
4170    scalars.primalTolerance = primalTolerance_;
4171    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
4172    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
4173    scalars.infeasibilityCost = infeasibilityCost_;
4174    scalars.numberRows = numberRows_;
4175    scalars.numberColumns = numberColumns_;
4176    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
4177    scalars.numberIterations = numberIterations_;
4178    scalars.problemStatus = problemStatus_;
4179    scalars.maximumIterations = maximumIterations();
4180    scalars.lengthNames = lengthNames_;
4181    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
4182    scalars.numberDualInfeasibilitiesWithoutFree
4183      = numberDualInfeasibilitiesWithoutFree_;
4184    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
4185    scalars.numberRefinements = numberRefinements_;
4186    scalars.scalingFlag = scalingFlag_;
4187    scalars.algorithm = algorithm_;
4188    scalars.specialOptions = specialOptions_;
4189    scalars.dualPivotChoice = dualRowPivot_->type();
4190    scalars.primalPivotChoice = primalColumnPivot_->type();
4191    scalars.matrixStorageChoice = matrix_->type();
4192
4193    // put out scalars
4194    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
4195    if (numberWritten!=1)
4196      return 1;
4197    // strings
4198    CoinBigIndex length;
4199    for (i=0;i<ClpLastStrParam;i++) {
4200      length = strParam_[i].size();
4201      numberWritten = fwrite(&length,sizeof(int),1,fp);
4202      if (numberWritten!=1)
4203        return 1;
4204      if (length) {
4205        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
4206        if (numberWritten!=1)
4207          return 1;
4208      }
4209    }
4210    // arrays - in no particular order
4211    if (outDoubleArray(rowActivity_,numberRows_,fp))
4212        return 1;
4213    if (outDoubleArray(columnActivity_,numberColumns_,fp))
4214        return 1;
4215    if (outDoubleArray(dual_,numberRows_,fp))
4216        return 1;
4217    if (outDoubleArray(reducedCost_,numberColumns_,fp))
4218        return 1;
4219    if (outDoubleArray(rowLower_,numberRows_,fp))
4220        return 1;
4221    if (outDoubleArray(rowUpper_,numberRows_,fp))
4222        return 1;
4223    if (outDoubleArray(objective(),numberColumns_,fp))
4224        return 1;
4225    if (outDoubleArray(rowObjective_,numberRows_,fp))
4226        return 1;
4227    if (outDoubleArray(columnLower_,numberColumns_,fp))
4228        return 1;
4229    if (outDoubleArray(columnUpper_,numberColumns_,fp))
4230        return 1;
4231    if (ray_) {
4232      if (problemStatus_==1)
4233        if (outDoubleArray(ray_,numberRows_,fp))
4234          return 1;
4235      else if (problemStatus_==2)
4236        if (outDoubleArray(ray_,numberColumns_,fp))
4237          return 1;
4238      else
4239        if (outDoubleArray(NULL,0,fp))
4240          return 1;
4241    } else {
4242      if (outDoubleArray(NULL,0,fp))
4243        return 1;
4244    }
4245    if (status_&&(numberRows_+numberColumns_)>0) {
4246      length = numberRows_+numberColumns_;
4247      numberWritten = fwrite(&length,sizeof(int),1,fp);
4248      if (numberWritten!=1)
4249        return 1;
4250      numberWritten = fwrite(status_,sizeof(char),length, fp);
4251      if (numberWritten!=length)
4252        return 1;
4253    } else {
4254      length = 0;
4255      numberWritten = fwrite(&length,sizeof(int),1,fp);
4256      if (numberWritten!=1)
4257        return 1;
4258    }
4259    if (lengthNames_) {
4260      char * array = 
4261        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4262      char * put = array;
4263      CoinAssert (numberRows_ == (int) rowNames_.size());
4264      for (i=0;i<numberRows_;i++) {
4265        assert((int)rowNames_[i].size()<=lengthNames_);
4266        strcpy(put,rowNames_[i].c_str());
4267        put += lengthNames_+1;
4268      }
4269      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
4270      if (numberWritten!=numberRows_)
4271        return 1;
4272      put=array;
4273      CoinAssert (numberColumns_ == (int) columnNames_.size());
4274      for (i=0;i<numberColumns_;i++) {
4275        assert((int) columnNames_[i].size()<=lengthNames_);
4276        strcpy(put,columnNames_[i].c_str());
4277        put += lengthNames_+1;
4278      }
4279      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
4280      if (numberWritten!=numberColumns_)
4281        return 1;
4282      delete [] array;
4283    }
4284    // just standard type at present
4285    assert (matrix_->type()==1);
4286    CoinAssert (matrix_->getNumCols() == numberColumns_);
4287    CoinAssert (matrix_->getNumRows() == numberRows_);
4288    // we are going to save with gaps
4289    length = matrix_->getVectorStarts()[numberColumns_-1]
4290      + matrix_->getVectorLengths()[numberColumns_-1];
4291    numberWritten = fwrite(&length,sizeof(int),1,fp);
4292    if (numberWritten!=1)
4293      return 1;
4294    numberWritten = fwrite(matrix_->getElements(),
4295                           sizeof(double),length,fp);
4296    if (numberWritten!=length)
4297      return 1;
4298    numberWritten = fwrite(matrix_->getIndices(),
4299                           sizeof(int),length,fp);
4300    if (numberWritten!=length)
4301      return 1;
4302    numberWritten = fwrite(matrix_->getVectorStarts(),
4303                           sizeof(int),numberColumns_+1,fp);
4304    if (numberWritten!=numberColumns_+1)
4305      return 1;
4306    numberWritten = fwrite(matrix_->getVectorLengths(),
4307                           sizeof(int),numberColumns_,fp);
4308    if (numberWritten!=numberColumns_)
4309      return 1;
4310    // finished
4311    fclose(fp);
4312    return 0;
4313  } else {
4314    return -1;
4315  }
4316}
4317
4318int inDoubleArray(double * &array, int length, FILE * fp)
4319{
4320  int numberRead;
4321  int length2;
4322  numberRead = fread(&length2,sizeof(int),1,fp);
4323  if (numberRead!=1)
4324    return 1;
4325  if (length2) {
4326    // lengths must match
4327    if (length!=length2)
4328      return 2;
4329    array = new double[length];
4330    numberRead = fread(array,sizeof(double),length,fp);
4331    if (numberRead!=length)
4332      return 1;
4333  } 
4334  return 0;
4335}
4336/* Restore model from file, returns 0 if success,
4337   deletes current model */
4338int 
4339ClpSimplex::restoreModel(const char * fileName)
4340{
4341  FILE * fp = fopen(fileName,"rb");
4342  if (fp) {
4343    // Get rid of current model
4344    // save event handler in case already set
4345    ClpEventHandler * handler = eventHandler_->clone();
4346    ClpModel::gutsOfDelete();
4347    eventHandler_ = handler;
4348    gutsOfDelete(0);
4349    int i;
4350    for (i=0;i<6;i++) {
4351      rowArray_[i]=NULL;
4352      columnArray_[i]=NULL;
4353    }
4354    // get an empty factorization so we can set tolerances etc
4355    factorization_ = new ClpFactorization();
4356    // Say sparse
4357    factorization_->sparseThreshold(1);
4358    Clp_scalars scalars;
4359    CoinBigIndex numberRead;
4360
4361    // get scalars
4362    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
4363    if (numberRead!=1)
4364      return 1;
4365    // Fill in scalars
4366    optimizationDirection_ = scalars.optimizationDirection;
4367    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
4368    objectiveValue_ = scalars.objectiveValue;
4369    dualBound_ = scalars.dualBound;
4370    dualTolerance_ = scalars.dualTolerance;
4371    primalTolerance_ = scalars.primalTolerance;
4372    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
4373    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
4374    infeasibilityCost_ = scalars.infeasibilityCost;
4375    numberRows_ = scalars.numberRows;
4376    numberColumns_ = scalars.numberColumns;
4377    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
4378    numberIterations_ = scalars.numberIterations;
4379    problemStatus_ = scalars.problemStatus;
4380    setMaximumIterations(scalars.maximumIterations);
4381    lengthNames_ = scalars.lengthNames;
4382    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
4383    numberDualInfeasibilitiesWithoutFree_
4384      = scalars.numberDualInfeasibilitiesWithoutFree;
4385    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
4386    numberRefinements_ = scalars.numberRefinements;
4387    scalingFlag_ = scalars.scalingFlag;
4388    algorithm_ = scalars.algorithm;
4389    specialOptions_ = scalars.specialOptions;
4390    // strings
4391    CoinBigIndex length;
4392    for (i=0;i<ClpLastStrParam;i++) {
4393      numberRead = fread(&length,sizeof(int),1,fp);
4394      if (numberRead!=1)
4395        return 1;
4396      if (length) {
4397        char * array = new char[length+1];
4398        numberRead = fread(array,length,1,fp);
4399        if (numberRead!=1)
4400          return 1;
4401        array[length]='\0';
4402        strParam_[i]=array;
4403        delete [] array;
4404      }
4405    }
4406    // arrays - in no particular order
4407    if (inDoubleArray(rowActivity_,numberRows_,fp))
4408        return 1;
4409    if (inDoubleArray(columnActivity_,numberColumns_,fp))
4410        return 1;
4411    if (inDoubleArray(dual_,numberRows_,fp))
4412        return 1;
4413    if (inDoubleArray(reducedCost_,numberColumns_,fp))
4414        return 1;
4415    if (inDoubleArray(rowLower_,numberRows_,fp))
4416        return 1;
4417    if (inDoubleArray(rowUpper_,numberRows_,fp))
4418        return 1;
4419    double * objective;
4420    if (inDoubleArray(objective,numberColumns_,fp))
4421        return 1;
4422    delete objective_;
4423    objective_ = new ClpLinearObjective(objective,numberColumns_);
4424    delete [] objective;
4425    if (inDoubleArray(rowObjective_,numberRows_,fp))
4426        return 1;
4427    if (inDoubleArray(columnLower_,numberColumns_,fp))
4428        return 1;
4429    if (inDoubleArray(columnUpper_,numberColumns_,fp))
4430        return 1;
4431    if (problemStatus_==1) {
4432      if (inDoubleArray(ray_,numberRows_,fp))
4433        return 1;
4434    } else if (problemStatus_==2) {
4435      if (inDoubleArray(ray_,numberColumns_,fp))
4436        return 1;
4437    } else {
4438      // ray should be null
4439      numberRead = fread(&length,sizeof(int),1,fp);
4440      if (numberRead!=1)
4441        return 1;
4442      if (length)
4443        return 2;
4444    }
4445    delete [] status_;
4446    status_=NULL;
4447    // status region
4448    numberRead = fread(&length,sizeof(int),1,fp);
4449    if (numberRead!=1)
4450        return 1;
4451    if (length) {
4452      if (length!=numberRows_+numberColumns_)
4453        return 1;
4454      status_ = new char unsigned[length];
4455      numberRead = fread(status_,sizeof(char),length, fp);
4456      if (numberRead!=length)
4457        return 1;
4458    }
4459    if (lengthNames_) {
4460      char * array = 
4461        new char[CoinMax(numberRows_,numberColumns_)*(lengthNames_+1)];
4462      char * get = array;
4463      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
4464      if (numberRead!=numberRows_)
4465        return 1;
4466      rowNames_ = std::vector<std::string> ();
4467      rowNames_.resize(numberRows_);
4468      for (i=0;i<numberRows_;i++) {
4469        rowNames_.push_back(get);
4470        get += lengthNames_+1;
4471      }
4472      get = array;
4473      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
4474      if (numberRead!=numberColumns_)
4475        return 1;
4476      columnNames_ = std::vector<std::string> ();
4477      columnNames_.resize(numberColumns_);
4478      for (i=0;i<numberColumns_;i++) {
4479        columnNames_.push_back(get);
4480        get += lengthNames_+1;
4481      }
4482      delete [] array;
4483    }
4484    // Pivot choices
4485    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
4486    delete dualRowPivot_;
4487    switch ((scalars.dualPivotChoice&63)) {
4488    default:
4489      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
4490    case 1:
4491      // Dantzig
4492      dualRowPivot_ = new ClpDualRowDantzig();
4493      break;
4494    case 2:
4495      // Steepest - use mode
4496      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
4497      break;
4498    }
4499    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
4500    delete primalColumnPivot_;
4501    switch ((scalars.primalPivotChoice&63)) {
4502    default:
4503      printf("Need another primalPivot case %d\n",
4504             scalars.primalPivotChoice&63);
4505    case 1:
4506      // Dantzig
4507      primalColumnPivot_ = new ClpPrimalColumnDantzig();
4508      break;
4509    case 2:
4510      // Steepest - use mode
4511      primalColumnPivot_
4512        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
4513      break;
4514    }
4515    assert(scalars.matrixStorageChoice==1);
4516    delete matrix_;
4517    // get arrays
4518    numberRead = fread(&length,sizeof(int),1,fp);
4519    if (numberRead!=1)
4520      return 1;
4521    double * elements = new double[length];
4522    int * indices = new int[length];
4523    CoinBigIndex * starts = new CoinBigIndex[numberColumns_+1];
4524    int * lengths = new int[numberColumns_];
4525    numberRead = fread(elements, sizeof(double),length,fp);
4526    if (numberRead!=length)
4527      return 1;
4528    numberRead = fread(indices, sizeof(int),length,fp);
4529    if (numberRead!=length)
4530      return 1;
4531    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
4532    if (numberRead!=numberColumns_+1)
4533      return 1;
4534    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
4535    if (numberRead!=numberColumns_)
4536      return 1;
4537    // assign matrix
4538    CoinPackedMatrix * matrix = new CoinPackedMatrix();
4539    matrix->assignMatrix(true, numberRows_, numberColumns_,
4540                         length, elements, indices, starts, lengths);
4541    // and transfer to Clp
4542    matrix_ = new ClpPackedMatrix(matrix);
4543    // finished
4544    fclose(fp);
4545    return 0;
4546  } else {
4547    return -1;
4548  }
4549  return 0;
4550}
4551// value of incoming variable (in Dual)
4552double 
4553ClpSimplex::valueIncomingDual() const
4554{
4555  // Need value of incoming for list of infeasibilities as may be infeasible
4556  double valueIncoming = (dualOut_/alpha_)*directionOut_;
4557  if (directionIn_==-1)
4558    valueIncoming = upperIn_-valueIncoming;
4559  else
4560    valueIncoming = lowerIn_-valueIncoming;
4561  return valueIncoming;
4562}
4563// Sanity check on input data - returns true if okay
4564bool 
4565ClpSimplex::sanityCheck()
4566{
4567  // bad if empty
4568  if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
4569    int infeasNumber[2];
4570    double infeasSum[2];
4571    problemStatus_=emptyProblem(infeasNumber,infeasSum,false);
4572    numberDualInfeasibilities_=infeasNumber[0];
4573    sumDualInfeasibilities_=infeasSum[0];
4574    numberPrimalInfeasibilities_=infeasNumber[1];
4575    sumPrimalInfeasibilities_=infeasSum[1];
4576    return false;
4577  }
4578  int numberBad ;
4579  double largestBound, smallestBound, minimumGap;
4580  double smallestObj, largestObj;
4581  int firstBad;
4582  int modifiedBounds=0;
4583  int i;
4584  numberBad=0;
4585  firstBad=-1;
4586  minimumGap=1.0e100;
4587  smallestBound=1.0e100;
4588  largestBound=0.0;
4589  smallestObj=1.0e100;
4590  largestObj=0.0;
4591  // If bounds are too close - fix
4592  double fixTolerance = 10.0*primalTolerance_;
4593  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
4594    double value;
4595    value = fabs(cost_[i]);
4596    if (value>1.0e50) {
4597      numberBad++;
4598      if (firstBad<0)
4599        firstBad=i;
4600    } else if (value) {
4601      if (value>largestObj)
4602        largestObj=value;
4603      if (value<smallestObj)
4604        smallestObj=value;
4605    }
4606    value=upper_[i]-lower_[i];
4607    if (value<-primalTolerance_) {
4608      numberBad++;
4609      if (firstBad<0)
4610        firstBad=i;
4611    } else if (value<=fixTolerance) {
4612      if (value) {
4613        // modify
4614        upper_[i] = lower_[i];
4615        modifiedBounds++;
4616      }
4617    } else {
4618      if (value<minimumGap)
4619        minimumGap=value;
4620    }
4621    if (lower_[i]>-1.0e100&&lower_[i]) {
4622      value = fabs(lower_[i]);
4623      if (value>largestBound)
4624        largestBound=value;
4625      if (value<smallestBound)
4626        smallestBound=value;
4627    }
4628    if (upper_[i]<1.0e100&&upper_[i]) {
4629      value = fabs(upper_[i]);
4630      if (value>largestBound)
4631        largestBound=value;
4632      if (value<smallestBound)
4633        smallestBound=value;
4634    }
4635  }
4636  if (largestBound)
4637    handler_->message(CLP_RIMSTATISTICS3,messages_)
4638      <<smallestBound
4639      <<largestBound
4640      <<minimumGap
4641      <<CoinMessageEol;
4642  minimumGap=1.0e100;
4643  smallestBound=1.0e100;
4644  largestBound=0.0;
4645  for (i=0;i<numberColumns_;i++) {
4646    double value;
4647    value = fabs(cost_[i]);
4648    if (value>1.0e50) {
4649      numberBad++;
4650      if (firstBad<0)
4651        firstBad=i;
4652    } else if (value) {
4653      if (value>largestObj)
4654        largestObj=value;
4655      if (value<smallestObj)
4656        smallestObj=value;
4657    }
4658    value=upper_[i]-lower_[i];
4659    if (value<-primalTolerance_) {
4660      numberBad++;
4661      if (firstBad<0)
4662        firstBad=i;
4663    } else if (value<=fixTolerance) {
4664      if (value) {
4665        // modify
4666        upper_[i] = lower_[i];
4667        modifiedBounds++;
4668      }
4669    } else {
4670      if (value<minimumGap)
4671        minimumGap=value;
4672    }
4673    if (lower_[i]>-1.0e100&&lower_[i]) {
4674      value = fabs(lower_[i]);
4675      if (value>largestBound)
4676        largestBound=value;
4677      if (value<smallestBound)
4678        smallestBound=value;
4679    }
4680    if (upper_[i]<1.0e100&&upper_[i]) {
4681      value = fabs(upper_[i]);
4682      if (value>largestBound)
4683        largestBound=value;
4684      if (value<smallestBound)
4685        smallestBound=value;
4686    }
4687  }
4688  char rowcol[]={'R','C'};
4689  if (numberBad) {
4690    handler_->message(CLP_BAD_BOUNDS,messages_)
4691      <<numberBad
4692      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
4693      <<CoinMessageEol;
4694    problemStatus_=4;
4695    return false;
4696  }
4697  if (modifiedBounds)
4698    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
4699      <<modifiedBounds
4700      <<CoinMessageEol;
4701  handler_->message(CLP_RIMSTATISTICS1,messages_)
4702    <<smallestObj
4703    <<largestObj
4704    <<CoinMessageEol;  if (largestBound)
4705    handler_->message(CLP_RIMSTATISTICS2,messages_)
4706      <<smallestBound
4707      <<largestBound
4708      <<minimumGap
4709      <<CoinMessageEol;
4710  return true;
4711}
4712// Set up status array (for OsiClp)
4713void 
4714ClpSimplex::createStatus() 
4715{
4716  if(!status_)
4717    status_ = new unsigned char [numberColumns_+numberRows_];
4718  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
4719  int i;
4720  // set column status to one nearest zero
4721  for (i=0;i<numberColumns_;i++) {
4722#if 0
4723    if (columnLower_[i]>=0.0) {
4724      setColumnStatus(i,atLowerBound);
4725    } else if (columnUpper_[i]<=0.0) {
4726      setColumnStatus(i,atUpperBound);
4727    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4728      // free
4729      setColumnStatus(i,isFree);
4730    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4731      setColumnStatus(i,atLowerBound);
4732    } else {
4733      setColumnStatus(i,atUpperBound);
4734    }
4735#else
4736    setColumnStatus(i,atLowerBound);
4737#endif
4738  }
4739  for (i=0;i<numberRows_;i++) {
4740    setRowStatus(i,basic);
4741  }
4742}
4743/* Sets up all slack basis and resets solution to
4744   as it was after initial load or readMps */
4745void ClpSimplex::allSlackBasis(bool resetSolution)
4746{
4747  createStatus();
4748  if (resetSolution) {
4749    // put back to as it was originally
4750    int i;
4751    // set column status to one nearest zero
4752    // But set value to zero if lb <0.0 and ub>0.0
4753    for (i=0;i<numberColumns_;i++) {
4754      if (columnLower_[i]>=0.0) {
4755        columnActivity_[i]=columnLower_[i];
4756        setColumnStatus(i,atLowerBound);
4757      } else if (columnUpper_[i]<=0.0) {
4758        columnActivity_[i]=columnUpper_[i];
4759        setColumnStatus(i,atUpperBound);
4760      } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
4761        // free
4762        columnActivity_[i]=0.0;
4763        setColumnStatus(i,isFree);
4764      } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
4765        columnActivity_[i]=0.0;
4766        setColumnStatus(i,atLowerBound);
4767      } else {
4768        columnActivity_[i]=0.0;
4769        setColumnStatus(i,atUpperBound);
4770      }
4771    }
4772  }
4773}
4774/* Loads a problem (the constraints on the
4775   rows are given by lower and upper bounds). If a pointer is 0 then the
4776   following values are the default:
4777   <ul>
4778   <li> <code>colub</code>: all columns have upper bound infinity
4779   <li> <code>collb</code>: all columns have lower bound 0
4780   <li> <code>rowub</code>: all rows have upper bound infinity
4781   <li> <code>rowlb</code>: all rows have lower bound -infinity
4782   <li> <code>obj</code>: all variables have 0 objective coefficient
4783   </ul>
4784*/
4785void 
4786ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
4787                    const double* collb, const double* colub,   
4788                    const double* obj,
4789                    const double* rowlb, const double* rowub,
4790                    const double * rowObjective)
4791{
4792  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4793                        rowObjective);
4794  createStatus();
4795}
4796void 
4797ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
4798                    const double* collb, const double* colub,   
4799                    const double* obj,
4800                    const double* rowlb, const double* rowub,
4801                    const double * rowObjective)
4802{
4803  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
4804                        rowObjective);
4805  createStatus();
4806}
4807
4808/* Just like the other loadProblem() method except that the matrix is
4809   given in a standard column major ordered format (without gaps). */
4810void 
4811ClpSimplex::loadProblem (  const int numcols, const int numrows,
4812                    const CoinBigIndex* start, const int* index,
4813                    const double* value,
4814                    const double* collb, const double* colub,   
4815                    const double* obj,
4816                    const double* rowlb, const double* rowub,
4817                    const double * rowObjective)
4818{
4819  ClpModel::loadProblem(numcols, numrows, start, index, value,
4820                          collb, colub, obj, rowlb, rowub,
4821                          rowObjective);
4822  createStatus();
4823}
4824// This loads a model from a coinModel object - returns number of errors
4825int 
4826ClpSimplex::loadProblem (  CoinModel & modelObject, bool keepSolution)
4827{
4828  unsigned char * status = NULL;
4829  double * psol = NULL;
4830  double * dsol = NULL;
4831  if (status_&&numberRows_&&numberRows_==modelObject.numberRows()&&
4832      numberColumns_==modelObject.numberColumns()) {
4833    status = new unsigned char [numberRows_+numberColumns_];
4834    memcpy(status,status_,numberRows_+numberColumns_);
4835    psol = new double [numberRows_+numberColumns_];
4836    memcpy(psol,columnActivity_,numberColumns_*sizeof(double));
4837    memcpy(psol+numberColumns_,rowActivity_,numberRows_*sizeof(double));
4838    dsol = new double [numberRows_+numberColumns_];
4839    memcpy(dsol,reducedCost_,numberColumns_*sizeof(double));
4840    memcpy(dsol+numberColumns_,dual_,numberRows_*sizeof(double));
4841  }
4842  int returnCode = ClpModel::loadProblem(modelObject);
4843  createStatus();
4844  if (status) {
4845    // copy back
4846    memcpy(status_,status,numberRows_+numberColumns_);
4847    memcpy(columnActivity_,psol,numberColumns_*sizeof(double));
4848    memcpy(rowActivity_,psol+numberColumns_,numberRows_*sizeof(double));
4849    memcpy(reducedCost_,dsol,numberColumns_*sizeof(double));
4850    memcpy(dual_,dsol+numberColumns_,numberRows_*sizeof(double));
4851    delete [] status;
4852    delete [] psol;
4853    delete [] dsol;
4854  }
4855  return returnCode;
4856}
4857void 
4858ClpSimplex::loadProblem (  const int numcols, const int numrows,
4859                           const CoinBigIndex* start, const int* index,
4860                           const double* value,const int * length,
4861                           const double* collb, const double* colub,   
4862                           const double* obj,
4863                           const double* rowlb, const double* rowub,
4864                           const double * rowObjective)
4865{
4866  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
4867                          collb, colub, obj, rowlb, rowub,
4868                          rowObjective);
4869  createStatus();
4870}
4871// Read an mps file from the given filename
4872int 
4873ClpSimplex::readMps(const char *filename,
4874            bool keepNames,
4875            bool ignoreErrors)
4876{
4877  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
4878  createStatus();
4879  return status;
4880}
4881// Just check solution (for external use)
4882void 
4883ClpSimplex::checkSolution()
4884{
4885  // put in standard form
4886  createRim(7+8+16+32);
4887  dualTolerance_=dblParam_[ClpDualTolerance];
4888  primalTolerance_=dblParam_[ClpPrimalTolerance];
4889  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
4890  checkDualSolution();
4891  if (!numberDualInfeasibilities_&&
4892      !numberPrimalInfeasibilities_)
4893    problemStatus_=0;
4894  else
4895    problemStatus_=-1;
4896#ifdef CLP_DEBUG
4897  int i;
4898  double value=0.0;
4899  for (i=0;i<numberRows_+numberColumns_;i++)
4900    value += dj_[i]*solution_[i];
4901  printf("dual value %g, primal %g\n",value,objectiveValue());
4902#endif
4903  // release extra memory
4904  deleteRim(0);
4905}
4906/* Crash - at present just aimed at dual, returns
4907   -2 if dual preferred and crash basis created
4908   -1 if dual preferred and all slack basis preferred
4909   0 if basis going in was not all slack
4910   1 if primal preferred and all slack basis preferred
4911   2 if primal preferred and crash basis created.
4912   
4913   if gap between bounds <="gap" variables can be flipped
4914   ( If pivot -1 then can be made super basic!)
4915   
4916   If "pivot" is
4917   -1 No pivoting - always primal
4918   0 No pivoting (so will just be choice of algorithm)
4919   1 Simple pivoting e.g. gub
4920   2 Mini iterations
4921*/
4922int 
4923ClpSimplex::crash(double gap,int pivot)
4924{
4925  CoinAssert(!rowObjective_); // not coded
4926  int iColumn;
4927  int numberBad=0;
4928  int numberBasic=0;
4929  double dualTolerance=dblParam_[ClpDualTolerance];
4930  //double primalTolerance=dblParam_[ClpPrimalTolerance];
4931  int returnCode=0;
4932  // If no basis then make all slack one
4933  if (!status_)
4934    createStatus();
4935 
4936  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4937    if (getColumnStatus(iColumn)==basic)
4938      numberBasic++;
4939  }
4940  if (!numberBasic) {
4941    // all slack
4942    double * dj = new double [numberColumns_];
4943    double * solution = columnActivity_;
4944    const double * linearObjective = objective();
4945    //double objectiveValue=0.0;
4946    int iColumn;
4947    double direction = optimizationDirection_;
4948    // direction is actually scale out not scale in
4949    if (direction)
4950      direction = 1.0/direction;
4951    for (iColumn=0;iColumn<numberColumns_;iColumn++)
4952      dj[iColumn] = direction*linearObjective[iColumn];
4953    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4954      // assume natural place is closest to zero
4955      double lowerBound = columnLower_[iColumn];
4956      double upperBound = columnUpper_[iColumn];
4957      if (lowerBound>-1.0e20||upperBound<1.0e20) {
4958        bool atLower;
4959        if (fabs(upperBound)<fabs(lowerBound)) {
4960          atLower=false;
4961          setColumnStatus(iColumn,atUpperBound);
4962          solution[iColumn]=upperBound;
4963        } else {
4964          atLower=true;
4965          setColumnStatus(iColumn,atLowerBound);
4966          solution[iColumn]=lowerBound;
4967        }
4968        if (dj[iColumn]<-dualTolerance_) {
4969          // should be at upper bound
4970          if (atLower) {
4971            // can we flip
4972            if (upperBound-lowerBound<=gap) {
4973              columnActivity_[iColumn]=upperBound;
4974              setColumnStatus(iColumn,atUpperBound);
4975            } else if (pivot<0) {
4976              // set superbasic
4977              columnActivity_[iColumn]=lowerBound+gap;
4978              setColumnStatus(iColumn,superBasic);
4979            } else if (dj[iColumn]<-dualTolerance) {
4980              numberBad++;
4981            }
4982          }
4983        } else if (dj[iColumn]>dualTolerance_) {
4984          // should be at lower bound
4985          if (!atLower) {
4986            // can we flip
4987            if (upperBound-lowerBound<=gap) {
4988              columnActivity_[iColumn]=lowerBound;
4989              setColumnStatus(iColumn,atLowerBound);
4990            } else if (pivot<0) {
4991              // set superbasic
4992              columnActivity_[iColumn]=upperBound-gap;
4993              setColumnStatus(iColumn,superBasic);
4994            } else if (dj[iColumn]>dualTolerance) {
4995              numberBad++;
4996            }
4997          }
4998        }
4999      } else {
5000        // free
5001        setColumnStatus(iColumn,isFree);
5002        if (fabs(dj[iColumn])>dualTolerance) 
5003          numberBad++;
5004      }
5005    }
5006    if (numberBad||pivot) {
5007      if (pivot<=0) {
5008        delete [] dj;
5009        returnCode = 1;
5010      } else {
5011        // see if can be made dual feasible with gubs etc
5012        double * pi = new double[numberRows_];
5013        memset (pi,0,numberRows_*sizeof(double));
5014        int * way = new int[numberColumns_];
5015        int numberIn = 0;
5016
5017        // Get column copy
5018        CoinPackedMatrix * columnCopy = matrix();
5019        // Get a row copy in standard format
5020        CoinPackedMatrix copy;
5021        copy.reverseOrderedCopyOf(*columnCopy);
5022        // get matrix data pointers
5023        const int * column = copy.getIndices();
5024        const CoinBigIndex * rowStart = copy.getVectorStarts();
5025        const int * rowLength = copy.getVectorLengths(); 
5026        const double * elementByRow = copy.getElements();
5027        //const int * row = columnCopy->getIndices();
5028        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
5029        //const int * columnLength = columnCopy->getVectorLengths();
5030        //const double * element = columnCopy->getElements();
5031
5032
5033        // if equality row and bounds mean artificial in basis bad
5034        // then do anyway
5035
5036        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5037          // - if we want to reduce dj, + if we want to increase
5038          int thisWay = 100;
5039          double lowerBound = columnLower_[iColumn];
5040          double upperBound = columnUpper_[iColumn];
5041          if (upperBound>lowerBound) {
5042            switch(getColumnStatus(iColumn)) {
5043             
5044            case basic:
5045              thisWay=0;
5046            case ClpSimplex::isFixed:
5047              break;
5048            case isFree:
5049            case superBasic:
5050              if (dj[iColumn]<-dualTolerance) 
5051                thisWay = 1;
5052              else if (dj[iColumn]>dualTolerance) 
5053                thisWay = -1;
5054              else
5055                thisWay =0;
5056              break;
5057            case atUpperBound:
5058              if (dj[iColumn]>dualTolerance) 
5059                thisWay = -1;
5060              else if (dj[iColumn]<-dualTolerance) 
5061                thisWay = -3;
5062              else
5063                thisWay = -2;
5064              break;
5065            case atLowerBound:
5066              if (dj[iColumn]<-dualTolerance) 
5067                thisWay = 1;
5068              else if (dj[iColumn]>dualTolerance) 
5069                thisWay = 3;
5070              else
5071                thisWay = 2;
5072              break;
5073            }
5074          }
5075          way[iColumn] = thisWay;
5076        }
5077        /*if (!numberBad)
5078          printf("Was dual feasible before passes - rows %d\n",
5079          numberRows_);*/
5080        int lastNumberIn = -100000;
5081        int numberPasses=5;
5082        while (numberIn>lastNumberIn+numberRows_/100) {
5083          lastNumberIn = numberIn;
5084          // we need to maximize chance of doing good
5085          int iRow;
5086          for (iRow=0;iRow<numberRows_;iRow++) {
5087            double lowerBound = rowLower_[iRow];
5088            double upperBound = rowUpper_[iRow];
5089            if (getRowStatus(iRow)==basic) {
5090              // see if we can find a column to pivot on
5091              int j;
5092              // down is amount pi can go down
5093              double maximumDown = COIN_DBL_MAX;
5094              double maximumUp = COIN_DBL_MAX;
5095              double minimumDown =0.0;
5096              double minimumUp =0.0;
5097              int iUp=-1;
5098              int iDown=-1;
5099              int iUpB=-1;
5100              int iDownB=-1;
5101              if (lowerBound<-1.0e20)
5102                maximumUp = -1.0;
5103              if (upperBound>1.0e20)
5104                maximumDown = -1.0;
5105              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
5106                int iColumn = column[j];
5107                double value = elementByRow[j];
5108                double djValue = dj[iColumn];
5109                /* way -
5110                   -3 - okay at upper bound with negative dj
5111                   -2 - marginal at upper bound with zero dj - can only decrease
5112                   -1 - bad at upper bound
5113                   0 - we can never pivot on this row
5114                   1 - bad at lower bound
5115                   2 - marginal at lower bound with zero dj - can only increase
5116                   3 - okay at lower bound with positive dj
5117                   100 - fine we can just ignore
5118                */
5119                if (way[iColumn]!=100) {
5120                  switch(way[iColumn]) {
5121                   
5122                  case -3:
5123                    if (value>0.0) {
5124                      if (maximumDown*value>-djValue) {
5125                        maximumDown = -djValue/value;
5126                        iDown = iColumn;
5127                      }
5128                    } else {
5129                      if (-maximumUp*value>-djValue) {
5130                        maximumUp = djValue/value;
5131                        iUp = iColumn;
5132                      }
5133                    }
5134                    break;
5135                  case -2:
5136                    if (value>0.0) {
5137                      maximumDown = 0.0;
5138                    } else {
5139                      maximumUp = 0.0;
5140                    }
5141                    break;
5142                  case -1:
5143                    // see if could be satisfied
5144                    // dj value > 0
5145                    if (value>0.0) {
5146                      maximumDown=0.0;
5147                      if (maximumUp*value<djValue-dualTolerance) {
5148                        maximumUp = 0.0; // would improve but not enough
5149                      } else {
5150                        if (minimumUp*value<djValue) {
5151                          minimumUp = djValue/value;
5152                          iUpB = iColumn;
5153                        }
5154                      }
5155                    } else {
5156                      maximumUp=0.0;
5157                      if (-maximumDown*value<djValue-dualTolerance) {
5158                        maximumDown = 0.0; // would improve but not enough
5159                      } else {
5160                        if (-minimumDown*value<djValue) {
5161                          minimumDown = -djValue/value;
5162                          iDownB = iColumn;
5163                        }
5164                      }
5165                    }
5166                   
5167                    break;
5168                  case 0:
5169                    maximumDown = -1.0;
5170                    maximumUp=-1.0;
5171                    break;
5172                  case 1:
5173                    // see if could be satisfied
5174                    // dj value < 0
5175                    if (value>0.0) {
5176                      maximumUp=0.0;
5177                      if (maximumDown*value<-djValue-dualTolerance) {
5178                        maximumDown = 0.0; // would improve but not enough
5179                      } else {
5180                        if (minimumDown*value<-djValue) {
5181                          minimumDown = -djValue/value;
5182                          iDownB = iColumn;
5183                        }
5184                      }
5185                    } else {
5186                      maximumDown=0.0;
5187                      if (-maximumUp*value<-djValue-dualTolerance) {
5188                        maximumUp = 0.0; // would improve but not enough
5189                      } else {
5190                        if (-minimumUp*value<-djValue) {
5191                          minimumUp = djValue/value;
5192                          iUpB = iColumn;
5193                        }
5194                      }
5195                    }
5196                   
5197                    break;
5198                  case 2:
5199                    if (value>0.0) {
5200                      maximumUp = 0.0;
5201                    } else {
5202                      maximumDown = 0.0;
5203                    }
5204                   
5205                    break;
5206                  case 3:
5207                    if (value>0.0) {
5208                      if (maximumUp*value>djValue) {
5209                        maximumUp = djValue/value;
5210                        iUp = iColumn;
5211                      }
5212                    } else {
5213                      if (-maximumDown*value>djValue) {
5214                        maximumDown = -djValue/value;
5215                        iDown = iColumn;
5216                      }
5217                    }
5218                   
5219                    break;
5220                  default:
5221                    break;
5222                  }
5223                }
5224              }
5225              if (iUpB>=0)
5226                iUp=iUpB;
5227              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
5228                iUp=-1;
5229              if (iDownB>=0)
5230                iDown=iDownB;
5231              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
5232                iDown=-1;
5233              if (iUp>=0||iDown>=0) {
5234                // do something
5235                if (iUp>=0&&iDown>=0) {
5236                  if (maximumDown>maximumUp)
5237                    iUp=-1;
5238                }
5239                double change;
5240                int kColumn;
5241                if (iUp>=0) {
5242                  kColumn=iUp;
5243                  change=maximumUp;
5244                  // just do minimum if was dual infeasible
5245                  // ? only if maximum large?
5246                  if (minimumUp>0.0)
5247                    change=minimumUp;
5248                  setRowStatus(iRow,atUpperBound);
5249                } else {
5250                  kColumn=iDown;
5251                  change=-maximumDown;
5252                  // just do minimum if was dual infeasible
5253                  // ? only if maximum large?
5254                  if (minimumDown>0.0)
5255                    change=-minimumDown;
5256                  setRowStatus(iRow,atLowerBound);
5257                }
5258                assert (fabs(change)<1.0e20);
5259                setColumnStatus(kColumn,basic);
5260                numberIn++;
5261                pi[iRow]=change;
5262                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
5263                  int iColumn = column[j];
5264                  double value = elementByRow[j];
5265                  double djValue = dj[iColumn]-change*value;
5266                  dj[iColumn]=djValue;
5267                  if (abs(way[iColumn])==1) {
5268                    numberBad--;
5269                    /*if (!numberBad)
5270                      printf("Became dual feasible at row %d out of %d\n",
5271                      iRow, numberRows_);*/
5272                    lastNumberIn=-1000000;
5273                  }
5274                  int thisWay = 100;
5275                  double lowerBound = columnLower_[iColumn];
5276                  double upperBound = columnUpper_[iColumn];
5277                  if (upperBound>lowerBound) {
5278                    switch(getColumnStatus(iColumn)) {
5279                     
5280                    case basic:
5281                      thisWay=0;
5282                    case isFixed:
5283                      break;
5284                    case isFree:
5285                    case superBasic:
5286                      if (djValue<-dualTolerance) 
5287                        thisWay = 1;
5288                      else if (djValue>dualTolerance) 
5289                        thisWay = -1;
5290                      else
5291                        { thisWay =0; abort();}
5292                      break;
5293                    case atUpperBound:
5294                      if (djValue>dualTolerance) 
5295                        { thisWay =-1; abort();}
5296                      else if (djValue<-dualTolerance) 
5297                        thisWay = -3;
5298                      else
5299                        thisWay = -2;
5300                      break;
5301                    case atLowerBound:
5302                      if (djValue<-dualTolerance) 
5303                        { thisWay =1; abort();}
5304                      else if (djValue>dualTolerance) 
5305                        thisWay = 3;
5306                      else
5307                        thisWay = 2;
5308                      break;
5309                    }
5310                  }
5311                  way[iColumn] = thisWay;
5312                }
5313              }
5314            }
5315          }
5316          if (numberIn==lastNumberIn||numberBad||pivot<2)
5317            break;
5318          if (!(--numberPasses))
5319            break;
5320          //printf("%d put in so far\n",numberIn);
5321        }
5322        // last attempt to flip
5323        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
5324          double lowerBound = columnLower_[iColumn];
5325          double upperBound = columnUpper_[iColumn];
5326          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
5327            double djValue=dj[iColumn];
5328            switch(getColumnStatus(iColumn)) {
5329             
5330            case basic:
5331            case ClpSimplex::isFixed:
5332              break;
5333            case isFree:
5334            case superBasic:
5335              break;
5336            case atUpperBound:
5337              if (djValue>dualTolerance) {
5338                setColumnStatus(iColumn,atUpperBound);
5339                solution[iColumn]=upperBound;
5340              } 
5341              break;
5342            case atLowerBound:
5343              if (djValue<-dualTolerance) {
5344                setColumnStatus(iColumn,atUpperBound);
5345                solution[iColumn]=upperBound;
5346              }
5347              break;
5348            }
5349          }
5350        }
5351        delete [] pi;
5352        delete [] dj;
5353        delete [] way;
5354        handler_->message(CLP_CRASH,messages_)
5355          <<numberIn
5356          <<numberBad
5357          <<CoinMessageEol;
5358        returnCode =  -1;
5359      }
5360    } else {
5361      delete [] dj;
5362      returnCode =  -1;
5363    }
5364    //cleanStatus();
5365  }
5366  return returnCode;
5367}
5368/* Pivot in a variable and out a variable.  Returns 0 if okay,
5369   1 if inaccuracy forced re-factorization, -1 if would be singular.
5370   Also updates primal/dual infeasibilities.
5371   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
5372*/
5373int ClpSimplex::pivot()
5374{
5375  // scaling not allowed
5376  assert (!scalingFlag_);
5377  // assume In_ and Out_ are correct and directionOut_ set
5378  // (or In_ if flip
5379  lowerIn_ = lower_[sequenceIn_];
5380  valueIn_ = solution_[sequenceIn_];
5381  upperIn_ = upper_[sequenceIn_];
5382  dualIn_ = dj_[sequenceIn_];
5383  lowerOut_ = lower_[sequenceOut_];
5384  valueOut_ = solution_[sequenceOut_];
5385  upperOut_ = upper_[sequenceOut_];
5386  // for now assume primal is feasible (or in dual)
5387  dualOut_ = dj_[sequenceOut_];
5388  assert(fabs(dualOut_)<1.0e-6);
5389  bool roundAgain = true;
5390  int returnCode=0;
5391  while (roundAgain) {
5392    roundAgain=false;
5393    unpack(rowArray_[1]);
5394    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
5395    alpha_=0.0;
5396    int i;
5397    int * index = rowArray_[1]->getIndices();
5398    int number = rowArray_[1]->getNumElements();
5399    double * element = rowArray_[1]->denseVector();
5400    for (i=0;i<number;i++) {
5401      int ii = index[i];
5402      if ( pivotVariable_[ii]==sequenceOut_) {
5403        pivotRow_=ii;
5404        alpha_=element[pivotRow_];
5405        break;
5406      }
5407    }
5408    assert (fabs(alpha_)>1.0e-12);
5409    // we are going to subtract movement from current basic
5410    double movement;
5411    // see where incoming will go to
5412    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
5413      // flip so go to bound
5414      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
5415    } else {
5416      // get where outgoing needs to get to
5417      double outValue = (directionOut_<0) ? upperOut_ : lowerOut_;
5418      // solutionOut_ - movement*alpha_ == outValue
5419      movement = (valueOut_-outValue)/alpha_;
5420      // set directionIn_ correctly
5421      directionIn_ = (movement>0) ? 1 :-1;
5422    }
5423    // update primal solution
5424    for (i=0;i<number;i++) {
5425      int ii = index[i];
5426      // get column
5427      int ij = pivotVariable_[ii];
5428      solution_[ij] -= movement*element[ii];
5429      element[ii]=0.0;
5430    }
5431    rowArray_[1]->setNumElements(0);
5432    // see where something went to
5433    if (sequenceOut_<0) {
5434      if (directionIn_<0) {
5435        assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
5436        solution_[sequenceIn_]=upperIn_;
5437      } else {
5438        assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
5439        solution_[sequenceIn_]=lowerIn_;
5440      }
5441    } else {
5442      if (directionOut_<0) {
5443        assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
5444        solution_[sequenceOut_]=upperOut_;
5445      } else {
5446        assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
5447        solution_[sequenceOut_]=lowerOut_;
5448      }
5449      solution_[sequenceIn_]=valueIn_+movement;
5450    }
5451    double objectiveChange = dualIn_*movement;
5452    // update duals
5453    if (pivotRow_>=0) {
5454      assert (fabs(alpha_)>1.0e-8);
5455      double multiplier = dualIn_/alpha_;
5456      rowArray_[0]->insert(pivotRow_,multiplier);
5457      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
5458      // put row of tableau in rowArray[0] and columnArray[0]
5459      matrix_->transposeTimes(this,-1.0,
5460                              rowArray_[0],columnArray_[1],columnArray_[0]);
5461      // update column djs
5462      int i;
5463      int * index = columnArray_[0]->getIndices();
5464      int number = columnArray_[0]->getNumElements();
5465      double * element = columnArray_[0]->denseVector();
5466      for (i=0;i<number;i++) {
5467        int ii = index[i];
5468        dj_[ii] += element[ii];
5469        element[ii]=0.0;
5470      }
5471      columnArray_[0]->setNumElements(0);
5472      // and row djs
5473      index = rowArray_[0]->getIndices();
5474      number = rowArray_[0]->getNumElements();
5475      element = rowArray_[0]->denseVector();
5476      for (i=0;i<number;i++) {
5477        int ii = index[i];
5478        dj_[ii+numberColumns_] += element[ii];
5479        dual_[ii] = dj_[ii+numberColumns_];
5480        element[ii]=0.0;
5481      }
5482      rowArray_[0]->setNumElements(0);
5483      // check incoming
5484      assert (fabs(dj_[sequenceIn_])<1.0e-6);
5485    }
5486   
5487    // if stable replace in basis
5488    int updateStatus = factorization_->replaceColumn(this,
5489                                                     rowArray_[2],
5490                                                     rowArray_[1],
5491                                                     pivotRow_,
5492                                                     alpha_);
5493    bool takePivot=true;
5494    // if no pivots, bad update but reasonable alpha - take and invert
5495    if (updateStatus==2&&
5496        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
5497      updateStatus=4;
5498    if (updateStatus==1||updateStatus==4) {
5499      // slight error
5500      if (factorization_->pivots()>5||updateStatus==4) {
5501        returnCode=-1;
5502      }
5503    } else if (updateStatus==2) {
5504      // major error
5505      rowArray_[1]->clear();
5506      takePivot=false;
5507      if (factorization_->pivots()) {
5508        // refactorize here
5509        statusOfProblem();
5510        roundAgain=true;
5511      } else {
5512        returnCode=1;
5513      }
5514    } else if (updateStatus==3) {
5515      // out of memory
5516      // increase space if not many iterations
5517      if (factorization_->pivots()<
5518          0.5*factorization_->maximumPivots()&&
5519          factorization_->pivots()<200)
5520        factorization_->areaFactor(
5521                                   factorization_->areaFactor() * 1.1);
5522      returnCode =-1; // factorize now
5523    }
5524    if (takePivot) {
5525      int save = algorithm_;
5526      // make simple so always primal
5527      algorithm_=1;
5528      housekeeping(objectiveChange);
5529      algorithm_=save;
5530    }
5531  }
5532  if (returnCode == -1) {
5533    // refactorize here
5534    statusOfProblem();
5535  } else {
5536    // just for now - recompute anyway
5537    gutsOfSolution(NULL,NULL);
5538  }
5539  return returnCode;
5540}
5541
5542/* Pivot in a variable and choose an outgoing one.  Assumes primal
5543   feasible - will not go through a bound.  Returns step length in theta
5544   Returns ray in ray_ (or NULL if no pivot)
5545   Return codes as before but -1 means no acceptable pivot
5546*/
5547int ClpSimplex::primalPivotResult()
5548{
5549  assert (sequenceIn_>=0);
5550  valueIn_=solution_[sequenceIn_];
5551  lowerIn_=lower_[sequenceIn_];
5552  upperIn_=upper_[sequenceIn_];
5553  dualIn_=dj_[sequenceIn_];
5554
5555  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
5556  if (returnCode<0&&returnCode>-4) {
5557    return 0;
5558  } else {
5559    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
5560           returnCode);
5561    return -1;
5562  }
5563}
5564 
5565/* Pivot out a variable and choose an incoing one.  Assumes dual
5566   feasible - will not go through a reduced cost. 
5567   Returns step length in theta
5568   Returns ray in ray_ (or NULL if no pivot)
5569   Return codes as before but -1 means no acceptable pivot
5570*/
5571int 
5572ClpSimplex::dualPivotResult()
5573{
5574  return ((ClpSimplexDual *) this)->pivotResult();
5575}
5576// Factorization frequency
5577int 
5578ClpSimplex::factorizationFrequency() const
5579{
5580  if (factorization_)
5581    return factorization_->maximumPivots();
5582  else 
5583    return -1;
5584}
5585void 
5586ClpSimplex::setFactorizationFrequency(int value)
5587{
5588  if (factorization_)
5589    factorization_->maximumPivots(value);
5590}
5591// Common bits of coding for dual and primal
5592int 
5593ClpSimplex::startup(int ifValuesPass, int startFinishOptions)
5594{
5595  // Get rid of some arrays and empty factorization
5596  int useFactorization=false;
5597  if ((startFinishOptions&2)!=0&&(whatsChanged_&(2+512))==2+512)
5598    useFactorization=true; // Keep factorization if possible
5599  // sanity check
5600  // bad if empty (trap here to avoid using bad matrix_)
5601  if (!matrix_||(!matrix_->getNumElements()&&objective_->type()<2)) {
5602    int infeasNumber[2];
5603    double infeasSum[2];
5604    problemStatus_=emptyProblem(infeasNumber,infeasSum);
5605    numberDualInfeasibilities_=infeasNumber[0];
5606    sumDualInfeasibilities_=infeasSum[0];
5607    numberPrimalInfeasibilities_=infeasNumber[1];
5608    sumPrimalInfeasibilities_=infeasSum[1];
5609    return 2;
5610  }
5611  pivotRow_=-1;
5612  sequenceIn_=-1;
5613  sequenceOut_=-1;
5614  secondaryStatus_=0;
5615
5616  primalTolerance_=dblParam_[ClpPrimalTolerance];
5617  dualTolerance_=dblParam_[ClpDualTolerance];
5618  if (problemStatus_!=10)
5619    numberIterations_=0;
5620
5621  // put in standard form (and make row copy)
5622  // create modifiable copies of model rim and do optional scaling
5623  bool goodMatrix=createRim(7+8+16+32,