source: branches/pre/ClpInterior.cpp @ 212

Last change on this file since 212 was 212, checked in by forrest, 16 years ago

lots of stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 134.1 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6
7#include "CoinPragma.hpp"
8
9#include <math.h>
10
11#include "CoinHelperFunctions.hpp"
12#include "ClpInterior.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
30ClpInterior::ClpInterior () :
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
130ClpInterior::ClpInterior ( 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
232ClpInterior::~ClpInterior ()
233{
234  gutsOfDelete(0);
235  delete nonLinearCost_;
236}
237//#############################################################################
238void ClpInterior::setLargeValue( double value) 
239{
240  if (value>0.0&&value<COIN_DBL_MAX)
241    largeValue_=value;
242}
243int
244ClpInterior::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
359ClpInterior::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
492ClpInterior::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 ClpInterior::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 ClpInterior::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 ClpInterior::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 
697ClpInterior::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 ClpInterior::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 ClpInterior::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 
1142ClpInterior::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.
1238ClpInterior::ClpInterior(const ClpInterior &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
1335ClpInterior::ClpInterior(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
1431ClpInterior & 
1432ClpInterior::operator=(const ClpInterior & 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 
1445ClpInterior::gutsOfCopy(const ClpInterior & 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 ClpInteriorProgress(*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 
1563ClpInterior::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 
1631ClpInterior::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 
1704ClpInterior::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 
1886ClpInterior::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 
1898ClpInterior::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 
1913ClpInterior::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 
1930ClpInterior::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
1947ClpInterior::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
2189ClpInterior::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 
2237ClpInterior::setDualBound(double value)
2238{
2239  if (value>0.0)
2240    dualBound_=value;
2241}
2242void 
2243ClpInterior::setInfeasibilityCost(double value)
2244{
2245  if (value>0.0)
2246    infeasibilityCost_=value;
2247}
2248void ClpInterior::setNumberRefinements( int value) 
2249{
2250  if (value>=0&&value<10)
2251    numberRefinements_=value;
2252}
2253// Sets row pivot choice algorithm in dual
2254void 
2255ClpInterior::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
2256{
2257  delete dualRowPivot_;
2258  dualRowPivot_ = choice.clone(true);
2259}
2260// Sets row pivot choice algorithm in dual
2261void 
2262ClpInterior::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 
2269ClpInterior::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 
2283ClpInterior::setFactorization( ClpFactorization & factorization)
2284{
2285  delete factorization_;
2286  factorization_= new ClpFactorization(factorization);
2287}
2288void 
2289ClpInterior::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 
2298ClpInterior::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 
2314ClpInterior::setPerturbation(int value)
2315{
2316  if(value<=100&&value >=-1000) {
2317    perturbation_=value;
2318  } 
2319}
2320// Sparsity on or off
2321bool 
2322ClpInterior::sparseFactorization() const
2323{
2324  return factorization_->sparseThreshold()!=0;
2325}
2326void 
2327ClpInterior::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 
2346ClpInterior::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 "ClpInteriorDual.hpp"
2670#include "ClpInteriorPrimal.hpp"
2671int ClpInterior::dual (int ifValuesPass )
2672{
2673  int returnCode = ((ClpInteriorDual *) this)->dual(ifValuesPass);
2674  if (problemStatus_==10) {
2675    //printf("Cleaning up with primal\n");
2676    int savePerturbation = perturbation_;
2677    perturbation_=100;
2678    bool denseFactorization = initialDenseFactorization();
2679    // It will be safe to allow dense
2680    setInitialDenseFactorization(true);
2681    returnCode = ((ClpInteriorPrimal *) this)->primal(1);
2682    setInitialDenseFactorization(denseFactorization);
2683    perturbation_=savePerturbation;
2684  }
2685  return returnCode;
2686}
2687// primal
2688int ClpInterior::primal (int ifValuesPass )
2689{
2690  int returnCode = ((ClpInteriorPrimal *) this)->primal(ifValuesPass);
2691  if (problemStatus_==10) {
2692    //printf("Cleaning up with dual\n");
2693    int savePerturbation = perturbation_;
2694    perturbation_=100;
2695    bool denseFactorization = initialDenseFactorization();
2696    // It will be safe to allow dense
2697    setInitialDenseFactorization(true);
2698    returnCode = ((ClpInteriorDual *) this)->dual(0);
2699    setInitialDenseFactorization(denseFactorization);
2700    perturbation_=savePerturbation;
2701  }
2702  return returnCode;
2703}
2704#include "ClpInteriorPrimalQuadratic.hpp"
2705/* Solves quadratic problem using SLP - may be used as crash
2706   for other algorithms when number of iterations small
2707*/
2708int 
2709ClpInterior::quadraticSLP(int numberPasses, double deltaTolerance)
2710{
2711  return ((ClpInteriorPrimalQuadratic *) this)->primalSLP(numberPasses,deltaTolerance);
2712}
2713// Solves quadratic using Dantzig's algorithm - primal
2714int 
2715ClpInterior::quadraticPrimal(int phase)
2716{
2717  return ((ClpInteriorPrimalQuadratic *) this)->primalQuadratic(phase);
2718}
2719/* For strong branching.  On input lower and upper are new bounds
2720   while on output they are objective function values (>1.0e50 infeasible).
2721   Return code is 0 if nothing interesting, -1 if infeasible both
2722   ways and +1 if infeasible one way (check values to see which one(s))
2723*/
2724int ClpInterior::strongBranching(int numberVariables,const int * variables,
2725                                double * newLower, double * newUpper,
2726                                double ** outputSolution,
2727                                int * outputStatus, int * outputIterations,
2728                                bool stopOnFirstInfeasible,
2729                                bool alwaysFinish)
2730{
2731  return ((ClpInteriorDual *) this)->strongBranching(numberVariables,variables,
2732                                                    newLower,  newUpper,outputSolution,
2733                                                    outputStatus, outputIterations,
2734                                                    stopOnFirstInfeasible,
2735                                                    alwaysFinish);
2736}
2737/* Borrow model.  This is so we dont have to copy large amounts
2738   of data around.  It assumes a derived class wants to overwrite
2739   an empty model with a real one - while it does an algorithm.
2740   This is same as ClpModel one, but sets scaling on etc. */
2741void 
2742ClpInterior::borrowModel(ClpModel & otherModel) 
2743{
2744  ClpModel::borrowModel(otherModel);
2745  createStatus();
2746  ClpDualRowSteepest steep1;
2747  setDualRowPivotAlgorithm(steep1);
2748  ClpPrimalColumnSteepest steep2;
2749  setPrimalColumnPivotAlgorithm(steep2);
2750}
2751typedef struct {
2752  double optimizationDirection;
2753  double dblParam[ClpLastDblParam];
2754  double objectiveValue;
2755  double dualBound;
2756  double dualTolerance;
2757  double primalTolerance;
2758  double sumDualInfeasibilities;
2759  double sumPrimalInfeasibilities;
2760  double infeasibilityCost;
2761  int numberRows;
2762  int numberColumns;
2763  int intParam[ClpLastIntParam];
2764  int numberIterations;
2765  int problemStatus;
2766  int maximumIterations;
2767  int lengthNames;
2768  int numberDualInfeasibilities;
2769  int numberDualInfeasibilitiesWithoutFree;
2770  int numberPrimalInfeasibilities;
2771  int numberRefinements;
2772  int scalingFlag;
2773  int algorithm;
2774  unsigned int specialOptions;
2775  int dualPivotChoice;
2776  int primalPivotChoice;
2777  int matrixStorageChoice;
2778} Clp_scalars;
2779
2780int outDoubleArray(double * array, int length, FILE * fp)
2781{
2782  int numberWritten;
2783  if (array&&length) {
2784    numberWritten = fwrite(&length,sizeof(int),1,fp);
2785    if (numberWritten!=1)
2786      return 1;
2787    numberWritten = fwrite(array,sizeof(double),length,fp);
2788    if (numberWritten!=length)
2789      return 1;
2790  } else {
2791    length = 0;
2792    numberWritten = fwrite(&length,sizeof(int),1,fp);
2793    if (numberWritten!=1)
2794      return 1;
2795  }
2796  return 0;
2797}
2798// Save model to file, returns 0 if success
2799int
2800ClpInterior::saveModel(const char * fileName)
2801{
2802  FILE * fp = fopen(fileName,"wb");
2803  if (fp) {
2804    Clp_scalars scalars;
2805    int i;
2806    CoinBigIndex numberWritten;
2807    // Fill in scalars
2808    scalars.optimizationDirection = optimizationDirection_;
2809    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2810    scalars.objectiveValue = objectiveValue_;
2811    scalars.dualBound = dualBound_;
2812    scalars.dualTolerance = dualTolerance_;
2813    scalars.primalTolerance = primalTolerance_;
2814    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2815    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2816    scalars.infeasibilityCost = infeasibilityCost_;
2817    scalars.numberRows = numberRows_;
2818    scalars.numberColumns = numberColumns_;
2819    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2820    scalars.numberIterations = numberIterations_;
2821    scalars.problemStatus = problemStatus_;
2822    scalars.maximumIterations = maximumIterations();
2823    scalars.lengthNames = lengthNames_;
2824    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2825    scalars.numberDualInfeasibilitiesWithoutFree
2826      = numberDualInfeasibilitiesWithoutFree_;
2827    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2828    scalars.numberRefinements = numberRefinements_;
2829    scalars.scalingFlag = scalingFlag_;
2830    scalars.algorithm = algorithm_;
2831    scalars.specialOptions = specialOptions_;
2832    scalars.dualPivotChoice = dualRowPivot_->type();
2833    scalars.primalPivotChoice = primalColumnPivot_->type();
2834    scalars.matrixStorageChoice = matrix_->type();
2835
2836    // put out scalars
2837    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2838    if (numberWritten!=1)
2839      return 1;
2840    // strings
2841    CoinBigIndex length;
2842    for (i=0;i<ClpLastStrParam;i++) {
2843      length = strParam_[i].size();
2844      numberWritten = fwrite(&length,sizeof(int),1,fp);
2845      if (numberWritten!=1)
2846        return 1;
2847      if (length) {
2848        numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2849        if (numberWritten!=1)
2850          return 1;
2851      }
2852    }
2853    // arrays - in no particular order
2854    if (outDoubleArray(rowActivity_,numberRows_,fp))
2855        return 1;
2856    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2857        return 1;
2858    if (outDoubleArray(dual_,numberRows_,fp))
2859        return 1;
2860    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2861        return 1;
2862    if (outDoubleArray(rowLower_,numberRows_,fp))
2863        return 1;
2864    if (outDoubleArray(rowUpper_,numberRows_,fp))
2865        return 1;
2866    if (outDoubleArray(objective(),numberColumns_,fp))
2867        return 1;
2868    if (outDoubleArray(rowObjective_,numberRows_,fp))
2869        return 1;
2870    if (outDoubleArray(columnLower_,numberColumns_,fp))
2871        return 1;
2872    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2873        return 1;
2874    if (ray_) {
2875      if (problemStatus_==1)
2876        if (outDoubleArray(ray_,numberRows_,fp))
2877          return 1;
2878      else if (problemStatus_==2)
2879        if (outDoubleArray(ray_,numberColumns_,fp))
2880          return 1;
2881      else
2882        if (outDoubleArray(NULL,0,fp))
2883          return 1;
2884    } else {
2885      if (outDoubleArray(NULL,0,fp))
2886        return 1;
2887    }
2888    if (status_&&(numberRows_+numberColumns_)>0) {
2889      length = numberRows_+numberColumns_;
2890      numberWritten = fwrite(&length,sizeof(int),1,fp);
2891      if (numberWritten!=1)
2892        return 1;
2893      numberWritten = fwrite(status_,sizeof(char),length, fp);
2894      if (numberWritten!=length)
2895        return 1;
2896    } else {
2897      length = 0;
2898      numberWritten = fwrite(&length,sizeof(int),1,fp);
2899      if (numberWritten!=1)
2900        return 1;
2901    }
2902    if (lengthNames_) {
2903      char * array = 
2904        new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
2905      char * put = array;
2906      assert (numberRows_ == (int) rowNames_.size());
2907      for (i=0;i<numberRows_;i++) {
2908        assert((int)rowNames_[i].size()<=lengthNames_);
2909        strcpy(put,rowNames_[i].c_str());
2910        put += lengthNames_+1;
2911      }
2912      numberWritten = fwrite(array,lengthNames_+1,numberRows_,fp);
2913      if (numberWritten!=numberRows_)
2914        return 1;
2915      put=array;
2916      assert (numberColumns_ == (int) columnNames_.size());
2917      for (i=0;i<numberColumns_;i++) {
2918        assert((int) columnNames_[i].size()<=lengthNames_);
2919        strcpy(put,columnNames_[i].c_str());
2920        put += lengthNames_+1;
2921      }
2922      numberWritten = fwrite(array,lengthNames_+1,numberColumns_,fp);
2923      if (numberWritten!=numberColumns_)
2924        return 1;
2925      delete [] array;
2926    }
2927    // just standard type at present
2928    assert (matrix_->type()==1);
2929    assert (matrix_->getNumCols() == numberColumns_);
2930    assert (matrix_->getNumRows() == numberRows_);
2931    // we are going to save with gaps
2932    length = matrix_->getVectorStarts()[numberColumns_-1]
2933      + matrix_->getVectorLengths()[numberColumns_-1];
2934    numberWritten = fwrite(&length,sizeof(int),1,fp);
2935    if (numberWritten!=1)
2936      return 1;
2937    numberWritten = fwrite(matrix_->getElements(),
2938                           sizeof(double),length,fp);
2939    if (numberWritten!=length)
2940      return 1;
2941    numberWritten = fwrite(matrix_->getIndices(),
2942                           sizeof(int),length,fp);
2943    if (numberWritten!=length)
2944      return 1;
2945    numberWritten = fwrite(matrix_->getVectorStarts(),
2946                           sizeof(int),numberColumns_+1,fp);
2947    if (numberWritten!=numberColumns_+1)
2948      return 1;
2949    numberWritten = fwrite(matrix_->getVectorLengths(),
2950                           sizeof(int),numberColumns_,fp);
2951    if (numberWritten!=numberColumns_)
2952      return 1;
2953    // finished
2954    fclose(fp);
2955    return 0;
2956  } else {
2957    return -1;
2958  }
2959}
2960
2961int inDoubleArray(double * &array, int length, FILE * fp)
2962{
2963  int numberRead;
2964  int length2;
2965  numberRead = fread(&length2,sizeof(int),1,fp);
2966  if (numberRead!=1)
2967    return 1;
2968  if (length2) {
2969    // lengths must match
2970    if (length!=length2)
2971      return 2;
2972    array = new double[length];
2973    numberRead = fread(array,sizeof(double),length,fp);
2974    if (numberRead!=length)
2975      return 1;
2976  } 
2977  return 0;
2978}
2979/* Restore model from file, returns 0 if success,
2980   deletes current model */
2981int 
2982ClpInterior::restoreModel(const char * fileName)
2983{
2984  FILE * fp = fopen(fileName,"rb");
2985  if (fp) {
2986    // Get rid of current model
2987    ClpModel::gutsOfDelete();
2988    gutsOfDelete(0);
2989    int i;
2990    for (i=0;i<6;i++) {
2991      rowArray_[i]=NULL;
2992      columnArray_[i]=NULL;
2993    }
2994    // get an empty factorization so we can set tolerances etc
2995    factorization_ = new ClpFactorization();
2996    // Say sparse
2997    factorization_->sparseThreshold(1);
2998    Clp_scalars scalars;
2999    CoinBigIndex numberRead;
3000
3001    // get scalars
3002    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
3003    if (numberRead!=1)
3004      return 1;
3005    // Fill in scalars
3006    optimizationDirection_ = scalars.optimizationDirection;
3007    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
3008    objectiveValue_ = scalars.objectiveValue;
3009    dualBound_ = scalars.dualBound;
3010    dualTolerance_ = scalars.dualTolerance;
3011    primalTolerance_ = scalars.primalTolerance;
3012    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
3013    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
3014    infeasibilityCost_ = scalars.infeasibilityCost;
3015    numberRows_ = scalars.numberRows;
3016    numberColumns_ = scalars.numberColumns;
3017    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
3018    numberIterations_ = scalars.numberIterations;
3019    problemStatus_ = scalars.problemStatus;
3020    setMaximumIterations(scalars.maximumIterations);
3021    lengthNames_ = scalars.lengthNames;
3022    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
3023    numberDualInfeasibilitiesWithoutFree_
3024      = scalars.numberDualInfeasibilitiesWithoutFree;
3025    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
3026    numberRefinements_ = scalars.numberRefinements;
3027    scalingFlag_ = scalars.scalingFlag;
3028    algorithm_ = scalars.algorithm;
3029    specialOptions_ = scalars.specialOptions;
3030    // strings
3031    CoinBigIndex length;
3032    for (i=0;i<ClpLastStrParam;i++) {
3033      numberRead = fread(&length,sizeof(int),1,fp);
3034      if (numberRead!=1)
3035        return 1;
3036      if (length) {
3037        char * array = new char[length+1];
3038        numberRead = fread(array,length,1,fp);
3039        if (numberRead!=1)
3040          return 1;
3041        array[length]='\0';
3042        strParam_[i]=array;
3043        delete [] array;
3044      }
3045    }
3046    // arrays - in no particular order
3047    if (inDoubleArray(rowActivity_,numberRows_,fp))
3048        return 1;
3049    if (inDoubleArray(columnActivity_,numberColumns_,fp))
3050        return 1;
3051    if (inDoubleArray(dual_,numberRows_,fp))
3052        return 1;
3053    if (inDoubleArray(reducedCost_,numberColumns_,fp))
3054        return 1;
3055    if (inDoubleArray(rowLower_,numberRows_,fp))
3056        return 1;
3057    if (inDoubleArray(rowUpper_,numberRows_,fp))
3058        return 1;
3059    double * objective;
3060    if (inDoubleArray(objective,numberColumns_,fp))
3061        return 1;
3062    delete objective_;
3063    objective_ = new ClpLinearObjective(objective,numberColumns_);
3064    delete [] objective;
3065    if (inDoubleArray(rowObjective_,numberRows_,fp))
3066        return 1;
3067    if (inDoubleArray(columnLower_,numberColumns_,fp))
3068        return 1;
3069    if (inDoubleArray(columnUpper_,numberColumns_,fp))
3070        return 1;
3071    if (problemStatus_==1) {
3072      if (inDoubleArray(ray_,numberRows_,fp))
3073        return 1;
3074    } else if (problemStatus_==2) {
3075      if (inDoubleArray(ray_,numberColumns_,fp))
3076        return 1;
3077    } else {
3078      // ray should be null
3079      numberRead = fread(&length,sizeof(int),1,fp);
3080      if (numberRead!=1)
3081        return 1;
3082      if (length)
3083        return 2;
3084    }
3085    delete [] status_;
3086    status_=NULL;
3087    // status region
3088    numberRead = fread(&length,sizeof(int),1,fp);
3089    if (numberRead!=1)
3090        return 1;
3091    if (length) {
3092      if (length!=numberRows_+numberColumns_)
3093        return 1;
3094      status_ = new char unsigned[length];
3095      numberRead = fread(status_,sizeof(char),length, fp);
3096      if (numberRead!=length)
3097        return 1;
3098    }
3099    if (lengthNames_) {
3100      char * array = 
3101        new char[max(numberRows_,numberColumns_)*(lengthNames_+1)];
3102      char * get = array;
3103      numberRead = fread(array,lengthNames_+1,numberRows_,fp);
3104      if (numberRead!=numberRows_)
3105        return 1;
3106      rowNames_ = std::vector<std::string> ();
3107      rowNames_.resize(numberRows_);
3108      for (i=0;i<numberRows_;i++) {
3109        rowNames_[i]=get;
3110        get += lengthNames_+1;
3111      }
3112      get = array;
3113      numberRead = fread(array,lengthNames_+1,numberColumns_,fp);
3114      if (numberRead!=numberColumns_)
3115        return 1;
3116      columnNames_ = std::vector<std::string> ();
3117      columnNames_.resize(numberColumns_);
3118      for (i=0;i<numberColumns_;i++) {
3119        columnNames_[i]=get;
3120        get += lengthNames_+1;
3121      }
3122      delete [] array;
3123    }
3124    // Pivot choices
3125    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
3126    delete dualRowPivot_;
3127    switch ((scalars.dualPivotChoice&63)) {
3128    default:
3129      printf("Need another dualPivot case %d\n",scalars.dualPivotChoice&63);
3130    case 1:
3131      // Dantzig
3132      dualRowPivot_ = new ClpDualRowDantzig();
3133      break;
3134    case 2:
3135      // Steepest - use mode
3136      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
3137      break;
3138    }
3139    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
3140    delete primalColumnPivot_;
3141    switch ((scalars.primalPivotChoice&63)) {
3142    default:
3143      printf("Need another primalPivot case %d\n",
3144             scalars.primalPivotChoice&63);
3145    case 1:
3146      // Dantzig
3147      primalColumnPivot_ = new ClpPrimalColumnDantzig();
3148      break;
3149    case 2:
3150      // Steepest - use mode
3151      primalColumnPivot_
3152        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
3153      break;
3154    }
3155    assert(scalars.matrixStorageChoice==1);
3156    delete matrix_;
3157    // get arrays
3158    numberRead = fread(&length,sizeof(int),1,fp);
3159    if (numberRead!=1)
3160      return 1;
3161    double * elements = new double[length];
3162    int * indices = new int[length];
3163    CoinBigIndex * starts = new CoinBigIndex[numberColumns_];
3164    int * lengths = new int[numberColumns_];
3165    numberRead = fread(elements, sizeof(double),length,fp);
3166    if (numberRead!=length)
3167      return 1;
3168    numberRead = fread(indices, sizeof(int),length,fp);
3169    if (numberRead!=length)
3170      return 1;
3171    numberRead = fread(starts, sizeof(int),numberColumns_+1,fp);
3172    if (numberRead!=numberColumns_+1)
3173      return 1;
3174    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
3175    if (numberRead!=numberColumns_)
3176      return 1;
3177    // assign matrix
3178    CoinPackedMatrix * matrix = new CoinPackedMatrix();
3179    matrix->assignMatrix(true, numberRows_, numberColumns_,
3180                         length, elements, indices, starts, lengths);
3181    // and transfer to Clp
3182    matrix_ = new ClpPackedMatrix(matrix);
3183    // finished
3184    fclose(fp);
3185    return 0;
3186  } else {
3187    return -1;
3188  }
3189  return 0;
3190}
3191// value of incoming variable (in Dual)
3192double 
3193ClpInterior::valueIncomingDual() const
3194{
3195  // Need value of incoming for list of infeasibilities as may be infeasible
3196  double valueIncoming = (dualOut_/alpha_)*directionOut_;
3197  if (directionIn_==-1)
3198    valueIncoming = upperIn_-valueIncoming;
3199  else
3200    valueIncoming = lowerIn_-valueIncoming;
3201  return valueIncoming;
3202}
3203// Sanity check on input data - returns true if okay
3204bool 
3205ClpInterior::sanityCheck()
3206{
3207  // bad if empty
3208  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
3209    handler_->message(CLP_EMPTY_PROBLEM,messages_)
3210      <<numberRows_
3211      <<numberColumns_
3212      <<matrix_->getNumElements()
3213      <<CoinMessageEol;
3214    problemStatus_=4;
3215    return false;
3216  }
3217  int numberBad ;
3218  double largestBound, smallestBound, minimumGap;
3219  double smallestObj, largestObj;
3220  int firstBad;
3221  int modifiedBounds=0;
3222  int i;
3223  numberBad=0;
3224  firstBad=-1;
3225  minimumGap=1.0e100;
3226  smallestBound=1.0e100;
3227  largestBound=0.0;
3228  smallestObj=1.0e100;
3229  largestObj=0.0;
3230  // If bounds are too close - fix
3231  double fixTolerance = 1.1*primalTolerance_;
3232  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
3233    double value;
3234    value = fabs(cost_[i]);
3235    if (value>1.0e50) {
3236      numberBad++;
3237      if (firstBad<0)
3238        firstBad=i;
3239    } else if (value) {
3240      if (value>largestObj)
3241        largestObj=value;
3242      if (value<smallestObj)
3243        smallestObj=value;
3244    }
3245    value=upper_[i]-lower_[i];
3246    if (value<-primalTolerance_) {
3247      numberBad++;
3248      if (firstBad<0)
3249        firstBad=i;
3250    } else if (value<=fixTolerance) {
3251      if (value) {
3252        // modify
3253        upper_[i] = lower_[i];
3254        modifiedBounds++;
3255      }
3256    } else {
3257      if (value<minimumGap)
3258        minimumGap=value;
3259    }
3260    if (lower_[i]>-1.0e100&&lower_[i]) {
3261      value = fabs(lower_[i]);
3262      if (value>largestBound)
3263        largestBound=value;
3264      if (value<smallestBound)
3265        smallestBound=value;
3266    }
3267    if (upper_[i]<1.0e100&&upper_[i]) {
3268      value = fabs(upper_[i]);
3269      if (value>largestBound)
3270        largestBound=value;
3271      if (value<smallestBound)
3272        smallestBound=value;
3273    }
3274  }
3275  if (largestBound)
3276    handler_->message(CLP_RIMSTATISTICS3,messages_)
3277      <<smallestBound
3278      <<largestBound
3279      <<minimumGap
3280      <<CoinMessageEol;
3281  minimumGap=1.0e100;
3282  smallestBound=1.0e100;
3283  largestBound=0.0;
3284  for (i=0;i<numberColumns_;i++) {
3285    double value;
3286    value = fabs(cost_[i]);
3287    if (value>1.0e50) {
3288      numberBad++;
3289      if (firstBad<0)
3290        firstBad=i;
3291    } else if (value) {
3292      if (value>largestObj)
3293        largestObj=value;
3294      if (value<smallestObj)
3295        smallestObj=value;
3296    }
3297    value=upper_[i]-lower_[i];
3298    if (value<-primalTolerance_) {
3299      numberBad++;
3300      if (firstBad<0)
3301        firstBad=i;
3302    } else if (value<=fixTolerance) {
3303      if (value) {
3304        // modify
3305        upper_[i] = lower_[i];
3306        modifiedBounds++;
3307      }
3308    } else {
3309      if (value<minimumGap)
3310        minimumGap=value;
3311    }
3312    if (lower_[i]>-1.0e100&&lower_[i]) {
3313      value = fabs(lower_[i]);
3314      if (value>largestBound)
3315        largestBound=value;
3316      if (value<smallestBound)
3317        smallestBound=value;
3318    }
3319    if (upper_[i]<1.0e100&&upper_[i]) {
3320      value = fabs(upper_[i]);
3321      if (value>largestBound)
3322        largestBound=value;
3323      if (value<smallestBound)
3324        smallestBound=value;
3325    }
3326  }
3327  char rowcol[]={'R','C'};
3328  if (numberBad) {
3329    handler_->message(CLP_BAD_BOUNDS,messages_)
3330      <<numberBad
3331      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
3332      <<CoinMessageEol;
3333    problemStatus_=4;
3334    return false;
3335  }
3336  if (modifiedBounds)
3337    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
3338      <<modifiedBounds
3339      <<CoinMessageEol;
3340  handler_->message(CLP_RIMSTATISTICS1,messages_)
3341    <<smallestObj
3342    <<largestObj
3343    <<CoinMessageEol;  if (largestBound)
3344    handler_->message(CLP_RIMSTATISTICS2,messages_)
3345      <<smallestBound
3346      <<largestBound
3347      <<minimumGap
3348      <<CoinMessageEol;
3349  return true;
3350}
3351// Set up status array (for OsiClp)
3352void 
3353ClpInterior::createStatus() 
3354{
3355  if(!status_)
3356    status_ = new unsigned char [numberColumns_+numberRows_];
3357  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
3358  int i;
3359  // set column status to one nearest zero
3360  for (i=0;i<numberColumns_;i++) {
3361#if 0
3362    if (columnLower_[i]>=0.0) {
3363      setColumnStatus(i,atLowerBound);
3364    } else if (columnUpper_[i]<=0.0) {
3365      setColumnStatus(i,atUpperBound);
3366    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
3367      // free
3368      setColumnStatus(i,isFree);
3369    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
3370      setColumnStatus(i,atLowerBound);
3371    } else {
3372      setColumnStatus(i,atUpperBound);
3373    }
3374#else
3375    setColumnStatus(i,atLowerBound);
3376#endif
3377  }
3378  for (i=0;i<numberRows_;i++) {
3379    setRowStatus(i,basic);
3380  }
3381}
3382/* Loads a problem (the constraints on the
3383   rows are given by lower and upper bounds). If a pointer is 0 then the
3384   following values are the default:
3385   <ul>
3386   <li> <code>colub</code>: all columns have upper bound infinity
3387   <li> <code>collb</code>: all columns have lower bound 0
3388   <li> <code>rowub</code>: all rows have upper bound infinity
3389   <li> <code>rowlb</code>: all rows have lower bound -infinity
3390   <li> <code>obj</code>: all variables have 0 objective coefficient
3391   </ul>
3392*/
3393void 
3394ClpInterior::loadProblem (  const ClpMatrixBase& matrix,
3395                    const double* collb, const double* colub,   
3396                    const double* obj,
3397                    const double* rowlb, const double* rowub,
3398                    const double * rowObjective)
3399{
3400  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3401                        rowObjective);
3402  createStatus();
3403}
3404void 
3405ClpInterior::loadProblem (  const CoinPackedMatrix& matrix,
3406                    const double* collb, const double* colub,   
3407                    const double* obj,
3408                    const double* rowlb, const double* rowub,
3409                    const double * rowObjective)
3410{
3411  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3412                        rowObjective);
3413  createStatus();
3414}
3415
3416/* Just like the other loadProblem() method except that the matrix is
3417   given in a standard column major ordered format (without gaps). */
3418void 
3419ClpInterior::loadProblem (  const int numcols, const int numrows,
3420                    const CoinBigIndex* start, const int* index,
3421                    const double* value,
3422                    const double* collb, const double* colub,   
3423                    const double* obj,
3424                    const double* rowlb, const double* rowub,
3425                    const double * rowObjective)
3426{
3427  ClpModel::loadProblem(numcols, numrows, start, index, value,
3428                          collb, colub, obj, rowlb, rowub,
3429                          rowObjective);
3430  createStatus();
3431}
3432void 
3433ClpInterior::loadProblem (  const int numcols, const int numrows,
3434                           const CoinBigIndex* start, const int* index,
3435                           const double* value,const int * length,
3436                           const double* collb, const double* colub,   
3437                           const double* obj,
3438                           const double* rowlb, const double* rowub,
3439                           const double * rowObjective)
3440{
3441  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3442                          collb, colub, obj, rowlb, rowub,
3443                          rowObjective);
3444  createStatus();
3445}
3446// Read an mps file from the given filename
3447int 
3448ClpInterior::readMps(const char *filename,
3449            bool keepNames,
3450            bool ignoreErrors)
3451{
3452  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3453  createStatus();
3454  return status;
3455}
3456// Just check solution (for external use)
3457void 
3458ClpInterior::checkSolution()
3459{
3460  // put in standard form
3461  createRim(7+8+16);
3462  dualTolerance_=dblParam_[ClpDualTolerance];
3463  primalTolerance_=dblParam_[ClpPrimalTolerance];
3464  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
3465  checkDualSolution();
3466  if (!numberDualInfeasibilities_&&
3467      !numberPrimalInfeasibilities_)
3468    problemStatus_=0;
3469  else
3470    problemStatus_=-1;
3471#ifdef CLP_DEBUG
3472  int i;
3473  double value=0.0;
3474  for (i=0;i<numberRows_+numberColumns_;i++)
3475    value += dj_[i]*solution_[i];
3476  printf("dual value %g, primal %g\n",value,objectiveValue());
3477#endif
3478  // release extra memory
3479  deleteRim(0);
3480}
3481/* Crash - at present just aimed at dual, returns
3482   -2 if dual preferred and crash basis created
3483   -1 if dual preferred and all slack basis preferred
3484   0 if basis going in was not all slack
3485   1 if primal preferred and all slack basis preferred
3486   2 if primal preferred and crash basis created.
3487   
3488   if gap between bounds <="gap" variables can be flipped
3489   
3490   If "pivot" is
3491   0 No pivoting (so will just be choice of algorithm)
3492   1 Simple pivoting e.g. gub
3493   2 Mini iterations
3494*/
3495int 
3496ClpInterior::crash(double gap,int pivot)
3497{
3498  assert(!rowObjective_); // not coded
3499  int iColumn;
3500  int numberBad=0;
3501  int numberBasic=0;
3502  double dualTolerance=dblParam_[ClpDualTolerance];
3503  //double primalTolerance=dblParam_[ClpPrimalTolerance];
3504
3505  // If no basis then make all slack one
3506  if (!status_)
3507    createStatus();
3508 
3509  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
3510    if (getColumnStatus(iColumn)==basic)
3511      numberBasic++;
3512  }
3513  if (numberBasic) {
3514    // not all slack
3515    return 0;
3516  } else {
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        return 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 ClpInterior::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 ClpInterior::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        return -1;
3926      }
3927    } else {
3928      delete [] dj;
3929      return -1;
3930    }
3931  }
3932}
3933/* Pivot in a variable and out a variable.  Returns 0 if okay,
3934   1 if inaccuracy forced re-factorization, -1 if would be singular.
3935   Also updates primal/dual infeasibilities.
3936   Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
3937*/
3938int ClpInterior::pivot()
3939{
3940  // scaling not allowed
3941  assert (!scalingFlag_);
3942  // assume In_ and Out_ are correct and directionOut_ set
3943  // (or In_ if flip
3944  lowerIn_ = lower_[sequenceIn_];
3945  valueIn_ = solution_[sequenceIn_];
3946  upperIn_ = upper_[sequenceIn_];
3947  dualIn_ = dj_[sequenceIn_];
3948  if (sequenceOut_>=0&&sequenceIn_!=sequenceIn_) {
3949    assert (pivotRow_>=0&&pivotRow_<numberRows_);
3950    assert (pivotVariable_[pivotRow_]==sequenceOut_);
3951    lowerOut_ = lower_[sequenceOut_];
3952    valueOut_ = solution_[sequenceOut_];
3953    upperOut_ = upper_[sequenceOut_];
3954    // for now assume primal is feasible (or in dual)
3955    dualOut_ = dj_[sequenceOut_];
3956    assert(fabs(dualOut_)<1.0e-6);
3957  } else {
3958    assert (pivotRow_<0);
3959  }
3960  bool roundAgain = true;
3961  int returnCode=0;
3962  while (roundAgain) {
3963    roundAgain=false;
3964    unpack(rowArray_[1]);
3965    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
3966    // we are going to subtract movement from current basic
3967    double movement;
3968    // see where incoming will go to
3969    if (sequenceOut_<0||sequenceIn_==sequenceOut_) {
3970      // flip so go to bound
3971      movement  = ((directionIn_>0) ? upperIn_ : lowerIn_) - valueIn_;
3972    } else {
3973      // get where outgoing needs to get to
3974      double outValue = (directionOut_>0) ? upperOut_ : lowerOut_;
3975      // solutionOut_ - movement*alpha_ == outValue
3976      movement = (outValue-valueOut_)/alpha_;
3977      // set directionIn_ correctly
3978      directionIn_ = (movement>0) ? 1 :-1;
3979    }
3980    // update primal solution
3981    {
3982      int i;
3983      int * index = rowArray_[1]->getIndices();
3984      int number = rowArray_[1]->getNumElements();
3985      double * element = rowArray_[1]->denseVector();
3986      for (i=0;i<number;i++) {
3987        int ii = index[i];
3988        // get column
3989        ii = pivotVariable_[ii];
3990        solution_[ii] -= movement*element[i];
3991        element[i]=0.0;
3992      }
3993      // see where something went to
3994      if (sequenceOut_<0) {
3995        if (directionIn_<0) {
3996          assert (fabs(solution_[sequenceIn_]-upperIn_)<1.0e-7);
3997          solution_[sequenceIn_]=upperIn_;
3998        } else {
3999          assert (fabs(solution_[sequenceIn_]-lowerIn_)<1.0e-7);
4000          solution_[sequenceIn_]=lowerIn_;
4001        }
4002      } else {
4003        if (directionOut_<0) {
4004          assert (fabs(solution_[sequenceOut_]-upperOut_)<1.0e-7);
4005          solution_[sequenceOut_]=upperOut_;
4006        } else {
4007          assert (fabs(solution_[sequenceOut_]-lowerOut_)<1.0e-7);
4008          solution_[sequenceOut_]=lowerOut_;
4009        }
4010        solution_[sequenceIn_]=valueIn_+movement;
4011      }
4012    }   
4013    double objectiveChange = dualIn_*movement;
4014    // update duals
4015    if (pivotRow_>=0) {
4016      alpha_ = rowArray_[1]->denseVector()[pivotRow_];
4017      assert (fabs(alpha_)>1.0e-8);
4018      double multiplier = dualIn_/alpha_;
4019      rowArray_[0]->insert(pivotRow_,multiplier);
4020      factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
4021      // put row of tableau in rowArray[0] and columnArray[0]
4022      matrix_->transposeTimes(this,-1.0,
4023                              rowArray_[0],columnArray_[1],columnArray_[0]);
4024      // update column djs
4025      int i;
4026      int * index = columnArray_[0]->getIndices();
4027      int number = columnArray_[0]->getNumElements();
4028      double * element = columnArray_[0]->denseVector();
4029      for (i=0;i<number;i++) {
4030        int ii = index[i];
4031        dj_[ii] += element[ii];
4032        element[ii]=0.0;
4033      }
4034      columnArray_[0]->setNumElements(0);
4035      // and row djs
4036      index = rowArray_[0]->getIndices();
4037      number = rowArray_[0]->getNumElements();
4038      element = rowArray_[0]->denseVector();
4039      for (i=0;i<number;i++) {
4040        int ii = index[i];
4041        dj_[ii+numberColumns_] += element[ii];
4042        dual_[ii] = dj_[ii+numberColumns_];
4043        element[ii]=0.0;
4044      }
4045      rowArray_[0]->setNumElements(0);
4046      // check incoming
4047      assert (fabs(dj_[sequenceIn_])<1.0e-6);
4048    }
4049   
4050    // if stable replace in basis
4051    int updateStatus = factorization_->replaceColumn(rowArray_[2],
4052                                                   pivotRow_,
4053                                                     alpha_);
4054    bool takePivot=true;
4055    // if no pivots, bad update but reasonable alpha - take and invert
4056    if (updateStatus==2&&
4057        lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
4058      updateStatus=4;
4059    if (updateStatus==1||updateStatus==4) {
4060      // slight error
4061      if (factorization_->pivots()>5||updateStatus==4) {
4062        returnCode=-1;
4063      }
4064    } else if (updateStatus==2) {
4065      // major error
4066      rowArray_[1]->clear();
4067      takePivot=false;
4068      if (factorization_->pivots()) {
4069        // refactorize here
4070        statusOfProblem();
4071        roundAgain=true;
4072      } else {
4073        returnCode=1;
4074      }
4075    } else if (updateStatus==3) {
4076      // out of memory
4077      // increase space if not many iterations
4078      if (factorization_->pivots()<
4079          0.5*factorization_->maximumPivots()&&
4080          factorization_->pivots()<200)
4081        factorization_->areaFactor(
4082                                   factorization_->areaFactor() * 1.1);
4083      returnCode =-1; // factorize now
4084    }
4085    if (takePivot) {
4086      int save = algorithm_;
4087      // make simple so always primal
4088      algorithm_=1;
4089      housekeeping(objectiveChange);
4090      algorithm_=save;
4091    }
4092  }
4093  if (returnCode == -1) {
4094    // refactorize here
4095    statusOfProblem();
4096  } else {
4097    // just for now - recompute anyway
4098    gutsOfSolution(NULL,NULL);
4099  }
4100  return returnCode;
4101}
4102
4103/* Pivot in a variable and choose an outgoing one.  Assumes primal
4104   feasible - will not go through a bound.  Returns step length in theta
4105   Returns ray in ray_ (or NULL if no pivot)
4106   Return codes as before but -1 means no acceptable pivot
4107*/
4108int ClpInterior::primalPivotResult()
4109{
4110  assert (sequenceIn_>=0);
4111  valueIn_=solution_[sequenceIn_];
4112  lowerIn_=lower_[sequenceIn_];
4113  upperIn_=upper_[sequenceIn_];
4114  dualIn_=dj_[sequenceIn_];
4115
4116  int returnCode = ((ClpInteriorPrimal *) this)->pivotResult();
4117  if (returnCode<0&&returnCode>-4) {
4118    return 0;
4119  } else {
4120    printf("Return code of %d from ClpInteriorPrimal::pivotResult\n",
4121           returnCode);
4122    return -1;
4123  }
4124}
4125 
4126/* Pivot out a variable and choose an incoing one.  Assumes dual
4127   feasible - will not go through a reduced cost. 
4128   Returns step length in theta
4129   Returns ray in ray_ (or NULL if no pivot)
4130   Return codes as before but -1 means no acceptable pivot
4131*/
4132int 
4133ClpInterior::dualPivotResult()
4134{
4135  return ((ClpInteriorDual *) this)->pivotResult();
4136}
4137// Factorization frequency
4138int 
4139ClpInterior::factorizationFrequency() const
4140{
4141  if (factorization_)
4142    return factorization_->maximumPivots();
4143  else 
4144    return -1;
4145}
4146void 
4147ClpInterior::setFactorizationFrequency(int value)
4148{
4149  if (factorization_)
4150    factorization_->maximumPivots(value);
4151}
4152// Common bits of coding for dual and primal
4153int 
4154ClpInterior::startup(int ifValuesPass)
4155{
4156  // sanity check
4157  assert (numberRows_==matrix_->getNumRows());
4158  assert (numberColumns_==matrix_->getNumCols());
4159  // for moment all arrays must exist
4160  assert(columnLower_);
4161  assert(columnUpper_);
4162  assert(rowLower_);
4163  assert(rowUpper_);
4164  pivotRow_=-1;
4165  sequenceIn_=-1;
4166  sequenceOut_=-1;
4167  secondaryStatus_=0;
4168
4169  primalTolerance_=dblParam_[ClpPrimalTolerance];
4170  dualTolerance_=dblParam_[ClpDualTolerance];
4171  if (problemStatus_!=10)
4172    numberIterations_=0;
4173
4174  // put in standard form (and make row copy)
4175  // create modifiable copies of model rim and do optional scaling
4176  bool goodMatrix=createRim(7+8+16,true);
4177
4178  if (goodMatrix) {
4179    // Model looks okay
4180    // Do initial factorization
4181    // and set certain stuff
4182    // We can either set increasing rows so ...IsBasic gives pivot row
4183    // or we can just increment iBasic one by one
4184    // for now let ...iBasic give pivot row
4185    factorization_->increasingRows(2);
4186    // row activities have negative sign
4187    factorization_->slackValue(-1.0);
4188    factorization_->zeroTolerance(1.0e-13);
4189    // Switch off dense (unless special option set)
4190    int saveThreshold = factorization_->denseThreshold();
4191    factorization_->setDenseThreshold(0);
4192    // If values pass then perturb (otherwise may be optimal so leave a bit)
4193    if (ifValuesPass) {
4194      // do perturbation if asked for
4195     
4196      if (perturbation_<100) {
4197        if (algorithm_>0) {
4198          ((ClpInteriorPrimal *) this)->perturb(0);
4199        } else if (algorithm_<0) {
4200        ((ClpInteriorDual *) this)->perturb();
4201        }
4202      }
4203    }
4204    // for primal we will change bounds using infeasibilityCost_
4205    if (nonLinearCost_==NULL&&algorithm_>0) {
4206      // get a valid nonlinear cost function
4207      delete nonLinearCost_;
4208      nonLinearCost_= new ClpNonLinearCost(this);
4209    }
4210   
4211    // loop round to clean up solution if values pass
4212    int numberThrownOut = -1;
4213    int totalNumberThrownOut=0;
4214    while(numberThrownOut) {
4215      int status = internalFactorize(0+10*ifValuesPass);
4216      if (status<0)
4217        return 1; // some error
4218      else
4219        numberThrownOut = status;
4220     
4221      // for this we need clean basis so it is after factorize
4222      if (!numberThrownOut)
4223        numberThrownOut = gutsOfSolution(  NULL,NULL,
4224                                         ifValuesPass!=0);
4225      totalNumberThrownOut+= numberThrownOut;
4226     
4227    }
4228   
4229    if (totalNumberThrownOut)
4230      handler_->message(CLP_SINGULARITIES,messages_)
4231        <<totalNumberThrownOut
4232        <<CoinMessageEol;
4233    // Switch back dense
4234    factorization_->setDenseThreshold(saveThreshold);
4235   
4236    problemStatus_ = -1;
4237   
4238    // number of times we have declared optimality
4239    numberTimesOptimal_=0;
4240
4241    return 0;
4242  } else {
4243    // bad matrix
4244    return 2;
4245  }
4246   
4247}
4248
4249
4250void 
4251ClpInterior::finish()
4252{
4253  // Get rid of some arrays and empty factorization
4254  deleteRim();
4255  // Skip message if changing algorithms
4256  if (problemStatus_!=10) {
4257    assert(problemStatus_>=0&&problemStatus_<5);
4258    handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
4259      <<objectiveValue()
4260      <<CoinMessageEol;
4261  }
4262  factorization_->relaxAccuracyCheck(1.0);
4263  // get rid of any network stuff - could do more
4264  factorization_->cleanUp();
4265}
4266// Save data
4267ClpDataSave
4268ClpInterior::saveData() 
4269{
4270  ClpDataSave saved;
4271  saved.dualBound_ = dualBound_;
4272  saved.infeasibilityCost_ = infeasibilityCost_;
4273  saved.sparseThreshold_ = factorization_->sparseThreshold();
4274  saved.perturbation_ = perturbation_;
4275  // Progress indicator
4276  delete progress_;
4277  progress_ = new ClpInteriorProgress (this);
4278  return saved;
4279}
4280// Restore data
4281void 
4282ClpInterior::restoreData(ClpDataSave saved)
4283{
4284  factorization_->sparseThreshold(saved.sparseThreshold_);
4285  perturbation_ = saved.perturbation_;
4286  infeasibilityCost_ = saved.infeasibilityCost_;
4287  dualBound_ = saved.dualBound_;
4288  delete progress_;
4289  progress_=NULL;
4290}
4291/* Factorizes and returns true if optimal.  Used by user */
4292bool
4293ClpInterior::statusOfProblem()
4294{
4295  // is factorization okay?
4296  assert (internalFactorize(1)==0);
4297  // put back original costs and then check
4298  // also move to work arrays
4299  createRim(4+32);
4300  //memcpy(rowActivityWork_,rowActivity_,numberRows_*sizeof(double));
4301  //memcpy(columnActivityWork_,columnActivity_,numberColumns_*sizeof(double));
4302  gutsOfSolution(NULL,NULL);
4303  //memcpy(rowActivity_,rowActivityWork_,numberRows_*sizeof(double));
4304  //memcpy(columnActivity_,columnActivityWork_,numberColumns_*sizeof(double));
4305  //memcpy(reducedCost_,dj_,numberColumns_*sizeof(double));
4306  deleteRim(-1);
4307  return (primalFeasible()&&dualFeasible());
4308}
4309/* Return model - updates any scalars */
4310void 
4311ClpInterior::returnModel(ClpInterior & otherModel)
4312{
4313  ClpModel::returnModel(otherModel);
4314  otherModel.columnPrimalInfeasibility_ = columnPrimalInfeasibility_;
4315  otherModel.columnPrimalSequence_ = columnPrimalSequence_;
4316  otherModel.rowPrimalInfeasibility_ = rowPrimalInfeasibility_;
4317  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
4318  otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
4319  otherModel.columnDualSequence_ = columnDualSequence_;
4320  otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
4321  otherModel.rowDualSequence_ = rowDualSequence_;
4322  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
4323  otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
4324  otherModel.largestPrimalError_ = largestPrimalError_;
4325  otherModel.largestDualError_ = largestDualError_;
4326  otherModel.largestSolutionError_ = largestSolutionError_;
4327  otherModel.alpha_ = alpha_;
4328  otherModel.theta_ = theta_;
4329  otherModel.lowerIn_ = lowerIn_;
4330  otherModel.valueIn_ = valueIn_;
4331  otherModel.upperIn_ = upperIn_;
4332  otherModel.dualIn_ = dualIn_;
4333  otherModel.sequenceIn_ = sequenceIn_;
4334  otherModel.directionIn_ = directionIn_;
4335  otherModel.lowerOut_ = lowerOut_;
4336  otherModel.valueOut_ = valueOut_;
4337  otherModel.upperOut_ = upperOut_;
4338  otherModel.dualOut_ = dualOut_;
4339  otherModel.sequenceOut_ = sequenceOut_;
4340  otherModel.directionOut_ = directionOut_;
4341  otherModel.pivotRow_ = pivotRow_;
4342  otherModel.sumDualInfeasibilities_ = sumDualInfeasibilities_;
4343  otherModel.numberDualInfeasibilities_ = numberDualInfeasibilities_;
4344  otherModel.numberDualInfeasibilitiesWithoutFree_ = 
4345    numberDualInfeasibilitiesWithoutFree_;
4346  otherModel.sumPrimalInfeasibilities_ = sumPrimalInfeasibilities_;
4347  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
4348  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
4349  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
4350  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
4351}
4352/* Constructs a non linear cost from list of non-linearities (columns only)
4353   First lower of each column is taken as real lower
4354   Last lower is taken as real upper and cost ignored
4355   
4356   Returns nonzero if bad data e.g. lowers not monotonic
4357*/
4358int 
4359ClpInterior::createPiecewiseLinearCosts(const int * starts,
4360                                       const double * lower, const double * gradient)
4361{
4362  delete nonLinearCost_;
4363  // Set up feasible bounds and check monotonicity
4364  int iColumn;
4365  int returnCode=0;
4366
4367  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
4368    int iIndex = starts[iColumn];
4369    int end = starts[iColumn+1]-1;
4370    columnLower_[iColumn] = lower[iIndex];
4371    columnUpper_[iColumn] = lower[end];
4372    double value = columnLower_[iColumn];
4373    iIndex++;
4374    for (;iIndex<end;iIndex++) {
4375      if (lower[iIndex]<value)
4376        returnCode++; // not monotonic
4377      value=lower[iIndex];
4378    }
4379  }
4380  nonLinearCost_ = new ClpNonLinearCost(this,starts,lower,gradient);
4381  specialOptions_ |= 2; // say keep
4382  return returnCode;
4383}
4384/* For advanced use.  When doing iterative solves things can get
4385   nasty so on values pass if incoming solution has largest
4386   infeasibility < incomingInfeasibility throw out variables
4387   from basis until largest infeasibility < allowedInfeasibility
4388   or incoming largest infeasibility.
4389   If allowedInfeasibility>= incomingInfeasibility this is
4390   always possible altough you may end up with an all slack basis.
4391   
4392   Defaults are 1.0,10.0
4393*/
4394void 
4395ClpInterior::setValuesPassAction(float incomingInfeasibility,
4396                                float allowedInfeasibility)
4397{
4398  incomingInfeasibility_=incomingInfeasibility;
4399  allowedInfeasibility_=allowedInfeasibility;
4400  assert(incomingInfeasibility_>=0.0);
4401  assert(allowedInfeasibility_>=incomingInfeasibility_);
4402}
4403//#############################################################################
4404
4405ClpInteriorProgress::ClpInteriorProgress () 
4406{
4407  int i;
4408  for (i=0;i<CLP_PROGRESS;i++) {
4409    objective_[i] = COIN_DBL_MAX;
4410    infeasibility_[i] = -1.0; // set to an impossible value
4411    numberInfeasibilities_[i]=-1; 
4412    iterationNumber_[i]=-1;
4413  }
4414  for (i=0;i<CLP_CYCLE;i++) {
4415    in_[i]=-1;
4416    out_[i]=-1;
4417    way_[i]=0;
4418  }
4419  numberTimes_ = 0;
4420  numberBadTimes_ = 0;
4421  model_ = NULL;
4422}
4423
4424
4425//-----------------------------------------------------------------------------
4426
4427ClpInteriorProgress::~ClpInteriorProgress ()
4428{
4429}
4430// Copy constructor.
4431ClpInteriorProgress::ClpInteriorProgress(const ClpInteriorProgress &rhs) 
4432{
4433  int i;
4434  for (i=0;i<CLP_PROGRESS;i++) {
4435    objective_[i] = rhs.objective_[i];
4436    infeasibility_[i] = rhs.infeasibility_[i];
4437    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
4438    iterationNumber_[i]=rhs.iterationNumber_[i];
4439  }
4440  for (i=0;i<CLP_CYCLE;i++) {
4441    in_[i]=rhs.in_[i];
4442    out_[i]=rhs.out_[i];
4443    way_[i]=rhs.way_[i];
4444  }
4445  numberTimes_ = rhs.numberTimes_;
4446  numberBadTimes_ = rhs.numberBadTimes_;
4447  model_ = rhs.model_;
4448}
4449// Copy constructor.from model
4450ClpInteriorProgress::ClpInteriorProgress(ClpInterior * model) 
4451{
4452  model_ = model;
4453  int i;
4454  for (i=0;i<CLP_PROGRESS;i++) {
4455    if (model_->algorithm()>=0)
4456      objective_[i] = COIN_DBL_MAX;
4457    else
4458      objective_[i] = -COIN_DBL_MAX;
4459    infeasibility_[i] = -1.0; // set to an impossible value
4460    numberInfeasibilities_[i]=-1; 
4461    iterationNumber_[i]=-1;
4462  }
4463  for (i=0;i<CLP_CYCLE;i++) {
4464    in_[i]=-1;
4465    out_[i]=-1;
4466    way_[i]=0;
4467  }
4468  numberTimes_ = 0;
4469  numberBadTimes_ = 0;
4470}
4471// Assignment operator. This copies the data
4472ClpInteriorProgress & 
4473ClpInteriorProgress::operator=(const ClpInteriorProgress & rhs)
4474{
4475  if (this != &rhs) {
4476    int i;
4477    for (i=0;i<CLP_PROGRESS;i++) {
4478      objective_[i] = rhs.objective_[i];
4479      infeasibility_[i] = rhs.infeasibility_[i];
4480      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
4481      iterationNumber_[i]=rhs.iterationNumber_[i];
4482    }
4483    for (i=0;i<CLP_CYCLE;i++) {
4484      in_[i]=rhs.in_[i];
4485      out_[i]=rhs.out_[i];
4486      way_[i]=rhs.way_[i];
4487    }
4488    numberTimes_ = rhs.numberTimes_;
4489    numberBadTimes_ = rhs.numberBadTimes_;
4490    model_ = rhs.model_;
4491  }
4492  return *this;
4493}
4494// Seems to be something odd about exact comparison of doubles on linux
4495static bool equalDouble(double value1, double value2)
4496{
4497  unsigned int *i1 = (unsigned int *) &value1;
4498  unsigned int *i2 = (unsigned int *) &value2;
4499  if (sizeof(unsigned int)*2==sizeof(double)) 
4500    return (i1[0]==i2[0]&&i1[1]==i2[1]);
4501  else
4502    return (i1[0]==i2[0]);
4503}
4504int
4505ClpInteriorProgress::looping()
4506{
4507  if (!model_)
4508    return -1;
4509  double objective = model_->rawObjectiveValue();
4510  double infeasibility;
4511  int numberInfeasibilities;
4512  int iterationNumber = model_->numberIterations();
4513  if (model_->algorithm()<0) {
4514    // dual
4515    infeasibility = model_->sumPrimalInfeasibilities();
4516    numberInfeasibilities = model_->numberPrimalInfeasibilities();
4517  } else {
4518    //primal
4519    infeasibility = model_->sumDualInfeasibilities();
4520    numberInfeasibilities = model_->numberDualInfeasibilities();
4521  }
4522  int i;
4523  int numberMatched=0;
4524  int matched=0;
4525  int nsame=0;
4526  for (i=0;i<CLP_PROGRESS;i++) {
4527    bool matchedOnObjective = equalDouble(objective,objective_[i]);
4528    bool matchedOnInfeasibility = equalDouble(infeasibility,infeasibility_[i]);
4529    bool matchedOnInfeasibilities = 
4530      (numberInfeasibilities==numberInfeasibilities_[i]);
4531   
4532    if (matchedOnObjective&&matchedOnInfeasibility&&matchedOnInfeasibilities) {
4533      matched |= (1<<i);
4534      // Check not same iteration
4535      if (iterationNumber!=iterationNumber_[i]) {
4536        numberMatched++;
4537        // here mainly to get over compiler bug?
4538        if (model_->messageHandler()->logLevel()>10)
4539          printf("%d %d %d %d %d loop check\n",i,numberMatched,
4540                 matchedOnObjective, matchedOnInfeasibility, 
4541                 matchedOnInfeasibilities);
4542      } else {
4543        // stuck but code should notice
4544        nsame++;
4545      }
4546    }
4547    if (i) {
4548      objective_[i-1] = objective_[i];
4549      infeasibility_[i-1] = infeasibility_[i];
4550      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
4551      iterationNumber_[i-1]=iterationNumber_[i];
4552    }
4553  }
4554  objective_[CLP_PROGRESS-1] = objective;
4555  infeasibility_[CLP_PROGRESS-1] = infeasibility;
4556  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
4557  iterationNumber_[CLP_PROGRESS-1]=iterationNumber;
4558  if (nsame==CLP_PROGRESS)
4559    numberMatched=CLP_PROGRESS; // really stuck
4560  if (model_->progressFlag())
4561    numberMatched=0;
4562  numberTimes_++;
4563  if (numberTimes_<10)
4564    numberMatched=0;
4565  // skip if just last time as may be checking something
4566  if (matched==(1<<(CLP_PROGRESS-1)))
4567    numberMatched=0;
4568  if (numberMatched) {
4569    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
4570      <<numberMatched
4571      <<matched
4572      <<numberTimes_
4573      <<CoinMessageEol;
4574    numberBadTimes_++;
4575    if (numberBadTimes_<10) {
4576      // make factorize every iteration
4577      model_->forceFactorization(1);
4578      if (model_->algorithm()<0) {
4579        // dual - change tolerance
4580        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
4581        // if infeasible increase dual bound
4582        if (model_->dualBound()<1.0e19) {
4583          model_->setDualBound(model_->dualBound()*1.1);
4584        }
4585      } else {
4586        // primal - change tolerance    model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
4587        // if infeasible increase infeasibility cost
4588        if (model_->nonLinearCost()->numberInfeasibilities()&&
4589            model_->infeasibilityCost()<1.0e19) {
4590          model_->setInfeasibilityCost(model_->infeasibilityCost()*1.1);
4591        }
4592      }
4593      return -2;
4594    } else {
4595      // look at solution and maybe declare victory
4596      if (infeasibility<1.0e-4) {
4597        return 0;
4598      } else {
4599        model_->messageHandler()->message(CLP_LOOP,model_->messages())
4600          <<CoinMessageEol;
4601#ifndef NDEBUG
4602        abort();
4603#endif
4604        return 3;
4605      }
4606    }
4607  }
4608  return -1;
4609}
4610// Returns previous objective (if -1) - current if (0)
4611double 
4612ClpInteriorProgress::lastObjective(int back) const
4613{
4614  return objective_[CLP_PROGRESS-1-back];
4615}
4616// Modify objective e.g. if dual infeasible in dual
4617void 
4618ClpInteriorProgress::modifyObjective(double value)
4619{
4620  objective_[CLP_PROGRESS-1]=value;
4621}
4622// Returns previous iteration number (if -1) - current if (0)
4623int 
4624ClpInteriorProgress::lastIterationNumber(int back) const
4625{
4626  return iterationNumber_[CLP_PROGRESS-1-back];
4627}
4628// Start check at beginning of whileIterating
4629void 
4630ClpInteriorProgress::startCheck()
4631{
4632  int i;
4633  for (i=0;i<CLP_CYCLE;i++) {
4634    in_[i]=-1;
4635    out_[i]=-1;
4636    way_[i]=0;
4637  }
4638}
4639// Returns cycle length in whileIterating
4640int 
4641ClpInteriorProgress::cycle(int in, int out,int wayIn,int wayOut)
4642{
4643  int i;
4644  int matched=0;
4645  return 0;
4646  // first see if in matches any out
4647  for (i=1;i<CLP_CYCLE;i++) {
4648    if (in==out_[i]) {
4649      // even if flip then suspicious
4650      matched=-1;
4651      break;
4652    }
4653  }
4654  if (!matched) {
4655    // can't be cycle
4656    for (i=0;i<CLP_CYCLE-1;i++) {
4657      in_[i]=in_[i+1];
4658      out_[i]=out_[i+1];
4659      way_[i]=way_[i+1];
4660    }
4661  } else {
4662    // possible cycle
4663    matched=0;
4664    for (i=0;i<CLP_CYCLE-1;i++) {
4665      int k;
4666      char wayThis = way_[i];
4667      int inThis = in_[i];
4668      int outThis = out_[i];
4669      for(k=i+1;k<CLP_CYCLE;k++) {
4670        if (inThis==in_[k]&&outThis==out_[k]&&wayThis==way_[k]) {
4671          matched=k-i;
4672        }
4673      }
4674      in_[i]=in_[i+1];
4675      out_[i]=out_[i+1];
4676      way_[i]=way_[i+1];
4677    }
4678  }
4679  char way = 1-wayIn+4*(1-wayOut);
4680  in_[CLP_CYCLE-1]=in;
4681  out_[CLP_CYCLE-1]=out;
4682  way_[CLP_CYCLE-1]=way;
4683  return matched;
4684}
4685// Allow initial dense factorization
4686void 
4687ClpInterior::setInitialDenseFactorization(bool onOff)
4688{
4689  if (onOff)
4690    specialOptions_ |= 8;
4691  else
4692    specialOptions_ &= ~8;
4693}
4694bool 
4695ClpInterior::initialDenseFactorization() const
4696{
4697  return (specialOptions_&8)!=0;
4698}
Note: See TracBrowser for help on using the repository browser.