1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24
25 /**@file clufactor.h
26 * @brief Implementation of sparse LU factorization.
27 */
28 #ifndef _CLUFACTOR_H_
29 #define _CLUFACTOR_H_
30
31 #include "soplex/spxdefines.h"
32 #include "soplex/slinsolver.h"
33 #include "soplex/timer.h"
34 #include "soplex/svector.h"
35
36 #include "vector"
37
38 #define SOPLEX_WITH_L_ROWS 1
39
40 namespace soplex
41 {
42 /**@brief Implementation of sparse LU factorization.
43 * @ingroup Algo
44 *
45 * This class implements a sparse LU factorization with either
46 * FOREST-TOMLIN or ETA updates, using dynamic Markowitz pivoting.
47 */
48 template <class R>
49 class CLUFactor
50 {
51 public:
52
53 /// set tolerances
54 void setTolerances(std::shared_ptr<Tolerances> tolerances)
55 {
56 this->_tolerances = tolerances;
57 }
58 /// get tolerances
59 const std::shared_ptr<Tolerances> tolerances() const
60 {
61 return _tolerances;
62 }
63
64 //----------------------------------------
65 /**@name Public types */
66 ///@{
67 /** Doubly linked ring structure for garbage collection of column or
68 * row file in working matrix.
69 */
70 struct Dring
71 {
72 Dring* next;
73 Dring* prev;
74 int idx;
75 };
76
77 /// Pivot Ring
(1) Event rule_of_three_violation: |
Class "soplex::CLUFactor<double>::Pring" has a user definition for at least one special function (copy constructor, copy assignment, destructor) but not all. If one of these functions requires a user definition then the others likely do as well. |
(4) Event remediation: |
Add user-definition for a destructor. |
Also see events: |
[copy_ctor][copy_assign] |
78 class Pring
79 {
80 public:
81 Pring* next; ///<
82 Pring* prev; ///<
83 int idx; ///< index of pivot row
84 int pos; ///< position of pivot column in row
85 int mkwtz; ///< markowitz number of pivot
86
87 Pring() : next(0), prev(0) ///< constructor
88 {
89 mkwtz = -1;
90 idx = -1;
91 pos = -1;
92 }
93
94 private:
95 Pring(const Pring&); ///< blocked copy constructor
96 Pring& operator= (const Pring&); ///< blocked assignment operator
97 };
98 ///@}
99
100 protected:
101
102 //----------------------------------------
103 /**@name Protected types */
104 ///@{
105 /// Temporary data structures.
106 class Temp
107 {
108 public:
109 int* s_mark; ///< marker
110 std::vector<R> s_max; ///< maximum absolute value per row (or -1)
111 int* s_cact; ///< lengths of columns of active submatrix
112 int stage; ///< stage of the structure
113 Pring pivots; ///< ring of selected pivot rows
114 Pring* pivot_col; ///< column index handlers for R linked list
115 Pring* pivot_colNZ; ///< lists for columns to number of nonzeros
116 Pring* pivot_row; ///< row index handlers for R linked list
117 Pring* pivot_rowNZ; ///< lists for rows to number of nonzeros
118
119 Temp(); ///< constructor
120 ~Temp(); ///< destructor
121 void init(int p_dim); ///< initialization
122 void clear(); ///< clears the structure
123
124 private:
125 Temp(const Temp&); ///< blocked copy constructor
126 Temp& operator= (const Temp&); ///< blocked assignment operator
127 };
128
129 /// Data structures for saving the row and column permutations.
130 struct Perm
131 {
132 int* orig; ///< orig[p] original index from p
133 int* perm; ///< perm[i] permuted index from i
134 };
135
136 /// Data structures for saving the working matrix and U factor.
137 struct U
138 {
139 ///
140 struct Row
141 {
142 Dring list; /*!< \brief Double linked ringlist of VectorBase<R>
143 indices in the order they appear
144 in the row file */
145 Dring* elem; ///< %Array of ring elements.
146 int size; ///< size of arrays val and idx
147 int used; ///< used entries of arrays idx and val
148 std::vector<R> val; ///< hold nonzero values
149 int* idx; ///< hold column indices of nonzeros
150 int* start; ///< starting positions in val and idx
151 int* len; ///< used nonzeros per row vectors
152 int* max; /*!< \brief maximum available nonzeros per row:
153 start[i] + max[i] == start[elem[i].next->idx]
154 len[i] <= max[i]. */
155 } row;
156
157 ///
158 struct Col
159 {
160 Dring list; /*!< \brief Double linked ringlist of VectorBase<R>
161 indices in the order they appear
162 in the column file */
163 Dring* elem; ///< %Array of ring elements.
164 int size; ///< size of array idx
165 int used; ///< used entries of array idx
166 int* idx; ///< hold row indices of nonzeros
167 std::vector<R> val; /*!< \brief hold nonzero values: this is only initialized
168 in the end of the factorization with DEFAULT
169 updates. */
170 int* start; ///< starting positions in val and idx
171 int* len; ///< used nonzeros per column vector
172 int* max; /*!< \brief maximum available nonzeros per colunn:
173 start[i] + max[i] == start[elem[i].next->idx]
174 len[i] <= max[i]. */
175 } col;
176 };
177
178
179 /// Data structures for saving the working matrix and L factor.
180 struct L
181 {
182 int size; ///< size of arrays val and idx
183 std::vector<R> val; ///< values of L vectors
184 int* idx; ///< indices of L vectors
185 int startSize; ///< size of array start
186 int firstUpdate; ///< number of first update L vector
187 int firstUnused; ///< number of first unused L vector
188 int* start; ///< starting positions in val and idx
189 int* row; ///< column indices of L vectors
190 int updateType; ///< type of updates to be used.
191
192 /* The following arrays have length |firstUpdate|, since they keep
193 * rows of the L-vectors occuring during the factorization (without
194 * updates), only:
195 */
196 std::vector<R> rval; ///< values of rows of L
197 int* ridx; ///< indices of rows of L
198 int* rbeg; ///< start of rows in rval and ridx
199 int* rorig; ///< original row permutation
200 int* rperm; ///< original row permutation
201 };
202 ///@}
203
204 //----------------------------------------
205 /**@name Protected data */
206 ///@{
207 typename SLinSolver<R>::Status stat; ///< Status indicator.
208
209 int thedim; ///< dimension of factorized matrix
210 int nzCnt; ///< number of nonzeros in U
211 R initMaxabs; ///< maximum abs number in initail Matrix
212 R maxabs; ///< maximum abs number in L and U
213
214 R rowMemMult; ///< factor of minimum Memory * number of nonzeros
215 R colMemMult; ///< factor of minimum Memory * number of nonzeros
216 R lMemMult; ///< factor of minimum Memory * number of nonzeros
217
218 Perm row; ///< row permutation matrices
219 Perm col; ///< column permutation matrices
220
221 L l; ///< L matrix
222 std::vector<R> diag; ///< Array of pivot elements
223 U u; ///< U matrix
224
225 R* work; ///< Working array: must always be left as 0!
226
227 Timer* factorTime; ///< Time spent in factorizations
228 int factorCount; ///< Number of factorizations
229 int hugeValues; ///< number of times huge values occurred during solve (only used in debug mode)
230 std::shared_ptr<Tolerances> _tolerances; ///< Tolerances for the factorization
231 ///@}
232
233 private:
234
235 //----------------------------------------
236 /**@name Private data */
237 ///@{
238 Temp temp; ///< Temporary storage
239 ///@}
240
241 //----------------------------------------
242 /**@name Solving
243 These helper methods are used during the factorization process.
244 The solve*-methods solve lower and upper triangular systems from
245 the left or from the right, respectively The methods with '2' in
246 the end solve two systems at the same time. The methods with
247 "Eps" in the end consider elements smaller then the passed epsilon
248 as zero.
249 */
250 ///@{
251 ///
252 void solveUright(R* wrk, R* vec) const;
253 ///
254 int solveUrightEps(R* vec, int* nonz, R eps, R* rhs);
255 ///
256 void solveUright2(R* work1, R* vec1, R* work2, R* vec2);
257 ///
258 int solveUright2eps(R* work1, R* vec1, R* work2, R* vec2, int* nonz, R eps);
259 ///
260 void solveLright2(R* vec1, R* vec2);
261 ///
262 void solveUpdateRight(R* vec);
263 ///
264 void solveUpdateRight2(R* vec1, R* vec2);
265 ///
266 void solveUleft(R* work, R* vec);
267 ///
268 void solveUleft2(R* work1, R* vec1, R* work2, R* vec2);
269 ///
270 int solveLleft2forest(R* vec1, int* /* nonz */, R* vec2, R /* eps */);
271 ///
272 void solveLleft2(R* vec1, int* /* nonz */, R* vec2, R /* eps */);
273 ///
274 int solveLleftForest(R* vec, int* /* nonz */, R /* eps */);
275 ///
276 void solveLleft(R* vec) const;
277 ///
278 int solveLleftEps(R* vec, int* nonz, R eps);
279 ///
280 void solveUpdateLeft(R* vec);
281 ///
282 void solveUpdateLeft2(R* vec1, R* vec2);
283
284 void inline updateSolutionVectorLright(R change, int j, R& vec, int* idx, int& nnz);
285 ///
286 void vSolveLright(R* vec, int* ridx, int& rn, R eps);
287 ///
288 void vSolveLright2(R* vec, int* ridx, int& rn, R eps,
289 R* vec2, int* ridx2, int& rn2, R eps2);
290 ///
291 void vSolveLright3(R* vec, int* ridx, int& rn, R eps,
292 R* vec2, int* ridx2, int& rn2, R eps2,
293 R* vec3, int* ridx3, int& rn3, R eps3);
294 ///
295 int vSolveUright(R* vec, int* vidx, R* rhs, int* ridx, int rn, R eps);
296 ///
297 void vSolveUrightNoNZ(R* vec, R* rhs, int* ridx, int rn, R eps);
298 ///
299 int vSolveUright2(R* vec, int* vidx, R* rhs, int* ridx, int rn, R eps,
300 R* vec2, R* rhs2, int* ridx2, int rn2, R eps2);
301 ///
302 int vSolveUpdateRight(R* vec, int* ridx, int n, R eps);
303 ///
304 void vSolveUpdateRightNoNZ(R* vec, R /*eps*/);
305 ///
306 int solveUleft(R eps, R* vec, int* vecidx, R* rhs, int* rhsidx, int rhsn);
307 ///
308 void solveUleftNoNZ(R eps, R* vec, R* rhs, int* rhsidx, int rhsn);
309 ///
310 int solveLleftForest(R eps, R* vec, int* nonz, int n);
311 ///
312 void solveLleftForestNoNZ(R* vec);
313 ///
314 int solveLleft(R eps, R* vec, int* nonz, int rn);
315 ///
316 void solveLleftNoNZ(R* vec);
317 ///
318 int solveUpdateLeft(R eps, R* vec, int* nonz, int n);
319
320 ///
321 void forestPackColumns();
322 ///
323 void forestMinColMem(int size);
324 ///
325 void forestReMaxCol(int col, int len);
326
327 ///
328 void initPerm();
329 ///
330 void initFactorMatrix(const SVectorBase<R>** vec, const R eps);
331 ///
332 void minLMem(int size);
333 ///
334 void setPivot(const int p_stage, const int p_col, const int p_row, const R val);
335 ///
336 void colSingletons();
337 ///
338 void rowSingletons();
339
340 ///
341 void initFactorRings();
342 ///
343 void freeFactorRings();
344
345 ///
346 int setupColVals();
347 ///
348 void setupRowVals();
349
350 ///
351 void eliminateRowSingletons();
352 ///
353 void eliminateColSingletons();
354 ///
355 void selectPivots(R threshold);
356 ///
357 int updateRow(int r, int lv, int prow, int pcol, R pval, R eps);
358
359 ///
360 void eliminatePivot(int prow, int pos, R eps);
361 ///
362 void eliminateNucleus(const R eps, const R threshold);
363 ///
364 void minRowMem(int size);
365 ///
366 void minColMem(int size);
367 ///
368 void remaxCol(int p_col, int len);
369 ///
370 void packRows();
371 ///
372 void packColumns();
373 ///
374 void remaxRow(int p_row, int len);
375 ///
376 int makeLvec(int p_len, int p_row);
377 ///@}
378
379 protected:
380
381 //----------------------------------------
382 /**@name Solver methods */
383 ///@{
384 ///
385 void solveLright(R* vec);
386 ///
387 int solveRight4update(R* vec, int* nonz, R eps, R* rhs,
388 R* forest, int* forestNum, int* forestIdx);
389 ///
390 void solveRight(R* vec, R* rhs);
391 ///
392 int solveRight2update(R* vec1, R* vec2, R* rhs1,
393 R* rhs2, int* nonz, R eps, R* forest, int* forestNum, int* forestIdx);
394 ///
395 void solveRight2(R* vec1, R* vec2, R* rhs1, R* rhs2);
396 ///
397 void solveLeft(R* vec, R* rhs);
398 ///
399 int solveLeftEps(R* vec, R* rhs, int* nonz, R eps);
400 ///
401 int solveLeft2(R* vec1, int* nonz, R* vec2, R eps, R* rhs1, R* rhs2);
402
403 ///
404 int vSolveRight4update(R eps,
405 R* vec, int* idx, /* result */
406 R* rhs, int* ridx, int rn, /* rhs & Forest */
407 R* forest, int* forestNum, int* forestIdx);
408 ///
409 int vSolveRight4update2(R eps,
410 R* vec, int* idx, /* result1 */
411 R* rhs, int* ridx, int rn, /* rhs1 */
412 R* vec2, R eps2, /* result2 */
413 R* rhs2, int* ridx2, int rn2, /* rhs2 */
414 R* forest, int* forestNum, int* forestIdx);
415 /// sparse version of above method
416 void vSolveRight4update2sparse(
417 R eps, R* vec, int* idx, /* result1 */
418 R* rhs, int* ridx, int& rn, /* rhs1 */
419 R eps2, R* vec2, int* idx2, /* result2 */
420 R* rhs2, int* ridx2, int& rn2, /* rhs2 */
421 R* forest, int* forestNum, int* forestIdx);
422 ///
423 int vSolveRight4update3(R eps,
424 R* vec, int* idx, /* result1 */
425 R* rhs, int* ridx, int rn, /* rhs1 */
426 R* vec2, R eps2, /* result2 */
427 R* rhs2, int* ridx2, int rn2, /* rhs2 */
428 R* vec3, R eps3, /* result3 */
429 R* rhs3, int* ridx3, int rn3, /* rhs3 */
430 R* forest, int* forestNum, int* forestIdx);
431 /// sparse version of above method
432 void vSolveRight4update3sparse(
433 R eps, R* vec, int* idx, /* result1 */
434 R* rhs, int* ridx, int& rn, /* rhs1 */
435 R eps2, R* vec2, int* idx2, /* result2 */
436 R* rhs2, int* ridx2, int& rn2, /* rhs2 */
437 R eps3, R* vec3, int* idx3, /* result3 */
438 R* rhs3, int* ridx3, int& rn3, /* rhs3 */
439 R* forest, int* forestNum, int* forestIdx);
440 ///
441 void vSolveRightNoNZ(R* vec, R eps, /* result */
442 R* rhs, int* ridx, int rn); /* rhs */
443 ///
444 int vSolveLeft(R eps,
445 R* vec, int* idx, /* result */
446 R* rhs, int* ridx, int rn); /* rhs */
447 ///
448 void vSolveLeftNoNZ(R eps,
449 R* vec, /* result */
450 R* rhs, int* ridx, int rn); /* rhs */
451 ///
452 int vSolveLeft2(R eps,
453 R* vec, int* idx, /* result */
454 R* rhs, int* ridx, int rn, /* rhs */
455 R* vec2, /* result2 */
456 R* rhs2, int* ridx2, int rn2); /* rhs2 */
457 /// sparse version of solving 2 systems of equations
458 void vSolveLeft2sparse(R eps,
459 R* vec, int* idx, /* result */
460 R* rhs, int* ridx, int& rn, /* rhs */
461 R* vec2, int* idx2, /* result2 */
462 R* rhs2, int* ridx2, int& rn2); /* rhs2 */
463 ///
464 int vSolveLeft3(R eps,
465 R* vec, int* idx, /* result */
466 R* rhs, int* ridx, int rn, /* rhs */
467 R* vec2, /* result2 */
468 R* rhs2, int* ridx2, int rn2, /* rhs2 */
469 R* vec3, /* result3 */
470 R* rhs3, int* ridx3, int rn3); /* rhs3 */
471 /// sparse version of solving 3 systems of equations
472 void vSolveLeft3sparse(R eps,
473 R* vec, int* idx, /* result */
474 R* rhs, int* ridx, int& rn, /* rhs */
475 R* vec2, int* idx2, /* result2 */
476 R* rhs2, int* ridx2, int& rn2, /* rhs2 */
477 R* vec3, int* idx3, /* result2 */
478 R* rhs3, int* ridx3, int& rn3); /* rhs2 */
479
480 void forestUpdate(int col, R* work, int num, int* nonz);
481
482 void update(int p_col, R* p_work, const int* p_idx, int num);
483 void updateNoClear(int p_col, const R* p_work, const int* p_idx, int num);
484
485 ///
486 void factor(const SVectorBase<R>** vec, ///< Array of column VectorBase<R> pointers
487 R threshold, ///< pivoting threshold
488 R eps); ///< epsilon for zero detection
489 ///@}
490
491 //----------------------------------------
492 /**@name Debugging */
493 ///@{
494 ///
495 void dump() const;
496
497 ///
498 bool isConsistent() const;
499 ///@}
500 };
501
502 } // namespace soplex
503
504 // For general templated functions
505 #include "clufactor.hpp"
506
507 #endif // _CLUFACTOR_H_
508