source: trunk/include/ClpInterior.hpp @ 369

Last change on this file since 369 was 369, checked in by forrest, 16 years ago

improving interior point code

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