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

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

For presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 86.9 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  if ((what&16)!=0) {
1680    // move information to work arrays
1681    // row reduced costs
1682    if (!dj_) {
1683      dj_ = new double[numberRows_+numberColumns_];
1684      reducedCostWork_ = dj_;
1685      rowReducedCost_ = dj_+numberColumns_;
1686      memcpy(reducedCostWork_,reducedCost_,numberColumns_*sizeof(double));
1687      memcpy(rowReducedCost_,dual_,numberRows_*sizeof(double));
1688    }
1689    if (!solution_) {
1690      solution_ = new double[numberRows_+numberColumns_];
1691      columnActivityWork_ = solution_;
1692      rowActivityWork_ = solution_+numberColumns_;
1693      memcpy(columnActivityWork_,columnActivity_,
1694             numberColumns_*sizeof(double));
1695      memcpy(rowActivityWork_,rowActivity_,
1696             numberRows_*sizeof(double));
1697    }
1698    //check matrix
1699    if (!matrix_->allElementsInRange(this,1.0e-20,1.0e20)) {
1700      problemStatus_=4;
1701      goodMatrix= false;
1702    }
1703    if (makeRowCopy) {
1704      delete rowCopy_;
1705      // may return NULL if can't give row copy
1706      rowCopy_ = matrix_->reverseOrderedCopy();
1707    }
1708  }
1709  int i;
1710  if ((what&4)!=0) {
1711    delete [] cost_;
1712    cost_ = new double[numberColumns_+numberRows_];
1713    objectiveWork_ = cost_;
1714    rowObjectiveWork_ = cost_+numberColumns_;
1715    memcpy(objectiveWork_,objective_,numberColumns_*sizeof(double));
1716    if (rowObjective_)
1717      memcpy(rowObjectiveWork_,rowObjective_,numberRows_*sizeof(double));
1718    else
1719      memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
1720  }
1721  if ((what&1)!=0) {
1722    delete [] lower_;
1723    delete [] upper_;
1724    lower_ = new double[numberColumns_+numberRows_];
1725    upper_ = new double[numberColumns_+numberRows_];
1726    rowLowerWork_ = lower_+numberColumns_;
1727    columnLowerWork_ = lower_;
1728    rowUpperWork_ = upper_+numberColumns_;
1729    columnUpperWork_ = upper_;
1730    memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
1731    memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
1732    memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
1733    memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
1734    // clean up any mismatches on infinity
1735    for (i=0;i<numberColumns_;i++) {
1736      if (columnLowerWork_[i]<-1.0e30)
1737        columnLowerWork_[i] = -DBL_MAX;
1738      if (columnUpperWork_[i]>1.0e30)
1739        columnUpperWork_[i] = DBL_MAX;
1740    }
1741    // clean up any mismatches on infinity
1742    for (i=0;i<numberRows_;i++) {
1743      if (rowLowerWork_[i]<-1.0e30)
1744        rowLowerWork_[i] = -DBL_MAX;
1745      if (rowUpperWork_[i]>1.0e30)
1746        rowUpperWork_[i] = DBL_MAX;
1747    }
1748  }
1749  // do scaling if needed
1750  if (scalingFlag_>0&&!rowScale_) {
1751    if (matrix_->scale(this))
1752      scalingFlag_=-scalingFlag_; // not scaled after all
1753  }
1754  if ((what&4)!=0) {
1755    double direction = optimizationDirection_;
1756    // direction is actually scale out not scale in
1757    if (direction)
1758      direction = 1.0/direction;
1759    // but also scale by scale factors
1760    // not really sure about row scaling
1761    if (!rowScale_) {
1762      if (direction!=1.0) {
1763        for (i=0;i<numberRows_;i++)
1764          rowObjectiveWork_[i] *= direction;
1765        for (i=0;i<numberColumns_;i++)
1766          objectiveWork_[i] *= direction;
1767      }
1768    } else {
1769      for (i=0;i<numberRows_;i++)
1770        rowObjectiveWork_[i] *= direction/rowScale_[i];
1771      for (i=0;i<numberColumns_;i++)
1772        objectiveWork_[i] *= direction*columnScale_[i];
1773    }
1774  }
1775  if ((what&1)!=0&&rowScale_) {
1776    for (i=0;i<numberColumns_;i++) {
1777      double multiplier = 1.0/columnScale_[i];
1778      if (columnLowerWork_[i]>-1.0e50)
1779        columnLowerWork_[i] *= multiplier;
1780      if (columnUpperWork_[i]<1.0e50)
1781        columnUpperWork_[i] *= multiplier;
1782    }
1783    for (i=0;i<numberRows_;i++) {
1784      if (rowLowerWork_[i]>-1.0e50)
1785        rowLowerWork_[i] *= rowScale_[i];
1786      if (rowUpperWork_[i]<1.0e50)
1787        rowUpperWork_[i] *= rowScale_[i];
1788    }
1789  }
1790  if ((what&8)!=0&&rowScale_) {
1791    // on entry
1792    for (i=0;i<numberColumns_;i++) {
1793      columnActivityWork_[i] /= columnScale_[i];
1794      reducedCostWork_[i] *= columnScale_[i];
1795    }
1796    for (i=0;i<numberRows_;i++) {
1797      rowActivityWork_[i] *= rowScale_[i];
1798      dual_[i] /= rowScale_[i];
1799      rowReducedCost_[i] = dual_[i];
1800    }
1801  }
1802 
1803  if ((what&16)!=0) {
1804    // check rim of problem okay
1805    if (!sanityCheck())
1806      goodMatrix=false;
1807  } 
1808  // we need to treat matrix as if each element by rowScaleIn and columnScaleout??
1809  // maybe we need to move scales to SimplexModel for factorization?
1810  if ((what&8)!=0&&!pivotVariable_) {
1811    pivotVariable_=new int[numberRows_];
1812  }
1813
1814  if ((what&16)!=0&&!rowArray_[2]) {
1815    // get some arrays
1816    int iRow,iColumn;
1817    // these are "indexed" arrays so we always know where nonzeros are
1818    /**********************************************************
1819      rowArray_[3] is long enough for columns as well
1820    *********************************************************/
1821    for (iRow=0;iRow<4;iRow++) {
1822      if (!rowArray_[iRow]) {
1823        rowArray_[iRow]=new CoinIndexedVector();
1824        int length =numberRows_+factorization_->maximumPivots();
1825        if (iRow==3)
1826          length = max(length,numberColumns_);
1827        rowArray_[iRow]->reserve(length);
1828      }
1829    }
1830   
1831    for (iColumn=0;iColumn<2;iColumn++) {
1832      if (!columnArray_[iColumn]) {
1833        columnArray_[iColumn]=new CoinIndexedVector();
1834        columnArray_[iColumn]->reserve(numberColumns_);
1835      }
1836    }   
1837  }
1838  return goodMatrix;
1839}
1840void
1841ClpSimplex::deleteRim(bool getRidOfFactorizationData)
1842{
1843  int i;
1844  if (problemStatus_!=1&&problemStatus_!=2) {
1845    delete [] ray_;
1846    ray_=NULL;
1847  }
1848  if (rowScale_) {
1849    for (i=0;i<numberColumns_;i++) {
1850      columnActivity_[i] = columnActivityWork_[i]*columnScale_[i];
1851      reducedCost_[i] = reducedCostWork_[i]/columnScale_[i];
1852    }
1853    for (i=0;i<numberRows_;i++) {
1854      rowActivity_[i] = rowActivityWork_[i]/rowScale_[i];
1855      dual_[i] *= rowScale_[i];
1856    }
1857    if (problemStatus_==2) {
1858      for (i=0;i<numberColumns_;i++) {
1859        ray_[i] *= columnScale_[i];
1860      }
1861    } else if (problemStatus_==1) {
1862      for (i=0;i<numberRows_;i++) {
1863        ray_[i] *= rowScale_[i];
1864      }
1865    }
1866  } else {
1867    for (i=0;i<numberColumns_;i++) {
1868      columnActivity_[i] = columnActivityWork_[i];
1869      reducedCost_[i] = reducedCostWork_[i];
1870    }
1871    for (i=0;i<numberRows_;i++) {
1872      rowActivity_[i] = rowActivityWork_[i];
1873    }
1874  }
1875  // direction may have been modified by scaling - clean up
1876  if (optimizationDirection_>0.0)
1877    optimizationDirection_ = 1.0;
1878  else if (optimizationDirection_<0.0)
1879    optimizationDirection_ = -1.0;
1880  // scaling may have been turned off
1881  scalingFlag_ = abs(scalingFlag_);
1882  if(getRidOfFactorizationData)
1883    gutsOfDelete(2);
1884  else
1885    gutsOfDelete(1);
1886}
1887void 
1888ClpSimplex::setDualBound(double value)
1889{
1890  if (value>0.0)
1891    dualBound_=value;
1892}
1893void 
1894ClpSimplex::setInfeasibilityCost(double value)
1895{
1896  if (value>0.0)
1897    infeasibilityCost_=value;
1898}
1899void ClpSimplex::setnumberRefinements( int value) 
1900{
1901  if (value>=0&&value<10)
1902    numberRefinements_=value;
1903}
1904// Sets row pivot choice algorithm in dual
1905void 
1906ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
1907{
1908  delete dualRowPivot_;
1909  dualRowPivot_ = choice.clone(true);
1910}
1911// Sets row pivot choice algorithm in dual
1912void 
1913ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
1914{
1915  delete primalColumnPivot_;
1916  primalColumnPivot_ = choice.clone(true);
1917}
1918// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
1919void 
1920ClpSimplex::scaling(int mode)
1921{
1922  if (mode>0&&mode<3) {
1923    scalingFlag_=mode;
1924  } else if (!mode) {
1925    scalingFlag_=0;
1926    delete [] rowScale_;
1927    rowScale_ = NULL;
1928    delete [] columnScale_;
1929    columnScale_ = NULL;
1930  }
1931}
1932// Passes in factorization
1933void 
1934ClpSimplex::setFactorization( ClpFactorization & factorization)
1935{
1936  delete factorization_;
1937  factorization_= new ClpFactorization(factorization);
1938}
1939void 
1940ClpSimplex::times(double scalar,
1941                  const double * x, double * y) const
1942{
1943  if (rowScale_)
1944    matrix_->times(scalar,x,y,rowScale_,columnScale_);
1945  else
1946    matrix_->times(scalar,x,y);
1947}
1948void 
1949ClpSimplex::transposeTimes(double scalar,
1950                           const double * x, double * y) const 
1951{
1952  if (rowScale_)
1953    matrix_->transposeTimes(scalar,x,y,rowScale_,columnScale_);
1954  else
1955    matrix_->transposeTimes(scalar,x,y);
1956}
1957/* Perturbation:
1958   -50 to +50 - perturb by this power of ten (-6 sounds good)
1959   100 - auto perturb if takes too long (1.0e-6 largest nonzero)
1960   101 - we are perturbed
1961   102 - don't try perturbing again
1962   default is 100
1963*/
1964void 
1965ClpSimplex::setPerturbation(int value)
1966{
1967  if(value<=100&&value >=-50) {
1968    perturbation_=value;
1969  } 
1970}
1971// Sparsity on or off
1972bool 
1973ClpSimplex::sparseFactorization() const
1974{
1975  return factorization_->sparseThreshold();
1976}
1977void 
1978ClpSimplex::setSparseFactorization(bool value)
1979{
1980  if (value) {
1981    if (!factorization_->sparseThreshold())
1982      factorization_->goSparse();
1983  } else {
1984    factorization_->sparseThreshold(0);
1985  }
1986}
1987/* Tightens primal bounds to make dual faster.  Unless
1988   fixed, bounds are slightly looser than they could be.
1989   This is to make dual go faster and is probably not needed
1990   with a presolve.  Returns non-zero if problem infeasible
1991*/
1992int 
1993ClpSimplex::tightenPrimalBounds()
1994{
1995 
1996  // Get a row copy in standard format
1997  CoinPackedMatrix copy;
1998  copy.reverseOrderedCopyOf(*matrix());
1999  // get matrix data pointers
2000  const int * column = copy.getIndices();
2001  const int * rowStart = copy.getVectorStarts();
2002  const int * rowLength = copy.getVectorLengths(); 
2003  const double * element = copy.getElements();
2004  int numberChanged=1,iPass=0;
2005  double large = largeValue(); // treat bounds > this as infinite
2006  int numberInfeasible=0;
2007  int totalTightened = 0;
2008
2009  double tolerance = primalTolerance();
2010
2011
2012  // Save column bounds
2013  double * saveLower = new double [numberColumns_];
2014  memcpy(saveLower,columnLower_,numberColumns_*sizeof(double));
2015  double * saveUpper = new double [numberColumns_];
2016  memcpy(saveUpper,columnUpper_,numberColumns_*sizeof(double));
2017#define MAXPASS 10 
2018
2019  // Loop round seeing if we can tighten bounds
2020  // Would be faster to have a stack of possible rows
2021  // and we put altered rows back on stack
2022
2023  int iRow, iColumn;
2024
2025  while(numberChanged) {
2026
2027    numberChanged = 0; // Bounds tightened this pass
2028   
2029    if (iPass==MAXPASS) break;
2030    iPass++;
2031   
2032    for (iRow = 0; iRow < numberRows_; iRow++) {
2033
2034      if (rowLower_[iRow]>-large||rowUpper_[iRow]<large) {
2035
2036        // possible row
2037        int infiniteUpper = 0;
2038        int infiniteLower = 0;
2039        double maximumUp = 0.0;
2040        double maximumDown = 0.0;
2041        double newBound;
2042        int rStart = rowStart[iRow];
2043        int rEnd = rowStart[iRow]+rowLength[iRow];
2044        int j;
2045
2046        // Compute possible lower and upper ranges
2047
2048        for (j = rStart; j < rEnd; ++j) {
2049          double value=element[j];
2050          iColumn = column[j];
2051          if (value > 0.0) {
2052            if (columnUpper_[iColumn] >= large) {
2053              maximumUp = 1e31;
2054              ++infiniteUpper;
2055            } else {
2056              maximumUp += columnUpper_[iColumn] * value;
2057            }
2058            if (columnLower_[iColumn] <= -large) {
2059              maximumDown = -1e31;
2060              ++infiniteLower;
2061            } else {
2062              maximumDown += columnLower_[iColumn] * value;
2063            }
2064          } else if (value<0.0) {
2065            if (columnUpper_[iColumn] >= large) {
2066              maximumDown = -1e31;
2067              ++infiniteLower;
2068            } else {
2069              maximumDown += columnUpper_[iColumn] * value;
2070            }
2071            if (columnLower_[iColumn] <= -large) {
2072              maximumUp = 1e31;
2073              ++infiniteUpper;
2074            } else {
2075              maximumUp += columnLower_[iColumn] * value;
2076            }
2077          }
2078        }
2079        if (maximumUp <= rowUpper_[iRow] + tolerance && 
2080            maximumDown >= rowLower_[iRow] - tolerance) {
2081
2082          // Row is redundant - make totally free
2083          rowLower_[iRow]=-DBL_MAX;
2084          rowUpper_[iRow]=DBL_MAX;
2085
2086        } else {
2087          if (maximumUp < rowLower_[iRow] -tolerance ||
2088              maximumDown > rowUpper_[iRow]+tolerance) {
2089            // problem is infeasible - exit at once
2090            numberInfeasible++;
2091            break;
2092          }
2093
2094          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
2095            for (j = rStart; j < rEnd; ++j) {
2096              double value=element[j];
2097              iColumn = column[j];
2098              if (value > 0.0) {
2099                if (columnUpper_[iColumn] < large) {
2100                  newBound = columnUpper_[iColumn] + 
2101                    (rowLower_[iRow] - maximumUp) / value;
2102                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2103                    // Tighten the lower bound
2104
2105                    columnLower_[iColumn] = newBound;
2106                    ++numberChanged;
2107
2108                    // check infeasible (relaxed)
2109                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2110                        -100.0*tolerance) 
2111                      numberInfeasible++;
2112                    infiniteLower=1; // skip looking at other bound
2113                  }
2114                }
2115              } else {
2116                if (columnLower_[iColumn] > -large) {
2117                  newBound = columnLower_[iColumn] + 
2118                    (rowLower_[iRow] - maximumUp) / value;
2119                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2120                    // Tighten the upper bound
2121
2122                    columnUpper_[iColumn] = newBound;
2123                    ++numberChanged;
2124
2125                    // check infeasible (relaxed)
2126                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2127                        -100.0*tolerance) 
2128                      numberInfeasible++;
2129                    infiniteLower=1; // skip looking at other bound
2130                  }
2131                }
2132              }
2133            }
2134          }
2135         
2136          // Try other way
2137          if (infiniteLower == 0 && rowUpper_[iRow] < large) {
2138            for (j = rStart; j < rEnd; ++j) {
2139              double value=element[j];
2140              iColumn = column[j];
2141              if (value < 0.0) {
2142                if (columnUpper_[iColumn] < large) {
2143                  newBound = columnUpper_[iColumn] + 
2144                    (rowUpper_[iRow] - maximumDown) / value;
2145                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
2146                    // Tighten the lower bound
2147
2148                    columnLower_[iColumn] = newBound;
2149                    ++numberChanged;
2150
2151                    // check infeasible (relaxed)
2152                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2153                        -100.0*tolerance) 
2154                      numberInfeasible++;
2155                  }
2156                } 
2157              } else {
2158                if (columnLower_[iColumn] > -large) {
2159                  newBound = columnLower_[iColumn] + 
2160                    (rowUpper_[iRow] - maximumDown) / value;
2161                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
2162                    // Tighten the upper bound
2163
2164                    columnUpper_[iColumn] = newBound;
2165                    ++numberChanged;
2166
2167                    // check infeasible (relaxed)
2168                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
2169                        -100.0*tolerance) 
2170                      numberInfeasible++;
2171                  }
2172                }
2173              }
2174            }
2175          }
2176        }
2177      }
2178    }
2179    totalTightened += numberChanged;
2180    if (numberInfeasible) break;
2181  }
2182  if (!numberInfeasible) {
2183    handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN,messages_)
2184      <<totalTightened
2185      <<CoinMessageEol;
2186    // Set bounds slightly loose
2187    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2188      if (columnUpper_[iColumn]-columnLower_[iColumn]<tolerance) {
2189        // fix
2190        if (fabs(columnUpper_[iColumn])<fabs(columnLower_[iColumn])) 
2191          columnLower_[iColumn]=columnUpper_[iColumn];
2192        else
2193          columnUpper_[iColumn]=columnLower_[iColumn];
2194      } else {
2195        if (columnUpper_[iColumn]<saveUpper[iColumn]) {
2196          // relax a bit
2197          columnUpper_[iColumn] = min(columnUpper_[iColumn]+100.0*tolerance,
2198                                      saveUpper[iColumn]);
2199        }
2200        if (columnLower_[iColumn]>saveLower[iColumn]) {
2201          // relax a bit
2202          columnLower_[iColumn] = max(columnLower_[iColumn]-100.0*tolerance,
2203                                      saveLower[iColumn]);
2204        }
2205      }
2206    }
2207  } else {
2208    handler_->message(CLP_SIMPLEX_INFEASIBILITIES,messages_)
2209      <<numberInfeasible
2210      <<CoinMessageEol;
2211    // restore column bounds
2212    memcpy(columnLower_,saveLower,numberColumns_*sizeof(double));
2213    memcpy(columnUpper_,saveUpper,numberColumns_*sizeof(double));
2214  }
2215  delete [] saveLower;
2216  delete [] saveUpper;
2217  return (numberInfeasible);
2218}
2219// dual
2220#include "ClpSimplexDual.hpp"
2221int ClpSimplex::dual ( )
2222{
2223  return ((ClpSimplexDual *) this)->dual();
2224}
2225// primal
2226#include "ClpSimplexPrimal.hpp"
2227int ClpSimplex::primal (int ifValuesPass )
2228{
2229  return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
2230}
2231/* For strong branching.  On input lower and upper are new bounds
2232   while on output they are objective function values (>1.0e50 infeasible).
2233   Return code is 0 if nothing interesting, -1 if infeasible both
2234   ways and +1 if infeasible one way (check values to see which one(s))
2235*/
2236int ClpSimplex::strongBranching(int numberVariables,const int * variables,
2237                                double * newLower, double * newUpper,
2238                                bool stopOnFirstInfeasible,
2239                                bool alwaysFinish)
2240{
2241  return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
2242                                                    newLower,  newUpper,
2243                                                    stopOnFirstInfeasible,
2244                                                    alwaysFinish);
2245}
2246/* Borrow model.  This is so we dont have to copy large amounts
2247   of data around.  It assumes a derived class wants to overwrite
2248   an empty model with a real one - while it does an algorithm.
2249   This is same as ClpModel one, but sets scaling on etc. */
2250void 
2251ClpSimplex::borrowModel(ClpModel & otherModel) 
2252{
2253  ClpModel::borrowModel(otherModel);
2254  createStatus();
2255  scaling();
2256  ClpDualRowSteepest steep1;
2257  setDualRowPivotAlgorithm(steep1);
2258  ClpPrimalColumnSteepest steep2;
2259  setPrimalColumnPivotAlgorithm(steep2);
2260}
2261typedef struct {
2262  double optimizationDirection;
2263  double dblParam[ClpLastDblParam];
2264  double objectiveValue;
2265  double dualBound;
2266  double dualTolerance;
2267  double primalTolerance;
2268  double sumDualInfeasibilities;
2269  double sumPrimalInfeasibilities;
2270  double infeasibilityCost;
2271  int numberRows;
2272  int numberColumns;
2273  int intParam[ClpLastIntParam];
2274  int numberIterations;
2275  int problemStatus;
2276  int maximumIterations;
2277  int lengthNames;
2278  int numberDualInfeasibilities;
2279  int numberDualInfeasibilitiesWithoutFree;
2280  int numberPrimalInfeasibilities;
2281  int numberRefinements;
2282  int scalingFlag;
2283  int algorithm;
2284  unsigned int specialOptions;
2285  int dualPivotChoice;
2286  int primalPivotChoice;
2287  int matrixStorageChoice;
2288} Clp_scalars;
2289
2290int outDoubleArray(double * array, int length, FILE * fp)
2291{
2292  int numberWritten;
2293  if (array&&length) {
2294    numberWritten = fwrite(&length,sizeof(int),1,fp);
2295    if (numberWritten!=1)
2296      return 1;
2297    numberWritten = fwrite(array,sizeof(double),length,fp);
2298    if (numberWritten!=length)
2299      return 1;
2300  } else {
2301    length = 0;
2302    numberWritten = fwrite(&length,sizeof(int),1,fp);
2303    if (numberWritten!=1)
2304      return 1;
2305  }
2306  return 0;
2307}
2308// Save model to file, returns 0 if success
2309int
2310ClpSimplex::saveModel(const char * fileName)
2311{
2312  FILE * fp = fopen(fileName,"wb");
2313  if (fp) {
2314    Clp_scalars scalars;
2315    int i, numberWritten;
2316    // Fill in scalars
2317    scalars.optimizationDirection = optimizationDirection_;
2318    memcpy(scalars.dblParam, dblParam_,ClpLastDblParam * sizeof(double));
2319    scalars.objectiveValue = objectiveValue_;
2320    scalars.dualBound = dualBound_;
2321    scalars.dualTolerance = dualTolerance_;
2322    scalars.primalTolerance = primalTolerance_;
2323    scalars.sumDualInfeasibilities = sumDualInfeasibilities_;
2324    scalars.sumPrimalInfeasibilities = sumPrimalInfeasibilities_;
2325    scalars.infeasibilityCost = infeasibilityCost_;
2326    scalars.numberRows = numberRows_;
2327    scalars.numberColumns = numberColumns_;
2328    memcpy(scalars.intParam, intParam_,ClpLastIntParam * sizeof(double));
2329    scalars.numberIterations = numberIterations_;
2330    scalars.problemStatus = problemStatus_;
2331    scalars.maximumIterations = maximumIterations_;
2332    scalars.lengthNames = lengthNames_;
2333    scalars.numberDualInfeasibilities = numberDualInfeasibilities_;
2334    scalars.numberDualInfeasibilitiesWithoutFree
2335      = numberDualInfeasibilitiesWithoutFree_;
2336    scalars.numberPrimalInfeasibilities = numberPrimalInfeasibilities_;
2337    scalars.numberRefinements = numberRefinements_;
2338    scalars.scalingFlag = scalingFlag_;
2339    scalars.algorithm = algorithm_;
2340    scalars.specialOptions = specialOptions_;
2341    scalars.dualPivotChoice = dualRowPivot_->type();
2342    scalars.primalPivotChoice = primalColumnPivot_->type();
2343    scalars.matrixStorageChoice = matrix_->type();
2344
2345    // put out scalars
2346    numberWritten = fwrite(&scalars,sizeof(Clp_scalars),1,fp);
2347    if (numberWritten!=1)
2348      return 1;
2349    // strings
2350    int length;
2351    for (i=0;i<ClpLastStrParam;i++) {
2352      length = strParam_[i].size();
2353      numberWritten = fwrite(&length,sizeof(int),1,fp);
2354      if (numberWritten!=1)
2355        return 1;
2356      numberWritten = fwrite(strParam_[i].c_str(),length,1,fp);
2357      if (numberWritten!=1)
2358        return 1;
2359    }
2360    // arrays - in no particular order
2361    if (outDoubleArray(rowActivity_,numberRows_,fp))
2362        return 1;
2363    if (outDoubleArray(columnActivity_,numberColumns_,fp))
2364        return 1;
2365    if (outDoubleArray(dual_,numberRows_,fp))
2366        return 1;
2367    if (outDoubleArray(reducedCost_,numberColumns_,fp))
2368        return 1;
2369    if (outDoubleArray(rowLower_,numberRows_,fp))
2370        return 1;
2371    if (outDoubleArray(rowUpper_,numberRows_,fp))
2372        return 1;
2373    if (outDoubleArray(objective_,numberColumns_,fp))
2374        return 1;
2375    if (outDoubleArray(rowObjective_,numberRows_,fp))
2376        return 1;
2377    if (outDoubleArray(columnLower_,numberColumns_,fp))
2378        return 1;
2379    if (outDoubleArray(columnUpper_,numberColumns_,fp))
2380        return 1;
2381    if (ray_) {
2382      if (problemStatus_==1)
2383        if (outDoubleArray(ray_,numberRows_,fp))
2384          return 1;
2385      else if (problemStatus_==2)
2386        if (outDoubleArray(ray_,numberColumns_,fp))
2387          return 1;
2388      else
2389        if (outDoubleArray(NULL,0,fp))
2390          return 1;
2391    } else {
2392      if (outDoubleArray(NULL,0,fp))
2393        return 1;
2394    }
2395    if (status_&&(numberRows_+numberColumns_)>0) {
2396      length = numberRows_+numberColumns_;
2397      numberWritten = fwrite(&length,sizeof(int),1,fp);
2398      if (numberWritten!=1)
2399        return 1;
2400      numberWritten = fwrite(status_,sizeof(char),length, fp);
2401      if (numberWritten!=length)
2402        return 1;
2403    } else {
2404      length = 0;
2405      numberWritten = fwrite(&length,sizeof(int),1,fp);
2406      if (numberWritten!=1)
2407        return 1;
2408    }
2409    if (lengthNames_) {
2410      assert (numberRows_ == (int) rowNames_.size());
2411      for (i=0;i<numberRows_;i++) {
2412        length = rowNames_[i].size();
2413        numberWritten = fwrite(&length,sizeof(int),1,fp);
2414        if (numberWritten!=1)
2415          return 1;
2416        numberWritten = fwrite(rowNames_[i].c_str(),length,1,fp);
2417        if (numberWritten!=1)
2418          return 1;
2419      }
2420      assert (numberColumns_ == (int) columnNames_.size());
2421      for (i=0;i<numberColumns_;i++) {
2422        length = columnNames_[i].size();
2423        numberWritten = fwrite(&length,sizeof(int),1,fp);
2424        if (numberWritten!=1)
2425          return 1;
2426        numberWritten = fwrite(columnNames_[i].c_str(),length,1,fp);
2427        if (numberWritten!=1)
2428          return 1;
2429      }
2430    }
2431    // just standard type at present
2432    assert (matrix_->type()==1);
2433    assert (matrix_->getNumCols() == numberColumns_);
2434    assert (matrix_->getNumRows() == numberRows_);
2435    // we are going to save with gaps
2436    length = matrix_->getVectorStarts()[numberColumns_-1]
2437      + matrix_->getVectorLengths()[numberColumns_-1];
2438    numberWritten = fwrite(&length,sizeof(int),1,fp);
2439    if (numberWritten!=1)
2440      return 1;
2441    numberWritten = fwrite(matrix_->getElements(),
2442                           sizeof(double),length,fp);
2443    if (numberWritten!=length)
2444      return 1;
2445    numberWritten = fwrite(matrix_->getIndices(),
2446                           sizeof(int),length,fp);
2447    if (numberWritten!=length)
2448      return 1;
2449    numberWritten = fwrite(matrix_->getVectorStarts(),
2450                           sizeof(int),numberColumns_,fp);
2451    if (numberWritten!=numberColumns_)
2452      return 1;
2453    numberWritten = fwrite(matrix_->getVectorLengths(),
2454                           sizeof(int),numberColumns_,fp);
2455    if (numberWritten!=numberColumns_)
2456      return 1;
2457    // finished
2458    fclose(fp);
2459    return 0;
2460  } else {
2461    return -1;
2462  }
2463}
2464
2465int inDoubleArray(double * &array, int length, FILE * fp)
2466{
2467  int numberRead;
2468  int length2;
2469  numberRead = fread(&length2,sizeof(int),1,fp);
2470  if (numberRead!=1)
2471    return 1;
2472  if (length2) {
2473    // lengths must match
2474    if (length!=length2)
2475      return 2;
2476    array = new double[length];
2477    numberRead = fread(array,sizeof(double),length,fp);
2478    if (numberRead!=length)
2479      return 1;
2480  } 
2481  return 0;
2482}
2483/* Restore model from file, returns 0 if success,
2484   deletes current model */
2485int 
2486ClpSimplex::restoreModel(const char * fileName)
2487{
2488  FILE * fp = fopen(fileName,"rb");
2489  if (fp) {
2490    // Get rid of current model
2491    ClpModel::gutsOfDelete();
2492    gutsOfDelete(0);
2493    int i;
2494    for (i=0;i<6;i++) {
2495      rowArray_[i]=NULL;
2496      columnArray_[i]=NULL;
2497    }
2498    // get an empty factorization so we can set tolerances etc
2499    factorization_ = new ClpFactorization();
2500    // Say sparse
2501    factorization_->sparseThreshold(1);
2502    Clp_scalars scalars;
2503    int numberRead;
2504
2505    // get scalars
2506    numberRead = fread(&scalars,sizeof(Clp_scalars),1,fp);
2507    if (numberRead!=1)
2508      return 1;
2509    // Fill in scalars
2510    optimizationDirection_ = scalars.optimizationDirection;
2511    memcpy(dblParam_, scalars.dblParam, ClpLastDblParam * sizeof(double));
2512    objectiveValue_ = scalars.objectiveValue;
2513    dualBound_ = scalars.dualBound;
2514    dualTolerance_ = scalars.dualTolerance;
2515    primalTolerance_ = scalars.primalTolerance;
2516    sumDualInfeasibilities_ = scalars.sumDualInfeasibilities;
2517    sumPrimalInfeasibilities_ = scalars.sumPrimalInfeasibilities;
2518    infeasibilityCost_ = scalars.infeasibilityCost;
2519    numberRows_ = scalars.numberRows;
2520    numberColumns_ = scalars.numberColumns;
2521    memcpy(intParam_, scalars.intParam,ClpLastIntParam * sizeof(double));
2522    numberIterations_ = scalars.numberIterations;
2523    problemStatus_ = scalars.problemStatus;
2524    maximumIterations_ = scalars.maximumIterations;
2525    lengthNames_ = scalars.lengthNames;
2526    numberDualInfeasibilities_ = scalars.numberDualInfeasibilities;
2527    numberDualInfeasibilitiesWithoutFree_
2528      = scalars.numberDualInfeasibilitiesWithoutFree;
2529    numberPrimalInfeasibilities_ = scalars.numberPrimalInfeasibilities;
2530    numberRefinements_ = scalars.numberRefinements;
2531    scalingFlag_ = scalars.scalingFlag;
2532    algorithm_ = scalars.algorithm;
2533    specialOptions_ = scalars.specialOptions;
2534    // strings
2535    int length;
2536    for (i=0;i<ClpLastStrParam;i++) {
2537      numberRead = fread(&length,sizeof(int),1,fp);
2538      if (numberRead!=1)
2539        return 1;
2540      char * array = new char[length+1];
2541      numberRead = fread(array,length,1,fp);
2542      if (numberRead!=1)
2543        return 1;
2544      array[length]='\0';
2545      strParam_[i]=array;
2546      delete [] array;
2547    }
2548    // arrays - in no particular order
2549    if (inDoubleArray(rowActivity_,numberRows_,fp))
2550        return 1;
2551    if (inDoubleArray(columnActivity_,numberColumns_,fp))
2552        return 1;
2553    if (inDoubleArray(dual_,numberRows_,fp))
2554        return 1;
2555    if (inDoubleArray(reducedCost_,numberColumns_,fp))
2556        return 1;
2557    if (inDoubleArray(rowLower_,numberRows_,fp))
2558        return 1;
2559    if (inDoubleArray(rowUpper_,numberRows_,fp))
2560        return 1;
2561    if (inDoubleArray(objective_,numberColumns_,fp))
2562        return 1;
2563    if (inDoubleArray(rowObjective_,numberRows_,fp))
2564        return 1;
2565    if (inDoubleArray(columnLower_,numberColumns_,fp))
2566        return 1;
2567    if (inDoubleArray(columnUpper_,numberColumns_,fp))
2568        return 1;
2569    if (problemStatus_==1) {
2570      if (inDoubleArray(ray_,numberRows_,fp))
2571        return 1;
2572    } else if (problemStatus_==2) {
2573      if (inDoubleArray(ray_,numberColumns_,fp))
2574        return 1;
2575    } else {
2576      // ray should be null
2577      numberRead = fread(&length,sizeof(int),1,fp);
2578      if (numberRead!=1)
2579        return 1;
2580      if (length)
2581        return 2;
2582    }
2583    delete [] status_;
2584    status_=NULL;
2585    // status region
2586    numberRead = fread(&length,sizeof(int),1,fp);
2587    if (numberRead!=1)
2588        return 1;
2589    if (length) {
2590      if (length!=numberRows_+numberColumns_)
2591        return 1;
2592      status_ = new char unsigned[length];
2593      numberRead = fread(status_,sizeof(char),length, fp);
2594      if (numberRead!=length)
2595        return 1;
2596    }
2597    if (lengthNames_) {
2598      char * array = new char[lengthNames_+1];
2599      rowNames_.resize(0);
2600      for (i=0;i<numberRows_;i++) {
2601        numberRead = fread(&length,sizeof(int),1,fp);
2602        if (numberRead!=1)
2603          return 1;
2604        numberRead = fread(array,length,1,fp);
2605        if (numberRead!=1)
2606          return 1;
2607        rowNames_[i]=array;
2608      }
2609      columnNames_.resize(0);
2610      for (i=0;i<numberColumns_;i++) {
2611        numberRead = fread(&length,sizeof(int),1,fp);
2612        if (numberRead!=1)
2613          return 1;
2614        numberRead = fread(array,length,1,fp);
2615        if (numberRead!=1)
2616          return 1;
2617        columnNames_[i]=array;
2618      }
2619      delete [] array;
2620    }
2621    // Pivot choices
2622    assert(scalars.dualPivotChoice>0&&(scalars.dualPivotChoice&63)<3);
2623    delete dualRowPivot_;
2624    switch ((scalars.dualPivotChoice&63)) {
2625    case 1:
2626      // Dantzig
2627      dualRowPivot_ = new ClpDualRowDantzig();
2628      break;
2629    case 2:
2630      // Steepest - use mode
2631      dualRowPivot_ = new ClpDualRowSteepest(scalars.dualPivotChoice>>6);
2632      break;
2633    default:
2634      abort();
2635    }
2636    assert(scalars.primalPivotChoice>0&&(scalars.primalPivotChoice&63)<3);
2637    delete primalColumnPivot_;
2638    switch ((scalars.primalPivotChoice&63)) {
2639    case 1:
2640      // Dantzig
2641      primalColumnPivot_ = new ClpPrimalColumnDantzig();
2642      break;
2643    case 2:
2644      // Steepest - use mode
2645      primalColumnPivot_
2646        = new ClpPrimalColumnSteepest(scalars.primalPivotChoice>>6);
2647      break;
2648    default:
2649      abort();
2650    }
2651    assert(scalars.matrixStorageChoice==1);
2652    delete matrix_;
2653    // get arrays
2654    numberRead = fread(&length,sizeof(int),1,fp);
2655    if (numberRead!=1)
2656      return 1;
2657    double * elements = new double[length];
2658    int * indices = new int[length];
2659    int * starts = new int[numberColumns_];
2660    int * lengths = new int[numberColumns_];
2661    numberRead = fread(elements, sizeof(double),length,fp);
2662    if (numberRead!=length)
2663      return 1;
2664    numberRead = fread(indices, sizeof(int),length,fp);
2665    if (numberRead!=length)
2666      return 1;
2667    numberRead = fread(starts, sizeof(int),numberColumns_,fp);
2668    if (numberRead!=numberColumns_)
2669      return 1;
2670    numberRead = fread(lengths, sizeof(int),numberColumns_,fp);
2671    if (numberRead!=numberColumns_)
2672      return 1;
2673    // assign matrix
2674    CoinPackedMatrix * matrix = new CoinPackedMatrix();
2675    matrix->assignMatrix(true, numberRows_, numberColumns_,
2676                         length, elements, indices, starts, lengths);
2677    // and transfer to Clp
2678    matrix_ = new ClpPackedMatrix(matrix);
2679    // finished
2680    fclose(fp);
2681    return 0;
2682  } else {
2683    return -1;
2684  }
2685  return 0;
2686}
2687// value of incoming variable (in Dual)
2688double 
2689ClpSimplex::valueIncomingDual() const
2690{
2691  // Need value of incoming for list of infeasibilities as may be infeasible
2692  double valueIncoming = (dualOut_/alpha_)*directionOut_;
2693  if (directionIn_==-1)
2694    valueIncoming = upperIn_-valueIncoming;
2695  else
2696    valueIncoming = lowerIn_-valueIncoming;
2697  return valueIncoming;
2698}
2699//#############################################################################
2700
2701ClpSimplexProgress::ClpSimplexProgress () 
2702{
2703  int i;
2704  for (i=0;i<CLP_PROGRESS;i++) {
2705    objective_[i] = 0.0;
2706    infeasibility_[i] = -1.0; // set to an impossible value
2707    numberInfeasibilities_[i]=-1; 
2708  }
2709  numberTimes_ = 0;
2710  numberBadTimes_ = 0;
2711  model_ = NULL;
2712}
2713
2714
2715//-----------------------------------------------------------------------------
2716
2717ClpSimplexProgress::~ClpSimplexProgress ()
2718{
2719}
2720// Copy constructor.
2721ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs) 
2722{
2723  int i;
2724  for (i=0;i<CLP_PROGRESS;i++) {
2725    objective_[i] = rhs.objective_[i];
2726    infeasibility_[i] = rhs.infeasibility_[i];
2727    numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2728  }
2729  numberTimes_ = rhs.numberTimes_;
2730  numberBadTimes_ = rhs.numberBadTimes_;
2731  model_ = rhs.model_;
2732}
2733// Copy constructor.from model
2734ClpSimplexProgress::ClpSimplexProgress(ClpSimplex * model) 
2735{
2736  int i;
2737  for (i=0;i<CLP_PROGRESS;i++) {
2738    objective_[i] = 0.0;
2739    infeasibility_[i] = -1.0; // set to an impossible value
2740    numberInfeasibilities_[i]=-1; 
2741  }
2742  numberTimes_ = 0;
2743  numberBadTimes_ = 0;
2744  model_ = model;
2745}
2746// Assignment operator. This copies the data
2747ClpSimplexProgress & 
2748ClpSimplexProgress::operator=(const ClpSimplexProgress & rhs)
2749{
2750  if (this != &rhs) {
2751    int i;
2752    for (i=0;i<CLP_PROGRESS;i++) {
2753      objective_[i] = rhs.objective_[i];
2754      infeasibility_[i] = rhs.infeasibility_[i];
2755      numberInfeasibilities_[i]=rhs.numberInfeasibilities_[i]; 
2756    }
2757    numberTimes_ = rhs.numberTimes_;
2758    numberBadTimes_ = rhs.numberBadTimes_;
2759    model_ = rhs.model_;
2760  }
2761  return *this;
2762}
2763int
2764ClpSimplexProgress::looping()
2765{
2766  assert(model_);
2767  double objective = model_->objectiveValue();
2768  double infeasibility;
2769  int numberInfeasibilities;
2770  if (model_->algorithm()<0) {
2771    // dual
2772    infeasibility = model_->sumPrimalInfeasibilities();
2773    numberInfeasibilities = model_->numberPrimalInfeasibilities();
2774  } else {
2775    //primal
2776    infeasibility = model_->sumDualInfeasibilities();
2777    numberInfeasibilities = model_->numberDualInfeasibilities();
2778  }
2779  int i;
2780  int numberMatched=0;
2781  int matched=0;
2782
2783  for (i=0;i<CLP_PROGRESS;i++) {
2784    if (objective==objective_[i]&&
2785        infeasibility==infeasibility_[i]&&
2786        numberInfeasibilities==numberInfeasibilities_[i]) {
2787      matched |= (1<<i);
2788      numberMatched++;
2789    }
2790    if (i) {
2791      objective_[i-1] = objective_[i];
2792      infeasibility_[i-1] = infeasibility_[i];
2793      numberInfeasibilities_[i-1]=numberInfeasibilities_[i]; 
2794    }
2795  }
2796  objective_[CLP_PROGRESS-1] = objective;
2797  infeasibility_[CLP_PROGRESS-1] = infeasibility;
2798  numberInfeasibilities_[CLP_PROGRESS-1]=numberInfeasibilities;
2799  if (model_->progressFlag())
2800    numberMatched=0;
2801  numberTimes_++;
2802  if (numberTimes_<10)
2803    numberMatched=0;
2804  // skip if just last time as may be checking something
2805  if (matched==(1<<(CLP_PROGRESS-1)))
2806    numberMatched=0;
2807  if (numberMatched) {
2808    model_->messageHandler()->message(CLP_POSSIBLELOOP,model_->messages())
2809      <<numberMatched
2810      <<matched
2811      <<numberTimes_
2812      <<CoinMessageEol;
2813    numberBadTimes_++;
2814    if (numberBadTimes_<10) {
2815      if (model_->algorithm()<0) {
2816        // dual - change tolerance
2817        model_->setCurrentDualTolerance(model_->currentDualTolerance()*1.05);
2818      } else {
2819        // primal - change tolerance
2820        model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance()*1.05);
2821      }
2822    } else {
2823      model_->messageHandler()->message(CLP_LOOP,model_->messages())
2824        <<CoinMessageEol;
2825      // look at solution and maybe declare victory
2826      if (infeasibility<1.0e-4)
2827        return 0;
2828      else
2829        return 3;
2830    }
2831  }
2832  return -1;
2833}
2834// Sanity check on input data - returns true if okay
2835bool 
2836ClpSimplex::sanityCheck()
2837{
2838  // bad if empty
2839  if (!numberRows_||!numberColumns_||!matrix_->getNumElements()) {
2840    handler_->message(CLP_EMPTY_PROBLEM,messages_)
2841      <<numberRows_
2842      <<numberColumns_
2843      <<matrix_->getNumElements()
2844      <<CoinMessageEol;
2845    problemStatus_=4;
2846    return false;
2847  }
2848  int numberBad , numberBadBounds;
2849  double largestBound, smallestBound, minimumGap;
2850  double smallestObj, largestObj;
2851  int firstBad;
2852  int modifiedBounds=0;
2853  int i;
2854  numberBad=0;
2855  numberBadBounds=0;
2856  firstBad=-1;
2857  minimumGap=1.0e100;
2858  smallestBound=1.0e100;
2859  largestBound=0.0;
2860  smallestObj=1.0e100;
2861  largestObj=0.0;
2862  for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
2863    double value;
2864    value = fabs(cost_[i]);
2865    if (value>1.0e50) {
2866      numberBad++;
2867      if (firstBad<0)
2868        firstBad=i;
2869    } else if (value) {
2870      if (value>largestObj)
2871        largestObj=value;
2872      if (value<smallestObj)
2873        smallestObj=value;
2874    }
2875    value=upper_[i]-lower_[i];
2876    if (value<-primalTolerance_) {
2877      numberBadBounds++;
2878      if (firstBad<0)
2879        firstBad=i;
2880    } else if (value<=primalTolerance_) {
2881      if (value) {
2882        // modify
2883        upper_[i] = lower_[i];
2884        modifiedBounds++;
2885      }
2886    } else {
2887      if (value<minimumGap)
2888        minimumGap=value;
2889    }
2890    if (lower_[i]>-1.0e100&&lower_[i]) {
2891      value = fabs(lower_[i]);
2892      if (value>largestBound)
2893        largestBound=value;
2894      if (value<smallestBound)
2895        smallestBound=value;
2896    }
2897    if (upper_[i]<1.0e100&&upper_[i]) {
2898      value = fabs(upper_[i]);
2899      if (value>largestBound)
2900        largestBound=value;
2901      if (value<smallestBound)
2902        smallestBound=value;
2903    }
2904  }
2905  if (largestBound)
2906    handler_->message(CLP_RIMSTATISTICS3,messages_)
2907      <<smallestBound
2908      <<largestBound
2909      <<minimumGap
2910      <<CoinMessageEol;
2911  minimumGap=1.0e100;
2912  smallestBound=1.0e100;
2913  largestBound=0.0;
2914  for (i=0;i<numberColumns_;i++) {
2915    double value;
2916    value = fabs(cost_[i]);
2917    if (value>1.0e50) {
2918      numberBad++;
2919      if (firstBad<0)
2920        firstBad=i;
2921    } else if (value) {
2922      if (value>largestObj)
2923        largestObj=value;
2924      if (value<smallestObj)
2925        smallestObj=value;
2926    }
2927    value=upper_[i]-lower_[i];
2928    if (value<-primalTolerance_) {
2929      numberBadBounds++;
2930      if (firstBad<0)
2931        firstBad=i;
2932    } else if (value<=primalTolerance_) {
2933      if (value) {
2934        // modify
2935        upper_[i] = lower_[i];
2936        modifiedBounds++;
2937      }
2938    } else {
2939      if (value<minimumGap)
2940        minimumGap=value;
2941    }
2942    if (lower_[i]>-1.0e100&&lower_[i]) {
2943      value = fabs(lower_[i]);
2944      if (value>largestBound)
2945        largestBound=value;
2946      if (value<smallestBound)
2947        smallestBound=value;
2948    }
2949    if (upper_[i]<1.0e100&&upper_[i]) {
2950      value = fabs(upper_[i]);
2951      if (value>largestBound)
2952        largestBound=value;
2953      if (value<smallestBound)
2954        smallestBound=value;
2955    }
2956  }
2957  char rowcol[]={'R','C'};
2958  if (numberBad) {
2959    handler_->message(CLP_BAD_BOUNDS,messages_)
2960      <<numberBad
2961      <<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
2962      <<CoinMessageEol;
2963    problemStatus_=4;
2964    return false;
2965  }
2966  if (modifiedBounds)
2967    handler_->message(CLP_MODIFIEDBOUNDS,messages_)
2968      <<modifiedBounds
2969      <<CoinMessageEol;
2970  handler_->message(CLP_RIMSTATISTICS1,messages_)
2971    <<smallestObj
2972    <<largestObj
2973    <<CoinMessageEol;
2974  if (largestBound)
2975    handler_->message(CLP_RIMSTATISTICS2,messages_)
2976      <<smallestBound
2977      <<largestBound
2978      <<minimumGap
2979      <<CoinMessageEol;
2980  return true;
2981}
2982// Set up status array (for OsiClp)
2983void 
2984ClpSimplex::createStatus() 
2985{
2986  delete [] status_;
2987  status_ = new unsigned char [numberColumns_+numberRows_];
2988  memset(status_,0,(numberColumns_+numberRows_)*sizeof(char));
2989  int i;
2990  // set column status to one nearest zero
2991  for (i=0;i<numberColumns_;i++) {
2992#if 0
2993    if (columnLower_[i]>=0.0) {
2994      setColumnStatus(i,atLowerBound);
2995    } else if (columnUpper_[i]<=0.0) {
2996      setColumnStatus(i,atUpperBound);
2997    } else if (columnLower_[i]<-1.0e20&&columnUpper_[i]>1.0e20) {
2998      // free
2999      setColumnStatus(i,isFree);
3000    } else if (fabs(columnLower_[i])<fabs(columnUpper_[i])) {
3001      setColumnStatus(i,atLowerBound);
3002    } else {
3003      setColumnStatus(i,atUpperBound);
3004    }
3005#else
3006    setColumnStatus(i,atLowerBound);
3007#endif
3008  }
3009  for (i=0;i<numberRows_;i++) {
3010    setRowStatus(i,basic);
3011  }
3012}
3013/* Loads a problem (the constraints on the
3014   rows are given by lower and upper bounds). If a pointer is 0 then the
3015   following values are the default:
3016   <ul>
3017   <li> <code>colub</code>: all columns have upper bound infinity
3018   <li> <code>collb</code>: all columns have lower bound 0
3019   <li> <code>rowub</code>: all rows have upper bound infinity
3020   <li> <code>rowlb</code>: all rows have lower bound -infinity
3021   <li> <code>obj</code>: all variables have 0 objective coefficient
3022   </ul>
3023*/
3024void 
3025ClpSimplex::loadProblem (  const ClpMatrixBase& matrix,
3026                    const double* collb, const double* colub,   
3027                    const double* obj,
3028                    const double* rowlb, const double* rowub,
3029                    const double * rowObjective)
3030{
3031  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3032                        rowObjective);
3033  createStatus();
3034}
3035void 
3036ClpSimplex::loadProblem (  const CoinPackedMatrix& matrix,
3037                    const double* collb, const double* colub,   
3038                    const double* obj,
3039                    const double* rowlb, const double* rowub,
3040                    const double * rowObjective)
3041{
3042  ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
3043                        rowObjective);
3044  createStatus();
3045}
3046
3047/* Just like the other loadProblem() method except that the matrix is
3048   given in a standard column major ordered format (without gaps). */
3049void 
3050ClpSimplex::loadProblem (  const int numcols, const int numrows,
3051                    const int* start, const int* index,
3052                    const double* value,
3053                    const double* collb, const double* colub,   
3054                    const double* obj,
3055                    const double* rowlb, const double* rowub,
3056                    const double * rowObjective)
3057{
3058  ClpModel::loadProblem(numcols, numrows, start, index, value,
3059                          collb, colub, obj, rowlb, rowub,
3060                          rowObjective);
3061  createStatus();
3062}
3063void 
3064ClpSimplex::loadProblem (  const int numcols, const int numrows,
3065                           const int* start, const int* index,
3066                           const double* value,const int * length,
3067                           const double* collb, const double* colub,   
3068                           const double* obj,
3069                           const double* rowlb, const double* rowub,
3070                           const double * rowObjective)
3071{
3072  ClpModel::loadProblem(numcols, numrows, start, index, value, length,
3073                          collb, colub, obj, rowlb, rowub,
3074                          rowObjective);
3075  createStatus();
3076}
3077// Read an mps file from the given filename
3078int 
3079ClpSimplex::readMps(const char *filename,
3080            bool keepNames,
3081            bool ignoreErrors)
3082{
3083  int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
3084  createStatus();
3085  return status;
3086}
3087// Just check solution (for external use)
3088void 
3089ClpSimplex::checkSolution()
3090{
3091  // put in standard form
3092  createRim(7+8+16);
3093  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
3094  checkDualSolution();
3095  // release extra memory
3096  deleteRim(false);
3097}
Note: See TracBrowser for help on using the repository browser.