Scippy

SoPlex

Sequential object-oriented simPlex

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