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

Last change on this file since 25 was 23, checked in by forrest, 18 years ago

Allowing for bad matrices

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