Scippy

SoPlex

Sequential object-oriented simPlex

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