source: trunk/include/ClpSimplex.hpp @ 427

Last change on this file since 427 was 427, checked in by forrest, 15 years ago

make some things protected

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.3 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6 
7   John Forrest
8
9 */
10#ifndef ClpSimplex_H
11#define ClpSimplex_H
12
13#include <iostream>
14#include <cfloat>
15#include "ClpModel.hpp"
16#include "ClpMatrixBase.hpp"
17#include "ClpSolve.hpp"
18class ClpDualRowPivot;
19class ClpPrimalColumnPivot;
20class ClpFactorization;
21class CoinIndexedVector;
22class ClpNonLinearCost;
23class ClpSimplexProgress;
24
25/** This solves LPs using the simplex method
26
27    It inherits from ClpModel and all its arrays are created at
28    algorithm time. Originally I tried to work with model arrays
29    but for simplicity of coding I changed to single arrays with
30    structural variables then row variables.  Some coding is still
31    based on old style and needs cleaning up.
32
33    For a description of algorithms:
34
35    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
36    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
37
38    There is an algorithm data member.  + for primal variations
39    and - for dual variations
40
41    This file also includes (at end) a very simple class ClpSimplexProgress
42    which is where anti-looping stuff should migrate to
43
44*/
45
46class ClpSimplex : public ClpModel {
47   friend void ClpSimplexUnitTest(const std::string & mpsDir,
48                                  const std::string & netlibDir);
49
50public:
51  /** enums for status of various sorts.
52      First 4 match CoinWarmStartBasis,
53      isFixed means fixed at lower bound and out of basis
54  */
55  enum Status {
56    isFree = 0x00,
57    basic = 0x01,
58    atUpperBound = 0x02,
59    atLowerBound = 0x03,
60    superBasic = 0x04,
61    isFixed = 0x05
62  };
63  // For Dual
64  enum FakeBound {
65    noFake = 0x00,
66    bothFake = 0x01,
67    upperFake = 0x02,
68    lowerFake = 0x03
69  };
70
71  /**@name Constructors and destructor and copy */
72  //@{
73  /// Default constructor
74    ClpSimplex (  );
75
76  /** Copy constructor. May scale depending on mode
77      -1 leave mode as is
78      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
79  */
80  ClpSimplex(const ClpSimplex & rhs, int scalingMode =-1);
81  /** Copy constructor from model. May scale depending on mode
82      -1 leave mode as is
83      0 -off, 1 equilibrium, 2 geometric, 3, auto, 4 dynamic(later)
84  */
85  ClpSimplex(const ClpModel & rhs, int scalingMode=-1);
86  /** Subproblem constructor.  A subset of whole model is created from the
87      row and column lists given.  The new order is given by list order and
88      duplicates are allowed.  Name and integer information can be dropped
89  */
90  ClpSimplex (const ClpModel * wholeModel,
91              int numberRows, const int * whichRows,
92              int numberColumns, const int * whichColumns,
93              bool dropNames=true, bool dropIntegers=true);
94  /** This constructor modifies original ClpSimplex and stores
95      original stuff in created ClpSimplex.  It is only to be used in
96      conjunction with originalModel */
97  ClpSimplex (ClpSimplex * wholeModel,
98              int numberColumns, const int * whichColumns);
99  /** This copies back stuff from miniModel and then deletes miniModel.
100      Only to be used with mini constructor */
101  void originalModel(ClpSimplex * miniModel);
102  /// Assignment operator. This copies the data
103    ClpSimplex & operator=(const ClpSimplex & rhs);
104  /// Destructor
105   ~ClpSimplex (  );
106  // Ones below are just ClpModel with some changes
107  /** Loads a problem (the constraints on the
108        rows are given by lower and upper bounds). If a pointer is 0 then the
109        following values are the default:
110        <ul>
111          <li> <code>colub</code>: all columns have upper bound infinity
112          <li> <code>collb</code>: all columns have lower bound 0
113          <li> <code>rowub</code>: all rows have upper bound infinity
114          <li> <code>rowlb</code>: all rows have lower bound -infinity
115          <li> <code>obj</code>: all variables have 0 objective coefficient
116        </ul>
117    */
118  void loadProblem (  const ClpMatrixBase& matrix,
119                     const double* collb, const double* colub,   
120                     const double* obj,
121                     const double* rowlb, const double* rowub,
122                      const double * rowObjective=NULL);
123  void loadProblem (  const CoinPackedMatrix& matrix,
124                     const double* collb, const double* colub,   
125                     const double* obj,
126                     const double* rowlb, const double* rowub,
127                      const double * rowObjective=NULL);
128
129  /** Just like the other loadProblem() method except that the matrix is
130        given in a standard column major ordered format (without gaps). */
131  void loadProblem (  const int numcols, const int numrows,
132                     const CoinBigIndex* start, const int* index,
133                     const double* value,
134                     const double* collb, const double* colub,   
135                     const double* obj,
136                      const double* rowlb, const double* rowub,
137                      const double * rowObjective=NULL);
138  /// This one is for after presolve to save memory
139  void loadProblem (  const int numcols, const int numrows,
140                     const CoinBigIndex* start, const int* index,
141                      const double* value,const int * length,
142                     const double* collb, const double* colub,   
143                     const double* obj,
144                      const double* rowlb, const double* rowub,
145                      const double * rowObjective=NULL);
146  /// Read an mps file from the given filename
147  int readMps(const char *filename,
148              bool keepNames=false,
149              bool ignoreErrors = false);
150  /** Borrow model.  This is so we dont have to copy large amounts
151      of data around.  It assumes a derived class wants to overwrite
152      an empty model with a real one - while it does an algorithm.
153      This is same as ClpModel one, but sets scaling on etc. */
154  void borrowModel(ClpModel & otherModel);
155  void borrowModel(ClpSimplex & otherModel);
156   /// Pass in Event handler (cloned and deleted at end)
157   void passInEventHandler(const ClpEventHandler * eventHandler);
158  //@}
159
160  /**@name Functions most useful to user */
161  //@{
162  /** General solve algorithm which can do presolve.
163      See  ClpSolve.hpp for options
164   */
165  int initialSolve(ClpSolve & options);
166  /// Default initial solve
167  int initialSolve();
168  /// Dual initial solve
169  int initialDualSolve();
170  /// Primal initial solve
171  int initialPrimalSolve();
172  /** Dual algorithm - see ClpSimplexDual.hpp for method.
173      ifValuesPass==2 just does values pass and then stops */
174  int dual(int ifValuesPass=0);
175  /** Primal algorithm - see ClpSimplexPrimal.hpp for method.
176      ifValuesPass==2 just does values pass and then stops */
177  int primal(int ifValuesPass=0);
178  /** Solves nonlinear problem using SLP - may be used as crash
179      for other algorithms when number of iterations small.
180      Also exits if all problematical variables are changing
181      less than deltaTolerance
182  */
183  int nonlinearSLP(int numberPasses,double deltaTolerance);
184  /** Solves using barrier (assumes you have good cholesky factor code).
185      Does crossover to simplex if asked*/
186  int barrier(bool crossover=true);
187  /** Solves non-linear using reduced gradient.  Phase = 0 get feasible,
188      =1 use solution */
189  int reducedGradient(int phase=0);
190  /** Dual ranging.
191      This computes increase/decrease in cost for each given variable and corresponding
192      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
193      and numberColumns.. for artificials/slacks.
194      For non-basic variables the sequence number will be that of the non-basic variables.
195      The increase/decrease value is always >= 0.0
196
197      Up to user to provide correct length arrays.
198
199      Returns non-zero if infeasible unbounded etc
200  */
201  int dualRanging(int numberCheck,const int * which,
202                  double * costIncrease, int * sequenceIncrease,
203                  double * costDecrease, int * sequenceDecrease);
204  /** Primal ranging.
205      This computes increase/decrease in value for each given variable and corresponding
206      sequence numbers which would change basis.  Sequence numbers are 0..numberColumns
207      and numberColumns.. for artificials/slacks.
208      For basic variables the sequence number will be that of the basic variables.
209
210      Up to user to providen correct length arrays.
211
212      Returns non-zero if infeasible unbounded etc
213  */
214  int primalRanging(int numberCheck,const int * which,
215                  double * valueIncrease, int * sequenceIncrease,
216                  double * valueDecrease, int * sequenceDecrease);
217  /// Passes in factorization
218  void setFactorization( ClpFactorization & factorization);
219  /** Tightens primal bounds to make dual faster.  Unless
220      fixed, bounds are slightly looser than they could be.
221      This is to make dual go faster and is probably not needed
222      with a presolve.  Returns non-zero if problem infeasible.
223
224      Fudge for branch and bound - put bounds on columns of factor *
225      largest value (at continuous) - should improve stability
226      in branch and bound on infeasible branches (0.0 is off)
227  */
228  int tightenPrimalBounds(double factor=0.0);
229  /** Crash - at present just aimed at dual, returns
230      -2 if dual preferred and crash basis created
231      -1 if dual preferred and all slack basis preferred
232       0 if basis going in was not all slack
233       1 if primal preferred and all slack basis preferred
234       2 if primal preferred and crash basis created.
235       
236       if gap between bounds <="gap" variables can be flipped
237
238       If "pivot" is
239       0 No pivoting (so will just be choice of algorithm)
240       1 Simple pivoting e.g. gub
241       2 Mini iterations
242  */
243  int crash(double gap,int pivot);
244  /// Sets row pivot choice algorithm in dual
245  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
246  /// Sets column pivot choice algorithm in primal
247  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
248  /** For strong branching.  On input lower and upper are new bounds
249      while on output they are change in objective function values
250      (>1.0e50 infeasible).
251      Return code is 0 if nothing interesting, -1 if infeasible both
252      ways and +1 if infeasible one way (check values to see which one(s))
253      Solutions are filled in as well - even down, odd up - also
254      status and number of iterations
255  */
256  int strongBranching(int numberVariables,const int * variables,
257                      double * newLower, double * newUpper,
258                      double ** outputSolution,
259                      int * outputStatus, int * outputIterations,
260                      bool stopOnFirstInfeasible=true,
261                      bool alwaysFinish=false);
262  //@}
263
264  /**@name Needed for functionality of OsiSimplexInterface */
265  //@{
266  /** Pivot in a variable and out a variable.  Returns 0 if okay,
267      1 if inaccuracy forced re-factorization, -1 if would be singular.
268      Also updates primal/dual infeasibilities.
269      Assumes sequenceIn_ and pivotRow_ set and also directionIn and Out.
270  */
271  int pivot();
272
273  /** Pivot in a variable and choose an outgoing one.  Assumes primal
274      feasible - will not go through a bound.  Returns step length in theta
275      Returns ray in ray_ (or NULL if no pivot)
276      Return codes as before but -1 means no acceptable pivot
277  */
278  int primalPivotResult();
279 
280  /** Pivot out a variable and choose an incoing one.  Assumes dual
281      feasible - will not go through a reduced cost. 
282      Returns step length in theta
283      Returns ray in ray_ (or NULL if no pivot)
284      Return codes as before but -1 means no acceptable pivot
285  */
286  int dualPivotResult();
287
288  /** Common bits of coding for dual and primal.  Return s0 if okay,
289      1 if bad matrix, 2 if very bad factorization
290  */
291  int startup(int ifValuesPass);
292  void finish();
293 
294  /** Factorizes and returns true if optimal.  Used by user */
295  bool statusOfProblem(bool initial=false);
296  //@}
297
298  /**@name most useful gets and sets */
299  //@{
300  /// If problem is primal feasible
301  inline bool primalFeasible() const
302         { return (numberPrimalInfeasibilities_==0);};
303  /// If problem is dual feasible
304  inline bool dualFeasible() const
305         { return (numberDualInfeasibilities_==0);};
306  /// factorization
307  inline ClpFactorization * factorization() const 
308          { return factorization_;};
309  /// Sparsity on or off
310  bool sparseFactorization() const;
311  void setSparseFactorization(bool value);
312  /// Factorization frequency
313  int factorizationFrequency() const;
314  void setFactorizationFrequency(int value);
315  /// Dual bound
316  inline double dualBound() const
317          { return dualBound_;};
318  void setDualBound(double value);
319  /// Infeasibility cost
320  inline double infeasibilityCost() const
321          { return infeasibilityCost_;};
322  void setInfeasibilityCost(double value);
323  /** Amount of print out:
324      0 - none
325      1 - just final
326      2 - just factorizations
327      3 - as 2 plus a bit more
328      4 - verbose
329      above that 8,16,32 etc just for selective debug
330  */
331  /** Perturbation:
332      50  - switch on perturbation
333      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
334      101 - we are perturbed
335      102 - don't try perturbing again
336      default is 100
337      others are for playing
338  */
339  inline int perturbation() const
340    { return perturbation_;};
341  void setPerturbation(int value);
342  /// Current (or last) algorithm
343  inline int algorithm() const 
344  {return algorithm_; } ;
345  /// Set algorithm
346  inline void setAlgorithm(int value)
347  {algorithm_=value; } ;
348  /// Sum of dual infeasibilities
349  inline double sumDualInfeasibilities() const 
350          { return sumDualInfeasibilities_;} ;
351  inline void setSumDualInfeasibilities(double value)
352          { sumDualInfeasibilities_=value;} ;
353  /// Sum of relaxed dual infeasibilities
354  inline double sumOfRelaxedDualInfeasibilities() const 
355          { return sumOfRelaxedDualInfeasibilities_;} ;
356  inline void setSumOfRelaxedDualInfeasibilities(double value)
357          { sumOfRelaxedDualInfeasibilities_=value;} ;
358  /// Number of dual infeasibilities
359  inline int numberDualInfeasibilities() const 
360          { return numberDualInfeasibilities_;} ;
361  inline void setNumberDualInfeasibilities(int value)
362          { numberDualInfeasibilities_=value;} ;
363  /// Sum of primal infeasibilities
364  inline double sumPrimalInfeasibilities() const 
365          { return sumPrimalInfeasibilities_;} ;
366  inline void setSumPrimalInfeasibilities(double value)
367          { sumPrimalInfeasibilities_=value;} ;
368  /// Sum of relaxed primal infeasibilities
369  inline double sumOfRelaxedPrimalInfeasibilities() const 
370          { return sumOfRelaxedPrimalInfeasibilities_;} ;
371  inline void setSumOfRelaxedPrimalInfeasibilities(double value)
372          { sumOfRelaxedPrimalInfeasibilities_=value;} ;
373  /// Number of primal infeasibilities
374  inline int numberPrimalInfeasibilities() const 
375          { return numberPrimalInfeasibilities_;} ;
376  inline void setNumberPrimalInfeasibilities(int value)
377          { numberPrimalInfeasibilities_=value;} ;
378  /** Save model to file, returns 0 if success.  This is designed for
379      use outside algorithms so does not save iterating arrays etc.
380  It does not save any messaging information.
381  Does not save scaling values.
382  It does not know about all types of virtual functions.
383  */
384  int saveModel(const char * fileName);
385  /** Restore model from file, returns 0 if success,
386      deletes current model */
387  int restoreModel(const char * fileName);
388 
389  /** Just check solution (for external use) - sets sum of
390      infeasibilities etc */
391  void checkSolution();
392  /// Useful row length arrays (0,1,2,3,4,5)
393  inline CoinIndexedVector * rowArray(int index) const
394  { return rowArray_[index];};
395  /// Useful column length arrays (0,1,2,3,4,5)
396  inline CoinIndexedVector * columnArray(int index) const
397  { return columnArray_[index];};
398  //@}
399
400  /******************** End of most useful part **************/
401  /**@name Functions less likely to be useful to casual user */
402  //@{
403  /** Given an existing factorization computes and checks
404      primal and dual solutions.  Uses input arrays for variables at
405      bounds.  Returns feasibility states */
406  int getSolution (  const double * rowActivities,
407                     const double * columnActivities);
408  /** Given an existing factorization computes and checks
409      primal and dual solutions.  Uses current problem arrays for
410      bounds.  Returns feasibility states */
411  int getSolution ();
412  /** Constructs a non linear cost from list of non-linearities (columns only)
413      First lower of each column is taken as real lower
414      Last lower is taken as real upper and cost ignored
415
416      Returns nonzero if bad data e.g. lowers not monotonic
417  */
418  int createPiecewiseLinearCosts(const int * starts,
419                   const double * lower, const double * gradient);
420  /** Return model - updates any scalars */
421  void returnModel(ClpSimplex & otherModel);
422  /** Factorizes using current basis. 
423      solveType - 1 iterating, 0 initial, -1 external
424      If 10 added then in primal values pass
425      Return codes are as from ClpFactorization unless initial factorization
426      when total number of singularities is returned
427  */
428  int internalFactorize(int solveType);
429  /// Save data
430  ClpDataSave saveData() ;
431  /// Restore data
432    void restoreData(ClpDataSave saved);
433  /// Clean up status
434  void cleanStatus();
435  /// Factorizes using current basis. For external use
436  int factorize();
437  /** Computes duals from scratch. If givenDjs then
438      allows for nonzero basic djs */
439  void computeDuals(double * givenDjs);
440  /// Computes primals from scratch
441  void computePrimals (  const double * rowActivities,
442                     const double * columnActivities);
443  /** Adds multiple of a column into an array */
444  void add(double * array,
445                   int column, double multiplier) const;
446  /**
447     Unpacks one column of the matrix into indexed array
448     Uses sequenceIn_
449     Also applies scaling if needed
450  */
451  void unpack(CoinIndexedVector * rowArray) const ;
452  /**
453     Unpacks one column of the matrix into indexed array
454     Slack if sequence>= numberColumns
455     Also applies scaling if needed
456  */
457  void unpack(CoinIndexedVector * rowArray,int sequence) const;
458  /**
459     Unpacks one column of the matrix into indexed array
460     ** as packed vector
461     Uses sequenceIn_
462     Also applies scaling if needed
463  */
464  void unpackPacked(CoinIndexedVector * rowArray) ;
465  /**
466     Unpacks one column of the matrix into indexed array
467     ** as packed vector
468     Slack if sequence>= numberColumns
469     Also applies scaling if needed
470  */
471  void unpackPacked(CoinIndexedVector * rowArray,int sequence);
472protected: 
473  /**
474      This does basis housekeeping and does values for in/out variables.
475      Can also decide to re-factorize
476  */
477  int housekeeping(double objectiveChange);
478  /** This sets largest infeasibility and most infeasible and sum
479      and number of infeasibilities (Primal) */
480  void checkPrimalSolution(const double * rowActivities=NULL,
481                           const double * columnActivies=NULL);
482  /** This sets largest infeasibility and most infeasible and sum
483      and number of infeasibilities (Dual) */
484  void checkDualSolution();
485public:
486  /** For advanced use.  When doing iterative solves things can get
487      nasty so on values pass if incoming solution has largest
488      infeasibility < incomingInfeasibility throw out variables
489      from basis until largest infeasibility < allowedInfeasibility
490      or incoming largest infeasibility.
491      If allowedInfeasibility>= incomingInfeasibility this is
492      always possible altough you may end up with an all slack basis.
493
494      Defaults are 1.0,10.0
495  */
496  void setValuesPassAction(float incomingInfeasibility,
497                           float allowedInfeasibility);
498  //@}
499  /**@name most useful gets and sets */
500  //@{
501  /// Worst column primal infeasibility
502  inline double columnPrimalInfeasibility() const 
503          { return columnPrimalInfeasibility_;} ;
504  /// Sequence of worst (-1 if feasible)
505  inline int columnPrimalSequence() const 
506          { return columnPrimalSequence_;} ;
507  /// Worst row primal infeasibility
508  inline double rowPrimalInfeasibility() const 
509          { return rowPrimalInfeasibility_;} ;
510  /// Sequence of worst (-1 if feasible)
511  inline int rowPrimalSequence() const 
512          { return rowPrimalSequence_;} ;
513  /** Worst column dual infeasibility (note - these may not be as meaningful
514      if the problem is primal infeasible */
515  inline double columnDualInfeasibility() const 
516          { return columnDualInfeasibility_;} ;
517  /// Sequence of worst (-1 if feasible)
518  inline int columnDualSequence() const 
519          { return columnDualSequence_;} ;
520  /// Worst row dual infeasibility
521  inline double rowDualInfeasibility() const 
522          { return rowDualInfeasibility_;} ;
523  /// Sequence of worst (-1 if feasible)
524  inline int rowDualSequence() const 
525          { return rowDualSequence_;} ;
526  /// Primal tolerance needed to make dual feasible (<largeTolerance)
527  inline double primalToleranceToGetOptimal() const 
528          { return primalToleranceToGetOptimal_;} ;
529  /// Remaining largest dual infeasibility
530  inline double remainingDualInfeasibility() const 
531          { return remainingDualInfeasibility_;} ;
532  /// Large bound value (for complementarity etc)
533  inline double largeValue() const 
534          { return largeValue_;} ;
535  void setLargeValue( double value) ;
536  /// Largest error on Ax-b
537  inline double largestPrimalError() const
538          { return largestPrimalError_;} ;
539  /// Largest error on basic duals
540  inline double largestDualError() const
541          { return largestDualError_;} ;
542  /// Largest difference between input primal solution and computed
543  inline double largestSolutionError() const
544          { return largestSolutionError_;} ;
545  /// Basic variables pivoting on which rows
546  inline int * pivotVariable() const
547          { return pivotVariable_;};
548  /// If automatic scaling on
549  inline bool automaticScaling() const
550  { return automaticScale_!=0;};
551  inline void setAutomaticScaling(bool onOff)
552  { automaticScale_ = onOff ? 1: 0;}; 
553  /// Current dual tolerance
554  inline double currentDualTolerance() const 
555          { return dualTolerance_;} ;
556  inline void setCurrentDualTolerance(double value)
557          { dualTolerance_ = value;} ;
558  /// Current primal tolerance
559  inline double currentPrimalTolerance() const 
560          { return primalTolerance_;} ;
561  inline void setCurrentPrimalTolerance(double value)
562          { primalTolerance_ = value;} ;
563  /// How many iterative refinements to do
564  inline int numberRefinements() const 
565          { return numberRefinements_;} ;
566  void setNumberRefinements( int value) ;
567  /// Alpha (pivot element) for use by classes e.g. steepestedge
568  inline double alpha() const { return alpha_;};
569  inline void setAlpha(double value) { alpha_ = value;};
570  /// Reduced cost of last incoming for use by classes e.g. steepestedge
571  inline double dualIn() const { return dualIn_;};
572  /// Pivot Row for use by classes e.g. steepestedge
573  inline int pivotRow() const{ return pivotRow_;};
574  inline void setPivotRow(int value) { pivotRow_=value;};
575  /// value of incoming variable (in Dual)
576  double valueIncomingDual() const;
577  //@}
578
579  protected:
580  /**@name protected methods */
581  //@{
582  /** May change basis and then returns number changed.
583      Computation of solutions may be overriden by given pi and solution
584  */
585  int gutsOfSolution ( double * givenDuals,
586                       const double * givenPrimals,
587                       bool valuesPass=false);
588  /// Does most of deletion (0 = all, 1 = most, 2 most + factorization)
589  void gutsOfDelete(int type);
590  /// Does most of copying
591  void gutsOfCopy(const ClpSimplex & rhs);
592  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
593      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
594      8 bit does solution scaling in
595      16 bit does rowArray and columnArray indexed vectors
596      and makes row copy if wanted, also sets columnStart_ etc
597      Also creates scaling arrays if needed.  It does scaling if needed.
598      16 also moves solutions etc in to work arrays
599      On 16 returns false if problem "bad" i.e. matrix or bounds bad
600  */
601  bool createRim(int what,bool makeRowCopy=false);
602  /** releases above arrays and does solution scaling out.  May also
603      get rid of factorization data -
604      0 get rid of nothing, 1 get rid of arrays, 2 also factorization
605  */
606  void deleteRim(int getRidOfFactorizationData=2);
607  /// Sanity check on input rim data (after scaling) - returns true if okay
608  bool sanityCheck();
609  //@}
610  public:
611  /**@name public methods */
612  //@{
613  /** Return row or column sections - not as much needed as it
614      once was.  These just map into single arrays */
615  inline double * solutionRegion(int section) const
616  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
617  inline double * djRegion(int section) const
618  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
619  inline double * lowerRegion(int section) const
620  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
621  inline double * upperRegion(int section) const
622  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
623  inline double * costRegion(int section) const
624  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
625  /// Return region as single array
626  inline double * solutionRegion() const
627  { return solution_;};
628  inline double * djRegion() const
629  { return dj_;};
630  inline double * lowerRegion() const
631  { return lower_;};
632  inline double * upperRegion() const
633  { return upper_;};
634  inline double * costRegion() const
635  { return cost_;};
636  inline Status getStatus(int sequence) const
637  {return static_cast<Status> (status_[sequence]&7);};
638  inline void setStatus(int sequence, Status status)
639  {
640    unsigned char & st_byte = status_[sequence];
641    st_byte &= ~7;
642    st_byte |= status;
643  };
644  /** Normally the first factorization does sparse coding because
645      the factorization could be singular.  This allows initial dense
646      factorization when it is known to be safe
647  */
648  void setInitialDenseFactorization(bool onOff);
649  bool  initialDenseFactorization() const;
650  /** Return sequence In or Out */
651  inline int sequenceIn() const
652  {return sequenceIn_;};
653  inline int sequenceOut() const
654  {return sequenceOut_;};
655  /** Set sequenceIn or Out */
656  inline void  setSequenceIn(int sequence)
657  { sequenceIn_=sequence;};
658  inline void  setSequenceOut(int sequence)
659  { sequenceOut_=sequence;};
660  /** Return direction In or Out */
661  inline int directionIn() const
662  {return directionIn_;};
663  inline int directionOut() const
664  {return directionOut_;};
665  /** Set directionIn or Out */
666  inline void  setDirectionIn(int direction)
667  { directionIn_=direction;};
668  inline void  setDirectionOut(int direction)
669  { directionOut_=direction;};
670  /// Value of Out variable
671  inline double valueOut() const
672  { return valueOut_;};
673  /// Returns 1 if sequence indicates column
674  inline int isColumn(int sequence) const
675  { return sequence<numberColumns_ ? 1 : 0;};
676  /// Returns sequence number within section
677  inline int sequenceWithin(int sequence) const
678  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
679  /// Return row or column values
680  inline double solution(int sequence)
681  { return solution_[sequence];};
682  /// Return address of row or column values
683  inline double & solutionAddress(int sequence)
684  { return solution_[sequence];};
685  inline double reducedCost(int sequence)
686   { return dj_[sequence];};
687  inline double & reducedCostAddress(int sequence)
688   { return dj_[sequence];};
689  inline double lower(int sequence)
690  { return lower_[sequence];};
691  /// Return address of row or column lower bound
692  inline double & lowerAddress(int sequence)
693  { return lower_[sequence];};
694  inline double upper(int sequence)
695  { return upper_[sequence];};
696  /// Return address of row or column upper bound
697  inline double & upperAddress(int sequence)
698  { return upper_[sequence];};
699  inline double cost(int sequence)
700  { return cost_[sequence];};
701  /// Return address of row or column cost
702  inline double & costAddress(int sequence)
703  { return cost_[sequence];};
704  /// Return original lower bound
705  inline double originalLower(int iSequence) const
706  { if (iSequence<numberColumns_) return columnLower_[iSequence]; else
707    return rowLower_[iSequence-numberColumns_];};
708  /// Return original lower bound
709  inline double originalUpper(int iSequence) const
710  { if (iSequence<numberColumns_) return columnUpper_[iSequence]; else
711    return rowUpper_[iSequence-numberColumns_];};
712  /// Theta (pivot change)
713  inline double theta() const
714  { return theta_;};
715  /// Return pointer to details of costs
716  inline ClpNonLinearCost * nonLinearCost() const
717  { return nonLinearCost_;};
718  //@}
719  /**@name status methods */
720  //@{
721  inline void setFakeBound(int sequence, FakeBound fakeBound)
722  {
723    unsigned char & st_byte = status_[sequence];
724    st_byte &= ~24;
725    st_byte |= fakeBound<<3;
726  };
727  inline FakeBound getFakeBound(int sequence) const
728  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
729  inline void setRowStatus(int sequence, Status status)
730  {
731    unsigned char & st_byte = status_[sequence+numberColumns_];
732    st_byte &= ~7;
733    st_byte |= status;
734  };
735  inline Status getRowStatus(int sequence) const
736  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
737  inline void setColumnStatus(int sequence, Status status)
738  {
739    unsigned char & st_byte = status_[sequence];
740    st_byte &= ~7;
741    st_byte |= status;
742  };
743  inline Status getColumnStatus(int sequence) const
744  {return static_cast<Status> (status_[sequence]&7);};
745  inline void setPivoted( int sequence)
746  { status_[sequence] |= 32;};
747  inline void clearPivoted( int sequence)
748  { status_[sequence] &= ~32; };
749  inline bool pivoted(int sequence) const
750  {return (((status_[sequence]>>5)&1)!=0);};
751  /// To flag a variable (not inline to allow for column generation)
752  void setFlagged( int sequence);
753  inline void clearFlagged( int sequence)
754  {
755    status_[sequence] &= ~64;
756  };
757  inline bool flagged(int sequence) const
758  {return ((status_[sequence]&64)!=0);};
759  /// To say row active in primal pivot row choice
760  inline void setActive( int iRow)
761  {
762    status_[iRow] |= 128;
763  };
764  inline void clearActive( int iRow)
765  {
766    status_[iRow] &= ~128;
767  };
768  inline bool active(int iRow) const
769  {return ((status_[iRow]&128)!=0);};
770  /** Set up status array (can be used by OsiClp).
771      Also can be used to set up all slack basis */
772  void createStatus() ;
773  inline void allSlackBasis()
774  { createStatus();};
775   
776  /// So we know when to be cautious
777  inline int lastBadIteration() const
778  {return lastBadIteration_;};
779  /// Progress flag - at present 0 bit says artificials out
780  inline int progressFlag() const
781  {return progressFlag_;};
782  /// Force re-factorization early
783  inline void forceFactorization(int value)
784  { forceFactorization_ = value;};
785  /// Raw objective value (so always minimize in primal)
786  inline double rawObjectiveValue() const
787  { return objectiveValue_;};
788  /** Number of extra rows.  These are ones which will be dynamically created
789      each iteration.  This is for GUB but may have other uses.
790  */
791  inline int numberExtraRows() const
792  { return numberExtraRows_;};
793  /** Maximum number of basic variables - can be more than number of rows if GUB
794  */
795  inline int maximumBasic() const
796  { return maximumBasic_;};
797  /** For advanced options
798      1 - Don't keep changing infeasibility weight
799      2 - Keep nonLinearCost round solves
800      4 - Force outgoing variables to exact bound (primal)
801      8 - Safe to use dense initial factorization
802      16 -Just use basic variables for operation if column generation
803      32 -Clean up with primal before strong branching
804  */
805  inline unsigned int specialOptions() const
806  { return specialOptions_;};
807  inline void setSpecialOptions(unsigned int value)
808  { specialOptions_=value;};
809  //@}
810
811////////////////// data //////////////////
812protected:
813
814  /**@name data.  Many arrays have a row part and a column part.
815   There is a single array with both - columns then rows and
816   then normally two arrays pointing to rows and columns.  The
817   single array is the owner of memory
818  */
819  //@{
820  /// Worst column primal infeasibility
821  double columnPrimalInfeasibility_;
822  /// Worst row primal infeasibility
823  double rowPrimalInfeasibility_;
824  /// Sequence of worst (-1 if feasible)
825  int columnPrimalSequence_;
826  /// Sequence of worst (-1 if feasible)
827  int rowPrimalSequence_;
828  /// Worst column dual infeasibility
829  double columnDualInfeasibility_;
830  /// Worst row dual infeasibility
831  double rowDualInfeasibility_;
832  /// Sequence of worst (-1 if feasible)
833  int columnDualSequence_;
834  /// Sequence of worst (-1 if feasible)
835  int rowDualSequence_;
836  /// Primal tolerance needed to make dual feasible (<largeTolerance)
837  double primalToleranceToGetOptimal_;
838  /// Remaining largest dual infeasibility
839  double remainingDualInfeasibility_;
840  /// Large bound value (for complementarity etc)
841  double largeValue_;
842  /// Largest error on Ax-b
843  double largestPrimalError_;
844  /// Largest error on basic duals
845  double largestDualError_;
846  /// Largest difference between input primal solution and computed
847  double largestSolutionError_;
848  /// Dual bound
849  double dualBound_;
850  /// Alpha (pivot element)
851  double alpha_;
852  /// Theta (pivot change)
853  double theta_;
854  /// Lower Bound on In variable
855  double lowerIn_;
856  /// Value of In variable
857  double valueIn_;
858  /// Upper Bound on In variable
859  double upperIn_;
860  /// Reduced cost of In variable
861  double dualIn_;
862  /// Lower Bound on Out variable
863  double lowerOut_;
864  /// Value of Out variable
865  double valueOut_;
866  /// Upper Bound on Out variable
867  double upperOut_;
868  /// Infeasibility (dual) or ? (primal) of Out variable
869  double dualOut_;
870  /// Current dual tolerance for algorithm
871  double dualTolerance_;
872  /// Current primal tolerance for algorithm
873  double primalTolerance_;
874  /// Sum of dual infeasibilities
875  double sumDualInfeasibilities_;
876  /// Sum of primal infeasibilities
877  double sumPrimalInfeasibilities_;
878  /// Weight assigned to being infeasible in primal
879  double infeasibilityCost_;
880  /// Sum of Dual infeasibilities using tolerance based on error in duals
881  double sumOfRelaxedDualInfeasibilities_;
882  /// Sum of Primal infeasibilities using tolerance based on error in primals
883  double sumOfRelaxedPrimalInfeasibilities_;
884  /// Working copy of lower bounds (Owner of arrays below)
885  double * lower_;
886  /// Row lower bounds - working copy
887  double * rowLowerWork_;
888  /// Column lower bounds - working copy
889  double * columnLowerWork_;
890  /// Working copy of upper bounds (Owner of arrays below)
891  double * upper_;
892  /// Row upper bounds - working copy
893  double * rowUpperWork_;
894  /// Column upper bounds - working copy
895  double * columnUpperWork_;
896  /// Working copy of objective (Owner of arrays below)
897  double * cost_;
898  /// Row objective - working copy
899  double * rowObjectiveWork_;
900  /// Column objective - working copy
901  double * objectiveWork_;
902  /// Useful row length arrays
903  CoinIndexedVector * rowArray_[6];
904  /// Useful column length arrays
905  CoinIndexedVector * columnArray_[6];
906  /// Sequence of In variable
907  int sequenceIn_;
908  /// Direction of In, 1 going up, -1 going down, 0 not a clude
909  int directionIn_;
910  /// Sequence of Out variable
911  int sequenceOut_;
912  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
913  int directionOut_;
914  /// Pivot Row
915  int pivotRow_;
916  /// Last good iteration (immediately after a re-factorization)
917  int lastGoodIteration_;
918  /// Working copy of reduced costs (Owner of arrays below)
919  double * dj_;
920  /// Reduced costs of slacks not same as duals (or - duals)
921  double * rowReducedCost_;
922  /// Possible scaled reduced costs
923  double * reducedCostWork_;
924  /// Working copy of primal solution (Owner of arrays below)
925  double * solution_;
926  /// Row activities - working copy
927  double * rowActivityWork_;
928  /// Column activities - working copy
929  double * columnActivityWork_;
930  /// Number of dual infeasibilities
931  int numberDualInfeasibilities_;
932  /// Number of dual infeasibilities (without free)
933  int numberDualInfeasibilitiesWithoutFree_;
934  /// Number of primal infeasibilities
935  int numberPrimalInfeasibilities_;
936  /// How many iterative refinements to do
937  int numberRefinements_;
938  /// dual row pivot choice
939  ClpDualRowPivot * dualRowPivot_;
940  /// primal column pivot choice
941  ClpPrimalColumnPivot * primalColumnPivot_;
942  /// Basic variables pivoting on which rows
943  int * pivotVariable_;
944  /// factorization
945  ClpFactorization * factorization_;
946  /// Saved version of solution
947  double * savedSolution_;
948  /// Number of times code has tentatively thought optimal
949  int numberTimesOptimal_;
950  /// If change has been made (first attempt at stopping looping)
951  int changeMade_;
952  /// Algorithm >0 == Primal, <0 == Dual
953  int algorithm_;
954  /** Now for some reliability aids
955      This forces re-factorization early */
956  int forceFactorization_;
957  /** Perturbation:
958      -50 to +50 - perturb by this power of ten (-6 sounds good)
959      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
960      101 - we are perturbed
961      102 - don't try perturbing again
962      default is 100
963  */
964  int perturbation_;
965  /// Saved status regions
966  unsigned char * saveStatus_;
967  /** Very wasteful way of dealing with infeasibilities in primal.
968      However it will allow non-linearities and use of dual
969      analysis.  If it doesn't work it can easily be replaced.
970  */
971  ClpNonLinearCost * nonLinearCost_;
972  /** For advanced options
973      1 - Don't keep changing infeasibility weight
974      2 - Keep nonLinearCost round solves
975      4 - Force outgoing variables to exact bound (primal)
976      8 - Safe to use dense initial factorization
977      16 -Just use basic variables for operation
978  */
979  unsigned int specialOptions_;
980  /// So we know when to be cautious
981  int lastBadIteration_;
982  /// So we know when to open up again
983  int lastFlaggedIteration_;
984  /// Can be used for count of fake bounds (dual) or fake costs (primal)
985  int numberFake_;
986  /// Progress flag - at present 0 bit says artificials out, 1 free in
987  int progressFlag_;
988  /// First free/super-basic variable (-1 if none)
989  int firstFree_;
990  /** Number of extra rows.  These are ones which will be dynamically created
991      each iteration.  This is for GUB but may have other uses.
992  */
993  int numberExtraRows_;
994  /** Maximum number of basic variables - can be more than number of rows if GUB
995  */
996  int maximumBasic_;
997  /** For advanced use.  When doing iterative solves things can get
998      nasty so on values pass if incoming solution has largest
999      infeasibility < incomingInfeasibility throw out variables
1000      from basis until largest infeasibility < allowedInfeasibility.
1001      if allowedInfeasibility>= incomingInfeasibility this is
1002      always possible altough you may end up with an all slack basis.
1003
1004      Defaults are 1.0,10.0
1005  */
1006  float incomingInfeasibility_;
1007  float allowedInfeasibility_;
1008  /// Automatic scaling of objective and rhs and bounds
1009  int automaticScale_;
1010  /// For dealing with all issues of cycling etc
1011  ClpSimplexProgress * progress_;
1012  //@}
1013};
1014//#############################################################################
1015/** A function that tests the methods in the ClpSimplex class. The
1016    only reason for it not to be a member method is that this way it doesn't
1017    have to be compiled into the library. And that's a gain, because the
1018    library should be compiled with optimization on, but this method should be
1019    compiled with debugging.
1020
1021    It also does some testing of ClpFactorization class
1022 */
1023void
1024ClpSimplexUnitTest(const std::string & mpsDir,
1025                   const std::string & netlibDir);
1026
1027
1028/// For saving extra information to see if looping.
1029class ClpSimplexProgress {
1030
1031public:
1032
1033
1034  /**@name Constructors and destructor and copy */
1035  //@{
1036  /// Default constructor
1037    ClpSimplexProgress (  );
1038
1039  /// Constructor from model
1040    ClpSimplexProgress ( ClpSimplex * model );
1041
1042  /// Copy constructor.
1043  ClpSimplexProgress(const ClpSimplexProgress &);
1044
1045  /// Assignment operator. This copies the data
1046    ClpSimplexProgress & operator=(const ClpSimplexProgress & rhs);
1047  /// Destructor
1048   ~ClpSimplexProgress (  );
1049  //@}
1050
1051  /**@name Check progress */
1052  //@{
1053  /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
1054      >=0 if give up and use as problem status
1055  */
1056    int looping (  );
1057  /// Start check at beginning of whileIterating
1058  void startCheck();
1059  /// Returns cycle length in whileIterating
1060  int cycle(int in, int out,int wayIn,int wayOut); 
1061
1062  /// Returns previous objective (if -1) - current if (0)
1063  double lastObjective(int back=1) const;
1064  /// Set real primal infeasibility and move back
1065  void setInfeasibility(double value);
1066  /// Returns real primal infeasibility (if -1) - current if (0)
1067  double lastInfeasibility(int back=1) const;
1068  /// Modify objective e.g. if dual infeasible in dual
1069  void modifyObjective(double value);
1070  /// Returns previous iteration number (if -1) - current if (0)
1071  int lastIterationNumber(int back=1) const;
1072  /// Odd state
1073  inline void newOddState()
1074  { oddState_= - oddState_-1;};
1075  inline void endOddState()
1076  { oddState_=abs(oddState_);};
1077  inline void clearOddState() 
1078  { oddState_=0;};
1079  inline int oddState() const
1080  { return oddState_;};
1081  /// number of bad times
1082  inline int badTimes() const
1083  { return numberBadTimes_;};
1084  inline void clearBadTimes()
1085  { numberBadTimes_=0;};
1086
1087  //@}
1088  /**@name Data  */
1089#define CLP_PROGRESS 5
1090  //@{
1091  /// Objective values
1092  double objective_[CLP_PROGRESS];
1093  /// Sum of infeasibilities for algorithm
1094  double infeasibility_[CLP_PROGRESS];
1095  /// Sum of real primal infeasibilities for primal
1096  double realInfeasibility_[CLP_PROGRESS];
1097#define CLP_CYCLE 12
1098  /// For cycle checking
1099  //double obj_[CLP_CYCLE];
1100  int in_[CLP_CYCLE];
1101  int out_[CLP_CYCLE];
1102  char way_[CLP_CYCLE];
1103  /// Pointer back to model so we can get information
1104  ClpSimplex * model_;
1105  /// Number of infeasibilities
1106  int numberInfeasibilities_[CLP_PROGRESS];
1107  /// Iteration number at which occurred
1108  int iterationNumber_[CLP_PROGRESS];
1109  /// Number of times checked (so won't stop too early)
1110  int numberTimes_;
1111  /// Number of times it looked like loop
1112  int numberBadTimes_;
1113  /// If things are in an odd state
1114  int oddState_;
1115  //@}
1116};
1117#endif
Note: See TracBrowser for help on using the repository browser.