source: branches/pre/ClpSimplex.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago

new stuff

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