source: trunk/ClpSimplex.cpp @ 2

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

Adding Clp to development branch

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