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

Last change on this file since 1402 was 1370, checked in by forrest, 10 years ago

add ids

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