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

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

Presolve in as option

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