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

Last change on this file since 19 was 19, checked in by ladanyi, 17 years ago

reordering Clp

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