1/* $Id: ClpInterior.hpp 1665 2011-01-04 17:55:54Z lou $ */
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 = nullptr);
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 = nullptr);
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 = nullptr);
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 = nullptr);
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
571