1/* $Id: ClpCholeskyBase.hpp 1753 2011-06-19 16:27:26Z stefan $ */
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#ifndef ClpCholeskyBase_H
7#define ClpCholeskyBase_H
8
9#include "CoinPragma.hpp"
10#include "CoinTypes.hpp"
11//#define CLP_LONG_CHOLESKY 0
12#ifndef CLP_LONG_CHOLESKY
13#define CLP_LONG_CHOLESKY 0
14#endif
15/* valid combinations are
16 CLP_LONG_CHOLESKY 0 and COIN_LONG_WORK 0
17 CLP_LONG_CHOLESKY 1 and COIN_LONG_WORK 1
18 CLP_LONG_CHOLESKY 2 and COIN_LONG_WORK 1
19*/
20#if COIN_LONG_WORK==0
21#if CLP_LONG_CHOLESKY>0
22#define CHOLESKY_BAD_COMBINATION
23#endif
24#else
25#if CLP_LONG_CHOLESKY==0
26#define CHOLESKY_BAD_COMBINATION
27#endif
28#endif
29#ifdef CHOLESKY_BAD_COMBINATION
30# warning("Bad combination of CLP_LONG_CHOLESKY and COIN_BIG_DOUBLE/COIN_LONG_WORK");
31"Bad combination of CLP_LONG_CHOLESKY and COIN_LONG_WORK"
32#endif
33#if CLP_LONG_CHOLESKY>1
34typedef long double longDouble;
35#define CHOL_SMALL_VALUE 1.0e-15
36#elif CLP_LONG_CHOLESKY==1
37typedef double longDouble;
38#define CHOL_SMALL_VALUE 1.0e-11
39#else
40typedef double longDouble;
41#define CHOL_SMALL_VALUE 1.0e-11
42#endif
43class ClpInterior;
44class ClpCholeskyDense;
45class ClpMatrixBase;
46
47/** Base class for Clp Cholesky factorization
48 Will do better factorization. very crude ordering
49
50 Derived classes may be using more sophisticated methods
51*/
52
53class ClpCholeskyBase {
54
55public:
56 /**@name Virtual methods that the derived classes may provide */
57 //@{
58 /** Orders rows and saves pointer to matrix.and model.
59 returns non-zero if not enough memory.
60 You can use preOrder to set up ADAT
61 If using default symbolic etc then must set sizeFactor_ to
62 size of input matrix to order (and to symbolic).
63 Also just permute_ and permuteInverse_ should be created */
64 virtual int order(ClpInterior * model);
65 /** Does Symbolic factorization given permutation.
66 This is called immediately after order. If user provides this then
67 user must provide factorize and solve. Otherwise the default factorization is used
68 returns non-zero if not enough memory */
69 virtual int symbolic();
70 /** Factorize - filling in rowsDropped and returning number dropped.
71 If return code negative then out of memory */
72 virtual int factorize(const CoinWorkDouble * diagonal, int * rowsDropped) ;
73 /** Uses factorization to solve. */
74 virtual void solve (CoinWorkDouble * region) ;
75 /** Uses factorization to solve. - given as if KKT.
76 region1 is rows+columns, region2 is rows */
77 virtual void solveKKT (CoinWorkDouble * region1, CoinWorkDouble * region2, const CoinWorkDouble * diagonal,
78 CoinWorkDouble diagonalScaleFactor);
79private:
80 /// AMD ordering
81 int orderAMD();
82public:
83 //@}
84
85 /**@name Gets */
86 //@{
87 /// status. Returns status
88 inline int status() const {
89 return status_;
90 }
91 /// numberRowsDropped. Number of rows gone
92 inline int numberRowsDropped() const {
93 return numberRowsDropped_;
94 }
95 /// reset numberRowsDropped and rowsDropped.
96 void resetRowsDropped();
97 /// rowsDropped - which rows are gone
98 inline char * rowsDropped() const {
99 return rowsDropped_;
100 }
101 /// choleskyCondition.
102 inline double choleskyCondition() const {
103 return choleskyCondition_;
104 }
105 /// goDense i.e. use dense factoriaztion if > this (default 0.7).
106 inline double goDense() const {
107 return goDense_;
108 }
109 /// goDense i.e. use dense factoriaztion if > this (default 0.7).
110 inline void setGoDense(double value) {
111 goDense_ = value;
112 }
113 /// rank. Returns rank
114 inline int rank() const {
115 return numberRows_ - numberRowsDropped_;
116 }
117 /// Return number of rows
118 inline int numberRows() const {
119 return numberRows_;
120 }
121 /// Return size
122 inline CoinBigIndex size() const {
123 return sizeFactor_;
124 }
125 /// Return sparseFactor
126 inline longDouble * sparseFactor() const {
127 return sparseFactor_;
128 }
129 /// Return diagonal
130 inline longDouble * diagonal() const {
131 return diagonal_;
132 }
133 /// Return workDouble
134 inline longDouble * workDouble() const {
135 return workDouble_;
136 }
137 /// If KKT on
138 inline bool kkt() const {
139 return doKKT_;
140 }
141 /// Set KKT
142 inline void setKKT(bool yesNo) {
143 doKKT_ = yesNo;
144 }
145 /// Set integer parameter
146 inline void setIntegerParameter(int i, int value) {
147 integerParameters_[i] = value;
148 }
149 /// get integer parameter
150 inline int getIntegerParameter(int i) {
151 return integerParameters_[i];
152 }
153 /// Set double parameter
154 inline void setDoubleParameter(int i, double value) {
155 doubleParameters_[i] = value;
156 }
157 /// get double parameter
158 inline double getDoubleParameter(int i) {
159 return doubleParameters_[i];
160 }
161 //@}
162
163
164public:
165
166 /**@name Constructors, destructor
167 */
168 //@{
169 /** Constructor which has dense columns activated.
170 Default is off. */
171 ClpCholeskyBase(int denseThreshold = -1);
172 /** Destructor (has to be public) */
173 virtual ~ClpCholeskyBase();
174 /// Copy
175 ClpCholeskyBase(const ClpCholeskyBase&);
176 /// Assignment
177 ClpCholeskyBase& operator=(const ClpCholeskyBase&);
178 //@}
179 //@{
180 ///@name Other
181 /// Clone
182 virtual ClpCholeskyBase * clone() const;
183
184 /// Returns type
185 inline int type() const {
186 if (doKKT_) return 100;
187 else return type_;
188 }
189protected:
190 /// Sets type
191 inline void setType(int type) {
192 type_ = type;
193 }
194 /// model.
195 inline void setModel(ClpInterior * model) {
196 model_ = model;
197 }
198 //@}
199
200 /**@name Symbolic, factor and solve */
201 //@{
202 /** Symbolic1 - works out size without clever stuff.
203 Uses upper triangular as much easier.
204 Returns size
205 */
206 int symbolic1(const CoinBigIndex * Astart, const int * Arow);
207 /** Symbolic2 - Fills in indices
208 Uses lower triangular so can do cliques etc
209 */
210 void symbolic2(const CoinBigIndex * Astart, const int * Arow);
211 /** Factorize - filling in rowsDropped and returning number dropped
212 in integerParam.
213 */
214 void factorizePart2(int * rowsDropped) ;
215 /** solve - 1 just first half, 2 just second half - 3 both.
216 If 1 and 2 then diagonal has sqrt of inverse otherwise inverse
217 */
218 void solve(CoinWorkDouble * region, int type);
219 /// Forms ADAT - returns nonzero if not enough memory
220 int preOrder(bool lowerTriangular, bool includeDiagonal, bool doKKT);
221 /// Updates dense part (broken out for profiling)
222 void updateDense(longDouble * d, /*longDouble * work,*/ int * first);
223 //@}
224
225protected:
226 /**@name Data members
227 The data members are protected to allow access for derived classes. */
228 //@{
229 /// type (may be useful) if > 20 do KKT
230 int type_;
231 /// Doing full KKT (only used if default symbolic and factorization)
232 bool doKKT_;
233 /// Go dense at this fraction
234 double goDense_;
235 /// choleskyCondition.
236 double choleskyCondition_;
237 /// model.
238 ClpInterior * model_;
239 /// numberTrials. Number of trials before rejection
240 int numberTrials_;
241 /// numberRows. Number of Rows in factorization
242 int numberRows_;
243 /// status. Status of factorization
244 int status_;
245 /// rowsDropped
246 char * rowsDropped_;
247 /// permute inverse.
248 int * permuteInverse_;
249 /// main permute.
250 int * permute_;
251 /// numberRowsDropped. Number of rows gone
252 int numberRowsDropped_;
253 /// sparseFactor.
254 longDouble * sparseFactor_;
255 /// choleskyStart - element starts
256 CoinBigIndex * choleskyStart_;
257 /// choleskyRow (can be shorter than sparsefactor)
258 int * choleskyRow_;
259 /// Index starts
260 CoinBigIndex * indexStart_;
261 /// Diagonal
262 longDouble * diagonal_;
263 /// double work array
264 longDouble * workDouble_;
265 /// link array
266 int * link_;
267 // Integer work array
268 CoinBigIndex * workInteger_;
269 // Clique information
270 int * clique_;
271 /// sizeFactor.
272 CoinBigIndex sizeFactor_;
273 /// Size of index array
274 CoinBigIndex sizeIndex_;
275 /// First dense row
276 int firstDense_;
277 /// integerParameters
278 int integerParameters_[64];
279 /// doubleParameters;
280 double doubleParameters_[64];
281 /// Row copy of matrix
282 ClpMatrixBase * rowCopy_;
283 /// Dense indicators
284 char * whichDense_;
285 /// Dense columns (updated)
286 longDouble * denseColumn_;
287 /// Dense cholesky
288 ClpCholeskyDense * dense_;
289 /// Dense threshold (for taking out of Cholesky)
290 int denseThreshold_;
291 //@}
292};
293
294#endif
295