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

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

Presolve nearly ready

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