source: stable/1.15/Clp/src/ClpInterior.hpp @ 1949

Last change on this file since 1949 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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