source: branches/pre/ClpSimplex.cpp @ 210

Last change on this file since 210 was 210, checked in by forrest, 17 years ago

Trying

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