source: branches/devel-1/include/ClpSimplex.hpp @ 15

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

Hope this works from wincvs

Fix error in values pass

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.2 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      Return codes are as from ClpFactorization unless initial factorization
180      when total number of singularities is returned
181  */
182  int internalFactorize(int solveType);
183  /// Factorizes using current basis. For external use
184  int factorize();
185  /// Computes duals from scratch
186  void computeDuals();
187  /// Computes primals from scratch
188  void computePrimals (  const double * rowActivities,
189                     const double * columnActivities);
190  /**
191     Unpacks one column of the matrix into indexed array
192     Uses sequenceIn_
193     Also applies scaling if needed
194  */
195  void unpack(OsiIndexedVector * rowArray);
196  /**
197     Unpacks one column of the matrix into indexed array
198     Slack if sequence>= numberColumns
199     Also applies scaling if needed
200  */
201  void unpack(OsiIndexedVector * rowArray,int sequence);
202 
203  /**
204      This does basis housekeeping and does values for in/out variables.
205      Can also decide to re-factorize
206  */
207  int housekeeping(double objectiveChange);
208  /** This sets largest infeasibility and most infeasible and sum
209      and number of infeasibilities (Primal) */
210  void checkPrimalSolution(const double * rowActivities=NULL,
211                           const double * columnActivies=NULL);
212  /** This sets largest infeasibility and most infeasible and sum
213      and number of infeasibilities (Dual) */
214  void checkDualSolution();
215  //@}
216  /**@name Matrix times vector methods
217     They can be faster if scalar is +- 1
218     These are covers so user need not worry about scaling
219     Also for simplex I am not using basic/non-basic split */
220  //@{
221    /** Return <code>y + A * x * scalar</code> in <code>y</code>.
222        @precond <code>x<code> must be of size <code>numColumns()</code>
223        @precond <code>y<code> must be of size <code>numRows()</code> */
224   void times(double scalar,
225                       const double * x, double * y) const;
226    /** Return <code>y + x * scalar * A</code> in <code>y</code>.
227        @precond <code>x<code> must be of size <code>numRows()</code>
228        @precond <code>y<code> must be of size <code>numColumns()</code> */
229    void transposeTimes(double scalar,
230                                const double * x, double * y) const ;
231  //@}
232
233  /**@name most useful gets and sets */
234  //@{
235  /// Worst column primal infeasibility
236  inline double columnPrimalInfeasibility() const 
237          { return columnPrimalInfeasibility_;} ;
238  /// Sequence of worst (-1 if feasible)
239  inline int columnPrimalSequence() const 
240          { return columnPrimalSequence_;} ;
241  /// Worst row primal infeasibility
242  inline double rowPrimalInfeasibility() const 
243          { return rowPrimalInfeasibility_;} ;
244  /// Sequence of worst (-1 if feasible)
245  inline int rowPrimalSequence() const 
246          { return rowPrimalSequence_;} ;
247  /** Worst column dual infeasibility (note - these may not be as meaningful
248      if the problem is primal infeasible */
249  inline double columnDualInfeasibility() const 
250          { return columnDualInfeasibility_;} ;
251  /// Sequence of worst (-1 if feasible)
252  inline int columnDualSequence() const 
253          { return columnDualSequence_;} ;
254  /// Worst row dual infeasibility
255  inline double rowDualInfeasibility() const 
256          { return rowDualInfeasibility_;} ;
257  /// Sequence of worst (-1 if feasible)
258  inline int rowDualSequence() const 
259          { return rowDualSequence_;} ;
260  /// Primal tolerance needed to make dual feasible (<largeTolerance)
261  inline double primalToleranceToGetOptimal() const 
262          { return primalToleranceToGetOptimal_;} ;
263  /// Remaining largest dual infeasibility
264  inline double remainingDualInfeasibility() const 
265          { return remainingDualInfeasibility_;} ;
266  /// Large bound value (for complementarity etc)
267  inline double largeValue() const 
268          { return largeValue_;} ;
269  void setLargeValue( double value) ;
270  /// Largest error on Ax-b
271  inline double largestPrimalError() const
272          { return largestPrimalError_;} ;
273  /// Largest error on basic duals
274  inline double largestDualError() const
275          { return largestDualError_;} ;
276  /// Largest difference between input primal solution and computed
277  inline double largestSolutionError() const
278          { return largestSolutionError_;} ;
279  /// Basic variables pivoting on which rows
280  inline const int * pivotVariable() const
281          { return pivotVariable_;};
282  /// Current dual tolerance
283  inline double currentDualTolerance() const 
284          { return dualTolerance_;} ;
285  /// Current primal tolerance
286  inline double currentPrimalTolerance() const 
287          { return primalTolerance_;} ;
288  /// How many iterative refinements to do
289  inline int numberRefinements() const 
290          { return numberRefinements_;} ;
291  void setnumberRefinements( int value) ;
292  /// Alpha (pivot element) for use by classes e.g. steepestedge
293  inline double alpha() const { return alpha_;};
294  /// Reduced cost of last incoming for use by classes e.g. steepestedge
295  inline double dualIn() const { return dualIn_;};
296  /// Pivot Row for use by classes e.g. steepestedge
297  inline int pivotRow() const{ return pivotRow_;};
298  //@}
299
300  protected:
301  /**@name protected methods */
302  //@{
303  /// May change basis and then returns number changed
304  int gutsOfSolution ( const double * rowActivities,
305                       const double * columnActivities,
306                       bool valuesPass=false);
307  /// Does most of deletion (0 = all, 1 = most)
308  void gutsOfDelete(int type);
309  /// Does most of copying
310  void gutsOfCopy(const ClpSimplex & rhs);
311  /** puts in format I like (rowLower,rowUpper) also see StandardMatrix
312      1 bit does rows, 2 bit does column bounds, 4 bit does objective(s).
313      8 bit does solution scaling in
314      16 bit does rowArray and columnArray indexed vectors
315      and makes row copy if wanted, also sets columnStart_ etc
316      Also creates scaling arrays if needed.  It does scaling if needed.
317      16 also moves solutions etc in to work arrays
318  */
319  void createRim(int what,bool makeRowCopy=false);
320  /// releases above arrays and does solution scaling out
321  void deleteRim();
322  //@}
323  public:
324  /**@name public methods */
325  //@{
326  /** Return row or column sections - not as much needed as it
327      once was.  These just map into single arrays */
328  inline double * solutionRegion(int section)
329  { if (!section) return rowActivityWork_; else return columnActivityWork_;};
330  inline double * djRegion(int section)
331  { if (!section) return rowReducedCost_; else return reducedCostWork_;};
332  inline double * lowerRegion(int section)
333  { if (!section) return rowLowerWork_; else return columnLowerWork_;};
334  inline double * upperRegion(int section)
335  { if (!section) return rowUpperWork_; else return columnUpperWork_;};
336  inline double * costRegion(int section)
337  { if (!section) return rowObjectiveWork_; else return objectiveWork_;};
338  /// Return region as single array
339  inline double * solutionRegion()
340  { return solution_;};
341  inline double * djRegion()
342  { return dj_;};
343  inline double * lowerRegion()
344  { return lower_;};
345  inline double * upperRegion()
346  { return upper_;};
347  inline double * costRegion()
348  { return cost_;};
349  inline Status getStatus(int sequence) const
350  {return static_cast<Status> (status_[sequence]&7);};
351  inline void setStatus(int sequence, Status status)
352  {
353    unsigned char & st_byte = status_[sequence];
354    st_byte &= ~7;
355    st_byte |= status;
356  };
357  /** Return sequence In or Out */
358  inline int sequenceIn() const
359  {return sequenceIn_;};
360  inline int sequenceOut() const
361  {return sequenceOut_;};
362  /** Set sequenceIn or Out */
363  inline void  setSequenceIn(int sequence)
364  { sequenceIn_=sequence;};
365  inline void  setSequenceOut(int sequence)
366  { sequenceOut_=sequence;};
367  /// Returns 1 if sequence indicates column
368  inline int isColumn(int sequence) const
369  { return sequence<numberColumns_ ? 1 : 0;};
370  /// Returns sequence number within section
371  inline int sequenceWithin(int sequence) const
372  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
373  /// Return row or column values
374  inline double solution(int sequence)
375  { return solution_[sequence];};
376  /// Return address of row or column values
377  inline double & solutionAddress(int sequence)
378  { return solution_[sequence];};
379  inline double reducedCost(int sequence)
380   { return dj_[sequence];};
381  inline double & reducedCostAddress(int sequence)
382   { return dj_[sequence];};
383  inline double lower(int sequence)
384  { return lower_[sequence];};
385  /// Return address of row or column lower bound
386  inline double & lowerAddress(int sequence)
387  { return lower_[sequence];};
388  inline double upper(int sequence)
389  { return upper_[sequence];};
390  /// Return address of row or column upper bound
391  inline double & upperAddress(int sequence)
392  { return upper_[sequence];};
393  inline double cost(int sequence)
394  { return cost_[sequence];};
395  /// Return address of row or column cost
396  inline double & costAddress(int sequence)
397  { return cost_[sequence];};
398  /// Scaling
399  const double * rowScale() const {return rowScale_;};
400  const double * columnScale() const {return columnScale_;};
401  void setRowScale(double * scale) { rowScale_ = scale;};
402  void setColumnScale(double * scale) { columnScale_ = scale;};
403  //@}
404  /**@name status methods */
405  //@{
406  inline void setFakeBound(int sequence, FakeBound fakeBound)
407  {
408    unsigned char & st_byte = status_[sequence];
409    st_byte &= ~24;
410    st_byte |= fakeBound<<3;
411  };
412  inline FakeBound getFakeBound(int sequence) const
413  {return static_cast<FakeBound> ((status_[sequence]>>3)&3);};
414  inline void setRowStatus(int sequence, Status status)
415  {
416    unsigned char & st_byte = status_[sequence+numberColumns_];
417    st_byte &= ~7;
418    st_byte |= status;
419  };
420  inline Status getRowStatus(int sequence) const
421  {return static_cast<Status> (status_[sequence+numberColumns_]&7);};
422  inline void setColumnStatus(int sequence, Status status)
423  {
424    unsigned char & st_byte = status_[sequence];
425    st_byte &= ~7;
426    st_byte |= status;
427  };
428  inline Status getColumnStatus(int sequence) const
429  {return static_cast<Status> (status_[sequence]&7);};
430  inline void setFixed( int sequence)
431  { status_[sequence] |= 32;};
432  inline void clearFixed( int sequence)
433  { status_[sequence] &= ~32; };
434  inline bool fixed(int sequence) const
435  {return (((status_[sequence]>>5)&1)!=0);};
436  inline void setFlagged( int sequence)
437  {
438    status_[sequence] |= 64;
439  };
440  inline void clearFlagged( int sequence)
441  {
442    status_[sequence] &= ~64;
443  };
444  inline bool flagged(int sequence) const
445  {return (((status_[sequence]>>6)&1)!=0);};
446  /// So we know when to be cautious
447  inline int lastBadIteration() const
448  {return lastBadIteration_;};
449  //@}
450
451////////////////// data //////////////////
452protected:
453
454  /**@name data.  Many arrays have a row part and a column part.
455   There is a single array with both - columns then rows and
456   then normally two arrays pointing to rows and columns.  The
457   single array is the owner of memory
458  */
459  //@{
460  /// Worst column primal infeasibility
461  double columnPrimalInfeasibility_;
462  /// Sequence of worst (-1 if feasible)
463  int columnPrimalSequence_;
464  /// Worst row primal infeasibility
465  double rowPrimalInfeasibility_;
466  /// Sequence of worst (-1 if feasible)
467  int rowPrimalSequence_;
468  /// Worst column dual infeasibility
469  double columnDualInfeasibility_;
470  /// Sequence of worst (-1 if feasible)
471  int columnDualSequence_;
472  /// Worst row dual infeasibility
473  double rowDualInfeasibility_;
474  /// Sequence of worst (-1 if feasible)
475  int rowDualSequence_;
476  /// Primal tolerance needed to make dual feasible (<largeTolerance)
477  double primalToleranceToGetOptimal_;
478  /// Remaining largest dual infeasibility
479  double remainingDualInfeasibility_;
480  /// Large bound value (for complementarity etc)
481  double largeValue_;
482  /// Largest error on Ax-b
483  double largestPrimalError_;
484  /// Largest error on basic duals
485  double largestDualError_;
486  /// Largest difference between input primal solution and computed
487  double largestSolutionError_;
488  /// Dual bound
489  double dualBound_;
490  /// Working copy of lower bounds (Owner of arrays below)
491  double * lower_;
492  /// Row lower bounds - working copy
493  double * rowLowerWork_;
494  /// Column lower bounds - working copy
495  double * columnLowerWork_;
496  /// Working copy of upper bounds (Owner of arrays below)
497  double * upper_;
498  /// Row upper bounds - working copy
499  double * rowUpperWork_;
500  /// Column upper bounds - working copy
501  double * columnUpperWork_;
502  /// Working copy of objective (Owner of arrays below)
503  double * cost_;
504  /// Row objective - working copy
505  double * rowObjectiveWork_;
506  /// Column objective - working copy
507  double * objectiveWork_;
508  /// Useful row length arrays
509  OsiIndexedVector * rowArray_[6];
510  /// Useful column length arrays
511  OsiIndexedVector * columnArray_[6];
512  /// Alpha (pivot element)
513  double alpha_;
514  /// Theta (pivot change)
515  double theta_;
516  /// Lower Bound on In variable
517  double lowerIn_;
518  /// Value of In variable
519  double valueIn_;
520  /// Upper Bound on In variable
521  double upperIn_;
522  /// Reduced cost of In variable
523  double dualIn_;
524  /// Sequence of In variable
525  int sequenceIn_;
526  /// Direction of In, 1 going up, -1 going down, 0 not a clude
527  int directionIn_;
528  /// Lower Bound on Out variable
529  double lowerOut_;
530  /// Value of Out variable
531  double valueOut_;
532  /// Upper Bound on Out variable
533  double upperOut_;
534  /// Infeasibility (dual) or ? (primal) of Out variable
535  double dualOut_;
536  /// Sequence of Out variable
537  int sequenceOut_;
538  /// Direction of Out, 1 to upper bound, -1 to lower bound, 0 - superbasic
539  int directionOut_;
540  /// Pivot Row
541  int pivotRow_;
542  /// Status Region (Owner of array below)
543  unsigned char * status_;
544  /// Working copy of reduced costs (Owner of arrays below)
545  double * dj_;
546  /// Reduced costs of slacks not same as duals (or - duals)
547  double * rowReducedCost_;
548  /// Possible scaled reduced costs
549  double * reducedCostWork_;
550  /// Working copy of primal solution (Owner of arrays below)
551  double * solution_;
552  /// Row activities - working copy
553  double * rowActivityWork_;
554  /// Column activities - working copy
555  double * columnActivityWork_;
556  /// Current dual tolerance for algorithm
557  double dualTolerance_;
558  /// Current primal tolerance for algorithm
559  double primalTolerance_;
560  /// Sum of dual infeasibilities
561  double sumDualInfeasibilities_;
562  /// Number of dual infeasibilities
563  int numberDualInfeasibilities_;
564  /// Number of dual infeasibilities (without free)
565  int numberDualInfeasibilitiesWithoutFree_;
566  /// Sum of primal infeasibilities
567  double sumPrimalInfeasibilities_;
568  /// Number of primal infeasibilities
569  int numberPrimalInfeasibilities_;
570  /// dual row pivot choice
571  ClpDualRowPivot * dualRowPivot_;
572  /// primal column pivot choice
573  ClpPrimalColumnPivot * primalColumnPivot_;
574  /// Basic variables pivoting on which rows
575  int * pivotVariable_;
576  /// factorization
577  ClpFactorization * factorization_;
578  /// How many iterative refinements to do
579  int numberRefinements_;
580  /// Row scale factors for matrix
581  // ****** get working simply then make coding more efficient
582  // on full matrix operations
583  double * rowScale_;
584  /// Saved version of solution
585  double * savedSolution_;
586  /// Column scale factors
587  double * columnScale_;
588  /// Scale flag, 0 none, 1 normal, 2 dynamic
589  int scalingFlag_;
590  /// Number of times code has tentatively thought optimal
591  int numberTimesOptimal_;
592  /// If change has been made (first attempt at stopping looping)
593  int changeMade_;
594  /// Algorithm >0 == Primal, <0 == Dual
595  int algorithm_;
596  /** Now for some reliability aids
597      This forces re-factorization early */
598  int forceFactorization_;
599  /// Saved status regions
600  unsigned char * saveStatus_;
601  /** Perturbation:
602      -50 to +50 - perturb by this power of ten (-6 sounds good)
603      100 - auto perturb if takes too long (1.0e-6 largest nonzero)
604      101 - we are perturbed
605      102 - don't try perturbing again
606      default is 100
607  */
608  int perturbation_;
609  /// Weight assigned to being infeasible in primal
610  double infeasibilityCost_;
611  /** Very wasteful way of dealing with infeasibilities in primal.
612      However it will allow non-linearities and use of dual
613      analysis.  If it doesn't work it can easily be replaced.
614  */
615  ClpNonLinearCost * nonLinearCost_;
616  /// For advanced options
617  unsigned int specialOptions_;
618  /// So we know when to be cautious
619  int lastBadIteration_;
620  //@}
621};
622//#############################################################################
623/** A function that tests the methods in the ClpSimplex class. The
624    only reason for it not to be a member method is that this way it doesn't
625    have to be compiled into the library. And that's a gain, because the
626    library should be compiled with optimization on, but this method should be
627    compiled with debugging.
628
629    It also does some testing of ClpFactorization class
630 */
631void
632ClpSimplexUnitTest(const std::string & mpsDir,
633                   const std::string & netlibDir);
634#endif
Note: See TracBrowser for help on using the repository browser.