source: branches/pre/ClpSimplex.cpp @ 1926

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

Stuff

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