source: trunk/ClpSimplex.cpp @ 115

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

OsiSimplexInterface?

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