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

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

Changes to make more reliable if problem size changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 84.7 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 "ClpPackedMatrix.hpp"
17#include "CoinIndexedVector.hpp"
18#include "ClpDualRowDantzig.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.0e6),
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_(1),
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  numberFake_(0),
130  progressFlag_(0),
131  sumOfRelaxedDualInfeasibilities_(0.0),
132  sumOfRelaxedPrimalInfeasibilities_(0.0)
133
134{
135  int i;
136  for (i=0;i<6;i++) {
137    rowArray_[i]=NULL;
138    columnArray_[i]=NULL;
139  }
140  saveStatus_=NULL;
141  // get an empty factorization so we can set tolerances etc
142  factorization_ = new ClpFactorization();
143  // Say sparse
144  factorization_->sparseThreshold(1);
145  // say Steepest pricing
146  dualRowPivot_ = new ClpDualRowSteepest();
147  // say Steepest pricing
148  primalColumnPivot_ = new ClpPrimalColumnSteepest();
149 
150}
151
152
153//-----------------------------------------------------------------------------
154
155ClpSimplex::~ClpSimplex ()
156{
157  gutsOfDelete(0);
158}
159//#############################################################################
160void ClpSimplex::setLargeValue( double value) 
161{
162  if (value>0.0&&value<DBL_MAX)
163    largeValue_=value;
164}
165int
166ClpSimplex::gutsOfSolution ( const double * rowActivities,
167                             const double * columnActivities,
168                             bool valuesPass)
169{
170
171
172  // if values pass, save values of basic variables
173  double * save = NULL;
174  double oldValue=0.0;
175  if (valuesPass) {
176    assert(algorithm_>0); // only primal at present
177    assert(nonLinearCost_);
178    int iRow;
179    checkPrimalSolution( rowActivities, columnActivities);
180    // get correct bounds on all variables
181    nonLinearCost_->checkInfeasibilities();
182    oldValue = nonLinearCost_->largestInfeasibility();
183    save = new double[numberRows_];
184    for (iRow=0;iRow<numberRows_;iRow++) {
185      int iPivot=pivotVariable_[iRow];
186      save[iRow] = solution_[iPivot];
187    }
188  }
189  // do work
190  computePrimals(rowActivities, columnActivities);
191  double objectiveModification = 0.0;
192  if (algorithm_>0&&nonLinearCost_!=NULL) {
193    // primal algorithm
194    // get correct bounds on all variables
195    nonLinearCost_->checkInfeasibilities();
196    objectiveModification += nonLinearCost_->changeInCost();
197    if (nonLinearCost_->numberInfeasibilities())
198      handler_->message(CLP_SIMPLEX_NONLINEAR,messages_)
199        <<nonLinearCost_->changeInCost()
200        <<nonLinearCost_->numberInfeasibilities()
201        <<CoinMessageEol;
202  }
203  if (valuesPass) {
204#ifdef CLP_DEBUG
205    std::cout<<"Largest given infeasibility "<<oldValue
206             <<" now "<<nonLinearCost_->largestInfeasibility()<<std::endl;
207#endif
208    int numberOut=0;
209    if (oldValue<1.0
210        &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
211        largestPrimalError_>1.0e-3) {
212      // throw out up to 1000 structurals
213      int iRow;
214      int * sort = new int[numberRows_];
215      // first put back solution and store difference
216      for (iRow=0;iRow<numberRows_;iRow++) {
217        int iPivot=pivotVariable_[iRow];
218        double difference = fabs(solution_[iPivot]=save[iRow]);
219        solution_[iPivot]=save[iRow];
220        save[iRow]=difference;
221      }
222      for (iRow=0;iRow<numberRows_;iRow++) {
223        int iPivot=pivotVariable_[iRow];
224
225        if (iPivot<numberColumns_) {
226          // column
227          double difference= save[iRow];
228          if (difference>1.0e-4) {
229            sort[numberOut]=iPivot;
230            save[numberOut++]=difference;
231          }
232        }
233      }
234      CoinSort_2(save, save + numberOut, sort,
235                 CoinFirstGreater_2<double, int>());
236      numberOut = min(1000,numberOut);
237      for (iRow=0;iRow<numberOut;iRow++) {
238        int iColumn=sort[iRow];
239        setColumnStatus(iColumn,superBasic);
240
241      }
242      delete [] sort;
243    }
244    delete [] save;
245    if (numberOut)
246      return numberOut;
247  }
248
249  computeDuals();
250
251  // now check solutions
252  checkPrimalSolution( rowActivities, columnActivities);
253  objectiveValue_ += objectiveModification;
254  checkDualSolution();
255  if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2||
256                               largestDualError_>1.0e-2)) 
257    handler_->message(CLP_SIMPLEX_ACCURACY,messages_)
258      <<largestPrimalError_
259      <<largestDualError_
260      <<CoinMessageEol;
261  return 0;
262}
263void
264ClpSimplex::computePrimals ( const double * rowActivities,
265                                     const double * columnActivities)
266{
267
268  //work space
269  CoinIndexedVector  * workSpace = rowArray_[0];
270
271  double * array = new double [numberRows_];
272  double * save = new double [numberRows_];
273  double * previous = new double [numberRows_];
274
275  // accumulate non basic stuff
276  CoinFillN(array,numberRows_,0.0);
277
278  int iRow;
279  // order is this way for scaling
280  // Use whole matrix every time to make it easier for ClpMatrixBase
281  // So zero out basic
282  if (columnActivities!=columnActivityWork_)
283    CoinDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_);
284  if (rowActivities!=rowActivityWork_)
285    CoinDisjointCopyN(rowActivities,numberRows_,rowActivityWork_);
286  for (iRow=0;iRow<numberRows_;iRow++) {
287    int iPivot=pivotVariable_[iRow];
288    solution_[iPivot] = 0.0;
289  }
290  times(-1.0,columnActivityWork_,array);
291
292  for (iRow=0;iRow<numberRows_;iRow++) {
293    array[iRow] += rowActivityWork_[iRow];
294  }
295
296  // Ftran adjusted RHS and iterate to improve accuracy
297  double lastError=DBL_MAX;
298  int iRefine;
299  double * work = workSpace->denseVector();
300  factorization_->updateColumn(workSpace,array);
301
302  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
303    // save in case we want to restore
304    CoinDisjointCopyN ( array, numberRows_ , save);
305
306    // put solution in correct place
307    for (iRow=0;iRow<numberRows_;iRow++) {
308      int iPivot=pivotVariable_[iRow];
309      solution_[iPivot] = array[iRow];
310    }
311    // check Ax == b  (for all)
312    times(-1.0,columnActivityWork_,work);
313    for (iRow=0;iRow<numberRows_;iRow++) {
314      work[iRow] += rowActivityWork_[iRow];
315    }
316
317    largestPrimalError_=0.0;
318    for (iRow=0;iRow<numberRows_;iRow++) {
319      if (fabs(work[iRow])>largestPrimalError_) {
320        largestPrimalError_=fabs(work[iRow]);
321      }
322      //work[iRow] -= save[iRow];
323    }
324    if (largestPrimalError_>=lastError) {
325      // restore
326      double * temp = array;
327      array = previous;
328      previous=temp;
329      break;
330    }
331    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
332      // try and make better
333      // save this
334      double * temp = array;
335      array = previous;
336      previous=temp;
337      double multiplier = 131072.0;
338      for (iRow=0;iRow<numberRows_;iRow++) {
339        array[iRow] = multiplier*work[iRow];
340        work[iRow]=0.0;
341      }
342      lastError=largestPrimalError_;
343      factorization_->updateColumn(workSpace,array);
344      multiplier = 1.0/multiplier;
345      for (iRow=0;iRow<numberRows_;iRow++) {
346        array[iRow] = previous[iRow] + multiplier*array[iRow];
347        work[iRow]=0.0;
348      }
349    }
350  }
351
352  // solution as accurate as we are going to get
353  CoinFillN(work,numberRows_,0.0);
354  // put solution in correct place
355  for (iRow=0;iRow<numberRows_;iRow++) {
356    int iPivot=pivotVariable_[iRow];
357    solution_[iPivot] = array[iRow];
358  }
359  delete [] previous;
360  delete [] array;
361  delete [] save;
362}
363// now dual side
364void
365ClpSimplex::computeDuals()
366{
367  double slackValue = factorization_->slackValue();
368  //work space
369  CoinIndexedVector  * workSpace = rowArray_[0];
370
371  double * array = new double [numberRows_];
372  double * save = new double [numberRows_];
373  double * previous = new double [numberRows_];
374
375  int iRow;
376#ifdef CLP_DEBUG
377  workSpace->checkClear();
378#endif
379  for (iRow=0;iRow<numberRows_;iRow++) {
380    int iPivot=pivotVariable_[iRow];
381    if (iPivot>=numberColumns_) {
382      // slack
383      array[iRow] = rowObjectiveWork_[iPivot-numberColumns_];
384    } else {
385      // column
386      array[iRow]=objectiveWork_[iPivot];
387    }
388  }
389  CoinDisjointCopyN ( array, numberRows_ , save);
390
391  // Btran basic costs and get as accurate as possible
392  double lastError=DBL_MAX;
393  int iRefine;
394  double * work = workSpace->denseVector();
395  factorization_->updateColumnTranspose(workSpace,array);
396  for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
397    // check basic reduced costs zero
398    largestDualError_=0.0;
399    // would be faster to do just for basic but this reduces code
400    CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
401    transposeTimes(-1.0,array,reducedCostWork_);
402    for (iRow=0;iRow<numberRows_;iRow++) {
403      int iPivot=pivotVariable_[iRow];
404      if (iPivot>=numberColumns_) {
405        // slack
406        //work[iRow] += slackValue*array[iPivot-numberColumns_];
407        //work[iRow] += rowObjectiveWork_[iPivot-numberColumns_];
408        work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
409        - slackValue*array[iPivot-numberColumns_];
410      } else {
411        // column
412        work[iRow] = reducedCostWork_[iPivot];
413      }
414      if (fabs(work[iRow])>largestDualError_) {
415        largestDualError_=fabs(work[iRow]);
416      }
417    }
418    if (largestDualError_>=lastError) {
419      // restore
420      double * temp = array;
421      array = previous;
422      previous=temp;
423      break;
424    }
425    if (iRefine<numberRefinements_&&largestDualError_>1.0e-20) {
426      // try and make better
427      // save this
428      double * temp = array;
429      array = previous;
430      previous=temp;
431      double multiplier = 131072.0;
432      for (iRow=0;iRow<numberRows_;iRow++) {
433        array[iRow] = multiplier*work[iRow];
434        work[iRow]=0.0;
435      }
436      lastError=largestDualError_;
437      factorization_->updateColumnTranspose(workSpace,array);
438      multiplier = 1.0/multiplier;
439      for (iRow=0;iRow<numberRows_;iRow++) {
440        array[iRow] = previous[iRow] + multiplier*array[iRow];
441        work[iRow]=0.0;
442      }
443    }
444  }
445  CoinFillN(work,numberRows_,0.0);
446  // now look at dual solution
447  for (iRow=0;iRow<numberRows_;iRow++) {
448    // slack
449    double value = array[iRow];
450    dual_[iRow]=value;
451    value = rowObjectiveWork_[iRow] - value*slackValue;
452    rowReducedCost_[iRow]=value;
453  }
454  CoinDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_);
455  transposeTimes(-1.0,dual_,reducedCostWork_);
456  delete [] previous;
457  delete [] array;
458  delete [] save;
459}
460/* Given an existing factorization computes and checks
461   primal and dual solutions.  Uses input arrays for variables at
462   bounds.  Returns feasibility states */
463int ClpSimplex::getSolution ( const double * rowActivities,
464                               const double * columnActivities)
465{
466  if (!factorization_->status()) {
467    // put in standard form
468    createRim(7+8+16);
469    // do work
470    gutsOfSolution ( rowActivities, columnActivities);
471    // release extra memory
472    deleteRim(false);
473  }
474  return factorization_->status();
475}
476/* Given an existing factorization computes and checks
477   primal and dual solutions.  Uses current problem arrays for
478   bounds.  Returns feasibility states */
479int ClpSimplex::getSolution ( )
480{
481  double * rowActivities = new double[numberRows_];
482  double * columnActivities = new double[numberColumns_];
483  CoinDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
484  CoinDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
485  int status = getSolution( rowActivities, columnActivities);
486  delete [] rowActivities;
487  delete [] columnActivities;
488  return status;
489}
490// Factorizes using current basis.  This is for external use
491// Return codes are as from ClpFactorization
492int ClpSimplex::factorize ()
493{
494  // put in standard form
495  createRim(7+8+16);
496  // do work
497  int status = internalFactorize(-1);
498  // release extra memory
499  deleteRim(false);
500
501  return status;
502}
503
504/* Factorizes using current basis. 
505      solveType - 1 iterating, 0 initial, -1 external
506      If 10 added then in primal values pass
507*/
508/* Return codes are as from ClpFactorization unless initial factorization
509   when total number of singularities is returned */
510int ClpSimplex::internalFactorize ( int solveType)
511{
512
513  int iRow,iColumn;
514  int totalSlacks=numberRows_;
515
516  bool valuesPass=false;
517  if (solveType>=10) {
518    valuesPass=true;
519    solveType -= 10;
520  }
521  if (!solveType) {
522    if (!valuesPass) {
523      // not values pass so set to bounds
524      bool allSlack=true;
525      if (status_) {
526        for (iRow=0;iRow<numberRows_;iRow++) {
527          if (getRowStatus(iRow)!=basic) {
528            allSlack=false;
529            break;
530          }
531        }
532      }
533      if (!allSlack) {
534        // set values from warm start (if sensible)
535        int numberBasic=0;
536        for (iRow=0;iRow<numberRows_;iRow++) {
537          switch(getRowStatus(iRow)) {
538           
539          case basic:
540            numberBasic++;
541            break;
542          case isFree:
543            assert(rowLowerWork_[iRow]<-largeValue_);
544            assert(rowUpperWork_[iRow]>largeValue_);
545            rowActivityWork_[iRow]=0.0;
546            break;
547          case atUpperBound:
548            rowActivityWork_[iRow]=rowUpperWork_[iRow];
549            if (rowActivityWork_[iRow]>largeValue_) {
550              if (rowLowerWork_[iRow]>-largeValue_) {
551                rowActivityWork_[iRow]=rowLowerWork_[iRow];
552                setRowStatus(iRow,atLowerBound);
553              } else {
554                // say free
555                setRowStatus(iRow,isFree);
556                rowActivityWork_[iRow]=0.0;
557              }
558            }
559            break;
560          case atLowerBound:
561            rowActivityWork_[iRow]=rowLowerWork_[iRow];
562            if (rowActivityWork_[iRow]<-largeValue_) {
563              if (rowUpperWork_[iRow]<largeValue_) {
564                rowActivityWork_[iRow]=rowUpperWork_[iRow];
565                setRowStatus(iRow,atUpperBound);
566              } else {
567                // say free
568                setRowStatus(iRow,isFree);
569                rowActivityWork_[iRow]=0.0;
570              }
571            }
572            break;
573          case superBasic:
574            abort();
575            break;
576          }
577        }
578        totalSlacks=numberBasic;
579        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
580          switch(getColumnStatus(iColumn)) {
581           
582          case basic:
583            if (numberBasic==numberRows_) {
584              // take out of basis
585              if (columnLowerWork_[iColumn]>-largeValue_) {
586                if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]<
587                    columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) {
588                  columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
589                  setColumnStatus(iColumn,atLowerBound);
590                } else {
591                  columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
592                  setColumnStatus(iColumn,atUpperBound);
593                }
594              } else if (columnUpperWork_[iColumn]<largeValue_) {
595                columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
596                setColumnStatus(iColumn,atUpperBound);
597              } else {
598                columnActivityWork_[iColumn]=0.0;
599                setColumnStatus(iColumn,isFree);
600              }
601            } else {
602              numberBasic++;
603            }
604            break;
605          case isFree:
606            columnActivityWork_[iColumn]=0.0;
607            break;
608          case atUpperBound:
609            columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
610            if (columnActivityWork_[iColumn]>largeValue_) {
611              assert(columnLowerWork_[iColumn]>-largeValue_);
612              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
613              setColumnStatus(iColumn,atLowerBound);
614            }
615            break;
616          case atLowerBound:
617            columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
618            if (columnActivityWork_[iColumn]<-largeValue_) {
619              assert(columnUpperWork_[iColumn]<largeValue_);
620              columnActivityWork_[iColumn]=columnUpperWork_[iColumn];
621              setColumnStatus(iColumn,atUpperBound);
622            }
623            break;
624          case superBasic:
625            abort();
626            break;
627          }
628        }
629      } else {
630        //#define TESTFREE
631        // all slack basis
632        int numberBasic=0;
633        // changed to put free variables in basis
634        if (!status_) {
635          createStatus();
636        }
637        for (iRow=0;iRow<numberRows_;iRow++) {
638          double lower=rowLowerWork_[iRow];
639          double upper=rowUpperWork_[iRow];
640          if (lower>-largeValue_||upper<largeValue_) {
641            if (fabs(lower)<=fabs(upper)) {
642              setRowStatus(iRow,atLowerBound);
643              rowActivityWork_[iRow]=lower;
644            } else {
645              setRowStatus(iRow,atUpperBound);
646              rowActivityWork_[iRow]=upper;
647            }
648          } else {
649            setRowStatus(iRow,isFree);
650            rowActivityWork_[iRow]=0.0;
651          }
652#ifdef TESTFREE
653          setRowStatus(iRow,basic);
654          numberBasic++;
655#endif
656        }
657        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
658          double lower=columnLowerWork_[iColumn];
659          double upper=columnUpperWork_[iColumn];
660          double big_bound = largeValue_;
661          if (lower>-big_bound||upper<big_bound) {
662            if (fabs(lower)<=fabs(upper)) {
663              setColumnStatus(iColumn,atLowerBound);
664              columnActivityWork_[iColumn]=lower;
665            } else {
666              setColumnStatus(iColumn,atUpperBound);
667              columnActivityWork_[iColumn]=upper;
668            }
669          } else {
670#ifndef TESTFREE
671            numberBasic++;
672            setColumnStatus(iColumn,basic);
673#else
674            setColumnStatus(iColumn,isFree);
675#endif
676            columnActivityWork_[iColumn]=0.0;
677          }
678        }
679        assert(numberBasic<=numberRows_); // problems if too many free
680        if (!numberBasic) {
681          // might as well do all slack basis
682          for (iRow=0;iRow<numberRows_;iRow++) {
683            setRowStatus(iRow,basic);
684          }
685        }
686      }
687    } else {
688      // values pass has less coding
689      // make row activities correct
690      memset(rowActivityWork_,0,numberRows_*sizeof(double));
691      times(1.0,columnActivityWork_,rowActivityWork_);
692      if (status_) {
693        // set values from warm start (if sensible)
694        int numberBasic=0;
695        for (iRow=0;iRow<numberRows_;iRow++) {
696          if (getRowStatus(iRow)==basic) 
697            numberBasic++;
698          else
699            setRowStatus(iRow,superBasic);
700        }
701        totalSlacks=numberBasic;
702        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
703          if (getColumnStatus(iColumn)==basic) {
704            if (numberBasic==numberRows_) 
705              // take out of basis
706              setColumnStatus(iColumn,superBasic);
707            else 
708              numberBasic++;
709          } else {
710            setColumnStatus(iColumn,superBasic);
711          }
712        }
713      } else {
714        // all slack basis
715        int numberBasic=0;
716        if (!status_) {
717          createStatus();
718        }
719        for (iRow=0;iRow<numberRows_;iRow++) {
720          setRowStatus(iRow,basic);
721          numberBasic++;
722        }
723        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
724          setColumnStatus(iColumn,superBasic);
725        }
726      }
727    }
728    numberRefinements_=1;
729  }
730  int status=-99;
731  int * rowIsBasic = new int[numberRows_];
732  int * columnIsBasic = new int[numberColumns_];
733  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
734  while (status<-98) {
735 
736    int i;
737    int numberBasic=0;
738    for (i=0;i<numberRows_;i++) {
739      if (getRowStatus(i) == basic) {
740        rowIsBasic[i]=1;
741        numberBasic++;
742      } else {
743        rowIsBasic[i]=-1;
744      }
745    }
746    for (i=0;i<numberColumns_;i++) {
747      if (getColumnStatus(i) == basic) {
748        columnIsBasic[i]=1;
749        numberBasic++;
750      } else {
751        columnIsBasic[i]=-1;
752      }
753    }
754    assert (numberBasic<=numberRows_);
755    while (status==-99) {
756      status =  factorization_->factorize(this,matrix_,
757                                          numberRows_,numberColumns_,
758                                          rowIsBasic, columnIsBasic,
759                                          0.0);
760      if (status==-99) {
761        // get more memory
762        factorization_->areaFactor(2.0*factorization_->areaFactor());
763      }
764    }
765    if (!status) {
766      // do pivot information
767      for (i=0;i<numberRows_;i++) {
768        if (getRowStatus(i) == basic) {
769          pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
770        }
771      }
772      for (i=0;i<numberColumns_;i++) {
773        if (getColumnStatus(i) == basic) {
774          pivotVariable_[columnIsBasic[i]]=i;
775        }
776      }
777    } else {
778      // leave pivotVariable_ in useful form for cleaning basis
779      for (i=0;i<numberRows_;i++) {
780        pivotVariable_[i]=-1;
781      }
782      for (i=0;i<numberRows_;i++) {
783        if (getRowStatus(i) == basic) {
784          int iPivot = rowIsBasic[i];
785          if (iPivot>=0) 
786            pivotVariable_[iPivot]=i+numberColumns_;
787        }
788      }
789      for (i=0;i<numberColumns_;i++) {
790        if (getColumnStatus(i) == basic) {
791          int iPivot = columnIsBasic[i];
792          if (iPivot>=0) 
793            pivotVariable_[iPivot]=i;
794        }
795      }
796    }
797    if (status==-1) {
798      if (!solveType) {
799        //redo basis - first take ALL columns out
800        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
801          if (getColumnStatus(iColumn)==
802              basic) {
803            // take out
804            if (!valuesPass) {
805              double lower=columnLowerWork_[iColumn];
806              double upper=columnUpperWork_[iColumn];
807              double value=columnActivityWork_[iColumn];
808              if (lower>-largeValue_||upper<largeValue_) {
809                if (fabs(value-lower)<fabs(value-upper)) {
810                  setColumnStatus(iColumn,atLowerBound);
811                  columnActivityWork_[iColumn]=lower;
812                } else {
813                  setColumnStatus(iColumn,atUpperBound);
814                  columnActivityWork_[iColumn]=upper;
815                }
816              } else {
817                setColumnStatus(iColumn,isFree);
818              }
819            } else {
820              setColumnStatus(iColumn,superBasic);
821            }
822          }
823        }
824        for (iRow=0;iRow<numberRows_;iRow++) {
825          int iSequence=pivotVariable_[iRow];
826          if (iSequence>=0) {
827            // basic
828            if (iSequence>=numberColumns_) {
829              // slack in - leave
830              assert (iSequence-numberColumns_==iRow);
831            } else {
832              // put back structural
833              setColumnStatus(iSequence,basic);
834            }
835          } else {
836            // put in slack
837            setRowStatus(iRow,basic);
838          }
839        }
840        // signal repeat
841        status=-99;
842      }
843    }
844  } 
845  delete [] rowIsBasic;
846  delete [] columnIsBasic;
847  if (status) {
848    handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
849      <<status
850      <<CoinMessageEol;
851    return -1;
852  } else if (!solveType) {
853    // Initial basis - return number of singularities
854    int numberSlacks=0;
855    for (iRow=0;iRow<numberRows_;iRow++) {
856      if (getRowStatus(iRow) == basic)
857        numberSlacks++;
858    }
859    status= max(numberSlacks-totalSlacks,0);
860  }
861
862  // sparse methods
863  if (factorization_->sparseThreshold()) {
864    // get default value
865    factorization_->sparseThreshold(0);
866    factorization_->goSparse();
867  }
868 
869  return status;
870}
871/*
872   This does basis housekeeping and does values for in/out variables.
873   Can also decide to re-factorize
874*/
875int 
876ClpSimplex::housekeeping(double objectiveChange)
877{
878  numberIterations_++;
879  changeMade_++; // something has happened
880  // incoming variable
881
882  handler_->message(CLP_SIMPLEX_HOUSE1,messages_)
883    <<directionOut_
884    <<directionIn_<<theta_
885    <<dualOut_<<dualIn_<<alpha_
886    <<CoinMessageEol;
887  if (getStatus(sequenceIn_)==isFree) {
888    handler_->message(CLP_SIMPLEX_FREEIN,messages_)
889      <<sequenceIn_
890      <<CoinMessageEol;
891  }
892  // change of incoming
893  if (algorithm_<0) {
894    dualOut_ /= alpha_;
895    dualOut_ *= -directionOut_;
896  }
897  char rowcol[]={'R','C'};
898  double cost = cost_[sequenceIn_];
899  double value=valueIn_;
900  if (pivotRow_>=0)
901    pivotVariable_[pivotRow_]=sequenceIn();
902  if (directionIn_==-1) {
903    // as if from upper bound
904    if (sequenceIn_!=sequenceOut_) {
905      // variable becoming basic
906      setStatus(sequenceIn_,basic);
907      if (algorithm_<0) {
908        value = upperIn_+dualOut_;
909        dj_[sequenceIn_]=0.0;
910      } else {
911        value = valueIn_-fabs(theta_);
912      }
913    } else {
914      value=lowerIn_;
915      setStatus(sequenceIn_, atLowerBound);
916    }
917  } else {
918    // as if from lower bound
919    if (sequenceIn_!=sequenceOut_) {
920      // variable becoming basic
921      setStatus(sequenceIn_,basic);
922      if (algorithm_<0) {
923        value = lowerIn_+dualOut_;
924        dj_[sequenceIn_]=0.0;
925      } else {
926        value = valueIn_+fabs(theta_);
927      }
928    } else {
929      value=upperIn_;
930      setStatus(sequenceIn_, atUpperBound);
931    }
932  }
933  if (algorithm_<0)
934    objectiveChange += cost*(value-valueIn_);
935  else
936    objectiveChange += dualIn_*(value-valueIn_);
937  solution_[sequenceIn_]=value;
938
939  // outgoing
940  if (sequenceIn_!=sequenceOut_) {
941    assert( getStatus(sequenceOut_)== basic);
942    if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12)
943      progressFlag_ |= 1; // making real progress
944    if (algorithm_<0) {
945      if (directionOut_>0) {
946        value = lowerOut_;
947        setStatus(sequenceOut_,atLowerBound);
948        dj_[sequenceOut_]=theta_;
949      } else {
950        value = upperOut_;
951        setStatus(sequenceOut_,atUpperBound);
952        dj_[sequenceOut_]=-theta_;
953      }
954      solution_[sequenceOut_]=value;
955    } else {
956      if (directionOut_>0) {
957        value = lowerOut_;
958      } else {
959        value = upperOut_;
960      }
961      nonLinearCost_->setOne(sequenceOut_,value);
962      double lowerValue = lower_[sequenceOut_];
963      double upperValue = upper_[sequenceOut_];
964      assert(value>=lowerValue-primalTolerance_&&
965             value<=upperValue+primalTolerance_);
966      // may not be exactly at bound and bounds may have changed
967      if (value<=lowerValue+primalTolerance_) {
968        setStatus(sequenceOut_,atLowerBound);
969      } else if (value>=upperValue-primalTolerance_) {
970        setStatus(sequenceOut_,atUpperBound);
971      } else {
972        printf("*** variable wandered off bound %g %g %g!\n",
973               lowerValue,value,upperValue);
974        if (value-lowerValue<=upperValue-value) {
975          setStatus(sequenceOut_,atLowerBound);
976          value=lowerValue;
977        } else {
978          setStatus(sequenceOut_,atUpperBound);
979          value=upperValue;
980        }
981      }
982      solution_[sequenceOut_]=value;
983    }
984  }
985  // change cost and bounds on incoming if primal
986  if (algorithm_>0)
987    nonLinearCost_->setOne(sequenceIn_,solution_[sequenceIn_]); 
988  objectiveValue_ += objectiveChange;
989  handler_->message(CLP_SIMPLEX_HOUSE2,messages_)
990    <<numberIterations_<<objectiveValue_
991    <<rowcol[isColumn(sequenceIn_)]<<sequenceWithin(sequenceIn_)
992    <<rowcol[isColumn(sequenceOut_)]<<sequenceWithin(sequenceOut_);
993  handler_->printing(algorithm_<0)<<theta_<<dualOut_;
994  handler_->printing(algorithm_>0)<<dualIn_<<theta_;
995  handler_->message()<<CoinMessageEol;
996
997  if (numberIterations_>=maximumIterations_)
998    return 2;
999  // only time to re-factorize if one before real time
1000  // this is so user won't be surprised that maximumPivots has exact meaning
1001  if (factorization_->pivots()==factorization_->maximumPivots()) {
1002    return 1;
1003  } else {
1004    if (forceFactorization_>0&&
1005        factorization_->pivots()==forceFactorization_) {
1006      // relax
1007      forceFactorization_ *= 2;
1008      if (forceFactorization_>factorization_->maximumPivots())
1009        forceFactorization_ = -1; //off
1010      return 1;
1011    } else {
1012      // carry on iterating
1013      return 0;
1014    }
1015  }
1016}
1017// Copy constructor.
1018ClpSimplex::ClpSimplex(const ClpSimplex &rhs) :
1019  ClpModel(rhs),
1020  columnPrimalInfeasibility_(0.0),
1021  columnPrimalSequence_(-2),
1022  rowPrimalInfeasibility_(0.0),
1023  rowPrimalSequence_(-2), 
1024  columnDualInfeasibility_(0.0),
1025  columnDualSequence_(-2),
1026  rowDualInfeasibility_(0.0),
1027  rowDualSequence_(-2),
1028  primalToleranceToGetOptimal_(-1.0),
1029  remainingDualInfeasibility_(0.0),
1030  largeValue_(1.0e15),
1031  largestPrimalError_(0.0),
1032  largestDualError_(0.0),
1033  largestSolutionError_(0.0),
1034  dualBound_(1.0e7),
1035  lower_(NULL),
1036  rowLowerWork_(NULL),
1037  columnLowerWork_(NULL),
1038  upper_(NULL),
1039  rowUpperWork_(NULL),
1040  columnUpperWork_(NULL),
1041  cost_(NULL),
1042  rowObjectiveWork_(NULL),
1043  objectiveWork_(NULL),
1044  alpha_(0.0),
1045  theta_(0.0),
1046  lowerIn_(0.0),
1047  valueIn_(0.0),
1048  upperIn_(0.0),
1049  dualIn_(0.0),
1050  sequenceIn_(-1),
1051  directionIn_(-1),
1052  lowerOut_(-1),
1053  valueOut_(-1),
1054  upperOut_(-1),
1055  dualOut_(-1),
1056  sequenceOut_(-1),
1057  directionOut_(-1),
1058  pivotRow_(-1),
1059  status_(NULL),
1060  dj_(NULL),
1061  rowReducedCost_(NULL),
1062  reducedCostWork_(NULL),
1063  solution_(NULL),
1064  rowActivityWork_(NULL),
1065  columnActivityWork_(NULL),
1066  dualTolerance_(0.0),
1067  primalTolerance_(0.0),
1068  sumDualInfeasibilities_(0.0),
1069  numberDualInfeasibilities_(0),
1070  numberDualInfeasibilitiesWithoutFree_(0),
1071  sumPrimalInfeasibilities_(0.0),
1072  numberPrimalInfeasibilities_(0),
1073  pivotVariable_(NULL),
1074  factorization_(NULL),
1075  numberRefinements_(0),
1076  rowScale_(NULL),
1077  savedSolution_(NULL),
1078  columnScale_(NULL),
1079  scalingFlag_(0),
1080  numberTimesOptimal_(0),
1081  changeMade_(1),
1082  algorithm_(0),
1083  forceFactorization_(-1),
1084  perturbation_(100),
1085  infeasibilityCost_(1.0e7),
1086  nonLinearCost_(NULL),
1087  specialOptions_(0),
1088  lastBadIteration_(-999999),
1089  numberFake_(0),
1090  progressFlag_(0),
1091  sumOfRelaxedDualInfeasibilities_(0.0),
1092  sumOfRelaxedPrimalInfeasibilities_(0.0)
1093{
1094  int i;
1095  for (i=0;i<6;i++) {
1096    rowArray_[i]=NULL;
1097    columnArray_[i]=NULL;
1098  }
1099  saveStatus_=NULL;
1100  factorization_ = NULL;
1101  dualRowPivot_ = NULL;
1102  primalColumnPivot_ = NULL;
1103  gutsOfDelete(0);
1104  gutsOfCopy(rhs);
1105}
1106// Copy constructor from model
1107ClpSimplex::ClpSimplex(const ClpModel &rhs) :
1108  ClpModel(rhs),
1109  columnPrimalInfeasibility_(0.0),
1110  columnPrimalSequence_(-2),
1111  rowPrimalInfeasibility_(0.0),
1112  rowPrimalSequence_(-2), 
1113  columnDualInfeasibility_(0.0),
1114  columnDualSequence_(-2),
1115  rowDualInfeasibility_(0.0),
1116  rowDualSequence_(-2),
1117  primalToleranceToGetOptimal_(-1.0),
1118  remainingDualInfeasibility_(0.0),
1119  largeValue_(1.0e15),
1120  largestPrimalError_(0.0),
1121  largestDualError_(0.0),
1122  largestSolutionError_(0.0),
1123  dualBound_(1.0e7),
1124  lower_(NULL),
1125  rowLowerWork_(NULL),
1126  columnLowerWork_(NULL),
1127  upper_(NULL),
1128  rowUpperWork_(NULL),
1129  columnUpperWork_(NULL),
1130  cost_(NULL),
1131  rowObjectiveWork_(NULL),
1132  objectiveWork_(NULL),
1133  alpha_(0.0),
1134  theta_(0.0),
1135  lowerIn_(0.0),
1136  valueIn_(0.0),
1137  upperIn_(0.0),
1138  dualIn_(0.0),
1139  sequenceIn_(-1),
1140  directionIn_(-1),
1141  lowerOut_(-1),
1142  valueOut_(-1),
1143  upperOut_(-1),
1144  dualOut_(-1),
1145  sequenceOut_(-1),
1146  directionOut_(-1),
1147  pivotRow_(-1),
1148  status_(NULL),
1149  dj_(NULL),
1150  rowReducedCost_(NULL),
1151  reducedCostWork_(NULL),
1152  solution_(NULL),
1153  rowActivityWork_(NULL),
1154  columnActivityWork_(NULL),
1155  dualTolerance_(0.0),
1156  primalTolerance_(0.0),
1157  sumDualInfeasibilities_(0.0),
1158  numberDualInfeasibilities_(0),
1159  numberDualInfeasibilitiesWithoutFree_(0),
1160  sumPrimalInfeasibilities_(0.0),
1161  numberPrimalInfeasibilities_(0),
1162  pivotVariable_(NULL),
1163  factorization_(NULL),
1164  numberRefinements_(0),
1165  rowScale_(NULL),
1166  savedSolution_(NULL),
1167  columnScale_(NULL),
1168  scalingFlag_(0),
1169  numberTimesOptimal_(0),
1170  changeMade_(1),
1171  algorithm_(0),
1172  forceFactorization_(-1),
1173  perturbation_(100),
1174  infeasibilityCost_(1.0e7),
1175  nonLinearCost_(NULL),
1176  specialOptions_(0),
1177  lastBadIteration_(-999999),
1178  numberFake_(0),
1179  progressFlag_(0),
1180  sumOfRelaxedDualInfeasibilities_(0.0),
1181  sumOfRelaxedPrimalInfeasibilities_(0.0)
1182{
1183  int i;
1184  for (i=0;i<6;i++) {
1185    rowArray_[i]=NULL;
1186    columnArray_[i]=NULL;
1187  }
1188  saveStatus_=NULL;
1189  // get an empty factorization so we can set tolerances etc
1190  factorization_ = new ClpFactorization();
1191  // say Steepest pricing
1192  dualRowPivot_ = new ClpDualRowSteepest();
1193  // say Dantzig pricing
1194  primalColumnPivot_ = new ClpPrimalColumnDantzig();
1195 
1196}
1197// Assignment operator. This copies the data
1198ClpSimplex & 
1199ClpSimplex::operator=(const ClpSimplex & rhs)
1200{
1201  if (this != &rhs) {
1202    gutsOfDelete(0);
1203    ClpModel::operator=(rhs);
1204    gutsOfCopy(rhs);
1205  }
1206  return *this;
1207}
1208void 
1209ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
1210{
1211  lower_ = copyOfArray(rhs.lower_,numberColumns_+numberRows_);
1212  rowLowerWork_ = lower_+numberColumns_;
1213  columnLowerWork_ = lower_;
1214  upper_ = copyOfArray(rhs.upper_,numberColumns_+numberRows_);
1215  rowUpperWork_ = upper_+numberColumns_;
1216  columnUpperWork_ = upper_;
1217  cost_ = copyOfArray(rhs.cost_,(numberColumns_+numberRows_));
1218  objectiveWork_ = cost_;
1219  rowObjectiveWork_ = cost_+numberColumns_;
1220  dj_ = copyOfArray(rhs.dj_,numberRows_+numberColumns_);
1221  if (dj_) {
1222    reducedCostWork_ = dj_;
1223    rowReducedCost_ = dj_+numberColumns_;
1224  }
1225  solution_ = copyOfArray(rhs.solution_,numberRows_+numberColumns_);
1226  if (solution_) {
1227    columnActivityWork_ = solution_;
1228    rowActivityWork_ = solution_+numberColumns_;
1229  }
1230  if (rhs.pivotVariable_) {
1231    pivotVariable_ = new int[numberRows_];
1232    CoinDisjointCopyN ( rhs.pivotVariable_, numberRows_ , pivotVariable_);
1233  } else {
1234    pivotVariable_=NULL;
1235  }
1236  if (rhs.factorization_) {
1237    factorization_ = new ClpFactorization(*rhs.factorization_);
1238  } else {
1239    factorization_=NULL;
1240  }
1241  rowScale_ = copyOfArray(rhs.rowScale_,numberRows_);
1242  savedSolution_ = copyOfArray(rhs.savedSolution_,numberColumns_+numberRows_);
1243  columnScale_ = copyOfArray(rhs.columnScale_,numberColumns_);
1244  int i;
1245  for (i=0;i<6;i++) {
1246    rowArray_[i]=NULL;
1247    if (rhs.rowArray_[i]) 
1248      rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
1249    columnArray_[i]=NULL;
1250    if (rhs.columnArray_[i]) 
1251      columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
1252  }
1253  if (rhs.status_) {
1254    status_ = copyOfArray( rhs.status_,numberColumns_+numberRows_);
1255  }
1256  if (rhs.saveStatus_) {
1257    saveStatus_ = copyOfArray( rhs.saveStatus_,numberColumns_+numberRows_);
1258  }
1259  columnPrimalInfeasibility_ = rhs.columnPrimalInfeasibility_;
1260  columnPrimalSequence_ = rhs.columnPrimalSequence_;
1261  rowPrimalInfeasibility_ = rhs.rowPrimalInfeasibility_;
1262  rowPrimalSequence_ = rhs.rowPrimalSequence_;
1263  columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
1264  columnDualSequence_ = rhs.columnDualSequence_;
1265  rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
1266  rowDualSequence_ = rhs.rowDualSequence_;
1267  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
1268  remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
1269  largeValue_ = rhs.largeValue_;
1270  largestPrimalError_ = rhs.largestPrimalError_;
1271  largestDualError_ = rhs.largestDualError_;
1272  largestSolutionError_ = rhs.largestSolutionError_;
1273  dualBound_ = rhs.dualBound_;
1274  alpha_ = rhs.alpha_;
1275  theta_ = rhs.theta_;
1276  lowerIn_ = rhs.lowerIn_;
1277  valueIn_ = rhs.valueIn_;
1278  upperIn_ = rhs.upperIn_;
1279  dualIn_ = rhs.dualIn_;
1280  sequenceIn_ = rhs.sequenceIn_;
1281  directionIn_ = rhs.directionIn_;
1282  lowerOut_ = rhs.lowerOut_;
1283  valueOut_ = rhs.valueOut_;
1284  upperOut_ = rhs.upperOut_;
1285  dualOut_ = rhs.dualOut_;
1286  sequenceOut_ = rhs.sequenceOut_;
1287  directionOut_ = rhs.directionOut_;
1288  pivotRow_ = rhs.pivotRow_;
1289  numberRefinements_ = rhs.numberRefinements_;
1290  dualTolerance_ = rhs.dualTolerance_;
1291  primalTolerance_ = rhs.primalTolerance_;
1292  sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
1293  numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
1294  numberDualInfeasibilitiesWithoutFree_ = 
1295    rhs.numberDualInfeasibilitiesWithoutFree_;
1296  sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
1297  numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
1298  dualRowPivot_ = rhs.dualRowPivot_->clone(true);
1299  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
1300  scalingFlag_ = rhs.scalingFlag_;
1301  numberTimesOptimal_ = rhs.numberTimesOptimal_;
1302  changeMade_ = rhs.changeMade_;
1303  algorithm_ = rhs.algorithm_;
1304  forceFactorization_ = rhs.forceFactorization_;
1305  perturbation_ = rhs.perturbation_;
1306  infeasibilityCost_ = rhs.infeasibilityCost_;
1307  specialOptions_ = rhs.specialOptions_;
1308  lastBadIteration_ = rhs.lastBadIteration_;
1309  numberFake_ = rhs.numberFake_;
1310  progressFlag_ = rhs.progressFlag_;
1311  sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
1312  sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
1313  if (rhs.nonLinearCost_!=NULL)
1314    nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
1315}
1316// type == 0 do everything, most + pivot data, 2 factorization data as well
1317void 
1318ClpSimplex::gutsOfDelete(int type)
1319{
1320  delete [] lower_;
1321  lower_=NULL;
1322  rowLowerWork_=NULL;
1323  columnLowerWork_=NULL;
1324  delete [] upper_;
1325  upper_=NULL;
1326  rowUpperWork_=NULL;
1327  columnUpperWork_=NULL;
1328  delete [] cost_;
1329  cost_=NULL;
1330  objectiveWork_=NULL;
1331  rowObjectiveWork_=NULL;
1332  delete [] dj_;
1333  dj_=NULL;
1334  reducedCostWork_=NULL;
1335  rowReducedCost_=NULL;
1336  delete [] solution_;
1337  solution_=NULL;
1338  rowActivityWork_=NULL;
1339  columnActivityWork_=NULL;
1340  delete [] rowScale_;
1341  rowScale_ = NULL;
1342  delete [] savedSolution_;
1343  savedSolution_ = NULL;
1344  delete [] columnScale_;
1345  columnScale_ = NULL;
1346  delete nonLinearCost_;
1347  nonLinearCost_ = NULL;
1348  int i;
1349  for (i=0;i<6;i++) {
1350    delete rowArray_[i];
1351    rowArray_[i]=NULL;
1352    delete columnArray_[i];
1353    columnArray_[i]=NULL;
1354  }
1355  delete rowCopy_;
1356  rowCopy_=NULL;
1357  delete [] saveStatus_;
1358  saveStatus_=NULL;
1359  if (!type) {
1360    // delete everything
1361    delete factorization_;
1362    factorization_ = NULL;
1363    delete [] pivotVariable_;
1364    pivotVariable_=NULL;
1365    delete dualRowPivot_;
1366    dualRowPivot_ = NULL;
1367    delete primalColumnPivot_;
1368    primalColumnPivot_ = NULL;
1369    delete status_;
1370    status_=NULL;
1371  } else {
1372    // delete any size information in methods
1373    if (type>1) {
1374      factorization_->clearArrays();
1375      delete [] pivotVariable_;
1376      pivotVariable_=NULL;
1377    }
1378    dualRowPivot_->clearArrays();
1379    primalColumnPivot_->clearArrays();
1380  }
1381}
1382// This sets largest infeasibility and most infeasible
1383void 
1384ClpSimplex::checkPrimalSolution(const double * rowActivities,
1385                                        const double * columnActivities)
1386{
1387  double * solution;
1388  int iRow,iColumn;
1389
1390  objectiveValue_ = 0.0;
1391  // now look at primal solution
1392  columnPrimalInfeasibility_=0.0;
1393  columnPrimalSequence_=-1;
1394  rowPrimalInfeasibility_=0.0;
1395  rowPrimalSequence_=-1;
1396  largestSolutionError_=0.0;
1397  solution = rowActivityWork_;
1398  sumPrimalInfeasibilities_=0.0;
1399  numberPrimalInfeasibilities_=0;
1400  double primalTolerance = primalTolerance_;
1401  double relaxedTolerance=dualTolerance_;
1402  // we can't really trust infeasibilities if there is primal error
1403  double error = min(1.0e-3,largestPrimalError_);
1404  relaxedTolerance = max(relaxedTolerance, error);
1405  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
1406
1407  for (iRow=0;iRow<numberRows_;iRow++) {
1408    double infeasibility=0.0;
1409    objectiveValue_ += solution[iRow]*rowObjectiveWork_[iRow];
1410    if (solution[iRow]>rowUpperWork_[iRow]) {
1411      infeasibility=solution[iRow]-rowUpperWork_[iRow];
1412    } else if (solution[iRow]<rowLowerWork_[iRow]) {
1413      infeasibility=rowLowerWork_[iRow]-solution[iRow];
1414    }
1415    if (infeasibility>primalTolerance) {
1416      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1417      if (infeasibility>relaxedTolerance) 
1418        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1419      numberPrimalInfeasibilities_ ++;
1420    }
1421    if (infeasibility>rowPrimalInfeasibility_) {
1422      rowPrimalInfeasibility_=infeasibility;
1423      rowPrimalSequence_=iRow;
1424    }
1425    infeasibility = fabs(rowActivities[iRow]-solution[iRow]);
1426    if (infeasibility>largestSolutionError_)
1427      largestSolutionError_=infeasibility;
1428    if (rowLowerWork_[iRow]!=rowUpperWork_[iRow])
1429      clearFixed(iRow+numberColumns_);
1430    else
1431      setFixed(iRow+numberColumns_);
1432  }
1433  solution = columnActivityWork_;
1434  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1435    double infeasibility=0.0;
1436    objectiveValue_ += objectiveWork_[iColumn]*solution[iColumn];
1437    if (solution[iColumn]>columnUpperWork_[iColumn]) {
1438      infeasibility=solution[iColumn]-columnUpperWork_[iColumn];
1439    } else if (solution[iColumn]<columnLowerWork_[iColumn]) {
1440      infeasibility=columnLowerWork_[iColumn]-solution[iColumn];
1441    }
1442    if (infeasibility>columnPrimalInfeasibility_) {
1443      columnPrimalInfeasibility_=infeasibility;
1444      columnPrimalSequence_=iColumn;
1445    }
1446    if (infeasibility>primalTolerance) {
1447      sumPrimalInfeasibilities_ += infeasibility-primalTolerance_;
1448      if (infeasibility>relaxedTolerance) 
1449        sumOfRelaxedPrimalInfeasibilities_ += infeasibility-relaxedTolerance;
1450      numberPrimalInfeasibilities_ ++;
1451    }
1452    infeasibility = fabs(columnActivities[iColumn]-solution[iColumn]);
1453    if (infeasibility>largestSolutionError_)
1454      largestSolutionError_=infeasibility;
1455    if (columnUpperWork_[iColumn]-columnLowerWork_[iColumn]>primalTolerance_)
1456      clearFixed(iColumn);
1457    else
1458      setFixed(iColumn);
1459  }
1460}
1461void 
1462ClpSimplex::checkDualSolution()
1463{
1464
1465  double * solution;
1466  int iRow,iColumn;
1467  sumDualInfeasibilities_=0.0;
1468  numberDualInfeasibilities_=0;
1469  numberDualInfeasibilitiesWithoutFree_=0;
1470  columnDualInfeasibility_=0.0;
1471  columnDualSequence_=-1;
1472  rowDualInfeasibility_=0.0;
1473  rowDualSequence_=-1;
1474  primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
1475                                   columnPrimalInfeasibility_);
1476  remainingDualInfeasibility_=0.0;
1477  solution = rowActivityWork_;
1478  double relaxedTolerance=dualTolerance_;
1479  // we can't really trust infeasibilities if there is dual error
1480  double error = min(1.0e-3,largestDualError_);
1481  relaxedTolerance = max(relaxedTolerance, error);
1482  sumOfRelaxedDualInfeasibilities_ = 0.0;
1483
1484  for (iRow=0;iRow<numberRows_;iRow++) {
1485    if (getRowStatus(iRow) != basic) {
1486      // not basic
1487      double value = rowReducedCost_[iRow];
1488      double distance;
1489      distance = rowUpperWork_[iRow]-solution[iRow];
1490      if (distance>primalTolerance_) {
1491        // should not be negative
1492        if (value<0.0) {
1493          value = - value;
1494          if (value>rowDualInfeasibility_) {
1495            rowDualInfeasibility_=value;
1496            rowDualSequence_=iRow;
1497          }
1498          if (value>dualTolerance_) {
1499            sumDualInfeasibilities_ += value-dualTolerance_;
1500            if (value>relaxedTolerance) 
1501              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1502            numberDualInfeasibilities_ ++;
1503            if (getRowStatus(iRow) != isFree) 
1504              numberDualInfeasibilitiesWithoutFree_ ++;
1505            // maybe we can make feasible by increasing tolerance
1506            if (distance<largeValue_) {
1507              if (distance>primalToleranceToGetOptimal_)
1508                primalToleranceToGetOptimal_=distance;
1509            } else {
1510              //gap too big for any tolerance
1511              remainingDualInfeasibility_=
1512                max(remainingDualInfeasibility_,value);
1513            }
1514          }
1515        }
1516      }
1517      distance = solution[iRow] -rowLowerWork_[iRow];
1518      if (distance>primalTolerance_) {
1519        // should not be positive
1520        if (value>0.0) {
1521          if (value>rowDualInfeasibility_) {
1522            rowDualInfeasibility_=value;
1523            rowDualSequence_=iRow;
1524          }
1525          if (value>dualTolerance_) {
1526            sumDualInfeasibilities_ += value-dualTolerance_;
1527            if (value>relaxedTolerance) 
1528              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1529            numberDualInfeasibilities_ ++;
1530            if (getRowStatus(iRow) != isFree) 
1531              numberDualInfeasibilitiesWithoutFree_ ++;
1532            // maybe we can make feasible by increasing tolerance
1533            if (distance<largeValue_&&
1534                distance>primalToleranceToGetOptimal_)
1535              primalToleranceToGetOptimal_=distance;
1536          }
1537        }
1538      }
1539    }
1540  }
1541  solution = columnActivityWork_;
1542  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1543    if (getColumnStatus(iColumn) != basic) {
1544      // not basic
1545      double value = reducedCostWork_[iColumn];
1546      double distance;
1547      distance = columnUpperWork_[iColumn]-solution[iColumn];
1548      if (distance>primalTolerance_) {
1549        // should not be negative
1550        if (value<0.0) {
1551          value = - value;
1552          if (value>columnDualInfeasibility_) {
1553            columnDualInfeasibility_=value;
1554            columnDualSequence_=iColumn;
1555          }
1556          if (value>dualTolerance_) {
1557            sumDualInfeasibilities_ += value-dualTolerance_;
1558            if (value>relaxedTolerance) 
1559              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1560            numberDualInfeasibilities_ ++;
1561            if (getColumnStatus(iColumn) != isFree) 
1562              numberDualInfeasibilitiesWithoutFree_ ++;
1563            // maybe we can make feasible by increasing tolerance
1564            if (distance<largeValue_) {
1565              if (distance>primalToleranceToGetOptimal_)
1566                primalToleranceToGetOptimal_=distance;
1567            } else {
1568              //gap too big for any tolerance
1569              remainingDualInfeasibility_=
1570                max(remainingDualInfeasibility_,value);
1571            }
1572          }
1573        }
1574      }
1575      distance = solution[iColumn] -columnLowerWork_[iColumn];
1576      if (distance>primalTolerance_) {
1577        // should not be positive
1578        if (value>0.0) {
1579          if (value>columnDualInfeasibility_) {
1580            columnDualInfeasibility_=value;
1581            columnDualSequence_=iColumn;
1582          }
1583          if (value>dualTolerance_) {
1584            sumDualInfeasibilities_ += value-dualTolerance_;
1585            if (value>relaxedTolerance) 
1586              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
1587            numberDualInfeasibilities_ ++;
1588            if (getColumnStatus(iColumn) != isFree) 
1589              numberDualInfeasibilitiesWithoutFree_ ++;
1590            // maybe we can make feasible by increasing tolerance
1591            if (distance<largeValue_&&
1592                distance>primalToleranceToGetOptimal_)
1593              primalToleranceToGetOptimal_=distance;
1594          }
1595        }
1596      }
1597    }
1598  }
1599}
1600/*
1601  Unpacks one column of the matrix into indexed array
1602*/
1603void 
1604ClpSimplex::unpack(CoinIndexedVector * rowArray)
1605{
1606  rowArray->clear();
1607  if (sequenceIn_>=numberColumns_) {
1608    //slack
1609    rowArray->insert(sequenceIn_-numberColumns_,-1.0);
1610  } else {
1611    // column
1612    matrix_->unpack(this,rowArray,sequenceIn_);
1613  }
1614}
1615void 
1616ClpSimplex::unpack(CoinIndexedVector * rowArray,int sequence)
1617{
1618  rowArray->clear();
1619  if (sequence>=numberColumns_) {
1620    //slack
1621    rowArray->insert(sequence-numberColumns_,-1.0);
1622  } else {
1623    // column
1624    matrix_->unpack(this,rowArray,sequence);
1625  }
1626}
1627bool
1628ClpSimplex::createRim(int what,bool makeRowCopy)
1629{
1630  bool goodMatrix=true;
1631  if ((what&16)!=0) {
1632    // move information to work arrays
1633    // row reduced costs
1634    if (!dj_) {
1635      dj_ = new double[numberRows_+numberColumns_];
1636      reducedCostWork_ = dj_;
1637      rowReducedCost_ = dj_+numberColumns_;
1638      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1639    }
1640    if (!solution_) {
1641      solution_ = new double[numberRows_+numberColumns_];
1642      columnActivityWork_ = solution_;
1643      rowActivityWork_ = solution_+numberColumns_;
1644      memcpy(columnActivityWork_,columnActivity_,
1645             numberColumns_*sizeof(double));
1646      memcpy(rowActivityWork_,rowActivity_,
1647             numberRows_*sizeof(double));
1648    }
1649    //check matrix
1650    if (!matrix_->allElementsInRange(this,1.0e-20,1.0e20)) {
1651      problemStatus_=4;
1652      goodMatrix= false;
1653    }
1654    if (makeRowCopy) {
1655      delete rowCopy_;
1656      // may return NULL if can't give row copy
1657      rowCopy_ = matrix_->reverseOrderedCopy();
1658    }
1659  }
1660  int i;
1661  if ((what&4)!=0) {
1662    delete [] cost_;
1663    cost_ = new double[numberColumns_+numberRows_];
1664    objectiveWork_ = cost_;
1665    rowObjectiveWork_ = cost_+numberColumns_;
1666    memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
1667    if (rowObjective_)
1668      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1669    else
1670      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1671  }
1672  if ((what&1)!=0) {
1673    delete [] lower_;
1674    delete [] upper_;
1675    lower_ = new double[numberColumns_+numberRows_];
1676    upper_ = new double[numberColumns_+numberRows_];
1677    rowLowerWork_ = lower_+numberColumns_;
1678    columnLowerWork_ = lower_;
1679    rowUpperWork_ = upper_+numberColumns_;
1680    columnUpperWork_ = upper_;
1681    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1682    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1683    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1684    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1685    // clean up any mismatches on infinity
1686    for (i=0;i<numberColumns_;i++) {
1687      if (columnLowerWork_[i]<-1.0e30)
1688        columnLowerWork_[i] = -DBL_MAX;
1689      if (columnUpperWork_[i]>1.0e30)
1690        columnUpperWork_[i] = DBL_MAX;
1691    }
1692    // clean up any mismatches on infinity
1693    for (i=0;i<numberRows_;i++) {
1694      if (rowLowerWork_[i]<-1.0e30)
1695        rowLowerWork_[i] = -DBL_MAX;
1696      if (rowUpperWork_[i]>1.0e30)
1697        rowUpperWork_[i] = DBL_MAX;
1698    }
1699  }
1700  // do scaling if needed
1701  if (scalingFlag_>0&&!rowScale_) {
1702    if (matrix_->scale(this))
1703      scalingFlag_=-scalingFlag_; // not scaled after all
1704  }
1705  if ((what&4)!=0) {
1706    double direction = optimizationDirection_;
1707    // direction is actually scale out not scale in
1708    if (direction)
1709      direction = 1.0/direction;
1710    // but also scale by scale factors
1711    // not really sure about row scaling
1712    if (!rowScale_) {
1713      if (direction!=1.0) {
1714        for (i=0;i<numberRows_;i++)
1715          rowObjectiveWork_[i] *= direction;
1716        for (i=0;i<numberColumns_;i++)
1717          objectiveWork_[i] *= direction;
1718      }
1719    } else {
1720      for (i=0;i<numberRows_;i++)
1721        rowObjectiveWork_[i] *= direction/rowScale_[i];
1722      for (i=0;i<numberColumns_;i++)
1723        objectiveWork_[i] *= direction*columnScale_[i];
1724    }
1725  }
1726  if ((what&1)!=0&&rowScale_) {
1727    for (i=0;i<numberColumns_;i++) {
1728      double multiplier = 1.0/columnScale_[i];
1729      if (columnLowerWork_[i]>-1.0e50)
1730        columnLowerWork_[i] *= multiplier;
1731      if (columnUpperWork_[i]<1.0e50)
1732        columnUpperWork_[i] *= multiplier;
1733    }
1734    for (i=0;i<numberRows_;i++) {
1735      if (rowLowerWork_[i]>-1.0e50)
1736        rowLowerWork_[i] *= rowScale_[i];
1737      if (rowUpperWork_[i]<1.0e50)
1738        rowUpperWork_[i] *= rowScale_[i];
1739    }
1740  }
1741  if ((what&8)!=0&&rowScale_) {
1742    // on entry
1743    for (i=0;i<numberColumns_;i++) {
1744      columnActivityWork_[i] /= columnScale_[i];
1745      reducedCostWork_[i] *= columnScale_[i];
1746    }
1747    for (i=0;i<numberRows_;i++) {
1748      rowActivityWork_[i] *= rowScale_[i];
1749      dual_[i] /= rowScale_[i];
1750    }
1751  }
1752 
1753  if ((what&16)!=0) {
1754    // check rim of problem okay
1755    if (!sanityCheck())
1756      goodMatrix=false;
1757  } 
1758  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1759  // maybe we need to move scales to SimplexModel for factorization?
1760  if ((what&8)!=0&&!pivotVariable_) {
1761    pivotVariable_=new int[numberRows_];
1762  }
1763
1764  if ((what&16)!=0&&!rowArray_[2]) {
1765    // get some arrays
1766    int iRow,iColumn;
1767    // these are "indexed" arrays so we always know where nonzeros are
1768    /**********************************************************
1769      rowArray_[3] is long enough for columns as well
1770    *********************************************************/
1771    for (iRow=0;iRow<4;iRow++) {
1772      if (!rowArray_[iRow]) {
1773        rowArray_[iRow]=new CoinIndexedVector();
1774        int length =numberRows_+factorization_->maximumPivots();
1775        if (iRow==3)
1776          length = max(length,numberColumns_);
1777        rowArray_[iRow]->reserve(length);
1778      }
1779    }
1780   
1781    for (iColumn=0;iColumn<2;iColumn++) {
1782      if (!columnArray_[iColumn]) {
1783        columnArray_[iColumn]=new CoinIndexedVector();
1784        columnArray_[iColumn]->reserve(numberColumns_);
1785      }
1786    }   
1787  }
1788  return goodMatrix;
1789}
1790void
1791ClpSimplex::deleteRim(bool getRidOfFactorizationData)
1792{
1793  int i;
1794  if (problemStatus_!=1&&problemStatus_!=2) {
1795    delete [] ray_;
1796    ray_=NULL;
1797  }
1798  if (rowScale_) {
1799    for (i=0;i<numberColumns_;i++) {
1800      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1801      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1802    }
1803    for (i=0;i<numberRows_;i++) {
1804      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1805      dual_[i] *= rowScale_[i];
1806    }
1807    if (problemStatus_==2) {
1808      for (i=0;i<numberColumns_;i++) {
1809        ray_[i] *= columnScale_[i];
1810      }
1811    } else if (problemStatus_==1) {
1812      for (i=0;i<numberRows_;i++) {
1813        ray_[i] *= rowScale_[i];
1814      }
1815    }
1816  } else {
1817    for (i=0;i<numberColumns_;i++) {
1818      columnActivity_[i] = columnActivityWork_[i];
1819      reducedCost_[i] = reducedCostWork_[i];
1820    }
1821    for (i=0;i<numberRows_;i++) {
1822      rowActivity_[i] = rowActivityWork_[i];
1823    }
1824  }
1825  // direction may have been modified by scaling - clean up
1826  if (optimizationDirection_>0.0)
1827    optimizationDirection_ = 1.0;
1828  else if (optimizationDirection_<0.0)
1829    optimizationDirection_ = -1.0;
1830  // scaling may have been turned off
1831  scalingFlag_ = abs(scalingFlag_);
1832  if(getRidOfFactorizationData)
1833    gutsOfDelete(2);
1834  else
1835    gutsOfDelete(1);
1836}
1837void 
1838ClpSimplex::setDualBound(double value)
1839{
1840  if (value>0.0)
1841    dualBound_=value;
1842}
1843void 
1844ClpSimplex::setInfeasibilityCost(double value)
1845{
1846  if (value>0.0)
1847    infeasibilityCost_=value;
1848}
1849void ClpSimplex::setnumberRefinements( int value) 
1850{
1851  if (value>=0&&value<10)
1852    numberRefinements_=value;
1853}
1854// Sets row pivot choice algorithm in dual
1855void 
1856ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
1857{
1858  delete dualRowPivot_;
1859  dualRowPivot_ = choice.clone(true);
1860}
1861// Sets row pivot choice algorithm in dual
1862void 
1863ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
1864{
1865  delete primalColumnPivot_;
1866  primalColumnPivot_ = choice.clone(true);
1867}
1868// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
1869void 
1870ClpSimplex::scaling(int mode)
1871{
1872  if (mode>0&&mode<3) {
1873    scalingFlag_=mode;
1874  } else if (!mode) {
1875    scalingFlag_=0;
1876    delete [] rowScale_;
1877    rowScale_ = NULL;
1878    delete [] columnScale_;
1879    columnScale_ = NULL;
1880  }
1881}
1882// Passes in factorization
1883void 
1884ClpSimplex::setFactorization( ClpFactorization & factorization)
1885{
1886  delete factorization_;
1887  factorization_= new ClpFactorization(factorization);
1888}
1889void 
1890ClpSimplex::times(double scalar,
1891                  const double * x, double * y) const
1892{
1893  if (rowScale_)
1894    matrix_->times(scalar,x,y,rowScale_,columnScale_);
1895  else
1896    matrix_->times(scalar,x,y);
1897}
1898void 
1899ClpSimplex::transposeTimes(double scalar,
1900                           const double * x, double * y) const 
1901{
1902  if (rowScale_)
1903    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
1904  else
1905    matrix_->transposeTimes(scalar,x,y);
1906}
1907/* Perturbation:
1908   -50 to +50 - perturb by this power of ten (-6 sounds good)
1909   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1910   101 - we are perturbed
1911   102 - don't try perturbing again
1912   default is 100
1913*/
1914void 
1915ClpSimplex::setPerturbation(int value)
1916{
1917  if(value<=100&&value >=-50) {
1918    perturbation_=value;
1919  } 
1920}
1921// Sparsity on or off
1922bool 
1923ClpSimplex::sparseFactorization() const
1924{
1925  return factorization_->sparseThreshold();
1926}
1927void 
1928ClpSimplex::setSparseFactorization(bool value)
1929{
1930  if (value) {
1931    if (!factorization_->sparseThreshold())
1932      factorization_->goSparse();
1933  } else {
1934    factorization_->sparseThreshold(0);
1935  }
1936}
1937/* Tightens primal bounds to make dual faster.  Unless
1938   fixed, bounds are slightly looser than they could be.
1939   This is to make dual go faster and is probably not needed
1940   with a presolve.  Returns non-zero if problem infeasible
1941*/
1942int 
1943ClpSimplex::tightenPrimalBounds()
1944{
1945 
1946  // Get a row copy in standard format
1947  CoinPackedMatrix copy;
1948  copy.reverseOrderedCopyOf(*matrix());
1949  // get matrix data pointers
1950  const int * column = copy.getIndices();
1951  const int * rowStart = copy.getVectorStarts();
1952  const int * rowLength = copy.getVectorLengths(); 
1953  const double * element = copy.getElements();
1954  int numberChanged=1,iPass=0;
1955  double large = largeValue(); // treat bounds > this as infinite
1956  int numberInfeasible=0;
1957  int totalTightened = 0;
1958
1959  double tolerance = primalTolerance();
1960
1961
1962  // Save column bounds
1963  double * saveLower = new double [numberColumns_];
1964  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
1965  double * saveUpper = new double [numberColumns_];
1966  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
1967#define MAXPASS 10 
1968
1969  // Loop round seeing if we can tighten bounds
1970  // Would be faster to have a stack of possible rows
1971  // and we put altered rows back on stack
1972
1973  int iRow, iColumn;
1974
1975  while(numberChanged) {
1976
1977    numberChanged = 0; // Bounds tightened this pass
1978   
1979    if (iPass==MAXPASS) break;
1980    iPass++;
1981   
1982    for (iRow = 0; iRow < numberRows_; iRow++) {
1983
1984      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
1985
1986        // possible row
1987        int infiniteUpper = 0;
1988        int infiniteLower = 0;
1989        double maximumUp = 0.0;
1990        double maximumDown = 0.0;
1991        double newBound;
1992        int rStart = rowStart[iRow];
1993        int rEnd = rowStart[iRow]+rowLength[iRow];
1994        int j;
1995
1996        // Compute possible lower and upper ranges
1997
1998        for (j = rStart; j < rEnd; ++j) {
1999          double value=element[j];
2000          iColumn = column[j];
2001          if (value > 0.0) {
2002            if (columnUpper_[iColumn] >= large) {
2003              maximumUp = 1e31;
2004              ++infiniteUpper;
2005            } else {
2006              maximumUp += columnUpper_[iColumn] * value;
2007            }
2008            if (columnLower_[iColumn] <= -large) {
2009              maximumDown = -1e31;
2010              ++infiniteLower;
2011            } else {
2012              maximumDown += columnLower_[iColumn] * value;
2013            }
2014          } else if (value<0.0) {
2015            if (columnUpper_[iColumn] >= large) {
2016              maximumDown = -1e31;
2017              ++infiniteLower;
2018            } else {
2019              maximumDown += columnUpper_[iColumn] * value;
2020            }
2021            if (columnLower_[iColumn] <= -large) {
2022              maximumUp = 1e31;
2023              ++infiniteUpper;
2024            } else {
2025              maximumUp += columnLower_[iColumn] * value;
2026            }
2027          }
2028        }
2029        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2030            maximumDown >= rowLower_[iRow] - tolerance) {
2031
2032          // Row is redundant - make totally free
2033          rowLower_[iRow]=-DBL_MAX;
2034          rowUpper_[iRow]=DBL_MAX;
2035
2036        } else {
2037          if (maximumUp < rowLower_[iRow] -tolerance ||
2038              maximumDown > rowUpper_[iRow]+tolerance) {
2039            // problem is infeasible - exit at once
2040            numberInfeasible++;
2041            break;
2042          }
2043
2044          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2045            for (j = rStart; j < rEnd; ++j) {
2046              double value=element[j];
2047              iColumn = column[j];
2048              if (value > 0.0) {
2049                if (columnUpper_[iColumn] < large) {
2050                  newBound = columnUpper_[iColumn] + 
2051                    (rowLower_[iRow] - maximumUp) / value;
2052                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2053                    // Tighten the lower bound
2054
2055                    columnLower_[iColumn] = newBound;
2056                    ++numberChanged;
2057
2058                    // check infeasible (relaxed)
2059                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2060                        -100.0*tolerance) 
2061                      numberInfeasible++;
2062                    infiniteLower=1; // skip looking at other bound
2063                  }
2064                }
2065              } else {
2066                if (columnLower_[iColumn] > -large) {
2067                  newBound = columnLower_[iColumn] + 
2068                    (rowLower_[iRow] - maximumUp) / value;
2069                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2070                    // Tighten the upper bound
2071
2072                    columnUpper_[iColumn] = newBound;
2073                    ++numberChanged;
2074
2075                    // check infeasible (relaxed)
2076                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2077                        -100.0*tolerance) 
2078                      numberInfeasible++;
2079                    infiniteLower=1; // skip looking at other bound
2080                  }
2081                }
2082              }
2083            }
2084          }
2085         
2086          // Try other way
2087          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2088            for (j = rStart; j < rEnd; ++j) {
2089              double value=element[j];
2090              iColumn = column[j];
2091              if (value < 0.0) {
2092                if (columnUpper_[iColumn] < large) {
2093                  newBound = columnUpper_[iColumn] + 
2094                    (rowUpper_[iRow] - maximumDown) / value;
2095                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2096                    // Tighten the lower bound
2097
2098                    columnLower_[iColumn] = newBound;
2099                    ++numberChanged;
2100
2101                    // check infeasible (relaxed)
2102                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2103                        -100.0*tolerance) 
2104                      numberInfeasible++;
2105                  }
2106                } 
2107              } else {
2108                if (columnLower_[iColumn] > -large) {
2109                  newBound = columnLower_[iColumn] + 
2110                    (rowUpper_[iRow] - maximumDown) / value;
2111                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2112                    // Tighten the upper bound
2113
2114                    columnUpper_[iColumn] = newBound;
2115                    ++numberChanged;
2116
2117                    // check infeasible (relaxed)
2118                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2119                        -100.0*tolerance) 
2120                      numberInfeasible++;
2121                  }
2122                }
2123              }
2124            }
2125          }
2126        }
2127      }
2128    }
2129    totalTightened += numberChanged;
2130    if (numberInfeasible) break;
2131  }
2132  if (!numberInfeasible) {
2133    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2134      <<totalTightened
2135      <<CoinMessageEol;
2136    // Set bounds slightly loose
2137    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2138      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2139        // fix
2140        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2141          columnLower_[iColumn]=columnUpper_[iColumn];
2142        else
2143          columnUpper_[iColumn]=columnLower_[iColumn];
2144      } else {
2145        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2146          // relax a bit
2147          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2148                                      saveUpper[iColumn]);
2149        }
2150        if (columnLower_[iColumn]>saveLower[iColumn]) {
2151          // relax a bit
2152          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2153                                      saveLower[iColumn]);
2154        }
2155      }
2156    }
2157  } else {
2158    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2159      <<numberInfeasible
2160      <<CoinMessageEol;
2161    // restore column bounds
2162    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2163    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2164  }
2165  delete [] saveLower;
2166  delete [] saveUpper;
2167  return (numberInfeasible);
2168}
2169// dual
2170#include "ClpSimplexDual.hpp"
2171int ClpSimplex::dual ( )
2172{
2173  return ((ClpSimplexDual *) this)->dual();
2174}
2175// primal
2176#include "ClpSimplexPrimal.hpp"
2177int ClpSimplex::primal (int ifValuesPass )
2178{
2179  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2180}
2181/* For strong branching.  On input lower and upper are new bounds
2182   while on output they are objective function values (>1.0e50 infeasible).
2183   Return code is 0 if nothing interesting, -1 if infeasible both
2184   ways and +1 if infeasible one way (check values to see which one(s))
2185*/
2186int ClpSimplex::strongBranching(int numberVariables,const int * variables,
2187                                double * newLower, double * newUpper,
2188                                bool stopOnFirstInfeasible,
2189                                bool alwaysFinish)
2190{
2191  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
2192                                                    newLower,  newUpper,
2193                                                    stopOnFirstInfeasible,
2194                                                    alwaysFinish);
2195}
2196/* Borrow model.  This is so we dont have to copy large amounts
2197   of data around.  It assumes a derived class wants to overwrite
2198   an empty model with a real one - while it does an algorithm.
2199   This is same as ClpModel one, but sets scaling on etc. */
2200void 
2201ClpSimplex::borrowModel(ClpModel & otherModel) 
2202{
2203  ClpModel::borrowModel(otherModel);
2204  createStatus();
2205  scaling();
2206  ClpDualRowSteepest steep1;
2207  setDualRowPivotAlgorithm(steep1);
2208  ClpPrimalColumnSteepest steep2;
2209  setPrimalColumnPivotAlgorithm(steep2);
2210}
2211typedef struct {
2212  double optimizationDirection;
2213  double dblParam[ClpLastDblParam];
2214  double objectiveValue;
2215  double dualBound;
2216  double dualTolerance;
2217  double primalTolerance;
2218  double sumDualInfeasibilities;
2219  double sumPrimalInfeasibilities;
2220  double infeasibilityCost;
2221  int numberRows;
2222  int numberColumns;
2223  int intParam[ClpLastIntParam];
2224  int numberIterations;
2225  int problemStatus;
2226  int maximumIterations;
2227  int lengthNames;
2228  int numberDualInfeasibilities;
2229  int numberDualInfeasibilitiesWithoutFree;
2230  int numberPrimalInfeasibilities;
2231  int numberRefinements;
2232  int scalingFlag;
2233  int algorithm;
2234  unsigned int specialOptions;
2235  int dualPivotChoice;
2236  int primalPivotChoice;
2237  int matrixStorageChoice;
2238} Clp_scalars;
2239
2240int outDoubleArray(double * array, int length, FILE * fp)
2241{
2242  int numberWritten;
2243  if (array&&length) {
2244    numberWritten = fwrite(&length,sizeof(int),1,fp);
2245    if (numberWritten!=1)
2246      return 1;
2247    numberWritten = fwrite(array,sizeof(double),length,fp);
2248    if (numberWritten!=length)
2249      return 1;
2250  } else {
2251    length = 0;
2252    numberWritten = fwrite(&length,sizeof(int),1,fp);
2253    if (numberWritten!=1)
2254      return 1;
2255  }
2256  return 0;
2257}
2258// Save model to file, returns 0 if success
2259int
2260ClpSimplex::saveModel(const char * fileName)
2261{
2262  FILE * fp = fopen(fileName,"wb");
2263  if (fp) {
2264    Clp_scalars scalars;
2265    int i, numberWritten;
2266    // Fill in scalars
2267    scalars.optimizationDirection = optimizationDirection_;
2268    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2269    scalars.objectiveValue = objectiveValue_;
2270    scalars.dualBound = dualBound_;
2271    scalars.dualTolerance = dualTolerance_;
2272    scalars.primalTolerance = primalTolerance_;
2273    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2274    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2275    scalars.infeasibilityCost = infeasibilityCost_;
2276    scalars.numberRows = numberRows_;
2277    scalars.numberColumns = numberColumns_;
2278    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2279    scalars.numberIterations = numberIterations_;
2280    scalars.problemStatus = problemStatus_;
2281    scalars.maximumIterations = maximumIterations_;
2282    scalars.lengthNames = lengthNames_;
2283    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2284    scalars.numberDualInfeasibilitiesWithoutFree
2285      = numberDualInfeasibilitiesWithoutFree_;
2286    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2287    scalars.numberRefinements = numberRefinements_;
2288    scalars.scalingFlag = scalingFlag_;
2289    scalars.algorithm = algorithm_;
2290    scalars.specialOptions = specialOptions_;
2291    scalars.dualPivotChoice = dualRowPivot_->type();
2292    scalars.primalPivotChoice = primalColumnPivot_->type();
2293    scalars.matrixStorageChoice = matrix_->type();
2294
2295    // put out scalars
2296    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2297    if (numberWritten!=1)
2298      return 1;
2299    // strings
2300    int length;
2301    for (i=0;i<ClpLastStrParam;i++) {
2302      length = strParam_[i].size();
2303      numberWritten = fwrite(&length,sizeof(int),1,fp);
2304      if (numberWritten!=1)
2305        return 1;
2306      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2307      if (numberWritten!=1)
2308        return 1;
2309    }
2310    // arrays - in no particular order
2311    if (outDoubleArray(rowActivity_,numberRows_,fp))
2312        return 1;
2313    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2314        return 1;
2315    if (outDoubleArray(dual_,numberRows_,fp))
2316        return 1;
2317    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2318        return 1;
2319    if (outDoubleArray(rowLower_,numberRows_,fp))
2320        return 1;
2321    if (outDoubleArray(rowUpper_,numberRows_,fp))
2322        return 1;
2323    if (outDoubleArray(objective_,numberColumns_,fp))
2324        return 1;
2325    if (outDoubleArray(rowObjective_,numberRows_,fp))
2326        return 1;
2327    if (outDoubleArray(columnLower_,numberColumns_,fp))
2328        return 1;
2329    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2330        return 1;
2331    if (ray_) {
2332      if (problemStatus_==1)
2333        if (outDoubleArray(ray_,numberRows_,fp))
2334          return 1;
2335      else if (problemStatus_==2)
2336        if (outDoubleArray(ray_,numberColumns_,fp))
2337          return 1;
2338      else
2339        if (outDoubleArray(NULL,0,fp))
2340          return 1;
2341    } else {
2342      if (outDoubleArray(NULL,0,fp))
2343        return 1;
2344    }
2345    if (status_&&(numberRows_+numberColumns_)>0) {
2346      length = numberRows_+numberColumns_;
2347      numberWritten = fwrite(&length,sizeof(int),1,fp);
2348      if (numberWritten!=1)
2349        return 1;
2350      numberWritten = fwrite(status_,sizeof(char),length, fp);
2351      if (numberWritten!=length)
2352        return 1;
2353    } else {
2354      length = 0;
2355      numberWritten = fwrite(&length,sizeof(int),1,fp);
2356      if (numberWritten!=1)
2357        return 1;
2358    }
2359    if (lengthNames_) {
2360      assert (numberRows_ == (int) rowNames_.size());
2361      for (i=0;i<numberRows_;i++) {
2362        length = rowNames_[i].size();
2363        numberWritten = fwrite(&length,sizeof(int),1,fp);
2364        if (numberWritten!=1)
2365          return 1;
2366        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
2367        if (numberWritten!=1)
2368          return 1;
2369      }
2370      assert (numberColumns_ == (int) columnNames_.size());
2371      for (i=0;i<numberColumns_;i++) {
2372        length = columnNames_[i].size();
2373        numberWritten = fwrite(&length,sizeof(int),1,fp);
2374        if (numberWritten!=1)
2375          return 1;
2376        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
2377        if (numberWritten!=1)
2378          return 1;
2379      }
2380    }
2381    // just standard type at present
2382    assert (matrix_->type()==1);
2383    assert (matrix_->getNumCols() == numberColumns_);
2384    assert (matrix_->getNumRows() == numberRows_);
2385    // we are going to save with gaps
2386    length = matrix_->getVectorStarts()[numberColumns_-1]
2387      + matrix_->getVectorLengths()[numberColumns_-1];
2388    numberWritten = fwrite(&length,sizeof(int),1,fp);
2389    if (numberWritten!=1)
2390      return 1;
2391    numberWritten = fwrite(matrix_->getElements(),
2392                           sizeof(double),length,fp);
2393    if (numberWritten!=length)
2394      return 1;
2395    numberWritten = fwrite(matrix_->getIndices(),
2396                           sizeof(int),length,fp);
2397    if (numberWritten!=length)
2398      return 1;
2399    numberWritten = fwrite(matrix_->getVectorStarts(),
2400                           sizeof(int),numberColumns_,fp);
2401    if (numberWritten!=numberColumns_)
2402      return 1;
2403    numberWritten = fwrite(matrix_->getVectorLengths(),
2404                           sizeof(int),numberColumns_,fp);
2405    if (numberWritten!=numberColumns_)
2406      return 1;
2407    // finished
2408    fclose(fp);
2409    return 0;
2410  } else {
2411    return -1;
2412  }
2413}
2414
2415int inDoubleArray(double * &array, int length, FILE * fp)
2416{
2417  int numberRead;
2418  int length2;
2419  numberRead = fread(&length2,sizeof(int),1,fp);
2420  if (numberRead!=1)
2421    return 1;
2422  if (length2) {
2423    // lengths must match
2424    if (length!=length2)
2425      return 2;
2426    array = new double[length];
2427    numberRead = fread(array,sizeof(double),length,fp);
2428    if (numberRead!=length)
2429      return 1;
2430  } 
2431  return 0;
2432}
2433/* Restore model from file, returns 0 if success,
2434   deletes current model */
2435int 
2436ClpSimplex::restoreModel(const char * fileName)
2437{
2438  FILE * fp = fopen(fileName,"rb");
2439  if (fp) {
2440    // Get rid of current model
2441    ClpModel::gutsOfDelete();
2442    gutsOfDelete(0);
2443    int i;
2444    for (i=0;i<6;i++) {
2445      rowArray_[i]=NULL;
2446      columnArray_[i]=NULL;
2447    }
2448    // get an empty factorization so we can set tolerances etc
2449    factorization_ = new ClpFactorization();
2450    // Say sparse
2451    factorization_->sparseThreshold(1);
2452    Clp_scalars scalars;
2453    int numberRead;
2454
2455    // get scalars
2456    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
2457    if (numberRead!=1)
2458      return 1;
2459    // Fill in scalars
2460    optimizationDirection_ = scalars.optimizationDirection;
2461    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
2462    objectiveValue_ = scalars.objectiveValue;
2463    dualBound_ = scalars.dualBound;
2464    dualTolerance_ = scalars.dualTolerance;
2465    primalTolerance_ = scalars.primalTolerance;
2466    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
2467    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
2468    infeasibilityCost_ = scalars.infeasibilityCost;
2469    numberRows_ = scalars.numberRows;
2470    numberColumns_ = scalars.numberColumns;
2471    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
2472    numberIterations_ = scalars.numberIterations;
2473    problemStatus_ = scalars.problemStatus;
2474    maximumIterations_ = scalars.maximumIterations;
2475    lengthNames_ = scalars.lengthNames;
2476    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
2477    numberDualInfeasibilitiesWithoutFree_
2478      = scalars.numberDualInfeasibilitiesWithoutFree;
2479    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
2480    numberRefinements_ = scalars.numberRefinements;
2481    scalingFlag_ = scalars.scalingFlag;
2482    algorithm_ = scalars.algorithm;
2483    specialOptions_ = scalars.specialOptions;
2484    // strings
2485    int length;
2486    for (i=0;i<ClpLastStrParam;i++) {
2487      numberRead = fread(&length,sizeof(int),1,fp);
2488      if (numberRead!=1)
2489        return 1;
2490      char * array = new char[length+1];
2491      numberRead = fread(array,length,1,fp);
2492      if (numberRead!=1)
2493        return 1;
2494      array[length]='\0';
2495      strParam_[i]=array;
2496      delete [] array;
2497    }
2498    // arrays - in no particular order
2499    if (inDoubleArray(rowActivity_,numberRows_,fp))
2500        return 1;
2501    if (inDoubleArray(columnActivity_,numberColumns_,fp))
2502        return 1;
2503    if (inDoubleArray(dual_,numberRows_,fp))
2504        return 1;
2505    if (inDoubleArray(reducedCost_,numberColumns_,fp))
2506        return 1;
2507    if (inDoubleArray(rowLower_,numberRows_,fp))
2508        return 1;
2509    if (inDoubleArray(rowUpper_,numberRows_,fp))
2510        return 1;
2511    if (inDoubleArray(objective_,numberColumns_,fp))
2512        return 1;
2513    if (inDoubleArray(rowObjective_,numberRows_,fp))
2514        return 1;
2515    if (inDoubleArray(columnLower_,numberColumns_,fp))
2516        return 1;
2517    if (inDoubleArray(columnUpper_,numberColumns_,fp))
2518        return 1;
2519    if (problemStatus_==1) {
2520      if (inDoubleArray(ray_,numberRows_,fp))
2521        return 1;
2522    } else if (problemStatus_==2) {
2523      if (inDoubleArray(ray_,numberColumns_,fp))
2524        return 1;
2525    } else {
2526      // ray should be null
2527      numberRead = fread(&length,sizeof(int),1,fp);
2528      if (numberRead!=1)
2529        return 1;
2530      if (length)
2531        return 2;
2532    }
2533    delete [] status_;
2534    status_=NULL;
2535    // status region
2536    numberRead = fread(&length,sizeof(int),1,fp);
2537    if (numberRead!=1)
2538        return 1;
2539    if (length) {
2540      if (length!=numberRows_+numberColumns_)
2541        return 1;
2542      status_ = new char unsigned[length];
2543      numberRead = fread(status_,sizeof(char),length, fp);
2544      if (numberRead!=length)
2545        return 1;
2546    }
2547    if (lengthNames_) {
2548      char * array = new char[lengthNames_+1];
2549      rowNames_.resize(0);
2550      for (i=0;i<numberRows_;i++) {
2551        numberRead = fread(&length,sizeof(int),1,fp);
2552        if (numberRead!=1)
2553          return 1;
2554        numberRead = fread(array,length,1,fp);
2555        if (numberRead!=1)
2556          return 1;
2557        rowNames_[i]=array;
2558      }
2559      columnNames_.resize(0);
2560      for (i=0;i<numberColumns_;i++) {
2561        numberRead = fread(&length,sizeof(int),1,fp);
2562        if (numberRead!=1)
2563          return 1;
2564        numberRead = fread(array,length,1,fp);
2565        if (numberRead!=1)
2566          return 1;
2567        columnNames_[i]=array;
2568      }
2569      delete [] array;
2570    }
2571    // Pivot choices
2572    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
2573    delete dualRowPivot_;
2574    switch ((scalars.dualPivotChoice&63)) {
2575    case 1:
2576      // Dantzig
2577      dualRowPivot_ = new ClpDualRowDantzig();
2578      break;
2579    case 2:
2580      // Steepest - use mode
2581      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
2582      break;
2583    default:
2584      abort();
2585    }
2586    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
2587    delete primalColumnPivot_;
2588    switch ((scalars.primalPivotChoice&63)) {
2589    case 1:
2590      // Dantzig
2591      primalColumnPivot_ = new ClpPrimalColumnDantzig();
2592      break;
2593    case 2:
2594      // Steepest - use mode
2595      primalColumnPivot_
2596        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
2597      break;
2598    default:
2599      abort();
2600    }
2601    assert(scalars.matrixStorageChoice==1);
2602    delete matrix_;
2603    // get arrays
2604    numberRead = fread(&length,sizeof(int),1,fp);
2605    if (numberRead!=1)
2606      return 1;
2607    double * elements = new double[length];
2608    int * indices = new int[length];
2609    int * starts = new int[numberColumns_];
2610    int * lengths = new int[numberColumns_];
2611    numberRead = fread(elements, sizeof(double),length,fp);
2612    if (numberRead!=length)
2613      return 1;
2614    numberRead = fread(indices, sizeof(int),length,fp);
2615    if (numberRead!=length)
2616      return 1;
2617    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
2618    if (numberRead!=numberColumns_)
2619      return 1;
2620    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
2621    if (numberRead!=numberColumns_)
2622      return 1;
2623    // assign matrix
2624    CoinPackedMatrix * matrix = new CoinPackedMatrix();
2625    matrix->assignMatrix(true, numberRows_, numberColumns_,
2626                         length, elements, indices, starts, lengths);
2627    // and transfer to Clp
2628    matrix_ = new ClpPackedMatrix(matrix);
2629    // finished
2630    fclose(fp);
2631    return 0;
2632  } else {
2633    return -1;
2634  }
2635  return 0;
2636}
2637// value of incoming variable (in Dual)
2638double 
2639ClpSimplex::valueIncomingDual() const
2640{
2641  // Need value of incoming for list of infeasibilities as may be infeasible
2642  double valueIncoming = (dualOut_/alpha_)*directionOut_;
2643  if (directionIn_==-1)
2644    valueIncoming = upperIn_-valueIncoming;
2645  else
2646    valueIncoming = lowerIn_-valueIncoming;
2647  return valueIncoming;
2648}
2649//#############################################################################
2650
2651ClpSimplexProgress::ClpSimplexProgress () 
2652{
2653  int i;
2654  for (i=0;i<CLP_PROGRESS;i++) {
2655    objective_[i] = 0.0;
2656    infeasibility_[i] = -1.0; // set to an impossible value
2657    numberInfeasibilities_[i]=-1; 
2658  }
2659  numberTimes_ = 0;
2660  numberBadTimes_ = 0;
2661  model_ = NULL;
2662}
2663
2664
2665//-----------------------------------------------------------------------------
2666
2667ClpSimplexProgress::~ClpSimplexProgress ()
2668{
2669}
2670// Copy constructor.
2671ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2672{
2673  int i;
2674  for (i=0;i<CLP_PROGRESS;i++) {
2675    objective_[i] = rhs.objective_[i];
2676    infeasibility_[i] = rhs.infeasibility_[i];
2677    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2678  }
2679  numberTimes_ = rhs.numberTimes_;
2680  numberBadTimes_ = rhs.numberBadTimes_;
2681  model_ = rhs.model_;
2682}
2683// Copy constructor.from model
2684ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2685{
2686  int i;
2687  for (i=0;i<CLP_PROGRESS;i++) {
2688    objective_[i] = 0.0;
2689    infeasibility_[i] = -1.0; // set to an impossible value
2690    numberInfeasibilities_[i]=-1; 
2691  }
2692  numberTimes_ = 0;
2693  numberBadTimes_ = 0;
2694  model_ = model;
2695}
2696// Assignment operator. This copies the data
2697ClpSimplexProgress & 
2698ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2699{
2700  if (this != &rhs) {
2701    int i;
2702    for (i=0;i<CLP_PROGRESS;i++) {
2703      objective_[i] = rhs.objective_[i];
2704      infeasibility_[i] = rhs.infeasibility_[i];
2705      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2706    }
2707    numberTimes_ = rhs.numberTimes_;
2708    numberBadTimes_ = rhs.numberBadTimes_;
2709    model_ = rhs.model_;
2710  }
2711  return *this;
2712}
2713int
2714ClpSimplexProgress::looping()
2715{
2716  assert(model_);
2717  double objective = model_->objectiveValue();
2718  double infeasibility;
2719  int numberInfeasibilities;
2720  if (model_->algorithm()<0) {
2721    // dual
2722    infeasibility = model_->sumPrimalInfeasibilities();
2723    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2724  } else {
2725    //primal
2726    infeasibility = model_->sumDualInfeasibilities();
2727    numberInfeasibilities = model_->numberDualInfeasibilities();
2728  }
2729  int i;
2730  int numberMatched=0;
2731  int matched=0;
2732
2733  for (i=0;i<CLP_PROGRESS;i++) {
2734    if (objective==objective_[i]&&
2735        infeasibility==infeasibility_[i]&&
2736        numberInfeasibilities==numberInfeasibilities_[i]) {
2737      matched |= (1<<i);
2738      numberMatched++;
2739    }
2740    if (i) {
2741      objective_[i-1] = objective_[i];
2742      infeasibility_[i-1] = infeasibility_[i];
2743      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2744    }
2745  }
2746  objective_[CLP_PROGRESS-1] = objective;
2747  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2748  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2749  if (model_->progressFlag())
2750    numberMatched=0;
2751  numberTimes_++;
2752  if (numberTimes_<10)
2753    numberMatched=0;
2754  // skip if just last time as may be checking something
2755  if (matched==(1<<(CLP_PROGRESS-1)))
2756    numberMatched=0;
2757  if (numberMatched) {
2758    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2759      <<numberMatched
2760      <<matched
2761      <<numberTimes_
2762      <<CoinMessageEol;
2763    numberBadTimes_++;
2764    if (numberBadTimes_<10) {
2765      if (model_->algorithm()<0) {
2766        // dual - change tolerance
2767        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2768      } else {
2769        // primal - change tolerance
2770        model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2771      }
2772    } else {
2773      model_->messageHandler()->message(CLP_LOOP,model_->messages())
2774        <<CoinMessageEol;
2775      // look at solution and maybe declare victory
2776      if (infeasibility<1.0e-4)
2777        return 0;
2778      else
2779        return 3;
2780    }
2781  }
2782  return -1;
2783}
2784// Sanity check on input data - returns true if okay
2785bool 
2786ClpSimplex::sanityCheck()
2787{
2788  int numberBad , numberBadBounds;
2789  double largestBound, smallestBound, minimumGap;
2790  double smallestObj, largestObj;
2791  int firstBad;
2792  int modifiedBounds=0;
2793  int i;
2794  numberBad=0;
2795  numberBadBounds=0;
2796  firstBad=-1;
2797  minimumGap=1.0e100;
2798  smallestBound=1.0e100;
2799  largestBound=0.0;
2800  smallestObj=1.0e100;
2801  largestObj=0.0;
2802  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
2803    double value;
2804    value = fabs(cost_[i]);
2805    if (value>1.0e50) {
2806      numberBad++;
2807      if (firstBad<0)
2808        firstBad=i;
2809    } else if (value) {
2810      if (value>largestObj)
2811        largestObj=value;
2812      if (value<smallestObj)
2813        smallestObj=value;
2814    }
2815    value=upper_[i]-lower_[i];
2816    if (value<-primalTolerance_) {
2817      numberBadBounds++;
2818      if (firstBad<0)
2819        firstBad=i;
2820    } else if (value<=primalTolerance_) {
2821      if (value) {
2822        // modify
2823        upper_[i] = lower_[i];
2824        modifiedBounds++;
2825      }
2826    } else {
2827      if (value<minimumGap)
2828        minimumGap=value;
2829    }
2830    if (lower_[i]>-1.0e100&&lower_[i]) {
2831      value = fabs(lower_[i]);
2832      if (value>largestBound)
2833        largestBound=value;
2834      if (value<smallestBound)
2835        smallestBound=value;
2836    }
2837    if (upper_[i]<1.0e100&&upper_[i]) {
2838      value = fabs(upper_[i]);
2839      if (value>largestBound)
2840        largestBound=value;
2841      if (value<smallestBound)
2842        smallestBound=value;
2843    }
2844  }
2845  if (largestBound)
2846    handler_->message(CLP_RIMSTATISTICS3,messages_)
2847      <<smallestBound
2848      <<largestBound
2849      <<minimumGap
2850      <<CoinMessageEol;
2851  minimumGap=1.0e100;
2852  smallestBound=1.0e100;
2853  largestBound=0.0;
2854  for (i=0;i<numberColumns_;i++) {
2855    double value;
2856    value = fabs(cost_[i]);
2857    if (value>1.0e50) {
2858      numberBad++;
2859      if (firstBad<0)
2860        firstBad=i;
2861    } else if (value) {
2862      if (value>largestObj)
2863        largestObj=value;
2864      if (value<smallestObj)
2865        smallestObj=value;
2866    }
2867    value=upper_[i]-lower_[i];
2868    if (value<-primalTolerance_) {
2869      numberBadBounds++;
2870      if (firstBad<0)
2871        firstBad=i;
2872    } else if (value<=primalTolerance_) {
2873      if (value) {
2874        // modify
2875        upper_[i] = lower_[i];
2876        modifiedBounds++;
2877      }
2878    } else {
2879      if (value<minimumGap)
2880        minimumGap=value;
2881    }
2882    if (lower_[i]>-1.0e100&&lower_[i]) {
2883      value = fabs(lower_[i]);
2884      if (value>largestBound)
2885        largestBound=value;
2886      if (value<smallestBound)
2887        smallestBound=value;
2888    }
2889    if (upper_[i]<1.0e100&&upper_[i]) {
2890      value = fabs(upper_[i]);
2891      if (value>largestBound)
2892        largestBound=value;
2893      if (value<smallestBound)
2894        smallestBound=value;
2895    }
2896  }
2897  char rowcol[]={'R','C'};
2898  if (numberBad) {
2899    handler_->message(CLP_BAD_BOUNDS,messages_)
2900      <<numberBad
2901      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
2902      <<CoinMessageEol;
2903    problemStatus_=4;
2904    return false;
2905  }
2906  if (modifiedBounds)
2907    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
2908      <<modifiedBounds
2909      <<CoinMessageEol;
2910  handler_->message(CLP_RIMSTATISTICS1,messages_)
2911    <<smallestObj
2912    <<largestObj
2913    <<CoinMessageEol;
2914  if (largestBound)
2915    handler_->message(CLP_RIMSTATISTICS2,messages_)
2916      <<smallestBound
2917      <<largestBound
2918      <<minimumGap
2919      <<CoinMessageEol;
2920  return true;
2921}
2922// Set up status array (for OsiClp)
2923void 
2924ClpSimplex::createStatus() 
2925{
2926  delete [] status_;
2927  status_ = new unsigned char [numberColumns_+numberRows_];
2928  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2929  int i;
2930  // set column status to one nearest zero
2931  for (i=0;i<numberColumns_;i++) {
2932#if 0
2933    if (columnLower_[i]>=0.0) {
2934      setColumnStatus(i,atLowerBound);
2935    } else if (columnUpper_[i]<=0.0) {
2936      setColumnStatus(i,atUpperBound);
2937    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
2938      // free
2939      setColumnStatus(i,isFree);
2940    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
2941      setColumnStatus(i,atLowerBound);
2942    } else {
2943      setColumnStatus(i,atUpperBound);
2944    }
2945#else
2946    setColumnStatus(i,atLowerBound);
2947#endif
2948  }
2949  for (i=0;i<numberRows_;i++) {
2950    setRowStatus(i,basic);
2951  }
2952}
2953/* Loads a problem (the constraints on the
2954   rows are given by lower and upper bounds). If a pointer is 0 then the
2955   following values are the default:
2956   <ul>
2957   <li> <code>colub</code>: all columns have upper bound infinity
2958   <li> <code>collb</code>: all columns have lower bound 0
2959   <li> <code>rowub</code>: all rows have upper bound infinity
2960   <li> <code>rowlb</code>: all rows have lower bound -infinity
2961   <li> <code>obj</code>: all variables have 0 objective coefficient
2962   </ul>
2963*/
2964void 
2965ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
2966                    const double* collb, const double* colub,   
2967                    const double* obj,
2968                    const double* rowlb, const double* rowub,
2969                    const double * rowObjective)
2970{
2971  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
2972                        rowObjective);
2973  createStatus();
2974}
2975void 
2976ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
2977                    const double* collb, const double* colub,   
2978                    const double* obj,
2979                    const double* rowlb, const double* rowub,
2980                    const double * rowObjective)
2981{
2982  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
2983                        rowObjective);
2984  createStatus();
2985}
2986
2987/* Just like the other loadProblem() method except that the matrix is
2988   given in a standard column major ordered format (without gaps). */
2989void 
2990ClpSimplex::loadProblem (  const int numcols, const int numrows,
2991                    const int* start, const int* index,
2992                    const double* value,
2993                    const double* collb, const double* colub,   
2994                    const double* obj,
2995                    const double* rowlb, const double* rowub,
2996                    const double * rowObjective)
2997{
2998  ClpModel::loadProblem(numcols, numrows, start, index, value,
2999                          collb, colub, obj, rowlb, rowub,
3000                          rowObjective);
3001  createStatus();
3002}
3003void 
3004ClpSimplex::loadProblem (  const int numcols, const int numrows,
3005                           const int* start, const int* index,
3006                           const double* value,const int * length,
3007                           const double* collb, const double* colub,   
3008                           const double* obj,
3009                           const double* rowlb, const double* rowub,
3010                           const double * rowObjective)
3011{
3012  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3013                          collb, colub, obj, rowlb, rowub,
3014                          rowObjective);
3015  createStatus();
3016}
3017// Read an mps file from the given filename
3018int 
3019ClpSimplex::readMps(const char *filename,
3020            bool keepNames,
3021            bool ignoreErrors)
3022{
3023  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3024  createStatus();
3025  return status;
3026}
Note: See TracBrowser for help on using the repository browser.