source: branches/devel-1/ClpSimplex.cpp @ 14

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

Improving reliability using IBM Burlington problems (Primal)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 59.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6#if defined(_MSC_VER)
7// Turn off compiler warning about long names
8#  pragma warning(disable:4786)
9#endif
10
11#include <math.h>
12
13#include "CoinHelperFunctions.hpp"
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16#include "OsiPackedMatrix.hpp"
17#include "OsiIndexedVector.hpp"
18#include "OsiWarmStartBasis.hpp"
19#include "ClpDualRowSteepest.hpp"
20#include "ClpPrimalColumnDantzig.hpp"
21#include "ClpPrimalColumnSteepest.hpp"
22#include "ClpNonLinearCost.hpp"
23#include "ClpMessage.hpp"
24#include <cfloat>
25
26#include <string>
27#include <stdio.h>
28#include <iostream>
29// This returns a non const array filled with input from scalar
30// or actual array
31template <class T> inline T*
32copyOfArray( const T * array, const int size, T value)
33{
34  T * arrayNew = new T[size];
35  if (array) 
36    CoinDisjointCopyN(array,size,arrayNew);
37  else
38    CoinFillN ( arrayNew, size,value);
39  return arrayNew;
40}
41
42// This returns a non const array filled with actual array (or NULL)
43template <class T> inline T*
44copyOfArray( const T * array, const int size)
45{
46  if (array) {
47    T * arrayNew = new T[size];
48    CoinDisjointCopyN(array,size,arrayNew);
49    return arrayNew;
50  } else {
51    return NULL;
52  }
53}
54
55//#############################################################################
56
57ClpSimplex::ClpSimplex () :
58
59  ClpModel(),
60  columnPrimalInfeasibility_(0.0),
61  columnPrimalSequence_(-2),
62  rowPrimalInfeasibility_(0.0),
63  rowPrimalSequence_(-2), 
64  columnDualInfeasibility_(0.0),
65  columnDualSequence_(-2),
66  rowDualInfeasibility_(0.0),
67  rowDualSequence_(-2),
68  primalToleranceToGetOptimal_(-1.0),
69  remainingDualInfeasibility_(0.0),
70  largeValue_(1.0e15),
71  largestPrimalError_(0.0),
72  largestDualError_(0.0),
73  largestSolutionError_(0.0),
74  dualBound_(1.0e7),
75  lower_(NULL),
76  rowLowerWork_(NULL),
77  columnLowerWork_(NULL),
78  upper_(NULL),
79  rowUpperWork_(NULL),
80  columnUpperWork_(NULL),
81  cost_(NULL),
82  rowObjectiveWork_(NULL),
83  objectiveWork_(NULL),
84  alpha_(0.0),
85  theta_(0.0),
86  lowerIn_(0.0),
87  valueIn_(0.0),
88  upperIn_(0.0),
89  dualIn_(0.0),
90  sequenceIn_(-1),
91  directionIn_(-1),
92  lowerOut_(-1),
93  valueOut_(-1),
94  upperOut_(-1),
95  dualOut_(-1),
96  sequenceOut_(-1),
97  directionOut_(-1),
98  pivotRow_(-1),
99  status_(NULL),
100  dj_(NULL),
101  rowReducedCost_(NULL),
102  reducedCostWork_(NULL),
103  solution_(NULL),
104  rowActivityWork_(NULL),
105  columnActivityWork_(NULL),
106  dualTolerance_(0.0),
107  primalTolerance_(0.0),
108  sumDualInfeasibilities_(0.0),
109  numberDualInfeasibilities_(0),
110  numberDualInfeasibilitiesWithoutFree_(0),
111  sumPrimalInfeasibilities_(0.0),
112  numberPrimalInfeasibilities_(0),
113  pivotVariable_(NULL),
114  factorization_(NULL),
115  numberRefinements_(0),
116  rowScale_(NULL),
117  savedSolution_(NULL),
118  columnScale_(NULL),
119  scalingFlag_(0),
120  numberTimesOptimal_(0),
121  changeMade_(1),
122  algorithm_(0),
123  forceFactorization_(-1),
124  perturbation_(100),
125  infeasibilityCost_(1.0e7),
126  nonLinearCost_(NULL),
127  specialOptions_(0),
128  lastBadIteration_(-999999)
129
130{
131  int i;
132  for (i=0;i<6;i++) {
133    rowArray_[i]=NULL;
134    columnArray_[i]=NULL;
135  }
136  saveStatus_=NULL;
137  // get an empty factorization so we can set tolerances etc
138  factorization_ = new ClpFactorization();
139  // Say sparse
140  factorization_->sparseThreshold(1);
141  // say Steepest pricing
142  dualRowPivot_ = new ClpDualRowSteepest();
143  // say Steepest pricing
144  primalColumnPivot_ = new ClpPrimalColumnSteepest();
145 
146}
147
148
149//-----------------------------------------------------------------------------
150
151ClpSimplex::~ClpSimplex ()
152{
153  gutsOfDelete(0);
154}
155//#############################################################################
156void ClpSimplex::setLargeValue( double value) 
157{
158  if (value>0.0&&value<DBL_MAX)
159    largeValue_=value;
160}
161int
162ClpSimplex::gutsOfSolution ( const double * rowActivities,
163                             const double * columnActivities,
164                             bool valuesPass)
165{
166
167
168  // if values pass, save values of basic variables
169  double * save = NULL;
170  double oldValue=0.0;
171  if (valuesPass) {
172    assert(algorithm_>0); // only primal at present
173    assert(nonLinearCost_);
174    int iRow;
175    checkPrimalSolution( rowActivities, columnActivities);
176    // get correct bounds on all variables
177    nonLinearCost_->checkInfeasibilities();
178    oldValue = nonLinearCost_->largestInfeasibility();
179    save = new double[numberRows_];
180    for (iRow=0;iRow<numberRows_;iRow++) {
181      int iPivot=pivotVariable_[iRow];
182      save[iRow] = solution_[iPivot];
183    }
184  }
185  // do work
186  computePrimals(rowActivities, columnActivities);
187  double objectiveModification = 0.0;
188  if (algorithm_>0&&nonLinearCost_!=NULL) {
189    // primal algorithm
190    // get correct bounds on all variables
191    nonLinearCost_->checkInfeasibilities();
192    objectiveModification += nonLinearCost_->changeInCost();
193    if (nonLinearCost_->numberInfeasibilities())
194      handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
195        <<nonLinearCost_->changeInCost()
196        <<nonLinearCost_->numberInfeasibilities()
197        <<OsiMessageEol;
198  }
199  if (valuesPass) {
200#ifdef CLP_DEBUG
201    std::cout<<"Largest given infeasibility "<<oldValue
202             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
203#endif
204    int numberOut=0;
205    if (oldValue<1.0
206        &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
207        largestPrimalError_>1.0e-3) {
208      // throw out up to 1000 structurals
209      int iRow;
210      int * sort = new int[numberRows_];
211      // first put back solution and store difference
212      for (iRow=0;iRow<numberRows_;iRow++) {
213        int iPivot=pivotVariable_[iRow];
214        double difference = fabs(solution_[iPivot]=save[iRow]);
215        solution_[iPivot]=save[iRow];
216        save[iRow]=difference;
217      }
218      for (iRow=0;iRow<numberRows_;iRow++) {
219        int iPivot=pivotVariable_[iRow];
220       
221        if (iPivot<numberColumns_) {
222          // column
223          double difference= save[iRow];
224          if (difference>1.0e-4) {
225            sort[numberOut]=iPivot;
226            save[numberOut++]=difference;
227          }
228        }
229      }
230      CoinSort_2(save, save + numberOut, sort,
231                 CoinFirstGreater_2<double, int>());
232      numberOut = min(1000,numberOut);
233      for (iRow=0;iRow<numberOut;iRow++) {
234        int iColumn=sort[iRow];
235        setColumnStatus(iColumn,ClpSimplex::superBasic);
236       
237      }
238      delete [] sort;
239    }
240    delete [] save;
241    if (numberOut)
242      return numberOut;
243  }
244
245  computeDuals();
246
247  // now check solutions
248  checkPrimalSolution( rowActivities, columnActivities);
249  objectiveValue_ += objectiveModification;
250  checkDualSolution();
251  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
252                               largestDualError_>1.0e-2)) 
253    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
254      <<largestPrimalError_
255      <<largestDualError_
256      <<OsiMessageEol;
257  return 0;
258}
259void
260ClpSimplex::computePrimals ( const double * rowActivities,
261                                     const double * columnActivities)
262{
263
264  //work space
265  OsiIndexedVector  * workSpace = rowArray_[0];
266
267  double * array = new double [numberRows_];
268  double * save = new double [numberRows_];
269  double * previous = new double [numberRows_];
270
271  // accumulate non basic stuff
272  CoinFillN(array,numberRows_,0.0);
273
274  int iRow;
275  // order is this way for scaling
276  // Use whole matrix every time to make it easier for ClpMatrixBase
277  // So zero out basic
278  if (columnActivities!=columnActivityWork_)
279    CoinDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
280  if (rowActivities!=rowActivityWork_)
281    CoinDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
282  for (iRow=0;iRow<numberRows_;iRow++) {
283    int iPivot=pivotVariable_[iRow];
284    solution_[iPivot] = 0.0;
285  }
286  times(-1.0,columnActivityWork_,array);
287
288  for (iRow=0;iRow<numberRows_;iRow++) {
289    array[iRow] += rowActivityWork_[iRow];
290  }
291
292  // Ftran adjusted RHS and iterate to improve accuracy
293  double lastError=DBL_MAX;
294  int iRefine;
295  double * work = workSpace->denseVector();
296  factorization_->updateColumn(workSpace,array);
297
298  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
299    // save in case we want to restore
300    CoinDisjointCopyN ( array, numberRows_ , save);
301
302    // put solution in correct place
303    for (iRow=0;iRow<numberRows_;iRow++) {
304      int iPivot=pivotVariable_[iRow];
305      solution_[iPivot] = array[iRow];
306    }
307    // check Ax == b  (for all)
308    times(-1.0,columnActivityWork_,work);
309    for (iRow=0;iRow<numberRows_;iRow++) {
310      work[iRow] += rowActivityWork_[iRow];
311    }
312
313    largestPrimalError_=0.0;
314    for (iRow=0;iRow<numberRows_;iRow++) {
315      if (fabs(work[iRow])>largestPrimalError_) {
316        largestPrimalError_=fabs(work[iRow]);
317      }
318      //work[iRow] -= save[iRow];
319    }
320    if (largestPrimalError_>=lastError) {
321      // restore
322      double * temp = array;
323      array = previous;
324      previous=temp;
325      break;
326    }
327    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
328      // try and make better
329      // save this
330      double * temp = array;
331      array = previous;
332      previous=temp;
333      double multiplier = 131072.0;
334      for (iRow=0;iRow<numberRows_;iRow++) {
335        array[iRow] = multiplier*work[iRow];
336        work[iRow]=0.0;
337      }
338      lastError=largestPrimalError_;
339      factorization_->updateColumn(workSpace,array);
340      multiplier = 1.0/multiplier;
341      for (iRow=0;iRow<numberRows_;iRow++) {
342        array[iRow] = previous[iRow] + multiplier*array[iRow];
343        work[iRow]=0.0;
344      }
345    }
346  }
347
348  // solution as accurate as we are going to get
349  CoinFillN(work,numberRows_,0.0);
350  // put solution in correct place
351  for (iRow=0;iRow<numberRows_;iRow++) {
352    int iPivot=pivotVariable_[iRow];
353    solution_[iPivot] = array[iRow];
354  }
355  delete [] previous;
356  delete [] array;
357  delete [] save;
358}
359// now dual side
360void
361ClpSimplex::computeDuals()
362{
363  double slackValue = factorization_->slackValue();
364  //work space
365  OsiIndexedVector  * workSpace = rowArray_[0];
366
367  double * array = new double [numberRows_];
368  double * save = new double [numberRows_];
369  double * previous = new double [numberRows_];
370
371  int iRow;
372#ifdef CLP_DEBUG
373  workSpace->checkClear();
374#endif
375  for (iRow=0;iRow<numberRows_;iRow++) {
376    int iPivot=pivotVariable_[iRow];
377    if (iPivot>=numberColumns_) {
378      // slack
379      array[iRow] = rowObjectiveWork_[iPivot-numberColumns_];
380    } else {
381      // column
382      array[iRow]=objectiveWork_[iPivot];
383    }
384  }
385  CoinDisjointCopyN ( array, numberRows_ , save);
386
387  // Btran basic costs and get as accurate as possible
388  double lastError=DBL_MAX;
389  int iRefine;
390  double * work = workSpace->denseVector();
391  factorization_->updateColumnTranspose(workSpace,array);
392  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
393    // check basic reduced costs zero
394    largestDualError_=0.0;
395    // would be faster to do just for basic but this reduces code
396    CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
397    transposeTimes(-1.0,array,reducedCostWork_);
398    for (iRow=0;iRow<numberRows_;iRow++) {
399      int iPivot=pivotVariable_[iRow];
400      if (iPivot>=numberColumns_) {
401        // slack
402        //work[iRow] += slackValue*array[iPivot-numberColumns_];
403        //work[iRow] += rowObjectiveWork_[iPivot-numberColumns_];
404        work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
405        - slackValue*array[iPivot-numberColumns_];
406      } else {
407        // column
408        work[iRow] = reducedCostWork_[iPivot];
409      }
410      if (fabs(work[iRow])>largestDualError_) {
411        largestDualError_=fabs(work[iRow]);
412      }
413    }
414    if (largestDualError_>=lastError) {
415      // restore
416      double * temp = array;
417      array = previous;
418      previous=temp;
419      break;
420    }
421    if (iRefine<numberRefinements_&&largestDualError_>1.0e-20) {
422      // try and make better
423      // save this
424      double * temp = array;
425      array = previous;
426      previous=temp;
427      double multiplier = 131072.0;
428      for (iRow=0;iRow<numberRows_;iRow++) {
429        array[iRow] = multiplier*work[iRow];
430        work[iRow]=0.0;
431      }
432      lastError=largestDualError_;
433      factorization_->updateColumnTranspose(workSpace,array);
434      multiplier = 1.0/multiplier;
435      for (iRow=0;iRow<numberRows_;iRow++) {
436        array[iRow] = previous[iRow] + multiplier*array[iRow];
437        work[iRow]=0.0;
438      }
439    }
440  }
441  CoinFillN(work,numberRows_,0.0);
442  // now look at dual solution
443  for (iRow=0;iRow<numberRows_;iRow++) {
444    // slack
445    double value = array[iRow];
446    dual_[iRow]=value;
447    value = rowObjectiveWork_[iRow] - value*slackValue;
448    rowReducedCost_[iRow]=value;
449  }
450  CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
451  transposeTimes(-1.0,dual_,reducedCostWork_);
452  delete [] previous;
453  delete [] array;
454  delete [] save;
455}
456/* Given an existing factorization computes and checks
457   primal and dual solutions.  Uses input arrays for variables at
458   bounds.  Returns feasibility states */
459int ClpSimplex::getSolution ( const double * rowActivities,
460                               const double * columnActivities)
461{
462  if (!factorization_->status()) {
463    // put in standard form
464    createRim(7+8+16);
465    // do work
466    gutsOfSolution ( rowActivities, columnActivities);
467    // release extra memory
468    deleteRim();
469  }
470  return factorization_->status();
471}
472/* Given an existing factorization computes and checks
473   primal and dual solutions.  Uses current problem arrays for
474   bounds.  Returns feasibility states */
475int ClpSimplex::getSolution ( )
476{
477  double * rowActivities = new double[numberRows_];
478  double * columnActivities = new double[numberColumns_];
479  CoinDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
480  CoinDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
481  int status = getSolution( rowActivities, columnActivities);
482  delete [] rowActivities;
483  delete [] columnActivities;
484  return status;
485}
486// Factorizes using current basis.  This is for external use
487// Return codes are as from ClpFactorization
488int ClpSimplex::factorize ()
489{
490  // put in standard form
491  createRim(7+8+16);
492  // do work
493  int status = internalFactorize(-1);
494  // release extra memory
495  deleteRim();
496
497  return status;
498}
499
500/* Factorizes using current basis. 
501      solveType - 1 iterating, 0 initial, -1 external
502      If 10 added then in primal values pass
503*/
504// Return codes are as from ClpFactorization
505int ClpSimplex::internalFactorize ( int solveType)
506{
507
508  int iRow,iColumn;
509
510  bool valuesPass=false;
511  if (solveType>=10) {
512    valuesPass=true;
513    solveType -= 10;
514  }
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)!=ClpSimplex::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 ClpSimplex::basic:
534            numberBasic++;
535            break;
536          case ClpSimplex::isFree:
537            assert(rowLowerWork_[iRow]<-largeValue_);
538            assert(rowUpperWork_[iRow]>largeValue_);
539            rowActivityWork_[iRow]=0.0;
540            break;
541          case ClpSimplex::atUpperBound:
542            rowActivityWork_[iRow]=rowUpperWork_[iRow];
543            if (rowActivityWork_[iRow]>largeValue_) {
544              assert(rowLowerWork_[iRow]>-largeValue_);
545              rowActivityWork_[iRow]=rowLowerWork_[iRow];
546              setRowStatus(iRow,ClpSimplex::atLowerBound);
547            }
548            break;
549          case ClpSimplex::atLowerBound:
550            rowActivityWork_[iRow]=rowLowerWork_[iRow];
551            if (rowActivityWork_[iRow]<-largeValue_) {
552              assert(rowUpperWork_[iRow]<largeValue_);
553              rowActivityWork_[iRow]=rowUpperWork_[iRow];
554              setRowStatus(iRow,ClpSimplex::atUpperBound);
555            }
556            break;
557          case ClpSimplex::superBasic:
558            abort();
559            break;
560          }
561        }
562        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
563          switch(getColumnStatus(iColumn)) {
564           
565          case ClpSimplex::basic:
566            if (numberBasic==numberRows_) {
567              // take out of basis
568              if (columnLowerWork_[iColumn]>-largeValue_) {
569                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
570                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
571                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
572                  setColumnStatus(iColumn,ClpSimplex::atLowerBound);
573                } else {
574                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
575                  setColumnStatus(iColumn,ClpSimplex::atUpperBound);
576                }
577              } else if (columnUpperWork_[iColumn]<largeValue_) {
578                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
579                setColumnStatus(iColumn,ClpSimplex::atUpperBound);
580              } else {
581                columnActivityWork_[iColumn]=0.0;
582                setColumnStatus(iColumn,ClpSimplex::isFree);
583              }
584            } else {
585              numberBasic++;
586            }
587            break;
588          case ClpSimplex::isFree:
589            columnActivityWork_[iColumn]=0.0;
590            break;
591          case ClpSimplex::atUpperBound:
592            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
593            if (columnActivityWork_[iColumn]>largeValue_) {
594              assert(columnLowerWork_[iColumn]>-largeValue_);
595              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
596              setColumnStatus(iColumn,ClpSimplex::atLowerBound);
597            }
598            break;
599          case ClpSimplex::atLowerBound:
600            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
601            if (columnActivityWork_[iColumn]<-largeValue_) {
602              assert(columnUpperWork_[iColumn]<largeValue_);
603              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
604              setColumnStatus(iColumn,ClpSimplex::atUpperBound);
605            }
606            break;
607          case ClpSimplex::superBasic:
608            abort();
609            break;
610          }
611        }
612      } else {
613        //#define TESTFREE
614        // all slack basis
615        int numberBasic=0;
616        // changed to put free variables in basis
617        if (!status_) {
618          status_ = new unsigned char [numberColumns_+numberRows_];
619          memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
620        }
621        for (iRow=0;iRow<numberRows_;iRow++) {
622          double lower=rowLowerWork_[iRow];
623          double upper=rowUpperWork_[iRow];
624          if (lower>-largeValue_||upper<largeValue_) {
625            if (fabs(lower)<=fabs(upper)) {
626              setRowStatus(iRow,ClpSimplex::atLowerBound);
627              rowActivityWork_[iRow]=lower;
628            } else {
629              setRowStatus(iRow,ClpSimplex::atUpperBound);
630              rowActivityWork_[iRow]=upper;
631            }
632          } else {
633            setRowStatus(iRow,ClpSimplex::isFree);
634            rowActivityWork_[iRow]=0.0;
635          }
636#ifdef TESTFREE
637          setRowStatus(iRow,ClpSimplex::basic);
638          numberBasic++;
639#endif
640        }
641        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
642          double lower=columnLowerWork_[iColumn];
643          double upper=columnUpperWork_[iColumn];
644          double big_bound = largeValue_;
645          if (lower>-big_bound||upper<big_bound) {
646            if (fabs(lower)<=fabs(upper)) {
647              setColumnStatus(iColumn,ClpSimplex::atLowerBound);
648              columnActivityWork_[iColumn]=lower;
649            } else {
650              setColumnStatus(iColumn,ClpSimplex::atUpperBound);
651              columnActivityWork_[iColumn]=upper;
652            }
653          } else {
654#ifndef TESTFREE
655            numberBasic++;
656            setColumnStatus(iColumn,ClpSimplex::basic);
657#else
658            setColumnStatus(iColumn,ClpSimplex::isFree);
659#endif
660            columnActivityWork_[iColumn]=0.0;
661          }
662        }
663        assert(numberBasic<=numberRows_); // problems if too many free
664        if (!numberBasic) {
665          // might as well do all slack basis
666          for (iRow=0;iRow<numberRows_;iRow++) {
667            setRowStatus(iRow,ClpSimplex::basic);
668          }
669        }
670      }
671    } else {
672      // values pass has less coding
673      // make row activities correct
674      memset(rowActivityWork_,0,numberRows_*sizeof(double));
675      times(1.0,columnActivityWork_,rowActivityWork_);
676      if (status_) {
677        // set values from warm start (if sensible)
678        int numberBasic=0;
679        for (iRow=0;iRow<numberRows_;iRow++) {
680          if (getRowStatus(iRow)==ClpSimplex::basic) 
681            numberBasic++;
682          else
683            setRowStatus(iRow,ClpSimplex::superBasic);
684        }
685        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
686          if (getColumnStatus(iColumn)==ClpSimplex::basic) {
687            if (numberBasic==numberRows_) 
688              // take out of basis
689              setColumnStatus(iColumn,ClpSimplex::superBasic);
690            else 
691              numberBasic++;
692          } else {
693            setColumnStatus(iColumn,ClpSimplex::superBasic);
694          }
695        }
696      } else {
697        // all slack basis
698        int numberBasic=0;
699        if (!status_) {
700          status_ = new unsigned char [numberColumns_+numberRows_];
701          memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
702        }
703        for (iRow=0;iRow<numberRows_;iRow++) {
704          setRowStatus(iRow,ClpSimplex::basic);
705          numberBasic++;
706        }
707        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
708          setColumnStatus(iColumn,ClpSimplex::superBasic);
709        }
710      }
711    }
712    numberRefinements_=1;
713  }
714  int status=-99;
715  int * rowIsBasic = new int[numberRows_];
716  int * columnIsBasic = new int[numberColumns_];
717  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
718  while (status<-98) {
719 
720    int i;
721    int numberBasic=0;
722    for (i=0;i<numberRows_;i++) {
723      if (getRowStatus(i) == ClpSimplex::basic) {
724        rowIsBasic[i]=1;
725        numberBasic++;
726      } else {
727        rowIsBasic[i]=-1;
728      }
729    }
730    for (i=0;i<numberColumns_;i++) {
731      if (getColumnStatus(i) == ClpSimplex::basic) {
732        columnIsBasic[i]=1;
733        numberBasic++;
734      } else {
735        columnIsBasic[i]=-1;
736      }
737    }
738    assert (numberBasic<=numberRows_);
739    while (status==-99) {
740      status =  factorization_->factorize(this,matrix_,
741                                          numberRows_,numberColumns_,
742                                          rowIsBasic, columnIsBasic,
743                                          0.0);
744      if (status==-99) {
745        // get more memory
746        factorization_->areaFactor(2.0*factorization_->areaFactor());
747      }
748    }
749    if (!status) {
750      // do pivot information
751      for (i=0;i<numberRows_;i++) {
752        if (getRowStatus(i) == ClpSimplex::basic) {
753          pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
754        }
755      }
756      for (i=0;i<numberColumns_;i++) {
757        if (getColumnStatus(i) == ClpSimplex::basic) {
758          pivotVariable_[columnIsBasic[i]]=i;
759        }
760      }
761    } else {
762      // leave pivotVariable_ in useful form for cleaning basis
763      for (i=0;i<numberRows_;i++) {
764        pivotVariable_[i]=-1;
765      }
766      for (i=0;i<numberRows_;i++) {
767        if (getRowStatus(i) == ClpSimplex::basic) {
768          int iPivot = rowIsBasic[i];
769          if (iPivot>=0) 
770            pivotVariable_[iPivot]=i+numberColumns_;
771        }
772      }
773      for (i=0;i<numberColumns_;i++) {
774        if (getColumnStatus(i) == ClpSimplex::basic) {
775          int iPivot = columnIsBasic[i];
776          if (iPivot>=0) 
777            pivotVariable_[iPivot]=i;
778        }
779      }
780    }
781    if (status==-1) {
782      if (!solveType) {
783        //redo basis - first take ALL columns out
784        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
785          if (getColumnStatus(iColumn)==
786              ClpSimplex::basic) {
787            // take out
788            if (!valuesPass) {
789              double lower=columnLowerWork_[iColumn];
790              double upper=columnUpperWork_[iColumn];
791              double value=columnActivityWork_[iColumn];
792              if (lower>-largeValue_||upper<largeValue_) {
793                if (fabs(value-lower)<fabs(value-upper)) {
794                  setColumnStatus(iColumn,ClpSimplex::atLowerBound);
795                  columnActivityWork_[iColumn]=lower;
796                } else {
797                  setColumnStatus(iColumn,ClpSimplex::atUpperBound);
798                  columnActivityWork_[iColumn]=upper;
799                }
800              } else {
801                setColumnStatus(iColumn,ClpSimplex::isFree);
802              }
803            } else {
804              setColumnStatus(iColumn,ClpSimplex::superBasic);
805            }
806          }
807        }
808        for (iRow=0;iRow<numberRows_;iRow++) {
809          int iSequence=pivotVariable_[iRow];
810          if (iSequence>=0) {
811            // basic
812            if (iSequence>=numberColumns_) {
813              // slack in - leave
814              assert (iSequence-numberColumns_==iRow);
815            } else {
816              // put back structural
817              setColumnStatus(iSequence,ClpSimplex::basic);
818            }
819          } else {
820            // put in slack
821            setRowStatus(iRow,ClpSimplex::basic);
822          }
823        }
824        // signal repeat
825        status=-99;
826      }
827    }
828  } 
829  delete [] rowIsBasic;
830  delete [] columnIsBasic;
831  if (status) {
832    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
833      <<status
834      <<OsiMessageEol;
835    return -1;
836  }
837
838  // sparse methods
839  if (factorization_->sparseThreshold()) {
840    // get default value
841    factorization_->sparseThreshold(0);
842    factorization_->goSparse();
843  }
844 
845  return status;
846}
847/*
848   This does basis housekeeping and does values for in/out variables.
849   Can also decide to re-factorize
850*/
851int 
852ClpSimplex::housekeeping(double objectiveChange)
853{
854  numberIterations_++;
855  changeMade_++; // something has happened
856  // incoming variable
857
858  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
859    <<directionOut_
860    <<directionIn_<<theta_
861    <<dualOut_<<dualIn_<<alpha_
862    <<OsiMessageEol;
863  if (getStatus(sequenceIn_)==ClpSimplex::isFree) {
864    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
865      <<sequenceIn_
866      <<OsiMessageEol;
867  }
868  // change of incoming
869  if (algorithm_<0) {
870    dualOut_ /= alpha_;
871    dualOut_ *= -directionOut_;
872  }
873  char rowcol[]={'R','C'};
874  double cost = cost_[sequenceIn_];
875  double value=valueIn_;
876  if (pivotRow_>=0)
877    pivotVariable_[pivotRow_]=sequenceIn();
878  if (directionIn_==-1) {
879    // as if from upper bound
880    if (sequenceIn_!=sequenceOut_) {
881      // variable becoming basic
882      setStatus(sequenceIn_,ClpSimplex::basic);
883      if (algorithm_<0) {
884        value = upperIn_+dualOut_;
885        dj_[sequenceIn_]=0.0;
886      } else {
887        value = valueIn_-fabs(theta_);
888      }
889    } else {
890      value=lowerIn_;
891      setStatus(sequenceIn_, ClpSimplex::atLowerBound);
892    }
893  } else {
894    // as if from lower bound
895    if (sequenceIn_!=sequenceOut_) {
896      // variable becoming basic
897      setStatus(sequenceIn_,ClpSimplex::basic);
898      if (algorithm_<0) {
899        value = lowerIn_+dualOut_;
900        dj_[sequenceIn_]=0.0;
901      } else {
902        value = valueIn_+fabs(theta_);
903      }
904    } else {
905      value=upperIn_;
906      setStatus(sequenceIn_, ClpSimplex::atUpperBound);
907    }
908  }
909  if (algorithm_<0)
910    objectiveChange += cost*(value-valueIn_);
911  else
912    objectiveChange += dualIn_*(value-valueIn_);
913  solution_[sequenceIn_]=value;
914
915  // outgoing
916  if (sequenceIn_!=sequenceOut_) {
917    assert( getStatus(sequenceOut_)== ClpSimplex::basic);
918    if (algorithm_<0) {
919      if (directionOut_>0) {
920        value = lowerOut_;
921        setStatus(sequenceOut_,ClpSimplex::atLowerBound);
922        dj_[sequenceOut_]=theta_;
923      } else {
924        value = upperOut_;
925        setStatus(sequenceOut_,ClpSimplex::atUpperBound);
926        dj_[sequenceOut_]=-theta_;
927      }
928      solution_[sequenceOut_]=value;
929    } else {
930      if (directionOut_>0) {
931        value = lowerOut_;
932      } else {
933        value = upperOut_;
934      }
935      nonLinearCost_->setOne(sequenceOut_,value);
936      double lowerValue = lower_[sequenceOut_];
937      double upperValue = upper_[sequenceOut_];
938      assert(value>=lowerValue-primalTolerance_&&
939             value<=upperValue+primalTolerance_);
940      // may not be exactly at bound and bounds may have changed
941      if (value<=lowerValue+primalTolerance_) {
942        setStatus(sequenceOut_,ClpSimplex::atLowerBound);
943      } else if (value>=upperValue-primalTolerance_) {
944        setStatus(sequenceOut_,ClpSimplex::atUpperBound);
945      } else {
946        printf("*** variable wandered off bound %g %g %g!\n",
947               lowerValue,value,upperValue);
948        if (value-lowerValue<=upperValue-value) {
949          setStatus(sequenceOut_,ClpSimplex::atLowerBound);
950          value=lowerValue;
951        } else {
952          setStatus(sequenceOut_,ClpSimplex::atUpperBound);
953          value=upperValue;
954        }
955      }
956      solution_[sequenceOut_]=value;
957    }
958  }
959  // change cost and bounds on incoming if primal
960  if (algorithm_>0)
961    nonLinearCost_->setOne(sequenceIn_,solution_[sequenceIn_]); 
962  objectiveValue_ += objectiveChange;
963  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
964    <<numberIterations_<<objectiveValue_
965    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
966    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
967  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
968  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
969  handler_->message()<<OsiMessageEol;
970
971  if (numberIterations_>=maximumIterations_)
972    return 2;
973  // only time to re-factorize if one before real time
974  // this is so user won't be surprised that maximumPivots has exact meaning
975  if (factorization_->pivots()==factorization_->maximumPivots()) {
976    return 1;
977  } else {
978    if (forceFactorization_>0&&
979        factorization_->pivots()==forceFactorization_) {
980      // relax
981      forceFactorization_ *= 2;
982      if (forceFactorization_>factorization_->maximumPivots())
983        forceFactorization_ = -1; //off
984      return 1;
985    } else {
986      // carry on iterating
987      return 0;
988    }
989  }
990}
991// Copy constructor.
992ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
993  ClpModel(rhs),
994  columnPrimalInfeasibility_(0.0),
995  columnPrimalSequence_(-2),
996  rowPrimalInfeasibility_(0.0),
997  rowPrimalSequence_(-2), 
998  columnDualInfeasibility_(0.0),
999  columnDualSequence_(-2),
1000  rowDualInfeasibility_(0.0),
1001  rowDualSequence_(-2),
1002  primalToleranceToGetOptimal_(-1.0),
1003  remainingDualInfeasibility_(0.0),
1004  largeValue_(1.0e15),
1005  largestPrimalError_(0.0),
1006  largestDualError_(0.0),
1007  largestSolutionError_(0.0),
1008  dualBound_(1.0e7),
1009  lower_(NULL),
1010  rowLowerWork_(NULL),
1011  columnLowerWork_(NULL),
1012  upper_(NULL),
1013  rowUpperWork_(NULL),
1014  columnUpperWork_(NULL),
1015  cost_(NULL),
1016  rowObjectiveWork_(NULL),
1017  objectiveWork_(NULL),
1018  alpha_(0.0),
1019  theta_(0.0),
1020  lowerIn_(0.0),
1021  valueIn_(0.0),
1022  upperIn_(0.0),
1023  dualIn_(0.0),
1024  sequenceIn_(-1),
1025  directionIn_(-1),
1026  lowerOut_(-1),
1027  valueOut_(-1),
1028  upperOut_(-1),
1029  dualOut_(-1),
1030  sequenceOut_(-1),
1031  directionOut_(-1),
1032  pivotRow_(-1),
1033  status_(NULL),
1034  dj_(NULL),
1035  rowReducedCost_(NULL),
1036  reducedCostWork_(NULL),
1037  solution_(NULL),
1038  rowActivityWork_(NULL),
1039  columnActivityWork_(NULL),
1040  dualTolerance_(0.0),
1041  primalTolerance_(0.0),
1042  sumDualInfeasibilities_(0.0),
1043  numberDualInfeasibilities_(0),
1044  numberDualInfeasibilitiesWithoutFree_(0),
1045  sumPrimalInfeasibilities_(0.0),
1046  numberPrimalInfeasibilities_(0),
1047  pivotVariable_(NULL),
1048  factorization_(NULL),
1049  numberRefinements_(0),
1050  rowScale_(NULL),
1051  savedSolution_(NULL),
1052  columnScale_(NULL),
1053  scalingFlag_(0),
1054  numberTimesOptimal_(0),
1055  changeMade_(1),
1056  algorithm_(0),
1057  forceFactorization_(-1),
1058  perturbation_(100),
1059  infeasibilityCost_(1.0e7),
1060  nonLinearCost_(NULL),
1061  specialOptions_(0),
1062  lastBadIteration_(-999999)
1063{
1064  int i;
1065  for (i=0;i<6;i++) {
1066    rowArray_[i]=NULL;
1067    columnArray_[i]=NULL;
1068  }
1069  saveStatus_=NULL;
1070  factorization_ = NULL;
1071  dualRowPivot_ = NULL;
1072  primalColumnPivot_ = NULL;
1073  gutsOfDelete(0);
1074  gutsOfCopy(rhs);
1075}
1076// Copy constructor from model
1077ClpSimplex::ClpSimplex(const ClpModel &rhs) :
1078  ClpModel(rhs),
1079  columnPrimalInfeasibility_(0.0),
1080  columnPrimalSequence_(-2),
1081  rowPrimalInfeasibility_(0.0),
1082  rowPrimalSequence_(-2), 
1083  columnDualInfeasibility_(0.0),
1084  columnDualSequence_(-2),
1085  rowDualInfeasibility_(0.0),
1086  rowDualSequence_(-2),
1087  primalToleranceToGetOptimal_(-1.0),
1088  remainingDualInfeasibility_(0.0),
1089  largeValue_(1.0e15),
1090  largestPrimalError_(0.0),
1091  largestDualError_(0.0),
1092  largestSolutionError_(0.0),
1093  dualBound_(1.0e7),
1094  lower_(NULL),
1095  rowLowerWork_(NULL),
1096  columnLowerWork_(NULL),
1097  upper_(NULL),
1098  rowUpperWork_(NULL),
1099  columnUpperWork_(NULL),
1100  cost_(NULL),
1101  rowObjectiveWork_(NULL),
1102  objectiveWork_(NULL),
1103  alpha_(0.0),
1104  theta_(0.0),
1105  lowerIn_(0.0),
1106  valueIn_(0.0),
1107  upperIn_(0.0),
1108  dualIn_(0.0),
1109  sequenceIn_(-1),
1110  directionIn_(-1),
1111  lowerOut_(-1),
1112  valueOut_(-1),
1113  upperOut_(-1),
1114  dualOut_(-1),
1115  sequenceOut_(-1),
1116  directionOut_(-1),
1117  pivotRow_(-1),
1118  status_(NULL),
1119  dj_(NULL),
1120  rowReducedCost_(NULL),
1121  reducedCostWork_(NULL),
1122  solution_(NULL),
1123  rowActivityWork_(NULL),
1124  columnActivityWork_(NULL),
1125  dualTolerance_(0.0),
1126  primalTolerance_(0.0),
1127  sumDualInfeasibilities_(0.0),
1128  numberDualInfeasibilities_(0),
1129  numberDualInfeasibilitiesWithoutFree_(0),
1130  sumPrimalInfeasibilities_(0.0),
1131  numberPrimalInfeasibilities_(0),
1132  pivotVariable_(NULL),
1133  factorization_(NULL),
1134  numberRefinements_(0),
1135  rowScale_(NULL),
1136  savedSolution_(NULL),
1137  columnScale_(NULL),
1138  scalingFlag_(0),
1139  numberTimesOptimal_(0),
1140  changeMade_(1),
1141  algorithm_(0),
1142  forceFactorization_(-1),
1143  perturbation_(100),
1144  infeasibilityCost_(1.0e7),
1145  nonLinearCost_(NULL),
1146  specialOptions_(0),
1147  lastBadIteration_(-999999)
1148{
1149  int i;
1150  for (i=0;i<6;i++) {
1151    rowArray_[i]=NULL;
1152    columnArray_[i]=NULL;
1153  }
1154  saveStatus_=NULL;
1155  // get an empty factorization so we can set tolerances etc
1156  factorization_ = new ClpFactorization();
1157  // say Steepest pricing
1158  dualRowPivot_ = new ClpDualRowSteepest();
1159  // say Dantzig pricing
1160  primalColumnPivot_ = new ClpPrimalColumnDantzig();
1161 
1162}
1163// Assignment operator. This copies the data
1164ClpSimplex & 
1165ClpSimplex::operator=(const ClpSimplex & rhs)
1166{
1167  if (this != &rhs) {
1168    gutsOfDelete(0);
1169    ClpModel::operator=(rhs);
1170    gutsOfCopy(rhs);
1171  }
1172  return *this;
1173}
1174void 
1175ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1176{
1177  lower_ = copyOfArray(rhs.lower_,numberColumns_+numberRows_);
1178  rowLowerWork_ = lower_+numberColumns_;
1179  columnLowerWork_ = lower_;
1180  upper_ = copyOfArray(rhs.upper_,numberColumns_+numberRows_);
1181  rowUpperWork_ = upper_+numberColumns_;
1182  columnUpperWork_ = upper_;
1183  cost_ = copyOfArray(rhs.cost_,(numberColumns_+numberRows_));
1184  objectiveWork_ = cost_;
1185  rowObjectiveWork_ = cost_+numberColumns_;
1186  dj_ = copyOfArray(rhs.dj_,numberRows_+numberColumns_);
1187  if (dj_) {
1188    reducedCostWork_ = dj_;
1189    rowReducedCost_ = dj_+numberColumns_;
1190  }
1191  solution_ = copyOfArray(rhs.solution_,numberRows_+numberColumns_);
1192  if (solution_) {
1193    columnActivityWork_ = solution_;
1194    rowActivityWork_ = solution_+numberColumns_;
1195  }
1196  if (rhs.pivotVariable_) {
1197    pivotVariable_ = new int[numberRows_];
1198    CoinDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
1199  } else {
1200    pivotVariable_=NULL;
1201  }
1202  if (rhs.factorization_) {
1203    factorization_ = new ClpFactorization(*rhs.factorization_);
1204  } else {
1205    factorization_=NULL;
1206  }
1207  rowScale_ = copyOfArray(rhs.rowScale_,numberRows_);
1208  savedSolution_ = copyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1209  columnScale_ = copyOfArray(rhs.columnScale_,numberColumns_);
1210  int i;
1211  for (i=0;i<6;i++) {
1212    rowArray_[i]=NULL;
1213    if (rhs.rowArray_[i]) 
1214      rowArray_[i] = new OsiIndexedVector(*rhs.rowArray_[i]);
1215    columnArray_[i]=NULL;
1216    if (rhs.columnArray_[i]) 
1217      columnArray_[i] = new OsiIndexedVector(*rhs.columnArray_[i]);
1218  }
1219  if (rhs.status_) {
1220    status_ = copyOfArray( rhs.status_,numberColumns_+numberRows_);
1221  }
1222  if (rhs.saveStatus_) {
1223    saveStatus_ = copyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1224  }
1225  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1226  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1227  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1228  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1229  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1230  columnDualSequence_ = rhs.columnDualSequence_;
1231  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1232  rowDualSequence_ = rhs.rowDualSequence_;
1233  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1234  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1235  largeValue_ = rhs.largeValue_;
1236  largestPrimalError_ = rhs.largestPrimalError_;
1237  largestDualError_ = rhs.largestDualError_;
1238  largestSolutionError_ = rhs.largestSolutionError_;
1239  dualBound_ = rhs.dualBound_;
1240  alpha_ = rhs.alpha_;
1241  theta_ = rhs.theta_;
1242  lowerIn_ = rhs.lowerIn_;
1243  valueIn_ = rhs.valueIn_;
1244  upperIn_ = rhs.upperIn_;
1245  dualIn_ = rhs.dualIn_;
1246  sequenceIn_ = rhs.sequenceIn_;
1247  directionIn_ = rhs.directionIn_;
1248  lowerOut_ = rhs.lowerOut_;
1249  valueOut_ = rhs.valueOut_;
1250  upperOut_ = rhs.upperOut_;
1251  dualOut_ = rhs.dualOut_;
1252  sequenceOut_ = rhs.sequenceOut_;
1253  directionOut_ = rhs.directionOut_;
1254  pivotRow_ = rhs.pivotRow_;
1255  numberRefinements_ = rhs.numberRefinements_;
1256  dualTolerance_ = rhs.dualTolerance_;
1257  primalTolerance_ = rhs.primalTolerance_;
1258  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1259  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1260  numberDualInfeasibilitiesWithoutFree_ = 
1261    rhs.numberDualInfeasibilitiesWithoutFree_;
1262  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1263  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1264  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1265  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1266  scalingFlag_ = rhs.scalingFlag_;
1267  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1268  changeMade_ = rhs.changeMade_;
1269  algorithm_ = rhs.algorithm_;
1270  forceFactorization_ = rhs.forceFactorization_;
1271  perturbation_ = rhs.perturbation_;
1272  infeasibilityCost_ = rhs.infeasibilityCost_;
1273  specialOptions_ = rhs.specialOptions_;
1274  lastBadIteration_ = rhs.lastBadIteration_;
1275  if (rhs.nonLinearCost_!=NULL)
1276    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1277}
1278// type == 0 do everything
1279void 
1280ClpSimplex::gutsOfDelete(int type)
1281{
1282  delete [] lower_;
1283  lower_=NULL;
1284  rowLowerWork_=NULL;
1285  columnLowerWork_=NULL;
1286  delete [] upper_;
1287  upper_=NULL;
1288  rowUpperWork_=NULL;
1289  columnUpperWork_=NULL;
1290  delete [] cost_;
1291  cost_=NULL;
1292  objectiveWork_=NULL;
1293  rowObjectiveWork_=NULL;
1294  delete [] dj_;
1295  dj_=NULL;
1296  reducedCostWork_=NULL;
1297  rowReducedCost_=NULL;
1298  delete [] solution_;
1299  solution_=NULL;
1300  rowActivityWork_=NULL;
1301  columnActivityWork_=NULL;
1302  delete [] rowScale_;
1303  rowScale_ = NULL;
1304  delete [] savedSolution_;
1305  savedSolution_ = NULL;
1306  delete [] columnScale_;
1307  columnScale_ = NULL;
1308  delete nonLinearCost_;
1309  nonLinearCost_ = NULL;
1310  int i;
1311  for (i=0;i<6;i++) {
1312    delete rowArray_[i];
1313    rowArray_[i]=NULL;
1314    delete columnArray_[i];
1315    columnArray_[i]=NULL;
1316  }
1317  delete rowCopy_;
1318  rowCopy_=NULL;
1319  if (!type) {
1320    // delete everything
1321    delete [] pivotVariable_;
1322    pivotVariable_=NULL;
1323    delete factorization_;
1324    factorization_ = NULL;
1325    delete dualRowPivot_;
1326    dualRowPivot_ = NULL;
1327    delete primalColumnPivot_;
1328    primalColumnPivot_ = NULL;
1329    delete [] saveStatus_;
1330    saveStatus_=NULL;
1331    delete status_;
1332    status_=NULL;
1333  }
1334}
1335// This sets largest infeasibility and most infeasible
1336void 
1337ClpSimplex::checkPrimalSolution(const double * rowActivities,
1338                                        const double * columnActivities)
1339{
1340  double * solution;
1341  int iRow,iColumn;
1342
1343  objectiveValue_ = 0.0;
1344  // now look at primal solution
1345  columnPrimalInfeasibility_=0.0;
1346  columnPrimalSequence_=-1;
1347  rowPrimalInfeasibility_=0.0;
1348  rowPrimalSequence_=-1;
1349  largestSolutionError_=0.0;
1350  solution = rowActivityWork_;
1351  sumPrimalInfeasibilities_=0.0;
1352  numberPrimalInfeasibilities_=0;
1353  for (iRow=0;iRow<numberRows_;iRow++) {
1354    double infeasibility=0.0;
1355    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1356    if (solution[iRow]>rowUpperWork_[iRow]) {
1357      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1358    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1359      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1360    }
1361    if (infeasibility>primalTolerance_) {
1362      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1363      numberPrimalInfeasibilities_ ++;
1364    }
1365    if (infeasibility>rowPrimalInfeasibility_) {
1366      rowPrimalInfeasibility_=infeasibility;
1367      rowPrimalSequence_=iRow;
1368    }
1369    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1370    if (infeasibility>largestSolutionError_)
1371      largestSolutionError_=infeasibility;
1372    if (rowLowerWork_[iRow]!=rowUpperWork_[iRow])
1373      clearFixed(iRow+numberColumns_);
1374    else
1375      setFixed(iRow+numberColumns_);
1376  }
1377  solution = columnActivityWork_;
1378  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1379    double infeasibility=0.0;
1380    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1381    if (solution[iColumn]>columnUpperWork_[iColumn]) {
1382      infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1383    } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1384      infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1385    }
1386    if (infeasibility>columnPrimalInfeasibility_) {
1387      columnPrimalInfeasibility_=infeasibility;
1388      columnPrimalSequence_=iColumn;
1389    }
1390    if (infeasibility>primalTolerance_) {
1391      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1392      numberPrimalInfeasibilities_ ++;
1393    }
1394    infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1395    if (infeasibility>largestSolutionError_)
1396      largestSolutionError_=infeasibility;
1397    if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
1398      clearFixed(iColumn);
1399    else
1400      setFixed(iColumn);
1401  }
1402}
1403void 
1404ClpSimplex::checkDualSolution()
1405{
1406
1407  double * solution;
1408  int iRow,iColumn;
1409  sumDualInfeasibilities_=0.0;
1410  numberDualInfeasibilities_=0;
1411  numberDualInfeasibilitiesWithoutFree_=0;
1412  columnDualInfeasibility_=0.0;
1413  columnDualSequence_=-1;
1414  rowDualInfeasibility_=0.0;
1415  rowDualSequence_=-1;
1416  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
1417                                   columnPrimalInfeasibility_);
1418  remainingDualInfeasibility_=0.0;
1419  solution = rowActivityWork_;
1420  double dualTolerance=dualTolerance_;
1421  if (algorithm_>0) {
1422    // primal
1423    // we can't really trust infeasibilities if there is dual error
1424    if (largestDualError_>1.0e-6)
1425      dualTolerance *= largestDualError_/1.0e-6;
1426  }
1427  for (iRow=0;iRow<numberRows_;iRow++) {
1428    if (getRowStatus(iRow) != ClpSimplex::basic) {
1429      // not basic
1430      double value = rowReducedCost_[iRow];
1431      double distance;
1432      distance = rowUpperWork_[iRow]-solution[iRow];
1433      if (distance>primalTolerance_) {
1434        // should not be negative
1435        if (value<0.0) {
1436          value = - value;
1437          if (value>rowDualInfeasibility_) {
1438            rowDualInfeasibility_=value;
1439            rowDualSequence_=iRow;
1440          }
1441          if (value>dualTolerance) {
1442            sumDualInfeasibilities_ += value-dualTolerance;
1443            numberDualInfeasibilities_ ++;
1444            if (getRowStatus(iRow) != ClpSimplex::isFree) 
1445              numberDualInfeasibilitiesWithoutFree_ ++;
1446            // maybe we can make feasible by increasing tolerance
1447            if (distance<largeValue_) {
1448              if (distance>primalToleranceToGetOptimal_)
1449                primalToleranceToGetOptimal_=distance;
1450            } else {
1451              //gap too big for any tolerance
1452              remainingDualInfeasibility_=
1453                max(remainingDualInfeasibility_,value);
1454            }
1455          }
1456        }
1457      }
1458      distance = solution[iRow] -rowLowerWork_[iRow];
1459      if (distance>primalTolerance_) {
1460        // should not be positive
1461        if (value>0.0) {
1462          if (value>rowDualInfeasibility_) {
1463            rowDualInfeasibility_=value;
1464            rowDualSequence_=iRow;
1465          }
1466          if (value>dualTolerance) {
1467            sumDualInfeasibilities_ += value-dualTolerance;
1468            numberDualInfeasibilities_ ++;
1469            if (getRowStatus(iRow) != ClpSimplex::isFree) 
1470              numberDualInfeasibilitiesWithoutFree_ ++;
1471            // maybe we can make feasible by increasing tolerance
1472            if (distance<largeValue_&&
1473                distance>primalToleranceToGetOptimal_)
1474              primalToleranceToGetOptimal_=distance;
1475          }
1476        }
1477      }
1478    }
1479  }
1480  solution = columnActivityWork_;
1481  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1482    if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1483      // not basic
1484      double value = reducedCostWork_[iColumn];
1485      double distance;
1486      distance = columnUpperWork_[iColumn]-solution[iColumn];
1487      if (distance>primalTolerance_) {
1488        // should not be negative
1489        if (value<0.0) {
1490          value = - value;
1491          if (value>columnDualInfeasibility_) {
1492            columnDualInfeasibility_=value;
1493            columnDualSequence_=iColumn;
1494          }
1495          if (value>dualTolerance) {
1496            sumDualInfeasibilities_ += value-dualTolerance;
1497            numberDualInfeasibilities_ ++;
1498            if (getColumnStatus(iColumn) != ClpSimplex::isFree) 
1499              numberDualInfeasibilitiesWithoutFree_ ++;
1500            // maybe we can make feasible by increasing tolerance
1501            if (distance<largeValue_) {
1502              if (distance>primalToleranceToGetOptimal_)
1503                primalToleranceToGetOptimal_=distance;
1504            } else {
1505              //gap too big for any tolerance
1506              remainingDualInfeasibility_=
1507                max(remainingDualInfeasibility_,value);
1508            }
1509          }
1510        }
1511      }
1512      distance = solution[iColumn] -columnLowerWork_[iColumn];
1513      if (distance>primalTolerance_) {
1514        // should not be positive
1515        if (value>0.0) {
1516          if (value>columnDualInfeasibility_) {
1517            columnDualInfeasibility_=value;
1518            columnDualSequence_=iColumn;
1519          }
1520          if (value>dualTolerance) {
1521            sumDualInfeasibilities_ += value-dualTolerance;
1522            numberDualInfeasibilities_ ++;
1523            if (getColumnStatus(iColumn) != ClpSimplex::isFree) 
1524              numberDualInfeasibilitiesWithoutFree_ ++;
1525            // maybe we can make feasible by increasing tolerance
1526            if (distance<largeValue_&&
1527                distance>primalToleranceToGetOptimal_)
1528              primalToleranceToGetOptimal_=distance;
1529          }
1530        }
1531      }
1532    }
1533  }
1534}
1535/*
1536  Unpacks one column of the matrix into indexed array
1537*/
1538void 
1539ClpSimplex::unpack(OsiIndexedVector * rowArray)
1540{
1541  rowArray->clear();
1542  if (sequenceIn_>=numberColumns_) {
1543    //slack
1544    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
1545  } else {
1546    // column
1547    matrix_->unpack(this,rowArray,sequenceIn_);
1548  }
1549}
1550void 
1551ClpSimplex::unpack(OsiIndexedVector * rowArray,int sequence)
1552{
1553  rowArray->clear();
1554  if (sequence>=numberColumns_) {
1555    //slack
1556    rowArray->insert(sequence-numberColumns_,-1.0);
1557  } else {
1558    // column
1559    matrix_->unpack(this,rowArray,sequence);
1560  }
1561}
1562void
1563ClpSimplex::createRim(int what,bool makeRowCopy)
1564{
1565  if ((what&16)!=0) {
1566    // move information to work arrays
1567    // row reduced costs
1568    if (!dj_) {
1569      dj_ = new double[numberRows_+numberColumns_];
1570      reducedCostWork_ = dj_;
1571      rowReducedCost_ = dj_+numberColumns_;
1572      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1573    }
1574    if (!solution_) {
1575      solution_ = new double[numberRows_+numberColumns_];
1576      columnActivityWork_ = solution_;
1577      rowActivityWork_ = solution_+numberColumns_;
1578      memcpy(columnActivityWork_,columnActivity_,
1579             numberColumns_*sizeof(double));
1580      memcpy(rowActivityWork_,rowActivity_,
1581             numberRows_*sizeof(double));
1582    }
1583   
1584    if (makeRowCopy) {
1585      delete rowCopy_;
1586      // may return NULL if can't give row copy
1587      rowCopy_ = matrix_->reverseOrderedCopy();
1588    }
1589  }
1590  int i;
1591  if ((what&4)!=0) {
1592    delete [] cost_;
1593    cost_ = new double[numberColumns_+numberRows_];
1594    objectiveWork_ = cost_;
1595    rowObjectiveWork_ = cost_+numberColumns_;
1596    memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
1597    if (rowObjective_)
1598      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1599    else
1600      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1601  }
1602  if ((what&1)!=0) {
1603    delete [] lower_;
1604    delete [] upper_;
1605    lower_ = new double[numberColumns_+numberRows_];
1606    upper_ = new double[numberColumns_+numberRows_];
1607    rowLowerWork_ = lower_+numberColumns_;
1608    columnLowerWork_ = lower_;
1609    rowUpperWork_ = upper_+numberColumns_;
1610    columnUpperWork_ = upper_;
1611    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1612    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1613    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1614    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1615    // clean up any mismatches on infinity
1616    for (i=0;i<numberColumns_;i++) {
1617      if (columnLowerWork_[i]<-1.0e30)
1618        columnLowerWork_[i] = -DBL_MAX;
1619      if (columnUpperWork_[i]>1.0e30)
1620        columnUpperWork_[i] = DBL_MAX;
1621    }
1622    // clean up any mismatches on infinity
1623    for (i=0;i<numberRows_;i++) {
1624      if (rowLowerWork_[i]<-1.0e30)
1625        rowLowerWork_[i] = -DBL_MAX;
1626      if (rowUpperWork_[i]>1.0e30)
1627        rowUpperWork_[i] = DBL_MAX;
1628    }
1629  }
1630  // do scaling if needed
1631  if (scalingFlag_>0&&!rowScale_) {
1632    if (matrix_->scale(this))
1633      scalingFlag_=-scalingFlag_; // not scaled after all
1634  }
1635  if ((what&4)!=0) {
1636    double direction = optimizationDirection_;
1637    // direction is actually scale out not scale in
1638    if (direction)
1639      direction = 1.0/direction;
1640    // but also scale by scale factors
1641    // not really sure about row scaling
1642    if (!rowScale_) {
1643      if (direction!=1.0) {
1644        for (i=0;i<numberRows_;i++)
1645          rowObjectiveWork_[i] *= direction;
1646        for (i=0;i<numberColumns_;i++)
1647          objectiveWork_[i] *= direction;
1648      }
1649    } else {
1650      for (i=0;i<numberRows_;i++)
1651        rowObjectiveWork_[i] *= direction/rowScale_[i];
1652      for (i=0;i<numberColumns_;i++)
1653        objectiveWork_[i] *= direction*columnScale_[i];
1654    }
1655  }
1656  if ((what&1)!=0&&rowScale_) {
1657    for (i=0;i<numberColumns_;i++) {
1658      double multiplier = 1.0/columnScale_[i];
1659      if (columnLowerWork_[i]>-1.0e50)
1660        columnLowerWork_[i] *= multiplier;
1661      if (columnUpperWork_[i]<1.0e50)
1662        columnUpperWork_[i] *= multiplier;
1663    }
1664    for (i=0;i<numberRows_;i++) {
1665      if (rowLowerWork_[i]>-1.0e50)
1666        rowLowerWork_[i] *= rowScale_[i];
1667      if (rowUpperWork_[i]<1.0e50)
1668        rowUpperWork_[i] *= rowScale_[i];
1669    }
1670  }
1671  if ((what&8)!=0&&rowScale_) {
1672    // on entry
1673    for (i=0;i<numberColumns_;i++) {
1674      columnActivityWork_[i] /= columnScale_[i];
1675      reducedCostWork_[i] *= columnScale_[i];
1676    }
1677    for (i=0;i<numberRows_;i++) {
1678      rowActivityWork_[i] *= rowScale_[i];
1679      dual_[i] /= rowScale_[i];
1680    }
1681  }
1682 
1683  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1684  // maybe we need to move scales to SimplexModel for factorization?
1685  if ((what&8)!=0&&!pivotVariable_) {
1686    pivotVariable_=new int[numberRows_];
1687  }
1688
1689  if ((what&16)!=0&&!rowArray_[2]) {
1690    // get some arrays
1691    int iRow,iColumn;
1692    // these are "indexed" arrays so we always know where nonzeros are
1693    /**********************************************************
1694      rowArray_[3] is long enough for columns as well
1695    *********************************************************/
1696    for (iRow=0;iRow<4;iRow++) {
1697      if (!rowArray_[iRow]) {
1698        rowArray_[iRow]=new OsiIndexedVector();
1699        int length =numberRows_+factorization_->maximumPivots();
1700        if (iRow==3)
1701          length = max(length,numberColumns_);
1702        rowArray_[iRow]->reserve(length);
1703      }
1704    }
1705   
1706    for (iColumn=0;iColumn<2;iColumn++) {
1707      if (!columnArray_[iColumn]) {
1708        columnArray_[iColumn]=new OsiIndexedVector();
1709        columnArray_[iColumn]->reserve(numberColumns_);
1710      }
1711    }   
1712  }
1713
1714}
1715void
1716ClpSimplex::deleteRim()
1717{
1718  int i;
1719  if (problemStatus_!=1&&problemStatus_!=2) {
1720    delete [] ray_;
1721    ray_=NULL;
1722  }
1723  if (rowScale_) {
1724    for (i=0;i<numberColumns_;i++) {
1725      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1726      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1727    }
1728    for (i=0;i<numberRows_;i++) {
1729      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1730      dual_[i] *= rowScale_[i];
1731    }
1732    if (problemStatus_==2) {
1733      for (i=0;i<numberColumns_;i++) {
1734        ray_[i] *= columnScale_[i];
1735      }
1736    } else if (problemStatus_==1) {
1737      for (i=0;i<numberRows_;i++) {
1738        ray_[i] *= rowScale_[i];
1739      }
1740    }
1741  } else {
1742    for (i=0;i<numberColumns_;i++) {
1743      columnActivity_[i] = columnActivityWork_[i];
1744      reducedCost_[i] = reducedCostWork_[i];
1745    }
1746    for (i=0;i<numberRows_;i++) {
1747      rowActivity_[i] = rowActivityWork_[i];
1748    }
1749  }
1750  // direction may have been modified by scaling - clean up
1751  if (optimizationDirection_>0.0)
1752    optimizationDirection_ = 1.0;
1753  else if (optimizationDirection_<0.0)
1754    optimizationDirection_ = -1.0;
1755  // scaling may have been turned off
1756  scalingFlag_ = abs(scalingFlag_);
1757  gutsOfDelete(1);
1758}
1759void 
1760ClpSimplex::setDualBound(double value)
1761{
1762  if (value>0.0)
1763    dualBound_=value;
1764}
1765void 
1766ClpSimplex::setInfeasibilityCost(double value)
1767{
1768  if (value>0.0)
1769    infeasibilityCost_=value;
1770}
1771void ClpSimplex::setnumberRefinements( int value) 
1772{
1773  if (value>=0&&value<10)
1774    numberRefinements_=value;
1775}
1776// Sets row pivot choice algorithm in dual
1777void 
1778ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
1779{
1780  delete dualRowPivot_;
1781  dualRowPivot_ = choice.clone(true);
1782}
1783// Sets row pivot choice algorithm in dual
1784void 
1785ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
1786{
1787  delete primalColumnPivot_;
1788  primalColumnPivot_ = choice.clone(true);
1789}
1790// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
1791void 
1792ClpSimplex::scaling(int mode)
1793{
1794  if (mode>0&&mode<3) {
1795    scalingFlag_=mode;
1796  } else if (!mode) {
1797    scalingFlag_=0;
1798    delete [] rowScale_;
1799    rowScale_ = NULL;
1800    delete [] columnScale_;
1801    columnScale_ = NULL;
1802  }
1803}
1804// Sets up basis
1805void 
1806ClpSimplex::setBasis ( const OsiWarmStartBasis & basis)
1807{
1808  // transform basis to status arrays
1809  int iRow,iColumn;
1810  if (!status_) {
1811    /*
1812      get status arrays
1813      OsiWarmStartBasis would seem to have overheads and we will need
1814      extra bits anyway.
1815    */
1816    status_ = new unsigned char [numberColumns_+numberRows_];
1817    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
1818  }
1819  OsiWarmStartBasis basis2 = basis;
1820  // resize if necessary
1821  basis2.resize(numberRows_,numberColumns_);
1822  // move status
1823  for (iRow=0;iRow<numberRows_;iRow++) {
1824    setRowStatus(iRow,
1825                 (Status) basis2.getArtifStatus(iRow));
1826  }
1827  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1828    setColumnStatus(iColumn,
1829                    (Status) basis2.getStructStatus(iColumn));
1830  }
1831}
1832// Passes in factorization
1833void 
1834ClpSimplex::setFactorization( ClpFactorization & factorization)
1835{
1836  delete factorization_;
1837  factorization_= new ClpFactorization(factorization);
1838}
1839// Warm start
1840OsiWarmStartBasis 
1841ClpSimplex::getBasis() const
1842{
1843  int iRow,iColumn;
1844  OsiWarmStartBasis basis;
1845  basis.setSize(numberColumns_,numberRows_);
1846
1847  if(status_) {
1848    for (iRow=0;iRow<numberRows_;iRow++) {
1849      basis.setArtifStatus(iRow,
1850                           (OsiWarmStartBasis::Status) getRowStatus(iRow));
1851    }
1852    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1853      basis.setStructStatus(iColumn,
1854                       (OsiWarmStartBasis::Status) getColumnStatus(iColumn));
1855    }
1856  }
1857  return basis;
1858}
1859void 
1860ClpSimplex::times(double scalar,
1861                  const double * x, double * y) const
1862{
1863  if (rowScale_)
1864    matrix_->times(scalar,x,y,rowScale_,columnScale_);
1865  else
1866    matrix_->times(scalar,x,y);
1867}
1868void 
1869ClpSimplex::transposeTimes(double scalar,
1870                           const double * x, double * y) const 
1871{
1872  if (rowScale_)
1873    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
1874  else
1875    matrix_->transposeTimes(scalar,x,y);
1876}
1877/* Perturbation:
1878   -50 to +50 - perturb by this power of ten (-6 sounds good)
1879   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1880   101 - we are perturbed
1881   102 - don't try perturbing again
1882   default is 100
1883*/
1884void 
1885ClpSimplex::setPerturbation(int value)
1886{
1887  if(value<=100&&value >=-50) {
1888    perturbation_=value;
1889  } 
1890}
1891// Sparsity on or off
1892bool 
1893ClpSimplex::sparseFactorization() const
1894{
1895  return factorization_->sparseThreshold();
1896}
1897void 
1898ClpSimplex::setSparseFactorization(bool value)
1899{
1900  if (value) {
1901    if (!factorization_->sparseThreshold())
1902      factorization_->goSparse();
1903  } else {
1904    factorization_->sparseThreshold(0);
1905  }
1906}
1907/* Tightens primal bounds to make dual faster.  Unless
1908   fixed, bounds are slightly looser than they could be.
1909   This is to make dual go faster and is probably not needed
1910   with a presolve.  Returns non-zero if problem infeasible
1911*/
1912int 
1913ClpSimplex::tightenPrimalBounds()
1914{
1915 
1916  // Get a row copy in standard format
1917  OsiPackedMatrix copy;
1918  copy.reverseOrderedCopyOf(*matrix());
1919  // get matrix data pointers
1920  const int * column = copy.getIndices();
1921  const int * rowStart = copy.getVectorStarts();
1922  const int * rowLength = copy.getVectorLengths(); 
1923  const double * element = copy.getElements();
1924  int numberChanged=1,iPass=0;
1925  double large = largeValue(); // treat bounds > this as infinite
1926  int numberInfeasible=0;
1927  int totalTightened = 0;
1928
1929  double tolerance = primalTolerance();
1930
1931
1932  // Save column bounds
1933  double * saveLower = new double [numberColumns_];
1934  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
1935  double * saveUpper = new double [numberColumns_];
1936  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
1937#define MAXPASS 10 
1938
1939  // Loop round seeing if we can tighten bounds
1940  // Would be faster to have a stack of possible rows
1941  // and we put altered rows back on stack
1942
1943  int iRow, iColumn;
1944
1945  while(numberChanged) {
1946
1947    numberChanged = 0; // Bounds tightened this pass
1948   
1949    if (iPass==MAXPASS) break;
1950    iPass++;
1951   
1952    for (iRow = 0; iRow < numberRows_; iRow++) {
1953
1954      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
1955
1956        // possible row
1957        int infiniteUpper = 0;
1958        int infiniteLower = 0;
1959        double maximumUp = 0.0;
1960        double maximumDown = 0.0;
1961        double newBound;
1962        int rStart = rowStart[iRow];
1963        int rEnd = rowStart[iRow]+rowLength[iRow];
1964        int j;
1965
1966        // Compute possible lower and upper ranges
1967
1968        for (j = rStart; j < rEnd; ++j) {
1969          double value=element[j];
1970          iColumn = column[j];
1971          if (value > 0.0) {
1972            if (columnUpper_[iColumn] >= large) {
1973              maximumUp = 1e31;
1974              ++infiniteUpper;
1975            } else {
1976              maximumUp += columnUpper_[iColumn] * value;
1977            }
1978            if (columnLower_[iColumn] <= -large) {
1979              maximumDown = -1e31;
1980              ++infiniteLower;
1981            } else {
1982              maximumDown += columnLower_[iColumn] * value;
1983            }
1984          } else if (value<0.0) {
1985            if (columnUpper_[iColumn] >= large) {
1986              maximumDown = -1e31;
1987              ++infiniteLower;
1988            } else {
1989              maximumDown += columnUpper_[iColumn] * value;
1990            }
1991            if (columnLower_[iColumn] <= -large) {
1992              maximumUp = 1e31;
1993              ++infiniteUpper;
1994            } else {
1995              maximumUp += columnLower_[iColumn] * value;
1996            }
1997          }
1998        }
1999        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2000            maximumDown >= rowLower_[iRow] - tolerance) {
2001
2002          // Row is redundant - make totally free
2003          rowLower_[iRow]=-DBL_MAX;
2004          rowUpper_[iRow]=DBL_MAX;
2005
2006        } else {
2007          if (maximumUp < rowLower_[iRow] -tolerance ||
2008              maximumDown > rowUpper_[iRow]+tolerance) {
2009            // problem is infeasible - exit at once
2010            numberInfeasible++;
2011            break;
2012          }
2013
2014          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2015            for (j = rStart; j < rEnd; ++j) {
2016              double value=element[j];
2017              iColumn = column[j];
2018              if (value > 0.0) {
2019                if (columnUpper_[iColumn] < large) {
2020                  newBound = columnUpper_[iColumn] + 
2021                    (rowLower_[iRow] - maximumUp) / value;
2022                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2023                    // Tighten the lower bound
2024
2025                    columnLower_[iColumn] = newBound;
2026                    ++numberChanged;
2027
2028                    // check infeasible (relaxed)
2029                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2030                        -100.0*tolerance) 
2031                      numberInfeasible++;
2032                    infiniteLower=1; // skip looking at other bound
2033                  }
2034                }
2035              } else {
2036                if (columnLower_[iColumn] > -large) {
2037                  newBound = columnLower_[iColumn] + 
2038                    (rowLower_[iRow] - maximumUp) / value;
2039                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2040                    // Tighten the upper bound
2041
2042                    columnUpper_[iColumn] = newBound;
2043                    ++numberChanged;
2044
2045                    // check infeasible (relaxed)
2046                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2047                        -100.0*tolerance) 
2048                      numberInfeasible++;
2049                    infiniteLower=1; // skip looking at other bound
2050                  }
2051                }
2052              }
2053            }
2054          }
2055         
2056          // Try other way
2057          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2058            for (j = rStart; j < rEnd; ++j) {
2059              double value=element[j];
2060              iColumn = column[j];
2061              if (value < 0.0) {
2062                if (columnUpper_[iColumn] < large) {
2063                  newBound = columnUpper_[iColumn] + 
2064                    (rowUpper_[iRow] - maximumDown) / value;
2065                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2066                    // Tighten the lower bound
2067
2068                    columnLower_[iColumn] = newBound;
2069                    ++numberChanged;
2070
2071                    // check infeasible (relaxed)
2072                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2073                        -100.0*tolerance) 
2074                      numberInfeasible++;
2075                  }
2076                } 
2077              } else {
2078                if (columnLower_[iColumn] > -large) {
2079                  newBound = columnLower_[iColumn] + 
2080                    (rowUpper_[iRow] - maximumDown) / value;
2081                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2082                    // Tighten the upper bound
2083
2084                    columnUpper_[iColumn] = newBound;
2085                    ++numberChanged;
2086
2087                    // check infeasible (relaxed)
2088                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2089                        -100.0*tolerance) 
2090                      numberInfeasible++;
2091                  }
2092                }
2093              }
2094            }
2095          }
2096        }
2097      }
2098    }
2099    totalTightened += numberChanged;
2100    if (numberInfeasible) break;
2101  }
2102  if (!numberInfeasible) {
2103    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2104      <<totalTightened
2105      <<OsiMessageEol;
2106    // Set bounds slightly loose
2107    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2108      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2109        // fix
2110        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2111          columnLower_[iColumn]=columnUpper_[iColumn];
2112        else
2113          columnUpper_[iColumn]=columnLower_[iColumn];
2114      } else {
2115        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2116          // relax a bit
2117          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2118                                      saveUpper[iColumn]);
2119        }
2120        if (columnLower_[iColumn]>saveLower[iColumn]) {
2121          // relax a bit
2122          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2123                                      saveLower[iColumn]);
2124        }
2125      }
2126    }
2127  } else {
2128    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2129      <<numberInfeasible
2130      <<OsiMessageEol;
2131    // restore column bounds
2132    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2133    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2134  }
2135  delete [] saveLower;
2136  delete [] saveUpper;
2137  return (numberInfeasible);
2138}
2139// dual
2140#include "ClpSimplexDual.hpp"
2141int ClpSimplex::dual ( )
2142{
2143  return ((ClpSimplexDual *) this)->dual();
2144}
2145// primal
2146#include "ClpSimplexPrimal.hpp"
2147int ClpSimplex::primal (int ifValuesPass )
2148{
2149  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2150}
Note: See TracBrowser for help on using the repository browser.