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

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

BSP changes from 1247 to 1259

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.0 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 double sumDualInfeasibilities() const 
177          { return sumDualInfeasibilities_;} 
178  /// Sum of primal infeasibilities
179  inline double sumPrimalInfeasibilities() const 
180          { return sumPrimalInfeasibilities_;} 
181  /// dualObjective.
182  inline double dualObjective() const
183  { return dualObjective_;}
184  /// primalObjective.
185  inline double primalObjective() const
186  { return primalObjective_;}
187  /// diagonalNorm
188  inline double diagonalNorm() const
189  { return diagonalNorm_;}
190  /// linearPerturbation
191  inline double linearPerturbation() const
192  { return linearPerturbation_;}
193  inline void setLinearPerturbation(double value)
194  { linearPerturbation_=value;}
195  /// diagonalPerturbation
196  inline double diagonalPerturbation() const
197  { return diagonalPerturbation_;}
198  inline void setDiagonalPerturbation(double value)
199  { diagonalPerturbation_=value;}
200  /// gamma
201  inline double gamma() const
202  { return gamma_;}
203  inline void setGamma(double value)
204  { gamma_=value;}
205  /// delta
206  inline double delta() const
207  { return delta_;}
208  inline void setDelta(double value)
209  { delta_=value;}
210  /// ComplementarityGap
211  inline double complementarityGap() const 
212          { return complementarityGap_;} 
213  //@}
214
215  /**@name most useful gets and sets */
216  //@{
217  /// Largest error on Ax-b
218  inline double largestPrimalError() const
219          { return largestPrimalError_;} 
220  /// Largest error on basic duals
221  inline double largestDualError() const
222          { return largestDualError_;} 
223  /// Maximum iterations
224  inline int maximumBarrierIterations() const
225  { return maximumBarrierIterations_;}
226  inline void setMaximumBarrierIterations(int value)
227  { maximumBarrierIterations_=value;}
228  /// Set cholesky (and delete present one)
229  void setCholesky(ClpCholeskyBase * cholesky);
230  /// Return number fixed to see if worth presolving
231  int numberFixed() const;
232  /** fix variables interior says should be.  If reallyFix false then just
233      set values to exact bounds */
234  void fixFixed(bool reallyFix=true);
235  /// Primal erturbation vector
236  inline double * primalR() const
237  { return primalR_;}
238  /// Dual erturbation vector
239  inline double * dualR() const
240  { return dualR_;}
241  //@}
242
243  protected:
244  /**@name protected methods */
245  //@{
246  /// Does most of deletion
247  void gutsOfDelete();
248  /// Does most of copying
249  void gutsOfCopy(const ClpInterior & rhs);
250  /// Returns true if data looks okay, false if not
251  bool createWorkingData();
252  void deleteWorkingData();
253  /// Sanity check on input rim data
254  bool sanityCheck();
255  ///  This does housekeeping
256  int housekeeping();
257  //@}
258  public:
259  /**@name public methods */
260  //@{
261  /// Raw objective value (so always minimize)
262  inline double rawObjectiveValue() const
263  { return objectiveValue_;}
264  /// Returns 1 if sequence indicates column
265  inline int isColumn(int sequence) const
266  { return sequence<numberColumns_ ? 1 : 0;}
267  /// Returns sequence number within section
268  inline int sequenceWithin(int sequence) const
269  { return sequence<numberColumns_ ? sequence : sequence-numberColumns_;}
270  /// Checks solution
271  void checkSolution();
272  /** Modifies djs to allow for quadratic.
273      returns quadratic offset */
274  double quadraticDjs(double * djRegion, const double * solution,
275                      double scaleFactor);
276
277  /// To say a variable is fixed
278  inline void setFixed( int sequence)
279  {
280    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 1) ;
281  }
282  inline void clearFixed( int sequence)
283  {
284    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~1) ;
285  }
286  inline bool fixed(int sequence) const
287  {return ((status_[sequence]&1)!=0);}
288
289  /// To flag a variable
290  inline void setFlagged( int sequence)
291  {
292    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 2) ;
293  }
294  inline void clearFlagged( int sequence)
295  {
296    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~2) ;
297  }
298  inline bool flagged(int sequence) const
299  {return ((status_[sequence]&2)!=0);}
300
301  /// To say a variable is fixed OR free
302  inline void setFixedOrFree( int sequence)
303  {
304    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 4) ;
305  }
306  inline void clearFixedOrFree( int sequence)
307  {
308    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~4) ;
309  }
310  inline bool fixedOrFree(int sequence) const
311  {return ((status_[sequence]&4)!=0);}
312
313  /// To say a variable has lower bound
314  inline void setLowerBound( int sequence)
315  {
316    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 8) ;
317  }
318  inline void clearLowerBound( int sequence)
319  {
320    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~8) ;
321  }
322  inline bool lowerBound(int sequence) const
323  {return ((status_[sequence]&8)!=0);}
324
325  /// To say a variable has upper bound
326  inline void setUpperBound( int sequence)
327  {
328    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 16) ;
329  }
330  inline void clearUpperBound( int sequence)
331  {
332    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~16) ;
333  }
334  inline bool upperBound(int sequence) const
335  {return ((status_[sequence]&16)!=0);}
336
337  /// To say a variable has fake lower bound
338  inline void setFakeLower( int sequence)
339  {
340    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 32) ;
341  }
342  inline void clearFakeLower( int sequence)
343  {
344    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~32) ;
345  }
346  inline bool fakeLower(int sequence) const
347  {return ((status_[sequence]&32)!=0);}
348
349  /// To say a variable has fake upper bound
350  inline void setFakeUpper( int sequence)
351  {
352    status_[sequence] = static_cast<unsigned char>(status_[sequence] | 64) ;
353  }
354  inline void clearFakeUpper( int sequence)
355  {
356    status_[sequence] = static_cast<unsigned char>(status_[sequence] & ~64) ;
357  }
358  inline bool fakeUpper(int sequence) const
359  {return ((status_[sequence]&64)!=0);}
360  //@}
361
362////////////////// data //////////////////
363protected:
364
365  /**@name data.  Many arrays have a row part and a column part.
366   There is a single array with both - columns then rows and
367   then normally two arrays pointing to rows and columns.  The
368   single array is the owner of memory
369  */
370  //@{
371  /// Largest error on Ax-b
372  double largestPrimalError_;
373  /// Largest error on basic duals
374  double largestDualError_;
375  /// Sum of dual infeasibilities
376  double sumDualInfeasibilities_;
377  /// Sum of primal infeasibilities
378  double sumPrimalInfeasibilities_;
379  /// Worst complementarity
380  double worstComplementarity_;
381  ///
382public:
383  double xsize_;
384  double zsize_;
385protected:
386  /// Working copy of lower bounds (Owner of arrays below)
387  double * lower_;
388  /// Row lower bounds - working copy
389  double * rowLowerWork_;
390  /// Column lower bounds - working copy
391  double * columnLowerWork_;
392  /// Working copy of upper bounds (Owner of arrays below)
393  double * upper_;
394  /// Row upper bounds - working copy
395  double * rowUpperWork_;
396  /// Column upper bounds - working copy
397  double * columnUpperWork_;
398  /// Working copy of objective
399  double * cost_;
400public:
401  /// Rhs
402  double * rhs_;
403  double * x_;
404  double * y_;
405  double * dj_;
406protected:
407  /// Pointer to Lsqr object
408  ClpLsqr * lsqrObject_;
409  /// Pointer to stuff
410  ClpPdcoBase * pdcoStuff_;
411  /// Below here is standard barrier stuff
412  /// mu.
413  double mu_;
414  /// objectiveNorm.
415  double objectiveNorm_;
416  /// rhsNorm.
417  double rhsNorm_;
418  /// solutionNorm.
419  double solutionNorm_;
420  /// dualObjective.
421  double dualObjective_;
422  /// primalObjective.
423  double primalObjective_;
424  /// diagonalNorm.
425  double diagonalNorm_;
426  /// stepLength
427  double stepLength_;
428  /// linearPerturbation
429  double linearPerturbation_;
430  /// diagonalPerturbation
431  double diagonalPerturbation_;
432  // gamma from Saunders and Tomlin regularized
433  double gamma_;
434  // delta from Saunders and Tomlin regularized
435  double delta_;
436  /// targetGap
437  double targetGap_;
438  /// projectionTolerance
439  double projectionTolerance_;
440  /// maximumRHSError.  maximum Ax
441  double maximumRHSError_;
442  /// maximumBoundInfeasibility.
443  double maximumBoundInfeasibility_;
444  /// maximumDualError.
445  double maximumDualError_;
446  /// diagonalScaleFactor.
447  double diagonalScaleFactor_;
448  /// scaleFactor.  For scaling objective
449  double scaleFactor_;
450  /// actualPrimalStep
451  double actualPrimalStep_;
452  /// actualDualStep
453  double actualDualStep_;
454  /// smallestInfeasibility
455  double smallestInfeasibility_;
456  /// historyInfeasibility.
457#define LENGTH_HISTORY 5
458  double historyInfeasibility_[LENGTH_HISTORY];
459  /// complementarityGap.
460  double complementarityGap_;
461  /// baseObjectiveNorm
462  double baseObjectiveNorm_;
463  /// worstDirectionAccuracy
464  double worstDirectionAccuracy_;
465  /// maximumRHSChange
466  double maximumRHSChange_;
467  /// errorRegion. i.e. Ax
468  double * errorRegion_;
469  /// rhsFixRegion.
470  double * rhsFixRegion_;
471  /// upperSlack
472  double * upperSlack_;
473  /// lowerSlack
474  double * lowerSlack_;
475  /// diagonal
476  double * diagonal_;
477  /// solution
478  double * solution_;
479  /// work array
480  double * workArray_;
481  /// delta X
482  double * deltaX_;
483  /// delta Y
484  double * deltaY_;
485  /// deltaZ.
486  double * deltaZ_;
487  /// deltaW.
488  double * deltaW_;
489  /// deltaS.
490  double * deltaSU_;
491  double * deltaSL_;
492  /// Primal regularization array
493  double * primalR_;
494  /// Dual regularization array
495  double * dualR_;
496  /// rhs B
497  double * rhsB_;
498  /// rhsU.
499  double * rhsU_;
500  /// rhsL.
501  double * rhsL_;
502  /// rhsZ.
503  double * rhsZ_;
504  /// rhsW.
505  double * rhsW_;
506  /// rhs C
507  double * rhsC_;
508  /// zVec
509  double * zVec_;
510  /// wVec
511  double * wVec_;
512  /// cholesky.
513  ClpCholeskyBase * cholesky_;
514  /// numberComplementarityPairs i.e. ones with lower and/or upper bounds (not fixed)
515  int numberComplementarityPairs_;
516  /// numberComplementarityItems_ i.e. number of active bounds
517  int numberComplementarityItems_;
518  /// Maximum iterations
519  int maximumBarrierIterations_;
520  /// gonePrimalFeasible.
521  bool gonePrimalFeasible_;
522  /// goneDualFeasible.
523  bool goneDualFeasible_;
524  /// Which algorithm being used
525  int algorithm_;
526  //@}
527};
528//#############################################################################
529/** A function that tests the methods in the ClpInterior class. The
530    only reason for it not to be a member method is that this way it doesn't
531    have to be compiled into the library. And that's a gain, because the
532    library should be compiled with optimization on, but this method should be
533    compiled with debugging.
534
535    It also does some testing of ClpFactorization class
536 */
537void
538ClpInteriorUnitTest(const std::string & mpsDir,
539                   const std::string & netlibDir);
540
541
542#endif
Note: See TracBrowser for help on using the repository browser.