source: trunk/Clp/src/ClpInterior.hpp @ 1367

Last change on this file since 1367 was 1367, checked in by forrest, 11 years ago

try and improve stability - also long double interior point

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