source: trunk/include/ClpInterior.hpp @ 268

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

put back

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.1 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4/*
5   Authors
6   
7   John Tomlin (pdco)
8   John Forrest (standard predictor-corrector)
9
10 */
11#ifndef ClpInterior_H
12#define ClpInterior_H
13
14#include <iostream>
15#include <cfloat>
16#include "ClpModel.hpp"
17#include "ClpMatrixBase.hpp"
18#include "ClpSolve.hpp"
19class ClpLsqr;
20class ClpPdcoBase;
21// ******** DATA to be moved into protected section of ClpInterior
22typedef struct{
23  double  atolmin;
24  double  r3norm;
25  double  LSdamp;
26  double* deltay;
27} Info;
28
29typedef struct{
30  double  atolold;
31  double  atolnew;
32  double  r3ratio;
33  int   istop;
34  int   itncg;
35} Outfo;
36 
37typedef struct{
38double  gamma;
39double  delta;
40int MaxIter;
41double  FeaTol;
42double  OptTol;
43double  StepTol;
44double  x0min;
45double  z0min;
46double  mu0;
47int   LSmethod;   // 1=Cholesky    2=QR    3=LSQR
48int   LSproblem;  // See below
49int LSQRMaxIter;
50double  LSQRatol1; // Initial  atol
51double  LSQRatol2; // Smallest atol (unless atol1 is smaller)
52double  LSQRconlim;
53int  wait;
54} Options;
55class Lsqr;
56class ClpCholeskyBase;
57// ***** END
58/** This solves LPs using interior point methods
59
60    It inherits from ClpModel and all its arrays are created at
61    algorithm time.
62
63*/
64
65class ClpInterior : public ClpModel {
66   friend void ClpInteriorUnitTest(const std::string & mpsDir,
67                                  const std::string & netlibDir);
68
69public:
70
71  /**@name Constructors and destructor and copy */
72  //@{
73  /// Default constructor
74    ClpInterior (  );
75
76  /// Copy constructor.
77  ClpInterior(const ClpInterior &);
78  /// Copy constructor from model.
79  ClpInterior(const ClpModel &);
80  /** Subproblem constructor.  A subset of whole model is created from the
81      row and column lists given.  The new order is given by list order and
82      duplicates are allowed.  Name and integer information can be dropped
83  */
84  ClpInterior (const ClpModel * wholeModel,
85              int numberRows, const int * whichRows,
86              int numberColumns, const int * whichColumns,
87              bool dropNames=true, bool dropIntegers=true);
88  /// Assignment operator. This copies the data
89    ClpInterior & operator=(const ClpInterior & rhs);
90  /// Destructor
91   ~ClpInterior (  );
92  // Ones below are just ClpModel with some changes
93  /** Loads a problem (the constraints on the
94        rows are given by lower and upper bounds). If a pointer is 0 then the
95        following values are the default:
96        <ul>
97          <li> <code>colub</code>: all columns have upper bound infinity
98          <li> <code>collb</code>: all columns have lower bound 0
99          <li> <code>rowub</code>: all rows have upper bound infinity
100          <li> <code>rowlb</code>: all rows have lower bound -infinity
101          <li> <code>obj</code>: all variables have 0 objective coefficient
102        </ul>
103    */
104  void loadProblem (  const ClpMatrixBase& matrix,
105                     const double* collb, const double* colub,   
106                     const double* obj,
107                     const double* rowlb, const double* rowub,
108                      const double * rowObjective=NULL);
109  void loadProblem (  const CoinPackedMatrix& matrix,
110                     const double* collb, const double* colub,   
111                     const double* obj,
112                     const double* rowlb, const double* rowub,
113                      const double * rowObjective=NULL);
114
115  /** Just like the other loadProblem() method except that the matrix is
116        given in a standard column major ordered format (without gaps). */
117  void loadProblem (  const int numcols, const int numrows,
118                     const CoinBigIndex* start, const int* index,
119                     const double* value,
120                     const double* collb, const double* colub,   
121                     const double* obj,
122                      const double* rowlb, const double* rowub,
123                      const double * rowObjective=NULL);
124  /// This one is for after presolve to save memory
125  void loadProblem (  const int numcols, const int numrows,
126                     const CoinBigIndex* start, const int* index,
127                      const double* value,const int * length,
128                     const double* collb, const double* colub,   
129                     const double* obj,
130                      const double* rowlb, const double* rowub,
131                      const double * rowObjective=NULL);
132  /// Read an mps file from the given filename
133  int readMps(const char *filename,
134              bool keepNames=false,
135              bool ignoreErrors = false);
136  //@}
137
138  /**@name Functions most useful to user */
139  //@{
140  /** Pdco algorithm - see ClpPdco.hpp for method */
141  int pdco();
142  // ** Temporary version
143  int  pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo);
144  /// Primal-Dual Predictor-Corrector barrier
145  int primalDual();
146  //@}
147
148  /**@name most useful gets and sets */
149  //@{
150  /// If problem is primal feasible
151  inline bool primalFeasible() const
152         { return (sumPrimalInfeasibilities_<=1.0e-5);};
153  /// If problem is dual feasible
154  inline bool dualFeasible() const
155         { return (sumDualInfeasibilities_<=1.0e-5);};
156  /// Current (or last) algorithm
157  inline int algorithm() const 
158  {return algorithm_; } ;
159  /// Set algorithm
160  inline void setAlgorithm(int value)
161  {algorithm_=value; } ;
162  /// Sum of dual infeasibilities
163  inline double sumDualInfeasibilities() const 
164          { return sumDualInfeasibilities_;} ;
165  /// Sum of primal infeasibilities
166  inline double sumPrimalInfeasibilities() const 
167          { return sumPrimalInfeasibilities_;} ;
168  /// diagonalNorm
169  inline double diagonalNorm() const
170  { return diagonalNorm_;};
171  /// linearPerturbation
172  inline double linearPerturbation() const
173  { return linearPerturbation_;};
174  inline void setLinearPerturbation(double value)
175  { linearPerturbation_=value;};
176  /// diagonalPerturbation
177  inline double diagonalPerturbation() const
178  { return diagonalPerturbation_;};
179  inline void setDiagonalPerturbation(double value)
180  { diagonalPerturbation_=value;};
181  //@}
182
183  /**@name most useful gets and sets */
184  //@{
185  /// Largest error on Ax-b
186  inline double largestPrimalError() const
187          { return largestPrimalError_;} ;
188  /// Largest error on basic duals
189  inline double largestDualError() const
190          { return largestDualError_;} ;
191  /// Maximum iterations
192  inline int maximumBarrierIterations() const
193  { return maximumBarrierIterations_;};
194  inline void setMaximumBarrierIterations(int value)
195  { maximumBarrierIterations_=value;};
196  /// Set cholesky (and delete present one)
197  void setCholesky(ClpCholeskyBase * cholesky);
198  //@}
199
200  protected:
201  /**@name protected methods */
202  //@{
203  /// Does most of deletion
204  void gutsOfDelete();
205  /// Does most of copying
206  void gutsOfCopy(const ClpInterior & rhs);
207  /// Returns true if data looks okay, false if not
208  bool createWorkingData();
209  void deleteWorkingData();
210  /// Sanity check on input rim data
211  bool sanityCheck();
212  ///  This does housekeeping
213  int housekeeping();
214  //@}
215  public:
216  /**@name public methods */
217  //@{
218  /// Raw objective value (so always minimize)
219  inline double rawObjectiveValue() const
220  { return objectiveValue_;};
221  /// Returns 1 if sequence indicates column
222  inline int isColumn(int sequence) const
223  { return sequence<numberColumns_ ? 1 : 0;};
224  /// Returns sequence number within section
225  inline int sequenceWithin(int sequence) const
226  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;};
227  /// Checks solution
228  void checkSolution();
229
230  /// To say a variable is fixed
231  inline void setFixed( int sequence)
232  {
233    status_[sequence] |= 1;
234  };
235  inline void clearFixed( int sequence)
236  {
237    status_[sequence] &= ~1;
238  };
239  inline bool fixed(int sequence) const
240  {return ((status_[sequence]&1)!=0);};
241
242  /// To flag a variable
243  inline void setFlagged( int sequence)
244  {
245    status_[sequence] |= 2;
246  };
247  inline void clearFlagged( int sequence)
248  {
249    status_[sequence] &= ~2;
250  };
251  inline bool flagged(int sequence) const
252  {return ((status_[sequence]&2)!=0);};
253
254  /// To say a variable is fixed OR free
255  inline void setFixedOrFree( int sequence)
256  {
257    status_[sequence] |= 4;
258  };
259  inline void clearFixedOrFree( int sequence)
260  {
261    status_[sequence] &= ~4;
262  };
263  inline bool fixedOrFree(int sequence) const
264  {return ((status_[sequence]&4)!=0);};
265
266  /// To say a variable has lower bound
267  inline void setLowerBound( int sequence)
268  {
269    status_[sequence] |= 8;
270  };
271  inline void clearLowerBound( int sequence)
272  {
273    status_[sequence] &= ~8;
274  };
275  inline bool lowerBound(int sequence) const
276  {return ((status_[sequence]&8)!=0);};
277
278  /// To say a variable has upper bound
279  inline void setUpperBound( int sequence)
280  {
281    status_[sequence] |= 16;
282  };
283  inline void clearUpperBound( int sequence)
284  {
285    status_[sequence] &= ~16;
286  };
287  inline bool upperBound(int sequence) const
288  {return ((status_[sequence]&16)!=0);};
289
290  /// To say a variable has fake lower bound
291  inline void setFakeLower( int sequence)
292  {
293    status_[sequence] |= 32;
294  };
295  inline void clearFakeLower( int sequence)
296  {
297    status_[sequence] &= ~32;
298  };
299  inline bool fakeLower(int sequence) const
300  {return ((status_[sequence]&32)!=0);};
301
302  /// To say a variable has fake upper bound
303  inline void setFakeUpper( int sequence)
304  {
305    status_[sequence] |= 64;
306  };
307  inline void clearFakeUpper( int sequence)
308  {
309    status_[sequence] &= ~64;
310  };
311  inline bool fakeUpper(int sequence) const
312  {return ((status_[sequence]&64)!=0);};
313  //@}
314
315////////////////// data //////////////////
316protected:
317
318  /**@name data.  Many arrays have a row part and a column part.
319   There is a single array with both - columns then rows and
320   then normally two arrays pointing to rows and columns.  The
321   single array is the owner of memory
322  */
323  //@{
324  /// Largest error on Ax-b
325  double largestPrimalError_;
326  /// Largest error on basic duals
327  double largestDualError_;
328  /// Sum of dual infeasibilities
329  double sumDualInfeasibilities_;
330  /// Sum of primal infeasibilities
331  double sumPrimalInfeasibilities_;
332  /// Worst complementarity
333  double worstComplementarity_;
334  ///
335public:
336  double xsize_;
337  double zsize_;
338protected:
339  /// Working copy of lower bounds (Owner of arrays below)
340  double * lower_;
341  /// Row lower bounds - working copy
342  double * rowLowerWork_;
343  /// Column lower bounds - working copy
344  double * columnLowerWork_;
345  /// Working copy of upper bounds (Owner of arrays below)
346  double * upper_;
347  /// Row upper bounds - working copy
348  double * rowUpperWork_;
349  /// Column upper bounds - working copy
350  double * columnUpperWork_;
351  /// Working copy of objective
352  double * cost_;
353public:
354  /// Rhs
355  double * rhs_;
356  double * x_;
357  double * y_;
358  double * dj_;
359protected:
360  /// Pointer to Lsqr object
361  ClpLsqr * lsqrObject_;
362  /// Pointer to stuff
363  ClpPdcoBase * pdcoStuff_;
364  /// Below here is standard barrier stuff
365  /// mu.
366  double mu_;
367  /// objectiveNorm.
368  double objectiveNorm_;
369  /// rhsNorm.
370  double rhsNorm_;
371  /// solutionNorm.
372  double solutionNorm_;
373  /// dualObjective.
374  double dualObjective_;
375  /// primalObjective.
376  double primalObjective_;
377  /// diagonalNorm.
378  double diagonalNorm_;
379  /// stepLength
380  double stepLength_;
381  /// linearPerturbation
382  double linearPerturbation_;
383  /// diagonalPerturbation
384  double diagonalPerturbation_;
385  /// targetGap
386  double targetGap_;
387  /// projectionTolerance
388  double projectionTolerance_;
389  /// maximumRHSError.  maximum Ax
390  double maximumRHSError_;
391  /// maximumBoundInfeasibility.
392  double maximumBoundInfeasibility_;
393  /// maximumDualError.
394  double maximumDualError_;
395  /// diagonalScaleFactor.
396  double diagonalScaleFactor_;
397  /// scaleFactor.  For scaling objective
398  double scaleFactor_;
399  /// actualPrimalStep
400  double actualPrimalStep_;
401  /// actualDualStep
402  double actualDualStep_;
403  /// smallestInfeasibility
404  double smallestInfeasibility_;
405  /// historyInfeasibility.
406#define LENGTH_HISTORY 5
407  double historyInfeasibility_[LENGTH_HISTORY];
408  /// complementarityGap.
409  double complementarityGap_;
410  /// baseObjectiveNorm
411  double baseObjectiveNorm_;
412  /// worstDirectionAccuracy
413  double worstDirectionAccuracy_;
414  /// maximumRHSChange
415  double maximumRHSChange_;
416  /// errorRegion. i.e. Ax
417  double * errorRegion_;
418  /// rhsFixRegion.
419  double * rhsFixRegion_;
420  /// updateRegion.
421  double * updateRegion_;
422  /// upperSlack
423  double * upperSlack_;
424  /// lowerSlack
425  double * lowerSlack_;
426  /// diagonal
427  double * diagonal_;
428  /// weights
429  double * weights_;
430  /// solution
431  double * solution_;
432  /// work1 or deltaZ.
433  double * deltaZ_;
434  /// work2 or deltaW.
435  double * deltaW_;
436  /// work3 or deltaS.
437  double * deltaS_;
438  /// work4 or deltaT.
439  double * deltaT_;
440  /// zVec
441  double * zVec_;
442  /// wVec
443  double * wVec_;
444  /// cholesky.
445  ClpCholeskyBase * cholesky_;
446  /// numberComplementarityPairs;
447  int numberComplementarityPairs_;
448  /// Maximum iterations
449  int maximumBarrierIterations_;
450  /// gonePrimalFeasible.
451  bool gonePrimalFeasible_;
452  /// goneDualFeasible.
453  bool goneDualFeasible_;
454  /// Which algorithm being used
455  int algorithm_;
456  //@}
457};
458//#############################################################################
459/** A function that tests the methods in the ClpInterior class. The
460    only reason for it not to be a member method is that this way it doesn't
461    have to be compiled into the library. And that's a gain, because the
462    library should be compiled with optimization on, but this method should be
463    compiled with debugging.
464
465    It also does some testing of ClpFactorization class
466 */
467void
468ClpInteriorUnitTest(const std::string & mpsDir,
469                   const std::string & netlibDir);
470
471
472#endif
Note: See TracBrowser for help on using the repository browser.