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_rational.h
26 * @brief Implementation of sparse LU factorization with Rational precision.
27 */
28 #ifndef _CLUFACTOR_RATIONAL_H_
29 #define _CLUFACTOR_RATIONAL_H_
30
31 #include "soplex/spxdefines.h"
32 #include "soplex/slinsolver_rational.h"
33 #include "soplex/timer.h"
34 #include "soplex/svector.h"
35 #include "soplex/rational.h"
36 #include "soplex/basevectors.h"
37
38 #define SOPLEX_WITH_L_ROWS 1
39
40 namespace soplex
41 {
42 /**@brief Implementation of sparse LU factorization with Rational precision.
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 class CLUFactorRational
49 {
50 public:
51
52 //----------------------------------------
53 /**@name Public types */
54 ///@{
55 /** Doubly linked ring structure for garbage collection of column or
56 * row file in working matrix.
57 */
58 struct Dring
59 {
60 Dring* next;
61 Dring* prev;
62 int idx;
63 };
64
65 /// Pivot Ring
(1) Event rule_of_three_violation: |
Class "soplex::CLUFactorRational::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] |
66 class Pring
67 {
68 public:
69 Pring* next; ///<
70 Pring* prev; ///<
71 int idx; ///< index of pivot row
72 int pos; ///< position of pivot column in row
73 int mkwtz; ///< markowitz number of pivot
74
75 Pring() ///< constructor
76 : next(0)
77 , prev(0)
78 , idx(0)
79 , pos(0)
80 , mkwtz(0)
81 {}
82
83 private:
84 Pring(const Pring&); ///< blocked copy constructor
85 Pring& operator= (const Pring&); ///< blocked assignment operator
86 };
87 ///@}
88
89 protected:
90
91 //----------------------------------------
92 /**@name Protected types */
93 ///@{
94 /// Temporary data structures.
95 class Temp
96 {
97 public:
98 int* s_mark; ///< marker
99 VectorRational s_max; ///< maximum absolute value per row (or -1)
100 int* s_cact; ///< lengths of columns of active submatrix
101 int stage; ///< stage of the structure
102 Pring pivots; ///< ring of selected pivot rows
103 Pring* pivot_col; ///< column index handlers for Real linked list
104 Pring* pivot_colNZ; ///< lists for columns to number of nonzeros
105 Pring* pivot_row; ///< row index handlers for Real linked list
106 Pring* pivot_rowNZ; ///< lists for rows to number of nonzeros
107
108 Temp(); ///< constructor
109 ~Temp(); ///< destructor
110 void init(int p_dim); ///< initialization
111 void clear(); ///< clears the structure
112
113 private:
114 Temp(const Temp&); ///< blocked copy constructor
115 Temp& operator= (const Temp&); ///< blocked assignment operator
116 };
117
118 /// Data structures for saving the row and column permutations.
119 struct Perm
120 {
121 int* orig; ///< orig[p] original index from p
122 int* perm; ///< perm[i] permuted index from i
123 };
124
125 /// Data structures for saving the working matrix and U factor.
126 struct U
127 {
128 ///
129 struct Row
130 {
131 Dring list; /*!< \brief Double linked ringlist of vector
132 indices in the order they appear
133 in the row file */
134 Dring* elem; ///< %Array of ring elements.
135 int used; ///< used entries of arrays idx and val
136 VectorRational val; ///< hold nonzero values
137 int* idx; ///< array of length val.dim() to hold column indices of nonzeros in val
138 int* start; ///< starting positions in val and idx
139 int* len; ///< used nonzeros per row vectors
140 int* max; /*!< \brief maximum available nonzeros per row:
141 start[i] + max[i] == start[elem[i].next->idx]
142 len[i] <= max[i]. */
143 } row;
144
145 ///
146 struct Col
147 {
148 Dring list; /*!< \brief Double linked ringlist of vector
149 indices in the order they appear
150 in the column file */
151 Dring* elem; ///< %Array of ring elements.
152 int size; ///< size of array idx
153 int used; ///< used entries of array idx
154 int* idx; ///< hold row indices of nonzeros
155 VectorRational val; /*!< \brief hold nonzero values: this is only initialized
156 in the end of the factorization with DEFAULT
157 updates. */
158 int* start; ///< starting positions in val and idx
159 int* len; ///< used nonzeros per column vector
160 int* max; /*!< \brief maximum available nonzeros per colunn:
161 start[i] + max[i] == start[elem[i].next->idx]
162 len[i] <= max[i]. */
163 } col;
164 };
165
166
167 /// Data structures for saving the working matrix and L factor.
168 struct L
169 {
170 VectorRational val; ///< values of L vectors
171 int* idx; ///< array of size val.dim() storing indices of L vectors
172 int startSize; ///< size of array start
173 int firstUpdate; ///< number of first update L vector
174 int firstUnused; ///< number of first unused L vector
175 int* start; ///< starting positions in val and idx
176 int* row; ///< column indices of L vectors
177 int updateType; ///< type of updates to be used.
178
179 /* The following arrays have length |firstUpdate|, since they keep
180 * rows of the L-vectors occuring during the factorization (without
181 * updates), only:
182 */
183 VectorRational rval; ///< values of rows of L
184 int* ridx; ///< indices of rows of L
185 int* rbeg; ///< start of rows in rval and ridx
186 int* rorig; ///< original row permutation
187 int* rperm; ///< original row permutation
188 };
189 ///@}
190
191 //----------------------------------------
192 /**@name Protected data */
193 ///@{
194 SLinSolverRational::Status stat; ///< Status indicator.
195
196 int thedim; ///< dimension of factorized matrix
197 int nzCnt; ///< number of nonzeros in U
198 Rational initMaxabs; ///< maximum abs number in initail Matrix
199 Rational maxabs; ///< maximum abs number in L and U
200
201 Real rowMemMult; ///< factor of minimum Memory * number of nonzeros
202 Real colMemMult; ///< factor of minimum Memory * number of nonzeros
203 Real lMemMult; ///< factor of minimum Memory * number of nonzeros
204
205 Perm row; ///< row permutation matrices
206 Perm col; ///< column permutation matrices
207
208 L l; ///< L matrix
209 VectorRational diag; ///< Array of pivot elements
210 U u; ///< U matrix
211
212 Rational* work; ///< Working array: must always be left as 0!
213
214 Timer* factorTime; ///< Time spent in factorizations
215 int factorCount; ///< Number of factorizations
216 Real timeLimit; ///< Time limit on factorization or solves
217 ///@}
218
219 private:
220
221 //----------------------------------------
222 /**@name Private data */
223 ///@{
224 Temp temp; ///< Temporary storage
225 ///@}
226
227 //----------------------------------------
228 /**@name Solving
229 These helper methods are used during the factorization process.
230 The solve*-methods solve lower and upper triangular systems from
231 the left or from the right, respectively The methods with '2' in
232 the end solve two systems at the same time. The methods with
233 "Eps" in the end consider elements smaller then the passed epsilon
234 as zero.
235 */
236 ///@{
237 ///
238 void solveUright(Rational* wrk, Rational* vec);
239 ///
240 int solveUrightEps(Rational* vec, int* nonz, Rational* rhs);
241 ///
242 void solveUright2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2);
243 ///
244 int solveUright2eps(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2, int* nonz);
245 ///
246 void solveLright2(Rational* vec1, Rational* vec2);
247 ///
248 void solveUpdateRight(Rational* vec);
249 ///
250 void solveUpdateRight2(Rational* vec1, Rational* vec2);
251 ///
252 void solveUleft(Rational* work, Rational* vec);
253 ///
254 void solveUleft2(Rational* work1, Rational* vec1, Rational* work2, Rational* vec2);
255 ///
256 int solveLleft2forest(Rational* vec1, int* /* nonz */, Rational* vec2);
257 ///
258 void solveLleft2(Rational* vec1, int* /* nonz */, Rational* vec2);
259 ///
260 int solveLleftForest(Rational* vec, int* /* nonz */);
261 ///
262 void solveLleft(Rational* vec);
263 ///
264 int solveLleftEps(Rational* vec, int* nonz);
265 ///
266 void solveUpdateLeft(Rational* vec);
267 ///
268 void solveUpdateLeft2(Rational* vec1, Rational* vec2);
269
270 ///
271 int vSolveLright(Rational* vec, int* ridx, int rn);
272 ///
273 void vSolveLright2(Rational* vec, int* ridx, int* rnptr,
274 Rational* vec2, int* ridx2, int* rn2ptr);
275 ///
276 void vSolveLright3(Rational* vec, int* ridx, int* rnptr,
277 Rational* vec2, int* ridx2, int* rn2ptr,
278 Rational* vec3, int* ridx3, int* rn3ptr);
279 ///
280 int vSolveUright(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn);
281 ///
282 void vSolveUrightNoNZ(Rational* vec, Rational* rhs, int* ridx, int rn);
283 ///
284 int vSolveUright2(Rational* vec, int* vidx, Rational* rhs, int* ridx, int rn,
285 Rational* vec2, Rational* rhs2, int* ridx2, int rn2);
286 ///
287 int vSolveUpdateRight(Rational* vec, int* ridx, int n);
288 ///
289 void vSolveUpdateRightNoNZ(Rational* vec);
290 ///
291 int solveUleft(Rational* vec, int* vecidx, Rational* rhs, int* rhsidx, int rhsn);
292 ///
293 void solveUleftNoNZ(Rational* vec, Rational* rhs, int* rhsidx, int rhsn);
294 ///
295 int solveLleftForest(Rational* vec, int* nonz, int n);
296 ///
297 void solveLleftForestNoNZ(Rational* vec);
298 ///
299 int solveLleft(Rational* vec, int* nonz, int rn);
300 ///
301 void solveLleftNoNZ(Rational* vec);
302 ///
303 int solveUpdateLeft(Rational* vec, int* nonz, int n);
304
305 ///
306 void forestPackColumns();
307 ///
308 void forestMinColMem(int size);
309 ///
310 void forestReMaxCol(int col, int len);
311
312 ///
313 void initPerm();
314 ///
315 void initFactorMatrix(const SVectorRational** vec);
316 ///
317 void minLMem(int size);
318 ///
319 void setPivot(const int p_stage, const int p_col, const int p_row, const Rational& val);
320 ///
321 void colSingletons();
322 ///
323 void rowSingletons();
324
325 ///
326 void initFactorRings();
327 ///
328 void freeFactorRings();
329
330 ///
331 int setupColVals();
332 ///
333 void setupRowVals();
334
335 ///
336 void eliminateRowSingletons();
337 ///
338 void eliminateColSingletons();
339 ///
340 void selectPivots(const Rational& threshold);
341 ///
342 int updateRow(int r, int lv, int prow, int pcol, const Rational& pval);
343
344 ///
345 void eliminatePivot(int prow, int pos);
346 ///
347 void eliminateNucleus(const Rational& threshold);
348 ///
349 void minRowMem(int size);
350 ///
351 void minColMem(int size);
352 ///
353 void remaxCol(int p_col, int len);
354 ///
355 void packRows();
356 ///
357 void packColumns();
358 ///
359 void remaxRow(int p_row, int len);
360 ///
361 int makeLvec(int p_len, int p_row);
362 ///
363 bool timeLimitReached()
364 {
365 if(timeLimit >= 0.0 && factorTime->time() >= timeLimit)
366 {
367 stat = SLinSolverRational::TIME;
368 return true;
369 }
370 else
371 return false;
372 }
373 ///@}
374
375 protected:
376
377 //----------------------------------------
378 /**@name Solver methods */
379 ///@{
380 ///
381 void solveLright(Rational* vec);
382 ///
383 int solveRight4update(Rational* vec, int* nonz, Rational* rhs,
384 Rational* forest, int* forestNum, int* forestIdx);
385 ///
386 void solveRight(Rational* vec, Rational* rhs);
387 ///
388 int solveRight2update(Rational* vec1, Rational* vec2, Rational* rhs1,
389 Rational* rhs2, int* nonz, Rational* forest, int* forestNum, int* forestIdx);
390 ///
391 void solveRight2(Rational* vec1, Rational* vec2, Rational* rhs1, Rational* rhs2);
392 ///
393 void solveLeft(Rational* vec, Rational* rhs);
394 ///
395 int solveLeftEps(Rational* vec, Rational* rhs, int* nonz);
396 ///
397 int solveLeft2(Rational* vec1, int* nonz, Rational* vec2, Rational* rhs1, Rational* rhs2);
398
399 ///
400 int vSolveRight4update(
401 Rational* vec, int* idx, /* result */
402 Rational* rhs, int* ridx, int rn, /* rhs & Forest */
403 Rational* forest, int* forestNum, int* forestIdx);
404 ///
405 int vSolveRight4update2(
406 Rational* vec, int* idx, /* result1 */
407 Rational* rhs, int* ridx, int rn, /* rhs1 */
408 Rational* vec2, /* result2 */
409 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
410 Rational* forest, int* forestNum, int* forestIdx);
411 ///
412 int vSolveRight4update3(
413 Rational* vec, int* idx, /* result1 */
414 Rational* rhs, int* ridx, int rn, /* rhs1 */
415 Rational* vec2, /* result2 */
416 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
417 Rational* vec3, /* result3 */
418 Rational* rhs3, int* ridx3, int rn3, /* rhs3 */
419 Rational* forest, int* forestNum, int* forestIdx);
420 ///
421 void vSolveRightNoNZ(Rational* vec2, /* result2 */
422 Rational* rhs2, int* ridx2, int rn2); /* rhs2 */
423 ///
424 int vSolveLeft(
425 Rational* vec, int* idx, /* result */
426 Rational* rhs, int* ridx, int rn); /* rhs */
427 ///
428 void vSolveLeftNoNZ(
429 Rational* vec, /* result */
430 Rational* rhs, int* ridx, int rn); /* rhs */
431 ///
432 int vSolveLeft2(
433 Rational* vec, int* idx, /* result */
434 Rational* rhs, int* ridx, int rn, /* rhs */
435 Rational* vec2, /* result2 */
436 Rational* rhs2, int* ridx2, int rn2); /* rhs2 */
437 ///
438 int vSolveLeft3(Rational* vec, int* idx, /* result */
439 Rational* rhs, int* ridx, int rn, /* rhs */
440 Rational* vec2, /* result2 */
441 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */
442 Rational* vec3, /* result3 */
443 Rational* rhs3, int* ridx3, int rn3); /* rhs3 */
444
445 void forestUpdate(int col, Rational* work, int num, int* nonz);
446
447 void update(int p_col, Rational* p_work, const int* p_idx, int num);
448 void updateNoClear(int p_col, const Rational* p_work, const int* p_idx, int num);
449
450 ///
451 void factor(const SVectorRational** vec, ///< Array of column vector pointers
452 const Rational& threshold); ///< pivoting threshold
453 ///@}
454
455 //----------------------------------------
456 /**@name Debugging */
457 ///@{
458 ///
459 void dump() const;
460
461 ///
462 bool isConsistent() const;
463 ///@}
464 };
465
466 } // namespace soplex
467 #include "clufactor_rational.hpp"
468 #endif // _CLUFACTOR_RATIONAL_H_
469