source: branches/pre/ClpSimplex.cpp @ 182

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

Finished for now

  • 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#include "ClpSimplexPrimalQuadratic.hpp"
2790/* Solves quadratic problem using SLP - may be used as crash
2791   for other algorithms when number of iterations small
2792*/
2793int 
2794ClpSimplex::quadraticSLP(int numberPasses, double deltaTolerance)
2795{
2796  return ((ClpSimplexPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
2797}
2798// Solves quadratic using Dantzig's algorithm - primal
2799int 
2800ClpSimplex::quadraticPrimal(int phase)
2801{
2802  return ((ClpSimplexPrimalQuadratic *) this)->primalQuadratic(phase);
2803}
2804/* For strong branching.  On input lower and upper are new bounds
2805   while on output they are objective function values (>1.0e50 infeasible).
2806   Return code is 0 if nothing interesting, -1 if infeasible both
2807   ways and +1 if infeasible one way (check values to see which one(s))
2808*/
2809int ClpSimplex::strongBranching(int numberVariables,const int * variables,
2810                                double * newLower, double * newUpper,
2811                                bool stopOnFirstInfeasible,
2812                                bool alwaysFinish)
2813{
2814  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
2815                                                    newLower,  newUpper,
2816                                                    stopOnFirstInfeasible,
2817                                                    alwaysFinish);
2818}
2819/* Borrow model.  This is so we dont have to copy large amounts
2820   of data around.  It assumes a derived class wants to overwrite
2821   an empty model with a real one - while it does an algorithm.
2822   This is same as ClpModel one, but sets scaling on etc. */
2823void 
2824ClpSimplex::borrowModel(ClpModel & otherModel) 
2825{
2826  ClpModel::borrowModel(otherModel);
2827  createStatus();
2828  ClpDualRowSteepest steep1;
2829  setDualRowPivotAlgorithm(steep1);
2830  ClpPrimalColumnSteepest steep2;
2831  setPrimalColumnPivotAlgorithm(steep2);
2832}
2833typedef struct {
2834  double optimizationDirection;
2835  double dblParam[ClpLastDblParam];
2836  double objectiveValue;
2837  double dualBound;
2838  double dualTolerance;
2839  double primalTolerance;
2840  double sumDualInfeasibilities;
2841  double sumPrimalInfeasibilities;
2842  double infeasibilityCost;
2843  int numberRows;
2844  int numberColumns;
2845  int intParam[ClpLastIntParam];
2846  int numberIterations;
2847  int problemStatus;
2848  int maximumIterations;
2849  int lengthNames;
2850  int numberDualInfeasibilities;
2851  int numberDualInfeasibilitiesWithoutFree;
2852  int numberPrimalInfeasibilities;
2853  int numberRefinements;
2854  int scalingFlag;
2855  int algorithm;
2856  unsigned int specialOptions;
2857  int dualPivotChoice;
2858  int primalPivotChoice;
2859  int matrixStorageChoice;
2860} Clp_scalars;
2861
2862int outDoubleArray(double * array, int length, FILE * fp)
2863{
2864  int numberWritten;
2865  if (array&&length) {
2866    numberWritten = fwrite(&length,sizeof(int),1,fp);
2867    if (numberWritten!=1)
2868      return 1;
2869    numberWritten = fwrite(array,sizeof(double),length,fp);
2870    if (numberWritten!=length)
2871      return 1;
2872  } else {
2873    length = 0;
2874    numberWritten = fwrite(&length,sizeof(int),1,fp);
2875    if (numberWritten!=1)
2876      return 1;
2877  }
2878  return 0;
2879}
2880// Save model to file, returns 0 if success
2881int
2882ClpSimplex::saveModel(const char * fileName)
2883{
2884  FILE * fp = fopen(fileName,"wb");
2885  if (fp) {
2886    Clp_scalars scalars;
2887    int i;
2888    CoinBigIndex numberWritten;
2889    // Fill in scalars
2890    scalars.optimizationDirection = optimizationDirection_;
2891    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2892    scalars.objectiveValue = objectiveValue_;
2893    scalars.dualBound = dualBound_;
2894    scalars.dualTolerance = dualTolerance_;
2895    scalars.primalTolerance = primalTolerance_;
2896    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2897    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2898    scalars.infeasibilityCost = infeasibilityCost_;
2899    scalars.numberRows = numberRows_;
2900    scalars.numberColumns = numberColumns_;
2901    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2902    scalars.numberIterations = numberIterations_;
2903    scalars.problemStatus = problemStatus_;
2904    scalars.maximumIterations = maximumIterations();
2905    scalars.lengthNames = lengthNames_;
2906    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2907    scalars.numberDualInfeasibilitiesWithoutFree
2908      = numberDualInfeasibilitiesWithoutFree_;
2909    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2910    scalars.numberRefinements = numberRefinements_;
2911    scalars.scalingFlag = scalingFlag_;
2912    scalars.algorithm = algorithm_;
2913    scalars.specialOptions = specialOptions_;
2914    scalars.dualPivotChoice = dualRowPivot_->type();
2915    scalars.primalPivotChoice = primalColumnPivot_->type();
2916    scalars.matrixStorageChoice = matrix_->type();
2917
2918    // put out scalars
2919    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2920    if (numberWritten!=1)
2921      return 1;
2922    // strings
2923    CoinBigIndex length;
2924    for (i=0;i<ClpLastStrParam;i++) {
2925      length = strParam_[i].size();
2926      numberWritten = fwrite(&length,sizeof(int),1,fp);
2927      if (numberWritten!=1)
2928        return 1;
2929      if (length) {
2930        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2931        if (numberWritten!=1)
2932          return 1;
2933      }
2934    }
2935    // arrays - in no particular order
2936    if (outDoubleArray(rowActivity_,numberRows_,fp))
2937        return 1;
2938    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2939        return 1;
2940    if (outDoubleArray(dual_,numberRows_,fp))
2941        return 1;
2942    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2943        return 1;
2944    if (outDoubleArray(rowLower_,numberRows_,fp))
2945        return 1;
2946    if (outDoubleArray(rowUpper_,numberRows_,fp))
2947        return 1;
2948    if (outDoubleArray(objective(),numberColumns_,fp))
2949        return 1;
2950    if (outDoubleArray(rowObjective_,numberRows_,fp))
2951        return 1;
2952    if (outDoubleArray(columnLower_,numberColumns_,fp))
2953        return 1;
2954    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2955        return 1;
2956    if (ray_) {
2957      if (problemStatus_==1)
2958        if (outDoubleArray(ray_,numberRows_,fp))
2959          return 1;
2960      else if (problemStatus_==2)
2961        if (outDoubleArray(ray_,numberColumns_,fp))
2962          return 1;
2963      else
2964        if (outDoubleArray(NULL,0,fp))
2965          return 1;
2966    } else {
2967      if (outDoubleArray(NULL,0,fp))
2968        return 1;
2969    }
2970    if (status_&&(numberRows_+numberColumns_)>0) {
2971      length = numberRows_+numberColumns_;
2972      numberWritten = fwrite(&length,sizeof(int),1,fp);
2973      if (numberWritten!=1)
2974        return 1;
2975      numberWritten = fwrite(status_,sizeof(char),length, fp);
2976      if (numberWritten!=length)
2977        return 1;
2978    } else {
2979      length = 0;
2980      numberWritten = fwrite(&length,sizeof(int),1,fp);
2981      if (numberWritten!=1)
2982        return 1;
2983    }
2984    if (lengthNames_) {
2985      char * array = 
2986        new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
2987      char * put = array;
2988      assert (numberRows_ == (int) rowNames_.size());
2989      for (i=0;i<numberRows_;i++) {
2990        assert((int)rowNames_[i].size()<=lengthNames_);
2991        strcpy(put,rowNames_[i].c_str());
2992        put += lengthNames_+1;
2993      }
2994      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
2995      if (numberWritten!=numberRows_)
2996        return 1;
2997      put=array;
2998      assert (numberColumns_ == (int) columnNames_.size());
2999      for (i=0;i<numberColumns_;i++) {
3000        assert((int) columnNames_[i].size()<=lengthNames_);
3001        strcpy(put,columnNames_[i].c_str());
3002        put += lengthNames_+1;
3003      }
3004      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
3005      if (numberWritten!=numberColumns_)
3006        return 1;
3007      delete [] array;
3008    }
3009    // just standard type at present
3010    assert (matrix_->type()==1);
3011    assert (matrix_->getNumCols() == numberColumns_);
3012    assert (matrix_->getNumRows() == numberRows_);
3013    // we are going to save with gaps
3014    length = matrix_->getVectorStarts()[numberColumns_-1]
3015      + matrix_->getVectorLengths()[numberColumns_-1];
3016    numberWritten = fwrite(&length,sizeof(int),1,fp);
3017    if (numberWritten!=1)
3018      return 1;
3019    numberWritten = fwrite(matrix_->getElements(),
3020                           sizeof(double),length,fp);
3021    if (numberWritten!=length)
3022      return 1;
3023    numberWritten = fwrite(matrix_->getIndices(),
3024                           sizeof(int),length,fp);
3025    if (numberWritten!=length)
3026      return 1;
3027    numberWritten = fwrite(matrix_->getVectorStarts(),
3028                           sizeof(int),numberColumns_+1,fp);
3029    if (numberWritten!=numberColumns_+1)
3030      return 1;
3031    numberWritten = fwrite(matrix_->getVectorLengths(),
3032                           sizeof(int),numberColumns_,fp);
3033    if (numberWritten!=numberColumns_)
3034      return 1;
3035    // finished
3036    fclose(fp);
3037    return 0;
3038  } else {
3039    return -1;
3040  }
3041}
3042
3043int inDoubleArray(double * &array, int length, FILE * fp)
3044{
3045  int numberRead;
3046  int length2;
3047  numberRead = fread(&length2,sizeof(int),1,fp);
3048  if (numberRead!=1)
3049    return 1;
3050  if (length2) {
3051    // lengths must match
3052    if (length!=length2)
3053      return 2;
3054    array = new double[length];
3055    numberRead = fread(array,sizeof(double),length,fp);
3056    if (numberRead!=length)
3057      return 1;
3058  } 
3059  return 0;
3060}
3061/* Restore model from file, returns 0 if success,
3062   deletes current model */
3063int 
3064ClpSimplex::restoreModel(const char * fileName)
3065{
3066  FILE * fp = fopen(fileName,"rb");
3067  if (fp) {
3068    // Get rid of current model
3069    ClpModel::gutsOfDelete();
3070    gutsOfDelete(0);
3071    int i;
3072    for (i=0;i<6;i++) {
3073      rowArray_[i]=NULL;
3074      columnArray_[i]=NULL;
3075    }
3076    // get an empty factorization so we can set tolerances etc
3077    factorization_ = new ClpFactorization();
3078    // Say sparse
3079    factorization_->sparseThreshold(1);
3080    Clp_scalars scalars;
3081    CoinBigIndex numberRead;
3082
3083    // get scalars
3084    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
3085    if (numberRead!=1)
3086      return 1;
3087    // Fill in scalars
3088    optimizationDirection_ = scalars.optimizationDirection;
3089    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
3090    objectiveValue_ = scalars.objectiveValue;
3091    dualBound_ = scalars.dualBound;
3092    dualTolerance_ = scalars.dualTolerance;
3093    primalTolerance_ = scalars.primalTolerance;
3094    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
3095    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
3096    infeasibilityCost_ = scalars.infeasibilityCost;
3097    numberRows_ = scalars.numberRows;
3098    numberColumns_ = scalars.numberColumns;
3099    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
3100    numberIterations_ = scalars.numberIterations;
3101    problemStatus_ = scalars.problemStatus;
3102    setMaximumIterations(scalars.maximumIterations);
3103    lengthNames_ = scalars.lengthNames;
3104    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
3105    numberDualInfeasibilitiesWithoutFree_
3106      = scalars.numberDualInfeasibilitiesWithoutFree;
3107    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
3108    numberRefinements_ = scalars.numberRefinements;
3109    scalingFlag_ = scalars.scalingFlag;
3110    algorithm_ = scalars.algorithm;
3111    specialOptions_ = scalars.specialOptions;
3112    // strings
3113    CoinBigIndex length;
3114    for (i=0;i<ClpLastStrParam;i++) {
3115      numberRead = fread(&length,sizeof(int),1,fp);
3116      if (numberRead!=1)
3117        return 1;
3118      if (length) {
3119        char * array = new char[length+1];
3120        numberRead = fread(array,length,1,fp);
3121        if (numberRead!=1)
3122          return 1;
3123        array[length]='\0';
3124        strParam_[i]=array;
3125        delete [] array;
3126      }
3127    }
3128    // arrays - in no particular order
3129    if (inDoubleArray(rowActivity_,numberRows_,fp))
3130        return 1;
3131    if (inDoubleArray(columnActivity_,numberColumns_,fp))
3132        return 1;
3133    if (inDoubleArray(dual_,numberRows_,fp))
3134        return 1;
3135    if (inDoubleArray(reducedCost_,numberColumns_,fp))
3136        return 1;
3137    if (inDoubleArray(rowLower_,numberRows_,fp))
3138        return 1;
3139    if (inDoubleArray(rowUpper_,numberRows_,fp))
3140        return 1;
3141    double * objective;
3142    if (inDoubleArray(objective,numberColumns_,fp))
3143        return 1;
3144    delete objective_;
3145    objective_ = new ClpLinearObjective(objective,numberColumns_);
3146    delete [] objective;
3147    if (inDoubleArray(rowObjective_,numberRows_,fp))
3148        return 1;
3149    if (inDoubleArray(columnLower_,numberColumns_,fp))
3150        return 1;
3151    if (inDoubleArray(columnUpper_,numberColumns_,fp))
3152        return 1;
3153    if (problemStatus_==1) {
3154      if (inDoubleArray(ray_,numberRows_,fp))
3155        return 1;
3156    } else if (problemStatus_==2) {
3157      if (inDoubleArray(ray_,numberColumns_,fp))
3158        return 1;
3159    } else {
3160      // ray should be null
3161      numberRead = fread(&length,sizeof(int),1,fp);
3162      if (numberRead!=1)
3163        return 1;
3164      if (length)
3165        return 2;
3166    }
3167    delete [] status_;
3168    status_=NULL;
3169    // status region
3170    numberRead = fread(&length,sizeof(int),1,fp);
3171    if (numberRead!=1)
3172        return 1;
3173    if (length) {
3174      if (length!=numberRows_+numberColumns_)
3175        return 1;
3176      status_ = new char unsigned[length];
3177      numberRead = fread(status_,sizeof(char),length, fp);
3178      if (numberRead!=length)
3179        return 1;
3180    }
3181    if (lengthNames_) {
3182      char * array = 
3183        new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
3184      char * get = array;
3185      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
3186      if (numberRead!=numberRows_)
3187        return 1;
3188      rowNames_ = std::vector<std::string> ();
3189      rowNames_.resize(numberRows_);
3190      for (i=0;i<numberRows_;i++) {
3191        rowNames_[i]=get;
3192        get += lengthNames_+1;
3193      }
3194      get = array;
3195      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
3196      if (numberRead!=numberColumns_)
3197        return 1;
3198      columnNames_ = std::vector<std::string> ();
3199      columnNames_.resize(numberColumns_);
3200      for (i=0;i<numberColumns_;i++) {
3201        columnNames_[i]=get;
3202        get += lengthNames_+1;
3203      }
3204      delete [] array;
3205    }
3206    // Pivot choices
3207    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
3208    delete dualRowPivot_;
3209    switch ((scalars.dualPivotChoice&63)) {
3210    default:
3211      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
3212    case 1:
3213      // Dantzig
3214      dualRowPivot_ = new ClpDualRowDantzig();
3215      break;
3216    case 2:
3217      // Steepest - use mode
3218      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
3219      break;
3220    }
3221    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
3222    delete primalColumnPivot_;
3223    switch ((scalars.primalPivotChoice&63)) {
3224    default:
3225      printf("Need another primalPivot case %d\n",
3226             scalars.primalPivotChoice&63);
3227    case 1:
3228      // Dantzig
3229      primalColumnPivot_ = new ClpPrimalColumnDantzig();
3230      break;
3231    case 2:
3232      // Steepest - use mode
3233      primalColumnPivot_
3234        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
3235      break;
3236    }
3237    assert(scalars.matrixStorageChoice==1);
3238    delete matrix_;
3239    // get arrays
3240    numberRead = fread(&length,sizeof(int),1,fp);
3241    if (numberRead!=1)
3242      return 1;
3243    double * elements = new double[length];
3244    int * indices = new int[length];
3245    CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
3246    int * lengths = new int[numberColumns_];
3247    numberRead = fread(elements, sizeof(double),length,fp);
3248    if (numberRead!=length)
3249      return 1;
3250    numberRead = fread(indices, sizeof(int),length,fp);
3251    if (numberRead!=length)
3252      return 1;
3253    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
3254    if (numberRead!=numberColumns_+1)
3255      return 1;
3256    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
3257    if (numberRead!=numberColumns_)
3258      return 1;
3259    // assign matrix
3260    CoinPackedMatrix * matrix = new CoinPackedMatrix();
3261    matrix->assignMatrix(true, numberRows_, numberColumns_,
3262                         length, elements, indices, starts, lengths);
3263    // and transfer to Clp
3264    matrix_ = new ClpPackedMatrix(matrix);
3265    // finished
3266    fclose(fp);
3267    return 0;
3268  } else {
3269    return -1;
3270  }
3271  return 0;
3272}
3273// value of incoming variable (in Dual)
3274double 
3275ClpSimplex::valueIncomingDual() const
3276{
3277  // Need value of incoming for list of infeasibilities as may be infeasible
3278  double valueIncoming = (dualOut_/alpha_)*directionOut_;
3279  if (directionIn_==-1)
3280    valueIncoming = upperIn_-valueIncoming;
3281  else
3282    valueIncoming = lowerIn_-valueIncoming;
3283  return valueIncoming;
3284}
3285// Sanity check on input data - returns true if okay
3286bool 
3287ClpSimplex::sanityCheck()
3288{
3289  // bad if empty
3290  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
3291    handler_->message(CLP_EMPTY_PROBLEM,messages_)
3292      <<numberRows_
3293      <<numberColumns_
3294      <<matrix_->getNumElements()
3295      <<CoinMessageEol;
3296    problemStatus_=4;
3297    return false;
3298  }
3299  int numberBad ;
3300  double largestBound, smallestBound, minimumGap;
3301  double smallestObj, largestObj;
3302  int firstBad;
3303  int modifiedBounds=0;
3304  int i;
3305  numberBad=0;
3306  firstBad=-1;
3307  minimumGap=1.0e100;
3308  smallestBound=1.0e100;
3309  largestBound=0.0;
3310  smallestObj=1.0e100;
3311  largestObj=0.0;
3312  // If bounds are too close - fix
3313  double fixTolerance = 1.1*primalTolerance_;
3314  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
3315    double value;
3316    value = fabs(cost_[i]);
3317    if (value>1.0e50) {
3318      numberBad++;
3319      if (firstBad<0)
3320        firstBad=i;
3321    } else if (value) {
3322      if (value>largestObj)
3323        largestObj=value;
3324      if (value<smallestObj)
3325        smallestObj=value;
3326    }
3327    value=upper_[i]-lower_[i];
3328    if (value<-primalTolerance_) {
3329      numberBad++;
3330      if (firstBad<0)
3331        firstBad=i;
3332    } else if (value<=fixTolerance) {
3333      if (value) {
3334        // modify
3335        upper_[i] = lower_[i];
3336        modifiedBounds++;
3337      }
3338    } else {
3339      if (value<minimumGap)
3340        minimumGap=value;
3341    }
3342    if (lower_[i]>-1.0e100&&lower_[i]) {
3343      value = fabs(lower_[i]);
3344      if (value>largestBound)
3345        largestBound=value;
3346      if (value<smallestBound)
3347        smallestBound=value;
3348    }
3349    if (upper_[i]<1.0e100&&upper_[i]) {
3350      value = fabs(upper_[i]);
3351      if (value>largestBound)
3352        largestBound=value;
3353      if (value<smallestBound)
3354        smallestBound=value;
3355    }
3356  }
3357  if (largestBound)
3358    handler_->message(CLP_RIMSTATISTICS3,messages_)
3359      <<smallestBound
3360      <<largestBound
3361      <<minimumGap
3362      <<CoinMessageEol;
3363  minimumGap=1.0e100;
3364  smallestBound=1.0e100;
3365  largestBound=0.0;
3366  for (i=0;i<numberColumns_;i++) {
3367    double value;
3368    value = fabs(cost_[i]);
3369    if (value>1.0e50) {
3370      numberBad++;
3371      if (firstBad<0)
3372        firstBad=i;
3373    } else if (value) {
3374      if (value>largestObj)
3375        largestObj=value;
3376      if (value<smallestObj)
3377        smallestObj=value;
3378    }
3379    value=upper_[i]-lower_[i];
3380    if (value<-primalTolerance_) {
3381      numberBad++;
3382      if (firstBad<0)
3383        firstBad=i;
3384    } else if (value<=fixTolerance) {
3385      if (value) {
3386        // modify
3387        upper_[i] = lower_[i];
3388        modifiedBounds++;
3389      }
3390    } else {
3391      if (value<minimumGap)
3392        minimumGap=value;
3393    }
3394    if (lower_[i]>-1.0e100&&lower_[i]) {
3395      value = fabs(lower_[i]);
3396      if (value>largestBound)
3397        largestBound=value;
3398      if (value<smallestBound)
3399        smallestBound=value;
3400    }
3401    if (upper_[i]<1.0e100&&upper_[i]) {
3402      value = fabs(upper_[i]);
3403      if (value>largestBound)
3404        largestBound=value;
3405      if (value<smallestBound)
3406        smallestBound=value;
3407    }
3408  }
3409  char rowcol[]={'R','C'};
3410  if (numberBad) {
3411    handler_->message(CLP_BAD_BOUNDS,messages_)
3412      <<numberBad
3413      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
3414      <<CoinMessageEol;
3415    problemStatus_=4;
3416    return false;
3417  }
3418  if (modifiedBounds)
3419    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
3420      <<modifiedBounds
3421      <<CoinMessageEol;
3422  handler_->message(CLP_RIMSTATISTICS1,messages_)
3423    <<smallestObj
3424    <<largestObj
3425    <<CoinMessageEol;
3426  if (largestBound)
3427    handler_->message(CLP_RIMSTATISTICS2,messages_)
3428      <<smallestBound
3429      <<largestBound
3430      <<minimumGap
3431      <<CoinMessageEol;
3432  return true;
3433}
3434// Set up status array (for OsiClp)
3435void 
3436ClpSimplex::createStatus() 
3437{
3438  if(!status_)
3439    status_ = new unsigned char [numberColumns_+numberRows_];
3440  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3441  int i;
3442  // set column status to one nearest zero
3443  for (i=0;i<numberColumns_;i++) {
3444#if 0
3445    if (columnLower_[i]>=0.0) {
3446      setColumnStatus(i,atLowerBound);
3447    } else if (columnUpper_[i]<=0.0) {
3448      setColumnStatus(i,atUpperBound);
3449    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
3450      // free
3451      setColumnStatus(i,isFree);
3452    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
3453      setColumnStatus(i,atLowerBound);
3454    } else {
3455      setColumnStatus(i,atUpperBound);
3456    }
3457#else
3458    setColumnStatus(i,atLowerBound);
3459#endif
3460  }
3461  for (i=0;i<numberRows_;i++) {
3462    setRowStatus(i,basic);
3463  }
3464}
3465/* Loads a problem (the constraints on the
3466   rows are given by lower and upper bounds). If a pointer is 0 then the
3467   following values are the default:
3468   <ul>
3469   <li> <code>colub</code>: all columns have upper bound infinity
3470   <li> <code>collb</code>: all columns have lower bound 0
3471   <li> <code>rowub</code>: all rows have upper bound infinity
3472   <li> <code>rowlb</code>: all rows have lower bound -infinity
3473   <li> <code>obj</code>: all variables have 0 objective coefficient
3474   </ul>
3475*/
3476void 
3477ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
3478                    const double* collb, const double* colub,   
3479                    const double* obj,
3480                    const double* rowlb, const double* rowub,
3481                    const double * rowObjective)
3482{
3483  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3484                        rowObjective);
3485  createStatus();
3486}
3487void 
3488ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
3489                    const double* collb, const double* colub,   
3490                    const double* obj,
3491                    const double* rowlb, const double* rowub,
3492                    const double * rowObjective)
3493{
3494  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3495                        rowObjective);
3496  createStatus();
3497}
3498
3499/* Just like the other loadProblem() method except that the matrix is
3500   given in a standard column major ordered format (without gaps). */
3501void 
3502ClpSimplex::loadProblem (  const int numcols, const int numrows,
3503                    const CoinBigIndex* start, const int* index,
3504                    const double* value,
3505                    const double* collb, const double* colub,   
3506                    const double* obj,
3507                    const double* rowlb, const double* rowub,
3508                    const double * rowObjective)
3509{
3510  ClpModel::loadProblem(numcols, numrows, start, index, value,
3511                          collb, colub, obj, rowlb, rowub,
3512                          rowObjective);
3513  createStatus();
3514}
3515void 
3516ClpSimplex::loadProblem (  const int numcols, const int numrows,
3517                           const CoinBigIndex* start, const int* index,
3518                           const double* value,const int * length,
3519                           const double* collb, const double* colub,   
3520                           const double* obj,
3521                           const double* rowlb, const double* rowub,
3522                           const double * rowObjective)
3523{
3524  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3525                          collb, colub, obj, rowlb, rowub,
3526                          rowObjective);
3527  createStatus();
3528}
3529// Read an mps file from the given filename
3530int 
3531ClpSimplex::readMps(const char *filename,
3532            bool keepNames,
3533            bool ignoreErrors)
3534{
3535  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3536  createStatus();
3537  return status;
3538}
3539// Just check solution (for external use)
3540void 
3541ClpSimplex::checkSolution()
3542{
3543  // put in standard form
3544  createRim(7+8+16);
3545  dualTolerance_=dblParam_[ClpDualTolerance];
3546  primalTolerance_=dblParam_[ClpPrimalTolerance];
3547  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
3548  checkDualSolution();
3549  if (!numberDualInfeasibilities_&&
3550      !numberPrimalInfeasibilities_)
3551    problemStatus_=0;
3552  else
3553    problemStatus_=-1;
3554#ifdef CLP_DEBUG
3555  int i;
3556  double value=0.0;
3557  for (i=0;i<numberRows_+numberColumns_;i++)
3558    value += dj_[i]*solution_[i];
3559  printf("dual value %g, primal %g\n",value,objectiveValue());
3560#endif
3561  // release extra memory
3562  deleteRim(0);
3563}
3564/* Crash - at present just aimed at dual, returns
3565   -2 if dual preferred and crash basis created
3566   -1 if dual preferred and all slack basis preferred
3567   0 if basis going in was not all slack
3568   1 if primal preferred and all slack basis preferred
3569   2 if primal preferred and crash basis created.
3570   
3571   if gap between bounds <="gap" variables can be flipped
3572   
3573   If "pivot" is
3574   0 No pivoting (so will just be choice of algorithm)
3575   1 Simple pivoting e.g. gub
3576   2 Mini iterations
3577*/
3578int 
3579ClpSimplex::crash(double gap,int pivot)
3580{
3581  assert(!rowObjective_); // not coded
3582  int iColumn;
3583  int numberBad=0;
3584  int numberBasic=0;
3585  double dualTolerance=dblParam_[ClpDualTolerance];
3586  //double primalTolerance=dblParam_[ClpPrimalTolerance];
3587
3588  // If no basis then make all slack one
3589  if (!status_)
3590    createStatus();
3591 
3592  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3593    if (getColumnStatus(iColumn)==basic)
3594      numberBasic++;
3595  }
3596  if (numberBasic) {
3597    // not all slack
3598    return 0;
3599  } else {
3600    double * dj = new double [numberColumns_];
3601    double * solution = columnActivity_;
3602    const double * linearObjective = objective();
3603    //double objectiveValue=0.0;
3604    int iColumn;
3605    for (iColumn=0;iColumn<numberColumns_;iColumn++)
3606      dj[iColumn] = optimizationDirection_*linearObjective[iColumn];
3607    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3608      // assume natural place is closest to zero
3609      double lowerBound = columnLower_[iColumn];
3610      double upperBound = columnUpper_[iColumn];
3611      if (lowerBound>-1.0e20||upperBound<1.0e20) {
3612        bool atLower;
3613        if (fabs(upperBound)<fabs(lowerBound)) {
3614          atLower=false;
3615          setColumnStatus(iColumn,atUpperBound);
3616          solution[iColumn]=upperBound;
3617        } else {
3618          atLower=true;
3619          setColumnStatus(iColumn,atLowerBound);
3620          solution[iColumn]=lowerBound;
3621        }
3622        if (dj[iColumn]<0.0) {
3623          // should be at upper bound
3624          if (atLower) {
3625            // can we flip
3626            if (upperBound-lowerBound<=gap) {
3627              columnActivity_[iColumn]=upperBound;
3628              setColumnStatus(iColumn,atUpperBound);
3629            } else if (dj[iColumn]<-dualTolerance) {
3630              numberBad++;
3631            }
3632          }
3633        } else if (dj[iColumn]>0.0) {
3634          // should be at lower bound
3635          if (!atLower) {
3636            // can we flip
3637            if (upperBound-lowerBound<=gap) {
3638              columnActivity_[iColumn]=lowerBound;
3639              setColumnStatus(iColumn,atLowerBound);
3640            } else if (dj[iColumn]>dualTolerance) {
3641              numberBad++;
3642            }
3643          }
3644        }
3645      } else {
3646        // free
3647        setColumnStatus(iColumn,isFree);
3648        if (fabs(dj[iColumn])>dualTolerance) 
3649          numberBad++;
3650      }
3651    }
3652    if (numberBad||pivot) {
3653      if (!pivot) {
3654        delete [] dj;
3655        return 1;
3656      } else {
3657        // see if can be made dual feasible with gubs etc
3658        double * pi = new double[numberRows_];
3659        memset (pi,0,numberRows_*sizeof(double));
3660        int * way = new int[numberColumns_];
3661        int numberIn = 0;
3662
3663        // Get column copy
3664        CoinPackedMatrix * columnCopy = matrix();
3665        // Get a row copy in standard format
3666        CoinPackedMatrix copy;
3667        copy.reverseOrderedCopyOf(*columnCopy);
3668        // get matrix data pointers
3669        const int * column = copy.getIndices();
3670        const CoinBigIndex * rowStart = copy.getVectorStarts();
3671        const int * rowLength = copy.getVectorLengths(); 
3672        const double * elementByRow = copy.getElements();
3673        //const int * row = columnCopy->getIndices();
3674        //const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3675        //const int * columnLength = columnCopy->getVectorLengths();
3676        //const double * element = columnCopy->getElements();
3677
3678
3679        // if equality row and bounds mean artificial in basis bad
3680        // then do anyway
3681
3682        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3683          // - if we want to reduce dj, + if we want to increase
3684          int thisWay = 100;
3685          double lowerBound = columnLower_[iColumn];
3686          double upperBound = columnUpper_[iColumn];
3687          if (upperBound>lowerBound) {
3688            switch(getColumnStatus(iColumn)) {
3689             
3690            case basic:
3691              thisWay=0;
3692            case ClpSimplex::isFixed:
3693              break;
3694            case isFree:
3695            case superBasic:
3696              if (dj[iColumn]<-dualTolerance) 
3697                thisWay = 1;
3698              else if (dj[iColumn]>dualTolerance) 
3699                thisWay = -1;
3700              else
3701                thisWay =0;
3702              break;
3703            case atUpperBound:
3704              if (dj[iColumn]>dualTolerance) 
3705                thisWay = -1;
3706              else if (dj[iColumn]<-dualTolerance) 
3707                thisWay = -3;
3708              else
3709                thisWay = -2;
3710              break;
3711            case atLowerBound:
3712              if (dj[iColumn]<-dualTolerance) 
3713                thisWay = 1;
3714              else if (dj[iColumn]>dualTolerance) 
3715                thisWay = 3;
3716              else
3717                thisWay = 2;
3718              break;
3719            }
3720          }
3721          way[iColumn] = thisWay;
3722        }
3723        /*if (!numberBad)
3724          printf("Was dual feasible before passes - rows %d\n",
3725          numberRows_);*/
3726        int lastNumberIn = -100000;
3727        int numberPasses=5;
3728        while (numberIn>lastNumberIn+numberRows_/100) {
3729          lastNumberIn = numberIn;
3730          // we need to maximize chance of doing good
3731          int iRow;
3732          for (iRow=0;iRow<numberRows_;iRow++) {
3733            double lowerBound = rowLower_[iRow];
3734            double upperBound = rowUpper_[iRow];
3735            if (getRowStatus(iRow)==basic) {
3736              // see if we can find a column to pivot on
3737              int j;
3738              // down is amount pi can go down
3739              double maximumDown = COIN_DBL_MAX;
3740              double maximumUp = COIN_DBL_MAX;
3741              double minimumDown =0.0;
3742              double minimumUp =0.0;
3743              int iUp=-1;
3744              int iDown=-1;
3745              int iUpB=-1;
3746              int iDownB=-1;
3747              if (lowerBound<-1.0e20)
3748                maximumUp = -1.0;
3749              if (upperBound>1.0e20)
3750                maximumDown = -1.0;
3751              for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3752                int iColumn = column[j];
3753                double value = elementByRow[j];
3754                double djValue = dj[iColumn];
3755                /* way -
3756                   -3 - okay at upper bound with negative dj
3757                   -2 - marginal at upper bound with zero dj - can only decrease
3758                   -1 - bad at upper bound
3759                   0 - we can never pivot on this row
3760                   1 - bad at lower bound
3761                   2 - marginal at lower bound with zero dj - can only increase
3762                   3 - okay at lower bound with positive dj
3763                   100 - fine we can just ignore
3764                */
3765                if (way[iColumn]!=100) {
3766                  switch(way[iColumn]) {
3767                   
3768                  case -3:
3769                    if (value>0.0) {
3770                      if (maximumDown*value>-djValue) {
3771                        maximumDown = -djValue/value;
3772                        iDown = iColumn;
3773                      }
3774                    } else {
3775                      if (-maximumUp*value>-djValue) {
3776                        maximumUp = djValue/value;
3777                        iUp = iColumn;
3778                      }
3779                    }
3780                    break;
3781                  case -2:
3782                    if (value>0.0) {
3783                      maximumDown = 0.0;
3784                    } else {
3785                      maximumUp = 0.0;
3786                    }
3787                    break;
3788                  case -1:
3789                    // see if could be satisfied
3790                    // dj value > 0
3791                    if (value>0.0) {
3792                      maximumDown=0.0;
3793                      if (maximumUp*value<djValue-dualTolerance) {
3794                        maximumUp = 0.0; // would improve but not enough
3795                      } else {
3796                        if (minimumUp*value<djValue) {
3797                          minimumUp = djValue/value;
3798                          iUpB = iColumn;
3799                        }
3800                      }
3801                    } else {
3802                      maximumUp=0.0;
3803                      if (-maximumDown*value<djValue-dualTolerance) {
3804                        maximumDown = 0.0; // would improve but not enough
3805                      } else {
3806                        if (-minimumDown*value<djValue) {
3807                          minimumDown = -djValue/value;
3808                          iDownB = iColumn;
3809                        }
3810                      }
3811                    }
3812                   
3813                    break;
3814                  case 0:
3815                    maximumDown = -1.0;
3816                    maximumUp=-1.0;
3817                    break;
3818                  case 1:
3819                    // see if could be satisfied
3820                    // dj value < 0
3821                    if (value>0.0) {
3822                      maximumUp=0.0;
3823                      if (maximumDown*value<-djValue-dualTolerance) {
3824                        maximumDown = 0.0; // would improve but not enough
3825                      } else {
3826                        if (minimumDown*value<-djValue) {
3827                          minimumDown = -djValue/value;
3828                          iDownB = iColumn;
3829                        }
3830                      }
3831                    } else {
3832                      maximumDown=0.0;
3833                      if (-maximumUp*value<-djValue-dualTolerance) {
3834                        maximumUp = 0.0; // would improve but not enough
3835                      } else {
3836                        if (-minimumUp*value<-djValue) {
3837                          minimumUp = djValue/value;
3838                          iUpB = iColumn;
3839                        }
3840                      }
3841                    }
3842                   
3843                    break;
3844                  case 2:
3845                    if (value>0.0) {
3846                      maximumUp = 0.0;
3847                    } else {
3848                      maximumDown = 0.0;
3849                    }
3850                   
3851                    break;
3852                  case 3:
3853                    if (value>0.0) {
3854                      if (maximumUp*value>djValue) {
3855                        maximumUp = djValue/value;
3856                        iUp = iColumn;
3857                      }
3858                    } else {
3859                      if (-maximumDown*value>djValue) {
3860                        maximumDown = -djValue/value;
3861                        iDown = iColumn;
3862                      }
3863                    }
3864                   
3865                    break;
3866                  default:
3867                    break;
3868                  }
3869                }
3870              }
3871              if (iUpB>=0)
3872                iUp=iUpB;
3873              if (maximumUp<=dualTolerance||maximumUp<minimumUp)
3874                iUp=-1;
3875              if (iDownB>=0)
3876                iDown=iDownB;
3877              if (maximumDown<=dualTolerance||maximumDown<minimumDown)
3878                iDown=-1;
3879              if (iUp>=0||iDown>=0) {
3880                // do something
3881                if (iUp>=0&&iDown>=0) {
3882                  if (maximumDown>maximumUp)
3883                    iUp=-1;
3884                }
3885                double change;
3886                int kColumn;
3887                if (iUp>=0) {
3888                  kColumn=iUp;
3889                  change=maximumUp;
3890                  // just do minimum if was dual infeasible
3891                  // ? only if maximum large?
3892                  if (minimumUp>0.0)
3893                    change=minimumUp;
3894                  setRowStatus(iRow,atUpperBound);
3895                } else {
3896                  kColumn=iDown;
3897                  change=-maximumDown;
3898                  // just do minimum if was dual infeasible
3899                  // ? only if maximum large?
3900                  if (minimumDown>0.0)
3901                    change=-minimumDown;
3902                  setRowStatus(iRow,atLowerBound);
3903                }
3904                assert (fabs(change)<1.0e20);
3905                setColumnStatus(kColumn,basic);
3906                numberIn++;
3907                pi[iRow]=change;
3908                for (j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
3909                  int iColumn = column[j];
3910                  double value = elementByRow[j];
3911                  double djValue = dj[iColumn]-change*value;
3912                  dj[iColumn]=djValue;
3913                  if (abs(way[iColumn])==1) {
3914                    numberBad--;
3915                    /*if (!numberBad)
3916                      printf("Became dual feasible at row %d out of %d\n",
3917                      iRow, numberRows_);*/
3918                    lastNumberIn=-1000000;
3919                  }
3920                  int thisWay = 100;
3921                  double lowerBound = columnLower_[iColumn];
3922                  double upperBound = columnUpper_[iColumn];
3923                  if (upperBound>lowerBound) {
3924                    switch(getColumnStatus(iColumn)) {
3925                     
3926                    case basic:
3927                      thisWay=0;
3928                    case isFixed:
3929                      break;
3930                    case isFree:
3931                    case superBasic:
3932                      if (djValue<-dualTolerance) 
3933                        thisWay = 1;
3934                      else if (djValue>dualTolerance) 
3935                        thisWay = -1;
3936                      else
3937                        { thisWay =0; abort();}
3938                      break;
3939                    case atUpperBound:
3940                      if (djValue>dualTolerance) 
3941                        { thisWay =-1; abort();}
3942                      else if (djValue<-dualTolerance) 
3943                        thisWay = -3;
3944                      else
3945                        thisWay = -2;
3946                      break;
3947                    case atLowerBound:
3948                      if (djValue<-dualTolerance) 
3949                        { thisWay =1; abort();}
3950                      else if (djValue>dualTolerance) 
3951                        thisWay = 3;
3952                      else
3953                        thisWay = 2;
3954                      break;
3955                    }
3956                  }
3957                  way[iColumn] = thisWay;
3958                }
3959              }
3960            }
3961          }
3962          if (numberIn==lastNumberIn||numberBad||pivot<2)
3963            break;
3964          if (!(--numberPasses))
3965            break;
3966          //printf("%d put in so far\n",numberIn);
3967        }
3968        // last attempt to flip
3969        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3970          double lowerBound = columnLower_[iColumn];
3971          double upperBound = columnUpper_[iColumn];
3972          if (upperBound-lowerBound<=gap&&upperBound>lowerBound) {
3973            double djValue=dj[iColumn];
3974            switch(getColumnStatus(iColumn)) {
3975             
3976            case basic:
3977            case ClpSimplex::isFixed:
3978              break;
3979            case isFree:
3980            case superBasic:
3981              break;
3982            case atUpperBound:
3983              if (djValue>dualTolerance) {
3984                setColumnStatus(iColumn,atUpperBound);
3985                solution[iColumn]=upperBound;
3986              } 
3987              break;
3988            case atLowerBound:
3989              if (djValue<-dualTolerance) {
3990                setColumnStatus(iColumn,atUpperBound);
3991                solution[iColumn]=upperBound;
3992              }
3993              break;
3994            }
3995          }
3996        }
3997        delete [] pi;
3998        delete [] dj;
3999        delete [] way;
4000        handler_->message(CLP_CRASH,messages_)
4001          <<numberIn
4002          <<numberBad
4003          <<CoinMessageEol;
4004        return -1;
4005      }
4006    } else {
4007      delete [] dj;
4008      return -1;
4009    }
4010  }
4011}
4012int
4013ClpSimplex::nextSuperBasic()
4014{
4015  if (firstFree_>=0) {
4016    int returnValue=firstFree_;
4017    int iColumn=firstFree_+1;
4018    if (algorithm_>0) {
4019      // primal
4020      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
4021        if (getStatus(iColumn)==superBasic) {
4022          if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
4023            solution_[iColumn]=lower_[iColumn];
4024            setStatus(iColumn,atLowerBound);
4025          } else if (fabs(solution_[iColumn]-upper_[iColumn])
4026                     <=primalTolerance_) {
4027            solution_[iColumn]=upper_[iColumn];
4028            setStatus(iColumn,atUpperBound);
4029          } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
4030            setStatus(iColumn,isFree);
4031            break;
4032          } else {
4033            break;
4034          }
4035        }
4036      }
4037    } else {
4038      // dual
4039      for (;iColumn<numberRows_+numberColumns_;iColumn++) {
4040        if (getStatus(iColumn)==isFree) 
4041          if (fabs(dj_[iColumn])>1.0e2*dualTolerance_) 
4042            break;
4043      }
4044    }
4045    firstFree_ = iColumn;
4046    if (firstFree_==numberRows_+numberColumns_)
4047      firstFree_=-1;
4048    return returnValue;
4049  } else {
4050    return -1;
4051  }
4052}
4053/* Pivot in a variable and out a variable.  Returns 0 if okay,
4054   1 if inaccuracy forced re-factorization, -1 if would be singular.
4055   Also updates primal/dual infeasibilities.
4056   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
4057*/
4058int ClpSimplex::pivot()
4059{
4060  // scaling not allowed
4061  assert (!scalingFlag_);
4062  // assume In_ and Out_ are correct and directionOut_ set
4063  // (or In_ if flip
4064  lowerIn_ = lower_[sequenceIn_];
4065  valueIn_ = solution_[sequenceIn_];
4066  upperIn_ = upper_[sequenceIn_];
4067  dualIn_ = dj_[sequenceIn_];
4068  if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) {
4069    assert (pivotRow_>=0&&pivotRow_<numberRows_);
4070    assert (pivotVariable_[pivotRow_]==sequenceOut_);
4071    lowerOut_ = lower_[sequenceOut_];
4072    valueOut_ = solution_[sequenceOut_];
4073    upperOut_ = upper_[sequenceOut_];
4074    // for now assume primal is feasible (or in dual)
4075    dualOut_ = dj_[sequenceOut_];
4076    assert(fabs(dualOut_)<1.0e-6);
4077  } else {
4078    assert (pivotRow_<0);
4079  }
4080  bool roundAgain = true;
4081  int returnCode=0;
4082  while (roundAgain) {
4083    roundAgain=false;
4084    unpack(rowArray_[1]);
4085    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
4086    // we are going to subtract movement from current basic
4087    double movement;
4088    // see where incoming will go to
4089    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
4090      // flip so go to bound
4091      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
4092    } else {
4093      // get where outgoing needs to get to
4094      double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
4095      // solutionOut_ - movement*alpha_ == outValue
4096      movement = (outValue-valueOut_)/alpha_;
4097      // set directionIn_ correctly
4098      directionIn_ = (movement>0) ? 1 :-1;
4099    }
4100    // update primal solution
4101    {
4102      int i;
4103      int * index = rowArray_[1]->getIndices();
4104      int number = rowArray_[1]->getNumElements();
4105      double * element = rowArray_[1]->denseVector();
4106      for (i=0;i<number;i++) {
4107        int ii = index[i];
4108        // get column
4109        ii = pivotVariable_[ii];
4110        solution_[ii] -= movement*element[i];
4111        element[i]=0.0;
4112      }
4113      // see where something went to
4114      if (sequenceOut_<0) {
4115        if (directionIn_<0) {
4116          assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
4117          solution_[sequenceIn_]=upperIn_;
4118        } else {
4119          assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
4120          solution_[sequenceIn_]=lowerIn_;
4121        }
4122      } else {
4123        if (directionOut_<0) {
4124          assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
4125          solution_[sequenceOut_]=upperOut_;
4126        } else {
4127          assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
4128          solution_[sequenceOut_]=lowerOut_;
4129        }
4130        solution_[sequenceIn_]=valueIn_+movement;
4131      }
4132    }   
4133    double objectiveChange = dualIn_*movement;
4134    // update duals
4135    if (pivotRow_>=0) {
4136      alpha_ = rowArray_[1]->denseVector()[pivotRow_];
4137      assert (fabs(alpha_)>1.0e-8);
4138      double multiplier = dualIn_/alpha_;
4139      rowArray_[0]->insert(pivotRow_,multiplier);
4140      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
4141      // put row of tableau in rowArray[0] and columnArray[0]
4142      matrix_->transposeTimes(this,-1.0,
4143                              rowArray_[0],columnArray_[1],columnArray_[0]);
4144      // update column djs
4145      int i;
4146      int * index = columnArray_[0]->getIndices();
4147      int number = columnArray_[0]->getNumElements();
4148      double * element = columnArray_[0]->denseVector();
4149      for (i=0;i<number;i++) {
4150        int ii = index[i];
4151        dj_[ii] += element[ii];
4152        element[ii]=0.0;
4153      }
4154      columnArray_[0]->setNumElements(0);
4155      // and row djs
4156      index = rowArray_[0]->getIndices();
4157      number = rowArray_[0]->getNumElements();
4158      element = rowArray_[0]->denseVector();
4159      for (i=0;i<number;i++) {
4160        int ii = index[i];
4161        dj_[ii+numberColumns_] += element[ii];
4162        dual_[ii] = dj_[ii+numberColumns_];
4163        element[ii]=0.0;
4164      }
4165      rowArray_[0]->setNumElements(0);
4166      // check incoming
4167      assert (fabs(dj_[sequenceIn_])<1.0e-6);
4168    }
4169   
4170    // if stable replace in basis
4171    int updateStatus = factorization_->replaceColumn(rowArray_[2],
4172                                                   pivotRow_,
4173                                                     alpha_);
4174    bool takePivot=true;
4175    // if no pivots, bad update but reasonable alpha - take and invert
4176    if (updateStatus==2&&
4177        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
4178      updateStatus=4;
4179    if (updateStatus==1||updateStatus==4) {
4180      // slight error
4181      if (factorization_->pivots()>5||updateStatus==4) {
4182        returnCode=-1;
4183      }
4184    } else if (updateStatus==2) {
4185      // major error
4186      rowArray_[1]->clear();
4187      takePivot=false;
4188      if (factorization_->pivots()) {
4189        // refactorize here
4190        statusOfProblem();
4191        roundAgain=true;
4192      } else {
4193        returnCode=1;
4194      }
4195    } else if (updateStatus==3) {
4196      // out of memory
4197      // increase space if not many iterations
4198      if (factorization_->pivots()<
4199          0.5*factorization_->maximumPivots()&&
4200          factorization_->pivots()<200)
4201        factorization_->areaFactor(
4202                                   factorization_->areaFactor() * 1.1);
4203      returnCode =-1; // factorize now
4204    }
4205    if (takePivot) {
4206      int save = algorithm_;
4207      // make simple so always primal
4208      algorithm_=1;
4209      housekeeping(objectiveChange);
4210      algorithm_=save;
4211    }
4212  }
4213  if (returnCode == -1) {
4214    // refactorize here
4215    statusOfProblem();
4216  } else {
4217    // just for now - recompute anyway
4218    gutsOfSolution(NULL,NULL);
4219  }
4220  return returnCode;
4221}
4222
4223/* Pivot in a variable and choose an outgoing one.  Assumes primal
4224   feasible - will not go through a bound.  Returns step length in theta
4225   Returns ray in ray_ (or NULL if no pivot)
4226   Return codes as before but -1 means no acceptable pivot
4227*/
4228int ClpSimplex::primalPivotResult()
4229{
4230  assert (sequenceIn_>=0);
4231  valueIn_=solution_[sequenceIn_];
4232  lowerIn_=lower_[sequenceIn_];
4233  upperIn_=upper_[sequenceIn_];
4234  dualIn_=dj_[sequenceIn_];
4235
4236  int returnCode = ((ClpSimplexPrimal *) this)->pivotResult();
4237  if (returnCode<0&&returnCode>-4) {
4238    return 0;
4239  } else {
4240    printf("Return code of %d from ClpSimplexPrimal::pivotResult\n",
4241           returnCode);
4242    return -1;
4243  }
4244}
4245 
4246/* Pivot out a variable and choose an incoing one.  Assumes dual
4247   feasible - will not go through a reduced cost. 
4248   Returns step length in theta
4249   Returns ray in ray_ (or NULL if no pivot)
4250   Return codes as before but -1 means no acceptable pivot
4251*/
4252int 
4253ClpSimplex::dualPivotResult()
4254{
4255  return ((ClpSimplexDual *) this)->pivotResult();
4256}
4257// Factorization frequency
4258int 
4259ClpSimplex::factorizationFrequency() const
4260{
4261  if (factorization_)
4262    return factorization_->maximumPivots();
4263  else 
4264    return -1;
4265}
4266void 
4267ClpSimplex::setFactorizationFrequency(int value)
4268{
4269  if (factorization_)
4270    factorization_->maximumPivots(value);
4271}
4272// Common bits of coding for dual and primal
4273int 
4274ClpSimplex::startup(int ifValuesPass)
4275{
4276  // sanity check
4277  assert (numberRows_==matrix_->getNumRows());
4278  assert (numberColumns_==matrix_->getNumCols());
4279  // for moment all arrays must exist
4280  assert(columnLower_);
4281  assert(columnUpper_);
4282  assert(rowLower_);
4283  assert(rowUpper_);
4284  pivotRow_=-1;
4285  sequenceIn_=-1;
4286  sequenceOut_=-1;
4287
4288  primalTolerance_=dblParam_[ClpPrimalTolerance];
4289  dualTolerance_=dblParam_[ClpDualTolerance];
4290  if (problemStatus_!=10)
4291    numberIterations_=0;
4292
4293  // put in standard form (and make row copy)
4294  // create modifiable copies of model rim and do optional scaling
4295  bool goodMatrix=createRim(7+8+16,true);
4296
4297  if (goodMatrix) {
4298    // Model looks okay
4299    // Do initial factorization
4300    // and set certain stuff
4301    // We can either set increasing rows so ...IsBasic gives pivot row
4302    // or we can just increment iBasic one by one
4303    // for now let ...iBasic give pivot row
4304    factorization_->increasingRows(2);
4305    // row activities have negative sign
4306    factorization_->slackValue(-1.0);
4307    factorization_->zeroTolerance(1.0e-13);
4308    // Switch off dense
4309    int saveThreshold = factorization_->denseThreshold();
4310    factorization_->setDenseThreshold(0);
4311
4312    // do perturbation if asked for
4313
4314    if (perturbation_<100) {
4315      if (algorithm_>0) {
4316        ((ClpSimplexPrimal *) this)->perturb();
4317      } else if (algorithm_<0) {
4318        ((ClpSimplexDual *) this)->perturb();
4319      }
4320    }
4321    // for primal we will change bounds using infeasibilityCost_
4322    if (nonLinearCost_==NULL&&algorithm_>0) {
4323      // get a valid nonlinear cost function
4324      delete nonLinearCost_;
4325      nonLinearCost_= new ClpNonLinearCost(this);
4326    }
4327   
4328    // loop round to clean up solution if values pass
4329    int numberThrownOut = -1;
4330    int totalNumberThrownOut=0;
4331    while(numberThrownOut) {
4332      int status = internalFactorize(0+10*ifValuesPass);
4333      if (status<0)
4334        return 1; // some error
4335      else
4336        numberThrownOut = status;
4337     
4338      // for this we need clean basis so it is after factorize
4339      if (!numberThrownOut)
4340        numberThrownOut = gutsOfSolution(  NULL,NULL,
4341                                         ifValuesPass!=0);
4342      totalNumberThrownOut+= numberThrownOut;
4343     
4344    }
4345   
4346    if (totalNumberThrownOut)
4347      handler_->message(CLP_SINGULARITIES,messages_)
4348        <<totalNumberThrownOut
4349        <<CoinMessageEol;
4350    // Switch back dense
4351    factorization_->setDenseThreshold(saveThreshold);
4352   
4353    problemStatus_ = -1;
4354   
4355    // number of times we have declared optimality
4356    numberTimesOptimal_=0;
4357
4358    return 0;
4359  } else {
4360    // bad matrix
4361    return 2;
4362  }
4363   
4364}
4365
4366
4367void 
4368ClpSimplex::finish()
4369{
4370  // Get rid of some arrays and empty factorization
4371  deleteRim();
4372  // Skip message if changing algorithms
4373  if (problemStatus_!=10) {
4374    assert(problemStatus_>=0&&problemStatus_<5);
4375    handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
4376      <<objectiveValue()
4377      <<CoinMessageEol;
4378  }
4379  factorization_->relaxAccuracyCheck(1.0);
4380  // get rid of any network stuff - could do more
4381  factorization_->cleanUp();
4382}
4383// Save data
4384ClpDataSave
4385ClpSimplex::saveData() 
4386{
4387  ClpDataSave saved;
4388  saved.dualBound_ = dualBound_;
4389  saved.infeasibilityCost_ = infeasibilityCost_;
4390  saved.sparseThreshold_ = factorization_->sparseThreshold();
4391  saved.perturbation_ = perturbation_;
4392  // Progress indicator
4393  delete progress_;
4394  progress_ = new ClpSimplexProgress (this);
4395  return saved;
4396}
4397// Restore data
4398void 
4399ClpSimplex::restoreData(ClpDataSave saved)
4400{
4401  factorization_->sparseThreshold(saved.sparseThreshold_);
4402  perturbation_ = saved.perturbation_;
4403  infeasibilityCost_ = saved.infeasibilityCost_;
4404  dualBound_ = saved.dualBound_;
4405  delete progress_;
4406  progress_=NULL;
4407}
4408/* Factorizes and returns true if optimal.  Used by user */
4409bool
4410ClpSimplex::statusOfProblem()
4411{
4412  // is factorization okay?
4413  assert (internalFactorize(1)==0);
4414  // put back original costs and then check
4415  // also move to work arrays
4416  createRim(4+32);
4417  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
4418  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
4419  gutsOfSolution(NULL,NULL);
4420  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
4421  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
4422  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
4423  deleteRim(-1);
4424  return (primalFeasible()&&dualFeasible());
4425}
4426/* Return model - updates any scalars */
4427void 
4428ClpSimplex::returnModel(ClpSimplex & otherModel)
4429{
4430  ClpModel::returnModel(otherModel);
4431  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
4432  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
4433  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
4434  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
4435  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
4436  otherModel.columnDualSequence_ = columnDualSequence_;
4437  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
4438  otherModel.rowDualSequence_ = rowDualSequence_;
4439  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
4440  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
4441  otherModel.largestPrimalError_ = largestPrimalError_;
4442  otherModel.largestDualError_ = largestDualError_;
4443  otherModel.largestSolutionError_ = largestSolutionError_;
4444  otherModel.alpha_ = alpha_;
4445  otherModel.theta_ = theta_;
4446  otherModel.lowerIn_ = lowerIn_;
4447  otherModel.valueIn_ = valueIn_;
4448  otherModel.upperIn_ = upperIn_;
4449  otherModel.dualIn_ = dualIn_;
4450  otherModel.sequenceIn_ = sequenceIn_;
4451  otherModel.directionIn_ = directionIn_;
4452  otherModel.lowerOut_ = lowerOut_;
4453  otherModel.valueOut_ = valueOut_;
4454  otherModel.upperOut_ = upperOut_;
4455  otherModel.dualOut_ = dualOut_;
4456  otherModel.sequenceOut_ = sequenceOut_;
4457  otherModel.directionOut_ = directionOut_;
4458  otherModel.pivotRow_ = pivotRow_;
4459  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
4460  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
4461  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
4462    numberDualInfeasibilitiesWithoutFree_;
4463  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
4464  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
4465  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
4466  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
4467  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
4468}
4469/* Constructs a non linear cost from list of non-linearities (columns only)
4470   First lower of each column is taken as real lower
4471   Last lower is taken as real upper and cost ignored
4472   
4473   Returns nonzero if bad data e.g. lowers not monotonic
4474*/
4475int 
4476ClpSimplex::createPiecewiseLinearCosts(const int * starts,
4477                                       const double * lower, const double * gradient)
4478{
4479  delete nonLinearCost_;
4480  // Set up feasible bounds and check monotonicity
4481  int iColumn;
4482  int returnCode=0;
4483
4484  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4485    int iIndex = starts[iColumn];
4486    int end = starts[iColumn+1]-1;
4487    columnLower_[iColumn] = lower[iIndex];
4488    columnUpper_[iColumn] = lower[end];
4489    double value = columnLower_[iColumn];
4490    iIndex++;
4491    for (;iIndex<end;iIndex++) {
4492      if (lower[iIndex]<value)
4493        returnCode++; // not monotonic
4494      value=lower[iIndex];
4495    }
4496  }
4497  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
4498  specialOptions_ |= 2; // say keep
4499  return returnCode;
4500}
4501/* For advanced use.  When doing iterative solves things can get
4502   nasty so on values pass if incoming solution has largest
4503   infeasibility < incomingInfeasibility throw out variables
4504   from basis until largest infeasibility < allowedInfeasibility
4505   or incoming largest infeasibility.
4506   If allowedInfeasibility>= incomingInfeasibility this is
4507   always possible altough you may end up with an all slack basis.
4508   
4509   Defaults are 1.0,10.0
4510*/
4511void 
4512ClpSimplex::setValuesPassAction(float incomingInfeasibility,
4513                                float allowedInfeasibility)
4514{
4515  incomingInfeasibility_=incomingInfeasibility;
4516  allowedInfeasibility_=allowedInfeasibility;
4517  assert(incomingInfeasibility_>=0.0);
4518  assert(allowedInfeasibility_>=incomingInfeasibility_);
4519}
4520//#############################################################################
4521
4522ClpSimplexProgress::ClpSimplexProgress () 
4523{
4524  int i;
4525  for (i=0;i<CLP_PROGRESS;i++) {
4526    objective_[i] = COIN_DBL_MAX;
4527    infeasibility_[i] = -1.0; // set to an impossible value
4528    numberInfeasibilities_[i]=-1; 
4529    iterationNumber_[i]=-1;
4530  }
4531  for (i=0;i<CLP_CYCLE;i++) {
4532    in_[i]=-1;
4533    out_[i]=-1;
4534    way_[i]=0;
4535  }
4536  numberTimes_ = 0;
4537  numberBadTimes_ = 0;
4538  model_ = NULL;
4539}
4540
4541
4542//-----------------------------------------------------------------------------
4543
4544ClpSimplexProgress::~ClpSimplexProgress ()
4545{
4546}
4547// Copy constructor.
4548ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
4549{
4550  int i;
4551  for (i=0;i<CLP_PROGRESS;i++) {
4552    objective_[i] = rhs.objective_[i];
4553    infeasibility_[i] = rhs.infeasibility_[i];
4554    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
4555    iterationNumber_[i]=rhs.iterationNumber_[i];
4556  }
4557  for (i=0;i<CLP_CYCLE;i++) {
4558    in_[i]=rhs.in_[i];
4559    out_[i]=rhs.out_[i];
4560    way_[i]=rhs.way_[i];
4561  }
4562  numberTimes_ = rhs.numberTimes_;
4563  numberBadTimes_ = rhs.numberBadTimes_;
4564  model_ = rhs.model_;
4565}
4566// Copy constructor.from model
4567ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
4568{
4569  model_ = model;
4570  int i;
4571  for (i=0;i<CLP_PROGRESS;i++) {
4572    if (model_->algorithm()>=0)
4573      objective_[i] = COIN_DBL_MAX;
4574    else
4575      objective_[i] = -COIN_DBL_MAX;
4576    infeasibility_[i] = -1.0; // set to an impossible value
4577    numberInfeasibilities_[i]=-1; 
4578    iterationNumber_[i]=-1;
4579  }
4580  for (i=0;i<CLP_CYCLE;i++) {
4581    in_[i]=-1;
4582    out_[i]=-1;
4583    way_[i]=0;
4584  }
4585  numberTimes_ = 0;
4586  numberBadTimes_ = 0;
4587}
4588// Assignment operator. This copies the data
4589ClpSimplexProgress & 
4590ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
4591{
4592  if (this != &rhs) {
4593    int i;
4594    for (i=0;i<CLP_PROGRESS;i++) {
4595      objective_[i] = rhs.objective_[i];
4596      infeasibility_[i] = rhs.infeasibility_[i];
4597      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
4598      iterationNumber_[i]=rhs.iterationNumber_[i];
4599    }
4600    for (i=0;i<CLP_CYCLE;i++) {
4601      in_[i]=rhs.in_[i];
4602      out_[i]=rhs.out_[i];
4603      way_[i]=rhs.way_[i];
4604    }
4605    numberTimes_ = rhs.numberTimes_;
4606    numberBadTimes_ = rhs.numberBadTimes_;
4607    model_ = rhs.model_;
4608  }
4609  return *this;
4610}
4611// Seems to be something odd about exact comparison of doubles on linux
4612static bool equalDouble(double value1, double value2)
4613{
4614  assert(sizeof(unsigned int)*2==sizeof(double));
4615  unsigned int *i1 = (unsigned int *) &value1;
4616  unsigned int *i2 = (unsigned int *) &value2;
4617  return (i1[0]==i2[0]&&i1[1]==i2[1]);
4618}
4619int
4620ClpSimplexProgress::looping()
4621{
4622  if (!model_)
4623    return -1;
4624  double objective = model_->rawObjectiveValue();
4625  double infeasibility;
4626  int numberInfeasibilities;
4627  int iterationNumber = model_->numberIterations();
4628  if (model_->algorithm()<0) {
4629    // dual
4630    infeasibility = model_->sumPrimalInfeasibilities();
4631    numberInfeasibilities = model_->numberPrimalInfeasibilities();
4632  } else {
4633    //primal
4634    infeasibility = model_->sumDualInfeasibilities();
4635    numberInfeasibilities = model_->numberDualInfeasibilities();
4636  }
4637  int i;
4638  int numberMatched=0;
4639  int matched=0;
4640  int nsame=0;
4641  for (i=0;i<CLP_PROGRESS;i++) {
4642    bool matchedOnObjective = equalDouble(objective,objective_[i]);
4643    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
4644    bool matchedOnInfeasibilities = 
4645      (numberInfeasibilities==numberInfeasibilities_[i]);
4646   
4647    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
4648      matched |= (1<<i);
4649      // Check not same iteration
4650      if (iterationNumber!=iterationNumber_[i]) {
4651        numberMatched++;
4652        // here mainly to get over compiler bug?
4653        if (model_->messageHandler()->logLevel()>10)
4654          printf("%d %d %d %d %d loop check\n",i,numberMatched,
4655                 matchedOnObjective, matchedOnInfeasibility, 
4656                 matchedOnInfeasibilities);
4657      } else {
4658        // stuck but code should notice
4659        nsame++;
4660      }
4661    }
4662    if (i) {
4663      objective_[i-1] = objective_[i];
4664      infeasibility_[i-1] = infeasibility_[i];
4665      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
4666      iterationNumber_[i-1]=iterationNumber_[i];
4667    }
4668  }
4669  objective_[CLP_PROGRESS-1] = objective;
4670  infeasibility_[CLP_PROGRESS-1] = infeasibility;
4671  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
4672  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
4673  if (nsame==CLP_PROGRESS)
4674    numberMatched=CLP_PROGRESS; // really stuck
4675  if (model_->progressFlag())
4676    numberMatched=0;
4677  numberTimes_++;
4678  if (numberTimes_<10)
4679    numberMatched=0;
4680  // skip if just last time as may be checking something
4681  if (matched==(1<<(CLP_PROGRESS-1)))
4682    numberMatched=0;
4683  if (numberMatched) {
4684    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
4685      <<numberMatched
4686      <<matched
4687      <<numberTimes_
4688      <<CoinMessageEol;
4689    numberBadTimes_++;
4690    if (numberBadTimes_<10) {
4691      // make factorize every iteration
4692      model_->forceFactorization(1);
4693      if (model_->algorithm()<0) {
4694        // dual - change tolerance
4695        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
4696        // if infeasible increase dual bound
4697        if (model_->dualBound()<1.0e19) {
4698          model_->setDualBound(model_->dualBound()*1.1);
4699        }
4700      } else {
4701        // primal - change tolerance    model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
4702        // if infeasible increase infeasibility cost
4703        if (model_->nonLinearCost()->numberInfeasibilities()&&
4704            model_->infeasibilityCost()<1.0e19) {
4705          model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
4706        }
4707      }
4708      return -2;
4709    } else {
4710      model_->messageHandler()->message(CLP_LOOP,model_->messages())
4711        <<CoinMessageEol;
4712#ifndef NDEBUG
4713      abort();
4714#endif
4715      // look at solution and maybe declare victory
4716      if (infeasibility<1.0e-4)
4717        return 0;
4718      else
4719        return 3;
4720    }
4721  }
4722  return -1;
4723}
4724// Returns previous objective (if -1) - current if (0)
4725double 
4726ClpSimplexProgress::lastObjective(int back) const
4727{
4728  return objective_[CLP_PROGRESS-1-back];
4729}
4730// Modify objective e.g. if dual infeasible in dual
4731void 
4732ClpSimplexProgress::modifyObjective(double value)
4733{
4734  objective_[CLP_PROGRESS-1]=value;
4735}
4736// Returns previous iteration number (if -1) - current if (0)
4737int 
4738ClpSimplexProgress::lastIterationNumber(int back) const
4739{
4740  return iterationNumber_[CLP_PROGRESS-1-back];
4741}
4742// Start check at beginning of whileIterating
4743void 
4744ClpSimplexProgress::startCheck()
4745{
4746  int i;
4747  for (i=0;i<CLP_CYCLE;i++) {
4748    in_[i]=-1;
4749    out_[i]=-1;
4750    way_[i]=0;
4751  }
4752}
4753// Returns cycle length in whileIterating
4754int 
4755ClpSimplexProgress::cycle(int in, int out,int wayIn,int wayOut)
4756{
4757  int i;
4758  int matched=0;
4759  return 0;
4760  // first see if in matches any out
4761  for (i=1;i<CLP_CYCLE;i++) {
4762    if (in==out_[i]) {
4763      // even if flip then suspicious
4764      matched=-1;
4765      break;
4766    }
4767  }
4768  if (!matched) {
4769    // can't be cycle
4770    for (i=0;i<CLP_CYCLE-1;i++) {
4771      in_[i]=in_[i+1];
4772      out_[i]=out_[i+1];
4773      way_[i]=way_[i+1];
4774    }
4775  } else {
4776    // possible cycle
4777    matched=0;
4778    for (i=0;i<CLP_CYCLE-1;i++) {
4779      int k;
4780      char wayThis = way_[i];
4781      int inThis = in_[i];
4782      int outThis = out_[i];
4783      for(k=i+1;k<CLP_CYCLE;k++) {
4784        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
4785          matched=k-i;
4786        }
4787      }
4788      in_[i]=in_[i+1];
4789      out_[i]=out_[i+1];
4790      way_[i]=way_[i+1];
4791    }
4792  }
4793  char way = 1-wayIn+4*(1-wayOut);
4794  in_[CLP_CYCLE-1]=in;
4795  out_[CLP_CYCLE-1]=out;
4796  way_[CLP_CYCLE-1]=way;
4797  return matched;
4798}
Note: See TracBrowser for help on using the repository browser.