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

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

Hope this works from wincvs

Fix error in values pass

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 60.2 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 unless initial factorization
505   when total number of singularities is returned */
506int ClpSimplex::internalFactorize ( int solveType)
507{
508
509  int iRow,iColumn;
510  int totalSlacks=numberRows_;
511
512  bool valuesPass=false;
513  if (solveType>=10) {
514    valuesPass=true;
515    solveType -= 10;
516  }
517  if (!solveType) {
518    if (!valuesPass) {
519      // not values pass so set to bounds
520      bool allSlack=true;
521      if (status_) {
522        for (iRow=0;iRow<numberRows_;iRow++) {
523          if (getRowStatus(iRow)!=ClpSimplex::basic) {
524            allSlack=false;
525            break;
526          }
527        }
528      }
529      if (!allSlack) {
530        // set values from warm start (if sensible)
531        int numberBasic=0;
532        for (iRow=0;iRow<numberRows_;iRow++) {
533          switch(getRowStatus(iRow)) {
534           
535          case ClpSimplex::basic:
536            numberBasic++;
537            break;
538          case ClpSimplex::isFree:
539            assert(rowLowerWork_[iRow]<-largeValue_);
540            assert(rowUpperWork_[iRow]>largeValue_);
541            rowActivityWork_[iRow]=0.0;
542            break;
543          case ClpSimplex::atUpperBound:
544            rowActivityWork_[iRow]=rowUpperWork_[iRow];
545            if (rowActivityWork_[iRow]>largeValue_) {
546              assert(rowLowerWork_[iRow]>-largeValue_);
547              rowActivityWork_[iRow]=rowLowerWork_[iRow];
548              setRowStatus(iRow,ClpSimplex::atLowerBound);
549            }
550            break;
551          case ClpSimplex::atLowerBound:
552            rowActivityWork_[iRow]=rowLowerWork_[iRow];
553            if (rowActivityWork_[iRow]<-largeValue_) {
554              assert(rowUpperWork_[iRow]<largeValue_);
555              rowActivityWork_[iRow]=rowUpperWork_[iRow];
556              setRowStatus(iRow,ClpSimplex::atUpperBound);
557            }
558            break;
559          case ClpSimplex::superBasic:
560            abort();
561            break;
562          }
563        }
564        totalSlacks=numberBasic;
565        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
566          switch(getColumnStatus(iColumn)) {
567           
568          case ClpSimplex::basic:
569            if (numberBasic==numberRows_) {
570              // take out of basis
571              if (columnLowerWork_[iColumn]>-largeValue_) {
572                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
573                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
574                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
575                  setColumnStatus(iColumn,ClpSimplex::atLowerBound);
576                } else {
577                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
578                  setColumnStatus(iColumn,ClpSimplex::atUpperBound);
579                }
580              } else if (columnUpperWork_[iColumn]<largeValue_) {
581                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
582                setColumnStatus(iColumn,ClpSimplex::atUpperBound);
583              } else {
584                columnActivityWork_[iColumn]=0.0;
585                setColumnStatus(iColumn,ClpSimplex::isFree);
586              }
587            } else {
588              numberBasic++;
589            }
590            break;
591          case ClpSimplex::isFree:
592            columnActivityWork_[iColumn]=0.0;
593            break;
594          case ClpSimplex::atUpperBound:
595            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
596            if (columnActivityWork_[iColumn]>largeValue_) {
597              assert(columnLowerWork_[iColumn]>-largeValue_);
598              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
599              setColumnStatus(iColumn,ClpSimplex::atLowerBound);
600            }
601            break;
602          case ClpSimplex::atLowerBound:
603            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
604            if (columnActivityWork_[iColumn]<-largeValue_) {
605              assert(columnUpperWork_[iColumn]<largeValue_);
606              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
607              setColumnStatus(iColumn,ClpSimplex::atUpperBound);
608            }
609            break;
610          case ClpSimplex::superBasic:
611            abort();
612            break;
613          }
614        }
615      } else {
616        //#define TESTFREE
617        // all slack basis
618        int numberBasic=0;
619        // changed to put free variables in basis
620        if (!status_) {
621          status_ = new unsigned char [numberColumns_+numberRows_];
622          memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
623        }
624        for (iRow=0;iRow<numberRows_;iRow++) {
625          double lower=rowLowerWork_[iRow];
626          double upper=rowUpperWork_[iRow];
627          if (lower>-largeValue_||upper<largeValue_) {
628            if (fabs(lower)<=fabs(upper)) {
629              setRowStatus(iRow,ClpSimplex::atLowerBound);
630              rowActivityWork_[iRow]=lower;
631            } else {
632              setRowStatus(iRow,ClpSimplex::atUpperBound);
633              rowActivityWork_[iRow]=upper;
634            }
635          } else {
636            setRowStatus(iRow,ClpSimplex::isFree);
637            rowActivityWork_[iRow]=0.0;
638          }
639#ifdef TESTFREE
640          setRowStatus(iRow,ClpSimplex::basic);
641          numberBasic++;
642#endif
643        }
644        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
645          double lower=columnLowerWork_[iColumn];
646          double upper=columnUpperWork_[iColumn];
647          double big_bound = largeValue_;
648          if (lower>-big_bound||upper<big_bound) {
649            if (fabs(lower)<=fabs(upper)) {
650              setColumnStatus(iColumn,ClpSimplex::atLowerBound);
651              columnActivityWork_[iColumn]=lower;
652            } else {
653              setColumnStatus(iColumn,ClpSimplex::atUpperBound);
654              columnActivityWork_[iColumn]=upper;
655            }
656          } else {
657#ifndef TESTFREE
658            numberBasic++;
659            setColumnStatus(iColumn,ClpSimplex::basic);
660#else
661            setColumnStatus(iColumn,ClpSimplex::isFree);
662#endif
663            columnActivityWork_[iColumn]=0.0;
664          }
665        }
666        assert(numberBasic<=numberRows_); // problems if too many free
667        if (!numberBasic) {
668          // might as well do all slack basis
669          for (iRow=0;iRow<numberRows_;iRow++) {
670            setRowStatus(iRow,ClpSimplex::basic);
671          }
672        }
673      }
674    } else {
675      // values pass has less coding
676      // make row activities correct
677      memset(rowActivityWork_,0,numberRows_*sizeof(double));
678      times(1.0,columnActivityWork_,rowActivityWork_);
679      if (status_) {
680        // set values from warm start (if sensible)
681        int numberBasic=0;
682        for (iRow=0;iRow<numberRows_;iRow++) {
683          if (getRowStatus(iRow)==ClpSimplex::basic) 
684            numberBasic++;
685          else
686            setRowStatus(iRow,ClpSimplex::superBasic);
687        }
688        totalSlacks=numberBasic;
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            if (!valuesPass) {
793              double lower=columnLowerWork_[iColumn];
794              double upper=columnUpperWork_[iColumn];
795              double value=columnActivityWork_[iColumn];
796              if (lower>-largeValue_||upper<largeValue_) {
797                if (fabs(value-lower)<fabs(value-upper)) {
798                  setColumnStatus(iColumn,ClpSimplex::atLowerBound);
799                  columnActivityWork_[iColumn]=lower;
800                } else {
801                  setColumnStatus(iColumn,ClpSimplex::atUpperBound);
802                  columnActivityWork_[iColumn]=upper;
803                }
804              } else {
805                setColumnStatus(iColumn,ClpSimplex::isFree);
806              }
807            } else {
808              setColumnStatus(iColumn,ClpSimplex::superBasic);
809            }
810          }
811        }
812        for (iRow=0;iRow<numberRows_;iRow++) {
813          int iSequence=pivotVariable_[iRow];
814          if (iSequence>=0) {
815            // basic
816            if (iSequence>=numberColumns_) {
817              // slack in - leave
818              assert (iSequence-numberColumns_==iRow);
819            } else {
820              // put back structural
821              setColumnStatus(iSequence,ClpSimplex::basic);
822            }
823          } else {
824            // put in slack
825            setRowStatus(iRow,ClpSimplex::basic);
826          }
827        }
828        // signal repeat
829        status=-99;
830      }
831    }
832  } 
833  delete [] rowIsBasic;
834  delete [] columnIsBasic;
835  if (status) {
836    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
837      <<status
838      <<OsiMessageEol;
839    return -1;
840  } else if (!solveType) {
841    // Initial basis - return number of singularities
842    int numberSlacks=0;
843    for (iRow=0;iRow<numberRows_;iRow++) {
844      if (getRowStatus(iRow) == ClpSimplex::basic)
845        numberSlacks++;
846    }
847    status= numberSlacks-totalSlacks;
848  }
849
850  // sparse methods
851  if (factorization_->sparseThreshold()) {
852    // get default value
853    factorization_->sparseThreshold(0);
854    factorization_->goSparse();
855  }
856 
857  return status;
858}
859/*
860   This does basis housekeeping and does values for in/out variables.
861   Can also decide to re-factorize
862*/
863int 
864ClpSimplex::housekeeping(double objectiveChange)
865{
866  numberIterations_++;
867  changeMade_++; // something has happened
868  // incoming variable
869
870  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
871    <<directionOut_
872    <<directionIn_<<theta_
873    <<dualOut_<<dualIn_<<alpha_
874    <<OsiMessageEol;
875  if (getStatus(sequenceIn_)==ClpSimplex::isFree) {
876    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
877      <<sequenceIn_
878      <<OsiMessageEol;
879  }
880  // change of incoming
881  if (algorithm_<0) {
882    dualOut_ /= alpha_;
883    dualOut_ *= -directionOut_;
884  }
885  char rowcol[]={'R','C'};
886  double cost = cost_[sequenceIn_];
887  double value=valueIn_;
888  if (pivotRow_>=0)
889    pivotVariable_[pivotRow_]=sequenceIn();
890  if (directionIn_==-1) {
891    // as if from upper bound
892    if (sequenceIn_!=sequenceOut_) {
893      // variable becoming basic
894      setStatus(sequenceIn_,ClpSimplex::basic);
895      if (algorithm_<0) {
896        value = upperIn_+dualOut_;
897        dj_[sequenceIn_]=0.0;
898      } else {
899        value = valueIn_-fabs(theta_);
900      }
901    } else {
902      value=lowerIn_;
903      setStatus(sequenceIn_, ClpSimplex::atLowerBound);
904    }
905  } else {
906    // as if from lower bound
907    if (sequenceIn_!=sequenceOut_) {
908      // variable becoming basic
909      setStatus(sequenceIn_,ClpSimplex::basic);
910      if (algorithm_<0) {
911        value = lowerIn_+dualOut_;
912        dj_[sequenceIn_]=0.0;
913      } else {
914        value = valueIn_+fabs(theta_);
915      }
916    } else {
917      value=upperIn_;
918      setStatus(sequenceIn_, ClpSimplex::atUpperBound);
919    }
920  }
921  if (algorithm_<0)
922    objectiveChange += cost*(value-valueIn_);
923  else
924    objectiveChange += dualIn_*(value-valueIn_);
925  solution_[sequenceIn_]=value;
926
927  // outgoing
928  if (sequenceIn_!=sequenceOut_) {
929    assert( getStatus(sequenceOut_)== ClpSimplex::basic);
930    if (algorithm_<0) {
931      if (directionOut_>0) {
932        value = lowerOut_;
933        setStatus(sequenceOut_,ClpSimplex::atLowerBound);
934        dj_[sequenceOut_]=theta_;
935      } else {
936        value = upperOut_;
937        setStatus(sequenceOut_,ClpSimplex::atUpperBound);
938        dj_[sequenceOut_]=-theta_;
939      }
940      solution_[sequenceOut_]=value;
941    } else {
942      if (directionOut_>0) {
943        value = lowerOut_;
944      } else {
945        value = upperOut_;
946      }
947      nonLinearCost_->setOne(sequenceOut_,value);
948      double lowerValue = lower_[sequenceOut_];
949      double upperValue = upper_[sequenceOut_];
950      assert(value>=lowerValue-primalTolerance_&&
951             value<=upperValue+primalTolerance_);
952      // may not be exactly at bound and bounds may have changed
953      if (value<=lowerValue+primalTolerance_) {
954        setStatus(sequenceOut_,ClpSimplex::atLowerBound);
955      } else if (value>=upperValue-primalTolerance_) {
956        setStatus(sequenceOut_,ClpSimplex::atUpperBound);
957      } else {
958        printf("*** variable wandered off bound %g %g %g!\n",
959               lowerValue,value,upperValue);
960        if (value-lowerValue<=upperValue-value) {
961          setStatus(sequenceOut_,ClpSimplex::atLowerBound);
962          value=lowerValue;
963        } else {
964          setStatus(sequenceOut_,ClpSimplex::atUpperBound);
965          value=upperValue;
966        }
967      }
968      solution_[sequenceOut_]=value;
969    }
970  }
971  // change cost and bounds on incoming if primal
972  if (algorithm_>0)
973    nonLinearCost_->setOne(sequenceIn_,solution_[sequenceIn_]); 
974  objectiveValue_ += objectiveChange;
975  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
976    <<numberIterations_<<objectiveValue_
977    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
978    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
979  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
980  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
981  handler_->message()<<OsiMessageEol;
982
983  if (numberIterations_>=maximumIterations_)
984    return 2;
985  // only time to re-factorize if one before real time
986  // this is so user won't be surprised that maximumPivots has exact meaning
987  if (factorization_->pivots()==factorization_->maximumPivots()) {
988    return 1;
989  } else {
990    if (forceFactorization_>0&&
991        factorization_->pivots()==forceFactorization_) {
992      // relax
993      forceFactorization_ *= 2;
994      if (forceFactorization_>factorization_->maximumPivots())
995        forceFactorization_ = -1; //off
996      return 1;
997    } else {
998      // carry on iterating
999      return 0;
1000    }
1001  }
1002}
1003// Copy constructor.
1004ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
1005  ClpModel(rhs),
1006  columnPrimalInfeasibility_(0.0),
1007  columnPrimalSequence_(-2),
1008  rowPrimalInfeasibility_(0.0),
1009  rowPrimalSequence_(-2), 
1010  columnDualInfeasibility_(0.0),
1011  columnDualSequence_(-2),
1012  rowDualInfeasibility_(0.0),
1013  rowDualSequence_(-2),
1014  primalToleranceToGetOptimal_(-1.0),
1015  remainingDualInfeasibility_(0.0),
1016  largeValue_(1.0e15),
1017  largestPrimalError_(0.0),
1018  largestDualError_(0.0),
1019  largestSolutionError_(0.0),
1020  dualBound_(1.0e7),
1021  lower_(NULL),
1022  rowLowerWork_(NULL),
1023  columnLowerWork_(NULL),
1024  upper_(NULL),
1025  rowUpperWork_(NULL),
1026  columnUpperWork_(NULL),
1027  cost_(NULL),
1028  rowObjectiveWork_(NULL),
1029  objectiveWork_(NULL),
1030  alpha_(0.0),
1031  theta_(0.0),
1032  lowerIn_(0.0),
1033  valueIn_(0.0),
1034  upperIn_(0.0),
1035  dualIn_(0.0),
1036  sequenceIn_(-1),
1037  directionIn_(-1),
1038  lowerOut_(-1),
1039  valueOut_(-1),
1040  upperOut_(-1),
1041  dualOut_(-1),
1042  sequenceOut_(-1),
1043  directionOut_(-1),
1044  pivotRow_(-1),
1045  status_(NULL),
1046  dj_(NULL),
1047  rowReducedCost_(NULL),
1048  reducedCostWork_(NULL),
1049  solution_(NULL),
1050  rowActivityWork_(NULL),
1051  columnActivityWork_(NULL),
1052  dualTolerance_(0.0),
1053  primalTolerance_(0.0),
1054  sumDualInfeasibilities_(0.0),
1055  numberDualInfeasibilities_(0),
1056  numberDualInfeasibilitiesWithoutFree_(0),
1057  sumPrimalInfeasibilities_(0.0),
1058  numberPrimalInfeasibilities_(0),
1059  pivotVariable_(NULL),
1060  factorization_(NULL),
1061  numberRefinements_(0),
1062  rowScale_(NULL),
1063  savedSolution_(NULL),
1064  columnScale_(NULL),
1065  scalingFlag_(0),
1066  numberTimesOptimal_(0),
1067  changeMade_(1),
1068  algorithm_(0),
1069  forceFactorization_(-1),
1070  perturbation_(100),
1071  infeasibilityCost_(1.0e7),
1072  nonLinearCost_(NULL),
1073  specialOptions_(0),
1074  lastBadIteration_(-999999)
1075{
1076  int i;
1077  for (i=0;i<6;i++) {
1078    rowArray_[i]=NULL;
1079    columnArray_[i]=NULL;
1080  }
1081  saveStatus_=NULL;
1082  factorization_ = NULL;
1083  dualRowPivot_ = NULL;
1084  primalColumnPivot_ = NULL;
1085  gutsOfDelete(0);
1086  gutsOfCopy(rhs);
1087}
1088// Copy constructor from model
1089ClpSimplex::ClpSimplex(const ClpModel &rhs) :
1090  ClpModel(rhs),
1091  columnPrimalInfeasibility_(0.0),
1092  columnPrimalSequence_(-2),
1093  rowPrimalInfeasibility_(0.0),
1094  rowPrimalSequence_(-2), 
1095  columnDualInfeasibility_(0.0),
1096  columnDualSequence_(-2),
1097  rowDualInfeasibility_(0.0),
1098  rowDualSequence_(-2),
1099  primalToleranceToGetOptimal_(-1.0),
1100  remainingDualInfeasibility_(0.0),
1101  largeValue_(1.0e15),
1102  largestPrimalError_(0.0),
1103  largestDualError_(0.0),
1104  largestSolutionError_(0.0),
1105  dualBound_(1.0e7),
1106  lower_(NULL),
1107  rowLowerWork_(NULL),
1108  columnLowerWork_(NULL),
1109  upper_(NULL),
1110  rowUpperWork_(NULL),
1111  columnUpperWork_(NULL),
1112  cost_(NULL),
1113  rowObjectiveWork_(NULL),
1114  objectiveWork_(NULL),
1115  alpha_(0.0),
1116  theta_(0.0),
1117  lowerIn_(0.0),
1118  valueIn_(0.0),
1119  upperIn_(0.0),
1120  dualIn_(0.0),
1121  sequenceIn_(-1),
1122  directionIn_(-1),
1123  lowerOut_(-1),
1124  valueOut_(-1),
1125  upperOut_(-1),
1126  dualOut_(-1),
1127  sequenceOut_(-1),
1128  directionOut_(-1),
1129  pivotRow_(-1),
1130  status_(NULL),
1131  dj_(NULL),
1132  rowReducedCost_(NULL),
1133  reducedCostWork_(NULL),
1134  solution_(NULL),
1135  rowActivityWork_(NULL),
1136  columnActivityWork_(NULL),
1137  dualTolerance_(0.0),
1138  primalTolerance_(0.0),
1139  sumDualInfeasibilities_(0.0),
1140  numberDualInfeasibilities_(0),
1141  numberDualInfeasibilitiesWithoutFree_(0),
1142  sumPrimalInfeasibilities_(0.0),
1143  numberPrimalInfeasibilities_(0),
1144  pivotVariable_(NULL),
1145  factorization_(NULL),
1146  numberRefinements_(0),
1147  rowScale_(NULL),
1148  savedSolution_(NULL),
1149  columnScale_(NULL),
1150  scalingFlag_(0),
1151  numberTimesOptimal_(0),
1152  changeMade_(1),
1153  algorithm_(0),
1154  forceFactorization_(-1),
1155  perturbation_(100),
1156  infeasibilityCost_(1.0e7),
1157  nonLinearCost_(NULL),
1158  specialOptions_(0),
1159  lastBadIteration_(-999999)
1160{
1161  int i;
1162  for (i=0;i<6;i++) {
1163    rowArray_[i]=NULL;
1164    columnArray_[i]=NULL;
1165  }
1166  saveStatus_=NULL;
1167  // get an empty factorization so we can set tolerances etc
1168  factorization_ = new ClpFactorization();
1169  // say Steepest pricing
1170  dualRowPivot_ = new ClpDualRowSteepest();
1171  // say Dantzig pricing
1172  primalColumnPivot_ = new ClpPrimalColumnDantzig();
1173 
1174}
1175// Assignment operator. This copies the data
1176ClpSimplex & 
1177ClpSimplex::operator=(const ClpSimplex & rhs)
1178{
1179  if (this != &rhs) {
1180    gutsOfDelete(0);
1181    ClpModel::operator=(rhs);
1182    gutsOfCopy(rhs);
1183  }
1184  return *this;
1185}
1186void 
1187ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1188{
1189  lower_ = copyOfArray(rhs.lower_,numberColumns_+numberRows_);
1190  rowLowerWork_ = lower_+numberColumns_;
1191  columnLowerWork_ = lower_;
1192  upper_ = copyOfArray(rhs.upper_,numberColumns_+numberRows_);
1193  rowUpperWork_ = upper_+numberColumns_;
1194  columnUpperWork_ = upper_;
1195  cost_ = copyOfArray(rhs.cost_,(numberColumns_+numberRows_));
1196  objectiveWork_ = cost_;
1197  rowObjectiveWork_ = cost_+numberColumns_;
1198  dj_ = copyOfArray(rhs.dj_,numberRows_+numberColumns_);
1199  if (dj_) {
1200    reducedCostWork_ = dj_;
1201    rowReducedCost_ = dj_+numberColumns_;
1202  }
1203  solution_ = copyOfArray(rhs.solution_,numberRows_+numberColumns_);
1204  if (solution_) {
1205    columnActivityWork_ = solution_;
1206    rowActivityWork_ = solution_+numberColumns_;
1207  }
1208  if (rhs.pivotVariable_) {
1209    pivotVariable_ = new int[numberRows_];
1210    CoinDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
1211  } else {
1212    pivotVariable_=NULL;
1213  }
1214  if (rhs.factorization_) {
1215    factorization_ = new ClpFactorization(*rhs.factorization_);
1216  } else {
1217    factorization_=NULL;
1218  }
1219  rowScale_ = copyOfArray(rhs.rowScale_,numberRows_);
1220  savedSolution_ = copyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1221  columnScale_ = copyOfArray(rhs.columnScale_,numberColumns_);
1222  int i;
1223  for (i=0;i<6;i++) {
1224    rowArray_[i]=NULL;
1225    if (rhs.rowArray_[i]) 
1226      rowArray_[i] = new OsiIndexedVector(*rhs.rowArray_[i]);
1227    columnArray_[i]=NULL;
1228    if (rhs.columnArray_[i]) 
1229      columnArray_[i] = new OsiIndexedVector(*rhs.columnArray_[i]);
1230  }
1231  if (rhs.status_) {
1232    status_ = copyOfArray( rhs.status_,numberColumns_+numberRows_);
1233  }
1234  if (rhs.saveStatus_) {
1235    saveStatus_ = copyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1236  }
1237  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1238  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1239  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1240  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1241  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1242  columnDualSequence_ = rhs.columnDualSequence_;
1243  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1244  rowDualSequence_ = rhs.rowDualSequence_;
1245  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1246  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1247  largeValue_ = rhs.largeValue_;
1248  largestPrimalError_ = rhs.largestPrimalError_;
1249  largestDualError_ = rhs.largestDualError_;
1250  largestSolutionError_ = rhs.largestSolutionError_;
1251  dualBound_ = rhs.dualBound_;
1252  alpha_ = rhs.alpha_;
1253  theta_ = rhs.theta_;
1254  lowerIn_ = rhs.lowerIn_;
1255  valueIn_ = rhs.valueIn_;
1256  upperIn_ = rhs.upperIn_;
1257  dualIn_ = rhs.dualIn_;
1258  sequenceIn_ = rhs.sequenceIn_;
1259  directionIn_ = rhs.directionIn_;
1260  lowerOut_ = rhs.lowerOut_;
1261  valueOut_ = rhs.valueOut_;
1262  upperOut_ = rhs.upperOut_;
1263  dualOut_ = rhs.dualOut_;
1264  sequenceOut_ = rhs.sequenceOut_;
1265  directionOut_ = rhs.directionOut_;
1266  pivotRow_ = rhs.pivotRow_;
1267  numberRefinements_ = rhs.numberRefinements_;
1268  dualTolerance_ = rhs.dualTolerance_;
1269  primalTolerance_ = rhs.primalTolerance_;
1270  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1271  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1272  numberDualInfeasibilitiesWithoutFree_ = 
1273    rhs.numberDualInfeasibilitiesWithoutFree_;
1274  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1275  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1276  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1277  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1278  scalingFlag_ = rhs.scalingFlag_;
1279  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1280  changeMade_ = rhs.changeMade_;
1281  algorithm_ = rhs.algorithm_;
1282  forceFactorization_ = rhs.forceFactorization_;
1283  perturbation_ = rhs.perturbation_;
1284  infeasibilityCost_ = rhs.infeasibilityCost_;
1285  specialOptions_ = rhs.specialOptions_;
1286  lastBadIteration_ = rhs.lastBadIteration_;
1287  if (rhs.nonLinearCost_!=NULL)
1288    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1289}
1290// type == 0 do everything
1291void 
1292ClpSimplex::gutsOfDelete(int type)
1293{
1294  delete [] lower_;
1295  lower_=NULL;
1296  rowLowerWork_=NULL;
1297  columnLowerWork_=NULL;
1298  delete [] upper_;
1299  upper_=NULL;
1300  rowUpperWork_=NULL;
1301  columnUpperWork_=NULL;
1302  delete [] cost_;
1303  cost_=NULL;
1304  objectiveWork_=NULL;
1305  rowObjectiveWork_=NULL;
1306  delete [] dj_;
1307  dj_=NULL;
1308  reducedCostWork_=NULL;
1309  rowReducedCost_=NULL;
1310  delete [] solution_;
1311  solution_=NULL;
1312  rowActivityWork_=NULL;
1313  columnActivityWork_=NULL;
1314  delete [] rowScale_;
1315  rowScale_ = NULL;
1316  delete [] savedSolution_;
1317  savedSolution_ = NULL;
1318  delete [] columnScale_;
1319  columnScale_ = NULL;
1320  delete nonLinearCost_;
1321  nonLinearCost_ = NULL;
1322  int i;
1323  for (i=0;i<6;i++) {
1324    delete rowArray_[i];
1325    rowArray_[i]=NULL;
1326    delete columnArray_[i];
1327    columnArray_[i]=NULL;
1328  }
1329  delete rowCopy_;
1330  rowCopy_=NULL;
1331  if (!type) {
1332    // delete everything
1333    delete [] pivotVariable_;
1334    pivotVariable_=NULL;
1335    delete factorization_;
1336    factorization_ = NULL;
1337    delete dualRowPivot_;
1338    dualRowPivot_ = NULL;
1339    delete primalColumnPivot_;
1340    primalColumnPivot_ = NULL;
1341    delete [] saveStatus_;
1342    saveStatus_=NULL;
1343    delete status_;
1344    status_=NULL;
1345  }
1346}
1347// This sets largest infeasibility and most infeasible
1348void 
1349ClpSimplex::checkPrimalSolution(const double * rowActivities,
1350                                        const double * columnActivities)
1351{
1352  double * solution;
1353  int iRow,iColumn;
1354
1355  objectiveValue_ = 0.0;
1356  // now look at primal solution
1357  columnPrimalInfeasibility_=0.0;
1358  columnPrimalSequence_=-1;
1359  rowPrimalInfeasibility_=0.0;
1360  rowPrimalSequence_=-1;
1361  largestSolutionError_=0.0;
1362  solution = rowActivityWork_;
1363  sumPrimalInfeasibilities_=0.0;
1364  numberPrimalInfeasibilities_=0;
1365  for (iRow=0;iRow<numberRows_;iRow++) {
1366    double infeasibility=0.0;
1367    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1368    if (solution[iRow]>rowUpperWork_[iRow]) {
1369      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1370    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1371      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1372    }
1373    if (infeasibility>primalTolerance_) {
1374      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1375      numberPrimalInfeasibilities_ ++;
1376    }
1377    if (infeasibility>rowPrimalInfeasibility_) {
1378      rowPrimalInfeasibility_=infeasibility;
1379      rowPrimalSequence_=iRow;
1380    }
1381    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1382    if (infeasibility>largestSolutionError_)
1383      largestSolutionError_=infeasibility;
1384    if (rowLowerWork_[iRow]!=rowUpperWork_[iRow])
1385      clearFixed(iRow+numberColumns_);
1386    else
1387      setFixed(iRow+numberColumns_);
1388  }
1389  solution = columnActivityWork_;
1390  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1391    double infeasibility=0.0;
1392    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1393    if (solution[iColumn]>columnUpperWork_[iColumn]) {
1394      infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1395    } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1396      infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1397    }
1398    if (infeasibility>columnPrimalInfeasibility_) {
1399      columnPrimalInfeasibility_=infeasibility;
1400      columnPrimalSequence_=iColumn;
1401    }
1402    if (infeasibility>primalTolerance_) {
1403      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1404      numberPrimalInfeasibilities_ ++;
1405    }
1406    infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1407    if (infeasibility>largestSolutionError_)
1408      largestSolutionError_=infeasibility;
1409    if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
1410      clearFixed(iColumn);
1411    else
1412      setFixed(iColumn);
1413  }
1414}
1415void 
1416ClpSimplex::checkDualSolution()
1417{
1418
1419  double * solution;
1420  int iRow,iColumn;
1421  sumDualInfeasibilities_=0.0;
1422  numberDualInfeasibilities_=0;
1423  numberDualInfeasibilitiesWithoutFree_=0;
1424  columnDualInfeasibility_=0.0;
1425  columnDualSequence_=-1;
1426  rowDualInfeasibility_=0.0;
1427  rowDualSequence_=-1;
1428  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
1429                                   columnPrimalInfeasibility_);
1430  remainingDualInfeasibility_=0.0;
1431  solution = rowActivityWork_;
1432  double dualTolerance=dualTolerance_;
1433  if (algorithm_>0) {
1434    // primal
1435    // we can't really trust infeasibilities if there is dual error
1436    if (largestDualError_>1.0e-6)
1437      dualTolerance *= largestDualError_/1.0e-6;
1438  }
1439  for (iRow=0;iRow<numberRows_;iRow++) {
1440    if (getRowStatus(iRow) != ClpSimplex::basic) {
1441      // not basic
1442      double value = rowReducedCost_[iRow];
1443      double distance;
1444      distance = rowUpperWork_[iRow]-solution[iRow];
1445      if (distance>primalTolerance_) {
1446        // should not be negative
1447        if (value<0.0) {
1448          value = - value;
1449          if (value>rowDualInfeasibility_) {
1450            rowDualInfeasibility_=value;
1451            rowDualSequence_=iRow;
1452          }
1453          if (value>dualTolerance) {
1454            sumDualInfeasibilities_ += value-dualTolerance;
1455            numberDualInfeasibilities_ ++;
1456            if (getRowStatus(iRow) != ClpSimplex::isFree) 
1457              numberDualInfeasibilitiesWithoutFree_ ++;
1458            // maybe we can make feasible by increasing tolerance
1459            if (distance<largeValue_) {
1460              if (distance>primalToleranceToGetOptimal_)
1461                primalToleranceToGetOptimal_=distance;
1462            } else {
1463              //gap too big for any tolerance
1464              remainingDualInfeasibility_=
1465                max(remainingDualInfeasibility_,value);
1466            }
1467          }
1468        }
1469      }
1470      distance = solution[iRow] -rowLowerWork_[iRow];
1471      if (distance>primalTolerance_) {
1472        // should not be positive
1473        if (value>0.0) {
1474          if (value>rowDualInfeasibility_) {
1475            rowDualInfeasibility_=value;
1476            rowDualSequence_=iRow;
1477          }
1478          if (value>dualTolerance) {
1479            sumDualInfeasibilities_ += value-dualTolerance;
1480            numberDualInfeasibilities_ ++;
1481            if (getRowStatus(iRow) != ClpSimplex::isFree) 
1482              numberDualInfeasibilitiesWithoutFree_ ++;
1483            // maybe we can make feasible by increasing tolerance
1484            if (distance<largeValue_&&
1485                distance>primalToleranceToGetOptimal_)
1486              primalToleranceToGetOptimal_=distance;
1487          }
1488        }
1489      }
1490    }
1491  }
1492  solution = columnActivityWork_;
1493  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1494    if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1495      // not basic
1496      double value = reducedCostWork_[iColumn];
1497      double distance;
1498      distance = columnUpperWork_[iColumn]-solution[iColumn];
1499      if (distance>primalTolerance_) {
1500        // should not be negative
1501        if (value<0.0) {
1502          value = - value;
1503          if (value>columnDualInfeasibility_) {
1504            columnDualInfeasibility_=value;
1505            columnDualSequence_=iColumn;
1506          }
1507          if (value>dualTolerance) {
1508            sumDualInfeasibilities_ += value-dualTolerance;
1509            numberDualInfeasibilities_ ++;
1510            if (getColumnStatus(iColumn) != ClpSimplex::isFree) 
1511              numberDualInfeasibilitiesWithoutFree_ ++;
1512            // maybe we can make feasible by increasing tolerance
1513            if (distance<largeValue_) {
1514              if (distance>primalToleranceToGetOptimal_)
1515                primalToleranceToGetOptimal_=distance;
1516            } else {
1517              //gap too big for any tolerance
1518              remainingDualInfeasibility_=
1519                max(remainingDualInfeasibility_,value);
1520            }
1521          }
1522        }
1523      }
1524      distance = solution[iColumn] -columnLowerWork_[iColumn];
1525      if (distance>primalTolerance_) {
1526        // should not be positive
1527        if (value>0.0) {
1528          if (value>columnDualInfeasibility_) {
1529            columnDualInfeasibility_=value;
1530            columnDualSequence_=iColumn;
1531          }
1532          if (value>dualTolerance) {
1533            sumDualInfeasibilities_ += value-dualTolerance;
1534            numberDualInfeasibilities_ ++;
1535            if (getColumnStatus(iColumn) != ClpSimplex::isFree) 
1536              numberDualInfeasibilitiesWithoutFree_ ++;
1537            // maybe we can make feasible by increasing tolerance
1538            if (distance<largeValue_&&
1539                distance>primalToleranceToGetOptimal_)
1540              primalToleranceToGetOptimal_=distance;
1541          }
1542        }
1543      }
1544    }
1545  }
1546}
1547/*
1548  Unpacks one column of the matrix into indexed array
1549*/
1550void 
1551ClpSimplex::unpack(OsiIndexedVector * rowArray)
1552{
1553  rowArray->clear();
1554  if (sequenceIn_>=numberColumns_) {
1555    //slack
1556    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
1557  } else {
1558    // column
1559    matrix_->unpack(this,rowArray,sequenceIn_);
1560  }
1561}
1562void 
1563ClpSimplex::unpack(OsiIndexedVector * rowArray,int sequence)
1564{
1565  rowArray->clear();
1566  if (sequence>=numberColumns_) {
1567    //slack
1568    rowArray->insert(sequence-numberColumns_,-1.0);
1569  } else {
1570    // column
1571    matrix_->unpack(this,rowArray,sequence);
1572  }
1573}
1574void
1575ClpSimplex::createRim(int what,bool makeRowCopy)
1576{
1577  if ((what&16)!=0) {
1578    // move information to work arrays
1579    // row reduced costs
1580    if (!dj_) {
1581      dj_ = new double[numberRows_+numberColumns_];
1582      reducedCostWork_ = dj_;
1583      rowReducedCost_ = dj_+numberColumns_;
1584      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1585    }
1586    if (!solution_) {
1587      solution_ = new double[numberRows_+numberColumns_];
1588      columnActivityWork_ = solution_;
1589      rowActivityWork_ = solution_+numberColumns_;
1590      memcpy(columnActivityWork_,columnActivity_,
1591             numberColumns_*sizeof(double));
1592      memcpy(rowActivityWork_,rowActivity_,
1593             numberRows_*sizeof(double));
1594    }
1595   
1596    if (makeRowCopy) {
1597      delete rowCopy_;
1598      // may return NULL if can't give row copy
1599      rowCopy_ = matrix_->reverseOrderedCopy();
1600    }
1601  }
1602  int i;
1603  if ((what&4)!=0) {
1604    delete [] cost_;
1605    cost_ = new double[numberColumns_+numberRows_];
1606    objectiveWork_ = cost_;
1607    rowObjectiveWork_ = cost_+numberColumns_;
1608    memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
1609    if (rowObjective_)
1610      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1611    else
1612      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1613  }
1614  if ((what&1)!=0) {
1615    delete [] lower_;
1616    delete [] upper_;
1617    lower_ = new double[numberColumns_+numberRows_];
1618    upper_ = new double[numberColumns_+numberRows_];
1619    rowLowerWork_ = lower_+numberColumns_;
1620    columnLowerWork_ = lower_;
1621    rowUpperWork_ = upper_+numberColumns_;
1622    columnUpperWork_ = upper_;
1623    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1624    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1625    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1626    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1627    // clean up any mismatches on infinity
1628    for (i=0;i<numberColumns_;i++) {
1629      if (columnLowerWork_[i]<-1.0e30)
1630        columnLowerWork_[i] = -DBL_MAX;
1631      if (columnUpperWork_[i]>1.0e30)
1632        columnUpperWork_[i] = DBL_MAX;
1633    }
1634    // clean up any mismatches on infinity
1635    for (i=0;i<numberRows_;i++) {
1636      if (rowLowerWork_[i]<-1.0e30)
1637        rowLowerWork_[i] = -DBL_MAX;
1638      if (rowUpperWork_[i]>1.0e30)
1639        rowUpperWork_[i] = DBL_MAX;
1640    }
1641  }
1642  // do scaling if needed
1643  if (scalingFlag_>0&&!rowScale_) {
1644    if (matrix_->scale(this))
1645      scalingFlag_=-scalingFlag_; // not scaled after all
1646  }
1647  if ((what&4)!=0) {
1648    double direction = optimizationDirection_;
1649    // direction is actually scale out not scale in
1650    if (direction)
1651      direction = 1.0/direction;
1652    // but also scale by scale factors
1653    // not really sure about row scaling
1654    if (!rowScale_) {
1655      if (direction!=1.0) {
1656        for (i=0;i<numberRows_;i++)
1657          rowObjectiveWork_[i] *= direction;
1658        for (i=0;i<numberColumns_;i++)
1659          objectiveWork_[i] *= direction;
1660      }
1661    } else {
1662      for (i=0;i<numberRows_;i++)
1663        rowObjectiveWork_[i] *= direction/rowScale_[i];
1664      for (i=0;i<numberColumns_;i++)
1665        objectiveWork_[i] *= direction*columnScale_[i];
1666    }
1667  }
1668  if ((what&1)!=0&&rowScale_) {
1669    for (i=0;i<numberColumns_;i++) {
1670      double multiplier = 1.0/columnScale_[i];
1671      if (columnLowerWork_[i]>-1.0e50)
1672        columnLowerWork_[i] *= multiplier;
1673      if (columnUpperWork_[i]<1.0e50)
1674        columnUpperWork_[i] *= multiplier;
1675    }
1676    for (i=0;i<numberRows_;i++) {
1677      if (rowLowerWork_[i]>-1.0e50)
1678        rowLowerWork_[i] *= rowScale_[i];
1679      if (rowUpperWork_[i]<1.0e50)
1680        rowUpperWork_[i] *= rowScale_[i];
1681    }
1682  }
1683  if ((what&8)!=0&&rowScale_) {
1684    // on entry
1685    for (i=0;i<numberColumns_;i++) {
1686      columnActivityWork_[i] /= columnScale_[i];
1687      reducedCostWork_[i] *= columnScale_[i];
1688    }
1689    for (i=0;i<numberRows_;i++) {
1690      rowActivityWork_[i] *= rowScale_[i];
1691      dual_[i] /= rowScale_[i];
1692    }
1693  }
1694 
1695  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1696  // maybe we need to move scales to SimplexModel for factorization?
1697  if ((what&8)!=0&&!pivotVariable_) {
1698    pivotVariable_=new int[numberRows_];
1699  }
1700
1701  if ((what&16)!=0&&!rowArray_[2]) {
1702    // get some arrays
1703    int iRow,iColumn;
1704    // these are "indexed" arrays so we always know where nonzeros are
1705    /**********************************************************
1706      rowArray_[3] is long enough for columns as well
1707    *********************************************************/
1708    for (iRow=0;iRow<4;iRow++) {
1709      if (!rowArray_[iRow]) {
1710        rowArray_[iRow]=new OsiIndexedVector();
1711        int length =numberRows_+factorization_->maximumPivots();
1712        if (iRow==3)
1713          length = max(length,numberColumns_);
1714        rowArray_[iRow]->reserve(length);
1715      }
1716    }
1717   
1718    for (iColumn=0;iColumn<2;iColumn++) {
1719      if (!columnArray_[iColumn]) {
1720        columnArray_[iColumn]=new OsiIndexedVector();
1721        columnArray_[iColumn]->reserve(numberColumns_);
1722      }
1723    }   
1724  }
1725
1726}
1727void
1728ClpSimplex::deleteRim()
1729{
1730  int i;
1731  if (problemStatus_!=1&&problemStatus_!=2) {
1732    delete [] ray_;
1733    ray_=NULL;
1734  }
1735  if (rowScale_) {
1736    for (i=0;i<numberColumns_;i++) {
1737      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1738      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1739    }
1740    for (i=0;i<numberRows_;i++) {
1741      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1742      dual_[i] *= rowScale_[i];
1743    }
1744    if (problemStatus_==2) {
1745      for (i=0;i<numberColumns_;i++) {
1746        ray_[i] *= columnScale_[i];
1747      }
1748    } else if (problemStatus_==1) {
1749      for (i=0;i<numberRows_;i++) {
1750        ray_[i] *= rowScale_[i];
1751      }
1752    }
1753  } else {
1754    for (i=0;i<numberColumns_;i++) {
1755      columnActivity_[i] = columnActivityWork_[i];
1756      reducedCost_[i] = reducedCostWork_[i];
1757    }
1758    for (i=0;i<numberRows_;i++) {
1759      rowActivity_[i] = rowActivityWork_[i];
1760    }
1761  }
1762  // direction may have been modified by scaling - clean up
1763  if (optimizationDirection_>0.0)
1764    optimizationDirection_ = 1.0;
1765  else if (optimizationDirection_<0.0)
1766    optimizationDirection_ = -1.0;
1767  // scaling may have been turned off
1768  scalingFlag_ = abs(scalingFlag_);
1769  gutsOfDelete(1);
1770}
1771void 
1772ClpSimplex::setDualBound(double value)
1773{
1774  if (value>0.0)
1775    dualBound_=value;
1776}
1777void 
1778ClpSimplex::setInfeasibilityCost(double value)
1779{
1780  if (value>0.0)
1781    infeasibilityCost_=value;
1782}
1783void ClpSimplex::setnumberRefinements( int value) 
1784{
1785  if (value>=0&&value<10)
1786    numberRefinements_=value;
1787}
1788// Sets row pivot choice algorithm in dual
1789void 
1790ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
1791{
1792  delete dualRowPivot_;
1793  dualRowPivot_ = choice.clone(true);
1794}
1795// Sets row pivot choice algorithm in dual
1796void 
1797ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
1798{
1799  delete primalColumnPivot_;
1800  primalColumnPivot_ = choice.clone(true);
1801}
1802// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
1803void 
1804ClpSimplex::scaling(int mode)
1805{
1806  if (mode>0&&mode<3) {
1807    scalingFlag_=mode;
1808  } else if (!mode) {
1809    scalingFlag_=0;
1810    delete [] rowScale_;
1811    rowScale_ = NULL;
1812    delete [] columnScale_;
1813    columnScale_ = NULL;
1814  }
1815}
1816// Sets up basis
1817void 
1818ClpSimplex::setBasis ( const OsiWarmStartBasis & basis)
1819{
1820  // transform basis to status arrays
1821  int iRow,iColumn;
1822  if (!status_) {
1823    /*
1824      get status arrays
1825      OsiWarmStartBasis would seem to have overheads and we will need
1826      extra bits anyway.
1827    */
1828    status_ = new unsigned char [numberColumns_+numberRows_];
1829    memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
1830  }
1831  OsiWarmStartBasis basis2 = basis;
1832  // resize if necessary
1833  basis2.resize(numberRows_,numberColumns_);
1834  // move status
1835  for (iRow=0;iRow<numberRows_;iRow++) {
1836    setRowStatus(iRow,
1837                 (Status) basis2.getArtifStatus(iRow));
1838  }
1839  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1840    setColumnStatus(iColumn,
1841                    (Status) basis2.getStructStatus(iColumn));
1842  }
1843}
1844// Passes in factorization
1845void 
1846ClpSimplex::setFactorization( ClpFactorization & factorization)
1847{
1848  delete factorization_;
1849  factorization_= new ClpFactorization(factorization);
1850}
1851// Warm start
1852OsiWarmStartBasis 
1853ClpSimplex::getBasis() const
1854{
1855  int iRow,iColumn;
1856  OsiWarmStartBasis basis;
1857  basis.setSize(numberColumns_,numberRows_);
1858
1859  if(status_) {
1860    for (iRow=0;iRow<numberRows_;iRow++) {
1861      basis.setArtifStatus(iRow,
1862                           (OsiWarmStartBasis::Status) getRowStatus(iRow));
1863    }
1864    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1865      basis.setStructStatus(iColumn,
1866                       (OsiWarmStartBasis::Status) getColumnStatus(iColumn));
1867    }
1868  }
1869  return basis;
1870}
1871void 
1872ClpSimplex::times(double scalar,
1873                  const double * x, double * y) const
1874{
1875  if (rowScale_)
1876    matrix_->times(scalar,x,y,rowScale_,columnScale_);
1877  else
1878    matrix_->times(scalar,x,y);
1879}
1880void 
1881ClpSimplex::transposeTimes(double scalar,
1882                           const double * x, double * y) const 
1883{
1884  if (rowScale_)
1885    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
1886  else
1887    matrix_->transposeTimes(scalar,x,y);
1888}
1889/* Perturbation:
1890   -50 to +50 - perturb by this power of ten (-6 sounds good)
1891   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1892   101 - we are perturbed
1893   102 - don't try perturbing again
1894   default is 100
1895*/
1896void 
1897ClpSimplex::setPerturbation(int value)
1898{
1899  if(value<=100&&value >=-50) {
1900    perturbation_=value;
1901  } 
1902}
1903// Sparsity on or off
1904bool 
1905ClpSimplex::sparseFactorization() const
1906{
1907  return factorization_->sparseThreshold();
1908}
1909void 
1910ClpSimplex::setSparseFactorization(bool value)
1911{
1912  if (value) {
1913    if (!factorization_->sparseThreshold())
1914      factorization_->goSparse();
1915  } else {
1916    factorization_->sparseThreshold(0);
1917  }
1918}
1919/* Tightens primal bounds to make dual faster.  Unless
1920   fixed, bounds are slightly looser than they could be.
1921   This is to make dual go faster and is probably not needed
1922   with a presolve.  Returns non-zero if problem infeasible
1923*/
1924int 
1925ClpSimplex::tightenPrimalBounds()
1926{
1927 
1928  // Get a row copy in standard format
1929  OsiPackedMatrix copy;
1930  copy.reverseOrderedCopyOf(*matrix());
1931  // get matrix data pointers
1932  const int * column = copy.getIndices();
1933  const int * rowStart = copy.getVectorStarts();
1934  const int * rowLength = copy.getVectorLengths(); 
1935  const double * element = copy.getElements();
1936  int numberChanged=1,iPass=0;
1937  double large = largeValue(); // treat bounds > this as infinite
1938  int numberInfeasible=0;
1939  int totalTightened = 0;
1940
1941  double tolerance = primalTolerance();
1942
1943
1944  // Save column bounds
1945  double * saveLower = new double [numberColumns_];
1946  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
1947  double * saveUpper = new double [numberColumns_];
1948  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
1949#define MAXPASS 10 
1950
1951  // Loop round seeing if we can tighten bounds
1952  // Would be faster to have a stack of possible rows
1953  // and we put altered rows back on stack
1954
1955  int iRow, iColumn;
1956
1957  while(numberChanged) {
1958
1959    numberChanged = 0; // Bounds tightened this pass
1960   
1961    if (iPass==MAXPASS) break;
1962    iPass++;
1963   
1964    for (iRow = 0; iRow < numberRows_; iRow++) {
1965
1966      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
1967
1968        // possible row
1969        int infiniteUpper = 0;
1970        int infiniteLower = 0;
1971        double maximumUp = 0.0;
1972        double maximumDown = 0.0;
1973        double newBound;
1974        int rStart = rowStart[iRow];
1975        int rEnd = rowStart[iRow]+rowLength[iRow];
1976        int j;
1977
1978        // Compute possible lower and upper ranges
1979
1980        for (j = rStart; j < rEnd; ++j) {
1981          double value=element[j];
1982          iColumn = column[j];
1983          if (value > 0.0) {
1984            if (columnUpper_[iColumn] >= large) {
1985              maximumUp = 1e31;
1986              ++infiniteUpper;
1987            } else {
1988              maximumUp += columnUpper_[iColumn] * value;
1989            }
1990            if (columnLower_[iColumn] <= -large) {
1991              maximumDown = -1e31;
1992              ++infiniteLower;
1993            } else {
1994              maximumDown += columnLower_[iColumn] * value;
1995            }
1996          } else if (value<0.0) {
1997            if (columnUpper_[iColumn] >= large) {
1998              maximumDown = -1e31;
1999              ++infiniteLower;
2000            } else {
2001              maximumDown += columnUpper_[iColumn] * value;
2002            }
2003            if (columnLower_[iColumn] <= -large) {
2004              maximumUp = 1e31;
2005              ++infiniteUpper;
2006            } else {
2007              maximumUp += columnLower_[iColumn] * value;
2008            }
2009          }
2010        }
2011        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2012            maximumDown >= rowLower_[iRow] - tolerance) {
2013
2014          // Row is redundant - make totally free
2015          rowLower_[iRow]=-DBL_MAX;
2016          rowUpper_[iRow]=DBL_MAX;
2017
2018        } else {
2019          if (maximumUp < rowLower_[iRow] -tolerance ||
2020              maximumDown > rowUpper_[iRow]+tolerance) {
2021            // problem is infeasible - exit at once
2022            numberInfeasible++;
2023            break;
2024          }
2025
2026          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2027            for (j = rStart; j < rEnd; ++j) {
2028              double value=element[j];
2029              iColumn = column[j];
2030              if (value > 0.0) {
2031                if (columnUpper_[iColumn] < large) {
2032                  newBound = columnUpper_[iColumn] + 
2033                    (rowLower_[iRow] - maximumUp) / value;
2034                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2035                    // Tighten the lower bound
2036
2037                    columnLower_[iColumn] = newBound;
2038                    ++numberChanged;
2039
2040                    // check infeasible (relaxed)
2041                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2042                        -100.0*tolerance) 
2043                      numberInfeasible++;
2044                    infiniteLower=1; // skip looking at other bound
2045                  }
2046                }
2047              } else {
2048                if (columnLower_[iColumn] > -large) {
2049                  newBound = columnLower_[iColumn] + 
2050                    (rowLower_[iRow] - maximumUp) / value;
2051                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2052                    // Tighten the upper bound
2053
2054                    columnUpper_[iColumn] = newBound;
2055                    ++numberChanged;
2056
2057                    // check infeasible (relaxed)
2058                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2059                        -100.0*tolerance) 
2060                      numberInfeasible++;
2061                    infiniteLower=1; // skip looking at other bound
2062                  }
2063                }
2064              }
2065            }
2066          }
2067         
2068          // Try other way
2069          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2070            for (j = rStart; j < rEnd; ++j) {
2071              double value=element[j];
2072              iColumn = column[j];
2073              if (value < 0.0) {
2074                if (columnUpper_[iColumn] < large) {
2075                  newBound = columnUpper_[iColumn] + 
2076                    (rowUpper_[iRow] - maximumDown) / value;
2077                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2078                    // Tighten the lower bound
2079
2080                    columnLower_[iColumn] = newBound;
2081                    ++numberChanged;
2082
2083                    // check infeasible (relaxed)
2084                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2085                        -100.0*tolerance) 
2086                      numberInfeasible++;
2087                  }
2088                } 
2089              } else {
2090                if (columnLower_[iColumn] > -large) {
2091                  newBound = columnLower_[iColumn] + 
2092                    (rowUpper_[iRow] - maximumDown) / value;
2093                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2094                    // Tighten the upper bound
2095
2096                    columnUpper_[iColumn] = newBound;
2097                    ++numberChanged;
2098
2099                    // check infeasible (relaxed)
2100                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2101                        -100.0*tolerance) 
2102                      numberInfeasible++;
2103                  }
2104                }
2105              }
2106            }
2107          }
2108        }
2109      }
2110    }
2111    totalTightened += numberChanged;
2112    if (numberInfeasible) break;
2113  }
2114  if (!numberInfeasible) {
2115    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2116      <<totalTightened
2117      <<OsiMessageEol;
2118    // Set bounds slightly loose
2119    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2120      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2121        // fix
2122        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2123          columnLower_[iColumn]=columnUpper_[iColumn];
2124        else
2125          columnUpper_[iColumn]=columnLower_[iColumn];
2126      } else {
2127        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2128          // relax a bit
2129          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2130                                      saveUpper[iColumn]);
2131        }
2132        if (columnLower_[iColumn]>saveLower[iColumn]) {
2133          // relax a bit
2134          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2135                                      saveLower[iColumn]);
2136        }
2137      }
2138    }
2139  } else {
2140    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2141      <<numberInfeasible
2142      <<OsiMessageEol;
2143    // restore column bounds
2144    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2145    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2146  }
2147  delete [] saveLower;
2148  delete [] saveUpper;
2149  return (numberInfeasible);
2150}
2151// dual
2152#include "ClpSimplexDual.hpp"
2153int ClpSimplex::dual ( )
2154{
2155  return ((ClpSimplexDual *) this)->dual();
2156}
2157// primal
2158#include "ClpSimplexPrimal.hpp"
2159int ClpSimplex::primal (int ifValuesPass )
2160{
2161  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2162}
Note: See TracBrowser for help on using the repository browser.