source: trunk/include/ClpSimplex.hpp @ 2

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

Adding Clp to development branch

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.9 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"
17class ClpDualRowPivot;
18class ClpPrimalColumnPivot;
19class ClpFactorization;
20class OsiWarmStartBasis;
21class OsiIndexedVector;
22class ClpNonLinearCost;
23
24/** This solves LPs using the simplex method
25
26    It inherits from ClpModel and all its arrays are created at
27    algorithm time. Originally I tried to work with model arrays
28    but for simplicity of coding I changed to single arrays with
29    structural variables then row variables.  Some coding is still
30    based on old style and needs cleaning up.
31
32    For a description of algorithms:
33
34    for dual see ClpSimplexDual.hpp and at top of ClpSimplexDual.cpp
35    for primal see ClpSimplexPrimal.hpp and at top of ClpSimplexPrimal.cpp
36
37    There is an algorithm data member.  + for primal variations
38    and - for dual variations
39
40*/
41
42class ClpSimplex : public ClpModel {
43   friend void ClpSimplexUnitTest(const std::string & mpsDir,
44                                  const std::string & netlibDir);
45
46public:
47
48  /// enums for status of various sorts (matches OsiWarmStartBasis)
49  enum Status {
50    isFree = 0x00,
51    basic = 0x01,
52    atUpperBound = 0x02,
53    atLowerBound = 0x03,
54    superBasic = 0x04
55  };
56
57  enum FakeBound {
58    noFake = 0x00,
59    bothFake = 0x01,
60    upperFake = 0x02,
61    lowerFake = 0x03
62  };
63
64  /**@name Constructors and destructor and copy */
65  //@{
66  /// Default constructor
67    ClpSimplex (  );
68
69  /// Copy constructor.
70  ClpSimplex(const ClpSimplex &);
71  /// Copy constructor from model.
72  ClpSimplex(const ClpModel &);
73  /// Assignment operator. This copies the data
74    ClpSimplex & operator=(const ClpSimplex & rhs);
75  /// Destructor
76   ~ClpSimplex (  );
77  //@}
78
79  /**@name Functions most useful to user */
80  //@{
81  /** Dual algorithm - see ClpSimplexDual.hpp for method */
82  int dual();
83  /** Primal algorithm - see ClpSimplexPrimal.hpp for method */
84  int primal(int ifValuesPass=0);
85  /// Sets up working basis as a copy of input
86  void setBasis( const OsiWarmStartBasis & basis);
87  /// Passes in factorization
88  void setFactorization( ClpFactorization & factorization);
89  /// Sets or unsets scaling, 0 -off, 1 on, 2 dynamic(later)
90  void scaling(int mode=1);
91  /// Gets scalingFlag
92  inline int scalingFlag() const {return scalingFlag_;};
93  /** Tightens primal bounds to make dual faster.  Unless
94      fixed, bounds are slightly looser than they could be.
95      This is to make dual go faster and is probably not needed
96      with a presolve.  Returns non-zero if problem infeasible
97  */
98  int tightenPrimalBounds();
99  /// Sets row pivot choice algorithm in dual
100  void setDualRowPivotAlgorithm(ClpDualRowPivot & choice);
101  /// Sets column pivot choice algorithm in primal
102  void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
103  //@}
104
105  /**@name most useful gets and sets */
106  //@{
107  /// If problem is primal feasible
108  inline bool primalFeasible() const
109         { return (numberPrimalInfeasibilities_==0);};
110  /// If problem is dual feasible
111  inline bool dualFeasible() const
112         { return (numberDualInfeasibilities_==0);};
113  /// factorization
114  inline ClpFactorization * factorization() const 
115          { return factorization_;};
116  /// Sparsity on or off
117  bool sparseFactorization() const;
118  void setSparseFactorization(bool value);
119  /// Dual bound
120  inline double dualBound() const
121          { return dualBound_;};
122  void setDualBound(double value);
123  /// Infeasibility cost
124  inline double infeasibilityCost() const
125          { return infeasibilityCost_;};
126  void setInfeasibilityCost(double value);
127  /** Amount of print out:
128      0 - none
129      1 - just final
130      2 - just factorizations
131      3 - as 2 plus a bit more
132      4 - verbose
133      above that 8,16,32 etc just for selective debug
134  */
135  /** Perturbation:
136      -50 to +50 - perturb by this power of ten (-6 sounds good)
137      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
138      101 - we are perturbed
139      102 - don't try perturbing again
140      default is 100
141  */
142  inline int perturbation() const
143    { return perturbation_;};
144  void setPerturbation(int value);
145  /// Current (or last) algorithm
146  inline int algorithm() const 
147  {return algorithm_; } ;
148  /// Sum of dual infeasibilities
149  inline double sumDualInfeasibilities() const 
150          { return sumDualInfeasibilities_;} ;
151  /// Number of dual infeasibilities
152  inline int numberDualInfeasibilities() const 
153          { return numberDualInfeasibilities_;} ;
154  /// Sum of primal infeasibilities
155  inline double sumPrimalInfeasibilities() const 
156          { return sumPrimalInfeasibilities_;} ;
157  /// Number of primal infeasibilities
158  inline int numberPrimalInfeasibilities() const 
159          { return numberPrimalInfeasibilities_;} ;
160  /// Warm start
161  OsiWarmStartBasis getBasis() const;
162  //@}
163
164  /******************** End of most useful part **************/
165  /**@name Functions less likely to be useful to casual user */
166  //@{
167  /** Given an existing factorization computes and checks
168      primal and dual solutions.  Uses input arrays for variables at
169      bounds.  Returns feasibility states */
170  int getSolution (  const double * rowActivities,
171                     const double * columnActivities);
172  /** Given an existing factorization computes and checks
173      primal and dual solutions.  Uses current problem arrays for
174      bounds.  Returns feasibility states */
175  int getSolution ();
176  /** Factorizes using current basis. 
177      solveType - 1 iterating, 0 initial, -1 external
178      If 10 added then in primal values pass
179  */
180  int internalFactorize(int solveType);
181  /// Factorizes using current basis. For external use
182  int factorize();
183  /// Computes duals from scratch
184  void computeDuals();
185  /// Computes primals from scratch
186  void computePrimals (  const double * rowActivities,
187                     const double * columnActivities);
188  /**
189     Unpacks one column of the matrix into indexed array
190     Uses sequenceIn_
191     Also applies scaling if needed
192  */
193  void unpack(OsiIndexedVector * rowArray);
194  /**
195     Unpacks one column of the matrix into indexed array
196     Slack if sequence>= numberColumns
197     Also applies scaling if needed
198  */
199  void unpack(OsiIndexedVector * rowArray,int sequence);
200 
201  /**
202      This does basis housekeeping and does values for in/out variables.
203      Can also decide to re-factorize
204  */
205  int housekeeping(double objectiveChange);
206  /** This sets largest infeasibility and most infeasible and sum
207      and number of infeasibilities (Primal) */
208  void checkPrimalSolution(const double * rowActivities=NULL,
209                           const double * columnActivies=NULL);
210  /** This sets largest infeasibility and most infeasible and sum
211      and number of infeasibilities (Dual) */
212  void checkDualSolution();
213  //@}
214  /**@name Matrix times vector methods
215     They can be faster if scalar is +- 1
216     These are covers so user need not worry about scaling
217     Also for simplex I am not using basic/non-basic split */
218  //@{
219    /** Return <code>y + A * x * scalar</code> in <code>y</code>.
220        @precond <code>x<code> must be of size <code>numColumns()</code>
221        @precond <code>y<code> must be of size <code>numRows()</code> */
222   void times(double scalar,
223                       const double * x, double * y) const;
224    /** Return <code>y + x * scalar * A</code> in <code>y</code>.
225        @precond <code>x<code> must be of size <code>numRows()</code>
226        @precond <code>y<code> must be of size <code>numColumns()</code> */
227    void transposeTimes(double scalar,
228                                const double * x, double * y) const ;
229  //@}
230
231  /**@name most useful gets and sets */
232  //@{
233  /// Worst column primal infeasibility
234  inline double columnPrimalInfeasibility() const 
235          { return columnPrimalInfeasibility_;} ;
236  /// Sequence of worst (-1 if feasible)
237  inline int columnPrimalSequence() const 
238          { return columnPrimalSequence_;} ;
239  /// Worst row primal infeasibility
240  inline double rowPrimalInfeasibility() const 
241          { return rowPrimalInfeasibility_;} ;
242  /// Sequence of worst (-1 if feasible)
243  inline int rowPrimalSequence() const 
244          { return rowPrimalSequence_;} ;
245  /** Worst column dual infeasibility (note - these may not be as meaningful
246      if the problem is primal infeasible */
247  inline double columnDualInfeasibility() const 
248          { return columnDualInfeasibility_;} ;
249  /// Sequence of worst (-1 if feasible)
250  inline int columnDualSequence() const 
251          { return columnDualSequence_;} ;
252  /// Worst row dual infeasibility
253  inline double rowDualInfeasibility() const 
254          { return rowDualInfeasibility_;} ;
255  /// Sequence of worst (-1 if feasible)
256  inline int rowDualSequence() const 
257          { return rowDualSequence_;} ;
258  /// Primal tolerance needed to make dual feasible (<largeTolerance)
259  inline double primalToleranceToGetOptimal() const 
260          { return primalToleranceToGetOptimal_;} ;
261  /// Remaining largest dual infeasibility
262  inline double remainingDualInfeasibility() const 
263          { return remainingDualInfeasibility_;} ;
264  /// Large bound value (for complementarity etc)
265  inline double largeValue() const 
266          { return largeValue_;} ;
267  void setLargeValue( double value) ;
268  /// Largest error on Ax-b
269  inline double largestPrimalError() const
270          { return largestPrimalError_;} ;
271  /// Largest error on basic duals
272  inline double largestDualError() const
273          { return largestDualError_;} ;
274  /// Largest difference between input primal solution and computed
275  inline double largestSolutionError() const
276          { return largestSolutionError_;} ;
277  /// Basic variables pivoting on which rows
278  inline const int * pivotVariable() const
279          { return pivotVariable_;};
280  /// Current dual tolerance
281  inline double currentDualTolerance() const 
282          { return dualTolerance_;} ;
283  /// Current primal tolerance
284  inline double currentPrimalTolerance() const 
285          { return primalTolerance_;} ;
286  /// How many iterative refinements to do
287  inline int numberRefinements() const 
288          { return numberRefinements_;} ;
289  void setnumberRefinements( int value) ;
290  /// Alpha (pivot element) for use by classes e.g. steepestedge
291  inline double alpha() const { return alpha_;};
292  /// Reduced cost of last incoming for use by classes e.g. steepestedge
293  inline double dualIn() const { return dualIn_;};
294  /// Pivot Row for use by classes e.g. steepestedge
295  inline int pivotRow() const{ return pivotRow_;};
296  //@}
297
298  protected:
299  /**@name protected methods */
300  //@{
301  /// May change basis and then returns number changed
302  int gutsOfSolution ( const double * rowActivities,
303                       const double * columnActivities,
304                       bool valuesPass=false);
305  /// Does most of deletion (0 = all, 1 = most)
306  void gutsOfDelete(int type);
307  /// Does most of copying
308  void gutsOfCopy(const ClpSimplex & rhs);
309  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
310      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
311      8 bit does solution scaling in
312      16 bit does rowArray and columnArray indexed vectors
313      and makes row copy if wanted, also sets columnStart_ etc
314      Also creates scaling arrays if needed.  It does scaling if needed.
315      16 also moves solutions etc in to work arrays
316  */
317  void createRim(int what,bool makeRowCopy=false);
318  /// releases above arrays and does solution scaling out
319  void deleteRim();
320  //@}
321  public:
322  /**@name public methods */
323  //@{
324  /** Return row or column sections - not as much needed as it
325      once was.  These just map into single arrays */
326  inline double * solutionRegion(int section)
327  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
328  inline double * djRegion(int section)
329  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
330  inline double * lowerRegion(int section)
331  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
332  inline double * upperRegion(int section)
333  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
334  inline double * costRegion(int section)
335  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
336  /// Return region as single array
337  inline double * solutionRegion()
338  { return solution_;};
339  inline double * djRegion()
340  { return dj_;};
341  inline double * lowerRegion()
342  { return lower_;};
343  inline double * upperRegion()
344  { return upper_;};
345  inline double * costRegion()
346  { return cost_;};
347  inline Status getStatus(int sequence) const
348  {return static_cast<Status> (status_[sequence]&7);};
349  inline void setStatus(int sequence, Status status)
350  {
351    unsigned char & st_byte = status_[sequence];
352    st_byte &= ~7;
353    st_byte |= status;
354  };
355  /** Return sequence In or Out */
356  inline int sequenceIn() const
357  {return sequenceIn_;};
358  inline int sequenceOut() const
359  {return sequenceOut_;};
360  /** Set sequenceIn or Out */
361  inline void  setSequenceIn(int sequence)
362  { sequenceIn_=sequence;};
363  inline void  setSequenceOut(int sequence)
364  { sequenceOut_=sequence;};
365  /// Returns 1 if sequence indicates column
366  inline int isColumn(int sequence) const
367  { return sequence<numberColumns_ ? 1 : 0;};
368  /// Returns sequence number within section
369  inline int sequenceWithin(int sequence) const
370  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
371  /// Return row or column values
372  inline double solution(int sequence)
373  { return solution_[sequence];};
374  /// Return address of row or column values
375  inline double & solutionAddress(int sequence)
376  { return solution_[sequence];};
377  inline double reducedCost(int sequence)
378   { return dj_[sequence];};
379  inline double & reducedCostAddress(int sequence)
380   { return dj_[sequence];};
381  inline double lower(int sequence)
382  { return lower_[sequence];};
383  /// Return address of row or column lower bound
384  inline double & lowerAddress(int sequence)
385  { return lower_[sequence];};
386  inline double upper(int sequence)
387  { return upper_[sequence];};
388  /// Return address of row or column upper bound
389  inline double & upperAddress(int sequence)
390  { return upper_[sequence];};
391  inline double cost(int sequence)
392  { return cost_[sequence];};
393  /// Return address of row or column cost
394  inline double & costAddress(int sequence)
395  { return cost_[sequence];};
396  /// Scaling
397  const double * rowScale() const {return rowScale_;};
398  const double * columnScale() const {return columnScale_;};
399  void setRowScale(double * scale) { rowScale_ = scale;};
400  void setColumnScale(double * scale) { columnScale_ = scale;};
401  //@}
402  /**@name status methods */
403  //@{
404  inline void setFakeBound(int sequence, FakeBound fakeBound)
405  {
406    unsigned char & st_byte = status_[sequence];
407    st_byte &= ~24;
408    st_byte |= fakeBound<<3;
409  };
410  inline FakeBound getFakeBound(int sequence) const
411  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
412  inline void setRowStatus(int sequence, Status status)
413  {
414    unsigned char & st_byte = status_[sequence+numberColumns_];
415    st_byte &= ~7;
416    st_byte |= status;
417  };
418  inline Status getRowStatus(int sequence) const
419  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
420  inline void setColumnStatus(int sequence, Status status)
421  {
422    unsigned char & st_byte = status_[sequence];
423    st_byte &= ~7;
424    st_byte |= status;
425  };
426  inline Status getColumnStatus(int sequence) const
427  {return static_cast<Status> (status_[sequence]&7);};
428  inline void setFixed( int sequence)
429  { status_[sequence] |= 32;};
430  inline void clearFixed( int sequence)
431  { status_[sequence] &= ~32; };
432  inline bool fixed(int sequence) const
433  {return (((status_[sequence]>>5)&1)!=0);};
434  inline void setFlagged( int sequence)
435  {
436    status_[sequence] |= 64;
437  };
438  inline void clearFlagged( int sequence)
439  {
440    status_[sequence] &= ~64;
441  };
442  inline bool flagged(int sequence) const
443  {return (((status_[sequence]>>6)&1)!=0);};
444  //@}
445
446////////////////// data //////////////////
447protected:
448
449  /**@name data.  Many arrays have a row part and a column part.
450   There is a single array with both - columns then rows and
451   then normally two arrays pointing to rows and columns.  The
452   single array is the owner of memory
453  */
454  //@{
455  /// Worst column primal infeasibility
456  double columnPrimalInfeasibility_;
457  /// Sequence of worst (-1 if feasible)
458  int columnPrimalSequence_;
459  /// Worst row primal infeasibility
460  double rowPrimalInfeasibility_;
461  /// Sequence of worst (-1 if feasible)
462  int rowPrimalSequence_;
463  /// Worst column dual infeasibility
464  double columnDualInfeasibility_;
465  /// Sequence of worst (-1 if feasible)
466  int columnDualSequence_;
467  /// Worst row dual infeasibility
468  double rowDualInfeasibility_;
469  /// Sequence of worst (-1 if feasible)
470  int rowDualSequence_;
471  /// Primal tolerance needed to make dual feasible (<largeTolerance)
472  double primalToleranceToGetOptimal_;
473  /// Remaining largest dual infeasibility
474  double remainingDualInfeasibility_;
475  /// Large bound value (for complementarity etc)
476  double largeValue_;
477  /// Largest error on Ax-b
478  double largestPrimalError_;
479  /// Largest error on basic duals
480  double largestDualError_;
481  /// Largest difference between input primal solution and computed
482  double largestSolutionError_;
483  /// Dual bound
484  double dualBound_;
485  /// Working copy of lower bounds (Owner of arrays below)
486  double * lower_;
487  /// Row lower bounds - working copy
488  double * rowLowerWork_;
489  /// Column lower bounds - working copy
490  double * columnLowerWork_;
491  /// Working copy of upper bounds (Owner of arrays below)
492  double * upper_;
493  /// Row upper bounds - working copy
494  double * rowUpperWork_;
495  /// Column upper bounds - working copy
496  double * columnUpperWork_;
497  /// Working copy of objective (Owner of arrays below)
498  double * cost_;
499  /// Row objective - working copy
500  double * rowObjectiveWork_;
501  /// Column objective - working copy
502  double * objectiveWork_;
503  /// Useful row length arrays
504  OsiIndexedVector * rowArray_[6];
505  /// Useful column length arrays
506  OsiIndexedVector * columnArray_[6];
507  /// Alpha (pivot element)
508  double alpha_;
509  /// Theta (pivot change)
510  double theta_;
511  /// Lower Bound on In variable
512  double lowerIn_;
513  /// Value of In variable
514  double valueIn_;
515  /// Upper Bound on In variable
516  double upperIn_;
517  /// Reduced cost of In variable
518  double dualIn_;
519  /// Sequence of In variable
520  int sequenceIn_;
521  /// Direction of In, 1 going up, -1 going down, 0 not a clude
522  int directionIn_;
523  /// Lower Bound on Out variable
524  double lowerOut_;
525  /// Value of Out variable
526  double valueOut_;
527  /// Upper Bound on Out variable
528  double upperOut_;
529  /// Infeasibility (dual) or ? (primal) of Out variable
530  double dualOut_;
531  /// Sequence of Out variable
532  int sequenceOut_;
533  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
534  int directionOut_;
535  /// Pivot Row
536  int pivotRow_;
537  /// Status Region (Owner of array below)
538  unsigned char * status_;
539  /// Working copy of reduced costs (Owner of arrays below)
540  double * dj_;
541  /// Reduced costs of slacks not same as duals (or - duals)
542  double * rowReducedCost_;
543  /// Possible scaled reduced costs
544  double * reducedCostWork_;
545  /// Working copy of primal solution (Owner of arrays below)
546  double * solution_;
547  /// Row activities - working copy
548  double * rowActivityWork_;
549  /// Column activities - working copy
550  double * columnActivityWork_;
551  /// Current dual tolerance for algorithm
552  double dualTolerance_;
553  /// Current primal tolerance for algorithm
554  double primalTolerance_;
555  /// Sum of dual infeasibilities
556  double sumDualInfeasibilities_;
557  /// Number of dual infeasibilities
558  int numberDualInfeasibilities_;
559  /// Number of dual infeasibilities (without free)
560  int numberDualInfeasibilitiesWithoutFree_;
561  /// Sum of primal infeasibilities
562  double sumPrimalInfeasibilities_;
563  /// Number of primal infeasibilities
564  int numberPrimalInfeasibilities_;
565  /// dual row pivot choice
566  ClpDualRowPivot * dualRowPivot_;
567  /// primal column pivot choice
568  ClpPrimalColumnPivot * primalColumnPivot_;
569  /// Basic variables pivoting on which rows
570  int * pivotVariable_;
571  /// factorization
572  ClpFactorization * factorization_;
573  /// How many iterative refinements to do
574  int numberRefinements_;
575  /// Row scale factors for matrix
576  // ****** get working simply then make coding more efficient
577  // on full matrix operations
578  double * rowScale_;
579  /// Saved version of solution
580  double * savedSolution_;
581  /// Column scale factors
582  double * columnScale_;
583  /// Scale flag, 0 none, 1 normal, 2 dynamic
584  int scalingFlag_;
585  /// Number of times code has tentatively thought optimal
586  int numberTimesOptimal_;
587  /// If change has been made (first attempt at stopping looping)
588  int changeMade_;
589  /// Algorithm >0 == Primal, <0 == Dual
590  int algorithm_;
591  /** Now for some reliability aids
592      This forces re-factorization early */
593  int forceFactorization_;
594  /// Saved status regions
595  unsigned char * saveStatus_;
596  /** Perturbation:
597      -50 to +50 - perturb by this power of ten (-6 sounds good)
598      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
599      101 - we are perturbed
600      102 - don't try perturbing again
601      default is 100
602  */
603  int perturbation_;
604  /// Weight assigned to being infeasible in primal
605  double infeasibilityCost_;
606  /** Very wasteful way of dealing with infeasibilities in primal.
607      However it will allow non-linearities and use of dual
608      analysis.  If it doesn't work it can easily be replaced.
609  */
610  ClpNonLinearCost * nonLinearCost_;
611  /// For advanced options
612  unsigned int specialOptions_;
613  //@}
614};
615//#############################################################################
616/** A function that tests the methods in the ClpSimplex class. The
617    only reason for it not to be a member method is that this way it doesn't
618    have to be compiled into the library. And that's a gain, because the
619    library should be compiled with optimization on, but this method should be
620    compiled with debugging.
621
622    It also does some testing of ClpFactorization class
623 */
624void
625ClpSimplexUnitTest(const std::string & mpsDir,
626                   const std::string & netlibDir);
627#endif
Note: See TracBrowser for help on using the repository browser.