Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ldl.hpp
1 /* ========================================================================== */
2 /* === ldl.c: sparse LDL' factorization and solve package =================== */
3 /* ========================================================================== */
4 
5 /* LDL: a simple set of routines for sparse LDL' factorization. These routines
6  * are not terrifically fast (they do not use dense matrix kernels), but the
7  * code is very short. The purpose is to illustrate the algorithms in a very
8  * concise manner, primarily for educational purposes. Although the code is
9  * very concise, this package is slightly faster than the built-in sparse
10  * Cholesky factorization in MATLAB 7.0 (chol), when using the same input
11  * permutation.
12  *
13  * The routines compute the LDL' factorization of a real sparse symmetric
14  * matrix A (or PAP' if a permutation P is supplied), and solve upper
15  * and lower triangular systems with the resulting L and D factors. If A is
16  * positive definite then the factorization will be accurate. A can be
17  * indefinite (with negative values on the diagonal D), but in this case no
18  * guarantee of accuracy is provided, since no numeric pivoting is performed.
19  *
20  * The n-by-n sparse matrix A is in compressed-column form. The nonzero values
21  * in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with corresponding row
22  * indices in Ai [Ap [j] ... Ap [j+1]-1]. Ap [0] = 0 is required, and thus
23  * nz = Ap [n] is the number of nonzeros in A. Ap is an int array of size n+1.
24  * The int array Ai and the double array Ax are of size nz. This data structure
25  * is identical to the one used by MATLAB, except for the following
26  * generalizations. The row indices in each column of A need not be in any
27  * particular order, although they must be in the range 0 to n-1. Duplicate
28  * entries can be present; any duplicates are summed. That is, if row index i
29  * appears twice in a column j, then the value of A (i,j) is the sum of the two
30  * entries. The data structure used here for the input matrix A is more
31  * flexible than MATLAB's, which requires sorted columns with no duplicate
32  * entries.
33  *
34  * Only the diagonal and upper triangular part of A (or PAP' if a permutation
35  * P is provided) is accessed. The lower triangular parts of the matrix A or
36  * PAP' can be present, but they are ignored.
37  *
38  * The optional input permutation is provided as an array P of length n. If
39  * P [k] = j, the row and column j of A is the kth row and column of PAP'.
40  * If P is present then the factorization is LDL' = PAP' or L*D*L' = A(P,P) in
41  * 0-based MATLAB notation. If P is not present (a null pointer) then no
42  * permutation is performed, and the factorization is LDL' = A.
43  *
44  * The lower triangular matrix L is stored in the same compressed-column
45  * form (an int Lp array of size n+1, an int Li array of size Lp [n], and a
46  * double array Lx of the same size as Li). It has a unit diagonal, which is
47  * not stored. The row indices in each column of L are always returned in
48  * ascending order, with no duplicate entries. This format is compatible with
49  * MATLAB, except that it would be more convenient for MATLAB to include the
50  * unit diagonal of L. Doing so here would add additional complexity to the
51  * code, and is thus omitted in the interest of keeping this code short and
52  * readable.
53  *
54  * The elimination tree is held in the Parent [0..n-1] array. It is normally
55  * not required by the user, but it is required by ldl::numeric. The diagonal
56  * matrix D is held as an array D [0..n-1] of size n.
57  *
58  * --------------------
59  * C-callable routines:
60  * --------------------
61  *
62  * ldl::symbolic: Given the pattern of A, computes the Lp and Parent arrays
63  * required by ldl::numeric. Takes time proportional to the number of
64  * nonzeros in L. Computes the inverse Pinv of P if P is provided.
65  * Also returns Lnz, the count of nonzeros in each column of L below
66  * the diagonal (this is not required by ldl::numeric).
67  * ldl::numeric: Given the pattern and numerical values of A, the Lp array,
68  * the Parent array, and P and Pinv if applicable, computes the
69  * pattern and numerical values of L and D.
70  * ldl::lsolve: Solves Lx=b for a dense vector b.
71  * ldl::dsolve: Solves Dx=b for a dense vector b.
72  * ldl::ltsolve: Solves L'x=b for a dense vector b.
73  * ldl::perm: Computes x=Pb for a dense vector b.
74  * ldl::permt: Computes x=P'b for a dense vector b.
75  * ldl::valid_perm: checks the validity of a permutation vector
76  * ldl::valid_matrix: checks the validity of the sparse matrix A
77  *
78  * ----------------------------
79  * Limitations of this package:
80  * ----------------------------
81  *
82  * In the interest of keeping this code simple and readable, ldl::symbolic and
83  * ldl::numeric assume their inputs are valid. You can check your own inputs
84  * prior to calling these routines with the ldl::valid_perm and ldl::valid_matrix
85  * routines. Except for the two ldl::valid_* routines, no routine checks to see
86  * if the array arguments are present (non-NULL). Like all C routines, no
87  * routine can determine if the arrays are long enough and don't overlap.
88  *
89  * The ldl::numeric does check the numerical factorization, however. It returns
90  * n if the factorization is successful. If D (k,k) is zero, then k is
91  * returned, and L is only partially computed.
92  *
93  * No pivoting to control fill-in is performed, which is often critical for
94  * obtaining good performance. I recommend that you compute the permutation P
95  * using AMD or SYMAMD (approximate minimum degree ordering routines), or an
96  * appropriate graph-partitioning based ordering. See the ldldemo.m routine for
97  * an example in MATLAB, and the ldlmain.c stand-alone C program for examples of
98  * how to find P. Routines for manipulating compressed-column matrices are
99  * available in UMFPACK. AMD, SYMAMD, UMFPACK, and this LDL package are all
100  * available at http://www.cise.ufl.edu/research/sparse.
101  *
102  * -------------------------
103  * Possible simplifications:
104  * -------------------------
105  *
106  * These routines could be made even simpler with a few additional assumptions.
107  * If no input permutation were performed, the caller would have to permute the
108  * matrix first, but the computation of Pinv, and the use of P and Pinv could be
109  * removed. If only the diagonal and upper triangular part of A or PAP' are
110  * present, then the tests in the "if (i < k)" statement in ldl::symbolic and
111  * "if (i <= k)" in ldl::numeric, are always true, and could be removed (i can
112  * equal k in ldl::symbolic, but then the body of the if statement would
113  * correctly do no work since Flag [k] == k). If we could assume that no
114  * duplicate entries are present, then the statement Y [i] += Ax [p] could be
115  * replaced with Y [i] = Ax [p] in ldl::numeric.
116  *
117  * --------------------------
118  * Description of the method:
119  * --------------------------
120  *
121  * LDL computes the symbolic factorization by finding the pattern of L one row
122  * at a time. It does this based on the following theory. Consider a sparse
123  * system Lx=b, where L, x, and b, are all sparse, and where L comes from a
124  * Cholesky (or LDL') factorization. The elimination tree (etree) of L is
125  * defined as follows. The parent of node j is the smallest k > j such that
126  * L (k,j) is nonzero. Node j has no parent if column j of L is completely zero
127  * below the diagonal (j is a root of the etree in this case). The nonzero
128  * pattern of x is the union of the paths from each node i to the root, for
129  * each nonzero b (i). To compute the numerical solution to Lx=b, we can
130  * traverse the columns of L corresponding to nonzero values of x. This
131  * traversal does not need to be done in the order 0 to n-1. It can be done in
132  * any "topological" order, such that x (i) is computed before x (j) if i is a
133  * descendant of j in the elimination tree.
134  *
135  * The row-form of the LDL' factorization is shown in the MATLAB function
136  * ldlrow.m in this LDL package. Note that row k of L is found via a sparse
137  * triangular solve of L (1:k-1, 1:k-1) \ A (1:k-1, k), to use 1-based MATLAB
138  * notation. Thus, we can start with the nonzero pattern of the kth column of
139  * A (above the diagonal), follow the paths up to the root of the etree of the
140  * (k-1)-by-(k-1) leading submatrix of L, and obtain the pattern of the kth row
141  * of L. Note that we only need the leading (k-1)-by-(k-1) submatrix of L to
142  * do this. The elimination tree can be constructed as we go.
143  *
144  * The symbolic factorization does the same thing, except that it discards the
145  * pattern of L as it is computed. It simply counts the number of nonzeros in
146  * each column of L and then constructs the Lp index array when it's done. The
147  * symbolic factorization does not need to do this in topological order.
148  * Compare ldl::symbolic with the first part of ldl::numeric, and note that the
149  * while (len > 0) loop is not present in ldl::symbolic.
150  *
151  * LDL Version 1.1 (Apr. 22, 2005), Copyright (c) 2005 by Timothy A Davis,
152  * University of Florida. All Rights Reserved. Developed while on sabbatical
153  * at Stanford University and Lawrence Berkeley National Laboratory. Refer to
154  * the README file for the License. Available at
155  * http://www.cise.ufl.edu/research/sparse.
156  */
157 #include <boost/numeric/ublas/matrix_sparse.hpp>
158 
159 #include <feel/feelcore/feel.hpp>
160 #include <feel/feelcore/traits.hpp>
161 
162 namespace Feel
163 {
164 namespace glas
165 {
166 namespace ublas= boost::numeric::ublas;
167 
168 template<typename T = double>
169 class ldl
170 {
171 public:
172  typedef T value_type;
173  typedef ublas::compressed_matrix<value_type, ublas::column_major> matrix_type;
174 
175  ldl( matrix_type& A )
176  :
177  n( A.size1() ),
178  Ap( A.index1_data().begin() ),
179  Ai( A.index2_data().begin() ),
180  Ax( A.value_data().begin() ),
181  Lp( n+1 ),
182  Parent( n ),
183  Lnz( n ),
184  Li(),
185  Lx(),
186  D( n ),
187  Y( n ),
188  Flag( n ),
189  Pattern( n ),
190  P( 0 ),
191  Pinv( 0 )
192  {}
193  void symbolic () ;
194 
195  size_type numeric();
196 
197  void lsolve ( ublas::vector<value_type>& X );
198 
199  void dsolve ( ublas::vector<value_type>& X );
200 
201  void ltsolve ( ublas::vector<value_type>& X );
202 
206  void solve ( ublas::vector<value_type>& X )
207  {
208  lsolve( X );
209  dsolve( X );
210  ltsolve( X );
211  }
212 
213  void perm ( ublas::vector<value_type>& X, ublas::vector<value_type> const& B );
214  void permt ( ublas::vector<value_type>& X, ublas::vector<value_type> const& B );
215 
216  bool valid_perm ();
217  bool valid_matrix ();
218 
219 private:
220 
221  size_type n;
222 
223  size_type* Ap;
224  size_type* Ai;
225  value_type* Ax;
226 
227  ublas::vector<size_type> Lp;
228 
229  ublas::vector<int64_type> Parent;
230 
231  /* output of size n, not defined on input */
232  ublas::vector<size_type> Lnz;
233  ublas::vector<size_type> Li; /* output of size lnz=Lp[n], not defined on input */
234  ublas::vector<value_type> Lx; /* output of size lnz=Lp[n], not defined on input */
235  ublas::vector<value_type> D; /* output of size n, not defined on input */
236  ublas::vector<value_type> Y; /* workspace of size n, not defn. on input or output */
237 
238  /* workspace of size n, not defn. on input or output */
239  ublas::vector<size_type> Flag;
240  ublas::vector<size_type> Pattern;
241 
242  /* optional input of size n */
243  size_type* P;
244  /* optional output of size n (used if P is not NULL) */
245  size_type* Pinv;
246 };
247 
248 
249 /* ========================================================================== */
250 /* === ldl::symbolic ========================================================= */
251 /* ========================================================================== */
252 
253 /* The input to this routine is a sparse matrix A, stored in column form, and
254  * an optional permutation P. The output is the elimination tree
255  * and the number of nonzeros in each column of L. Parent [i] = k if k is the
256  * parent of i in the tree. The Parent array is required by ldl::numeric.
257  * Lnz [k] gives the number of nonzeros in the kth column of L, excluding the
258  * diagonal.
259  *
260  * One workspace vector (Flag) of size n is required.
261  *
262  * If P is NULL, then it is ignored. The factorization will be LDL' = A.
263  * Pinv is not computed. In this case, neither P nor Pinv are required by
264  * ldl::numeric.
265  *
266  * If P is not NULL, then it is assumed to be a valid permutation. If
267  * row and column j of A is the kth pivot, the P [k] = j. The factorization
268  * will be LDL' = PAP', or A (p,p) in MATLAB notation. The inverse permutation
269  * Pinv is computed, where Pinv [j] = k if P [k] = j. In this case, both P
270  * and Pinv are required as inputs to ldl::numeric.
271  *
272  * The floating-point operation count of the subsequent call to ldl::numeric
273  * is not returned, but could be computed after ldl::symbolic is done. It is
274  * the sum of (Lnz [k]) * (Lnz [k] + 2) for k = 0 to n-1.
275  */
276 template<typename T>
277 void ldl<T>::symbolic()
278 {
279  size_type i, k, p, kk, p2 ;
280 
281  if ( P )
282  {
283  /* If P is present then compute Pinv, the inverse of P */
284  for ( k = 0 ; k < n ; k++ )
285  {
286  Pinv [P [k]] = k ;
287  }
288  }
289 
290  for ( k = 0 ; k < n ; k++ )
291  {
292  /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
293  Parent [k] = -1 ; /* parent of k is not yet known */
294  Flag [k] = k ; /* mark node k as visited */
295  Lnz [k] = 0 ; /* count of nonzeros in column k of L */
296  kk = ( P ) ? ( P [k] ) : ( k ) ; /* kth original, or permuted, column */
297  p2 = Ap [kk+1] ;
298 
299  for ( p = Ap [kk] ; p < p2 ; p++ )
300  {
301  /* A (i,k) is nonzero (original or permuted A) */
302  i = ( Pinv ) ? ( Pinv [Ai [p]] ) : ( Ai [p] ) ;
303 
304  if ( i < k )
305  {
306  /* follow path from i to root of etree, stop at flagged node */
307  for ( ; Flag [i] != k ; i = Parent [i] )
308  {
309  /* find parent of i if not yet determined */
310  if ( Parent [i] == -1 ) Parent [i] = k ;
311 
312  Lnz [i]++ ; /* L (k,i) is nonzero */
313  Flag [i] = k ; /* mark i as visited */
314  }
315  }
316  }
317  }
318 
319  /* construct Lp index array from Lnz column counts */
320  Lp [0] = 0 ;
321 
322  for ( k = 0 ; k < n ; k++ )
323  {
324  Lp [k+1] = Lp [k] + Lnz [k] ;
325  }
326 
327  Li.resize( Lp[n] );
328  Lx.resize( Lp[n] );
329 }
330 
331 
332 /* ========================================================================== */
333 /* === ldl::numeric ========================================================== */
334 /* ========================================================================== */
335 
336 /* Given a sparse matrix A (the arguments n, Ap, Ai, and Ax) and its symbolic
337  * analysis (Lp and Parent, and optionally P and Pinv), compute the numeric LDL'
338  * factorization of A or PAP'. The outputs of this routine are arguments Li,
339  * Lx, and D. It also requires three size-n workspaces (Y, Pattern, and Flag).
340  */
341 
342 /* returns n if successful, k if D (k,k) is zero */
343 template<typename T>
344 size_type ldl<T>::numeric()
345 {
346  value_type yi, l_ki ;
347  size_type i, k, p, kk, p2;
348  size_type len, top ;
349 
350  for ( k = 0 ; k < n ; k++ )
351  {
352  /* compute nonzero Pattern of kth row of L, in topological order */
353  Y [k] = 0.0 ; /* Y(0:k) is now all zero */
354  top = n ; /* stack for pattern is empty */
355  Flag [k] = k ; /* mark node k as visited */
356  Lnz [k] = 0 ; /* count of nonzeros in column k of L */
357  kk = ( P ) ? ( P [k] ) : ( k ) ; /* kth original, or permuted, column */
358  p2 = Ap [kk+1] ;
359 
360  for ( p = Ap [kk] ; p < p2 ; p++ )
361  {
362  i = ( Pinv ) ? ( Pinv [Ai [p]] ) : ( Ai [p] ) ; /* get A(i,k) */
363 
364  if ( i <= k )
365  {
366  Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */
367 
368  for ( len = 0 ; Flag [i] != k ; i = Parent [i] )
369  {
370  Pattern [len++] = i ; /* L(k,i) is nonzero */
371  Flag [i] = k ; /* mark i as visited */
372  }
373 
374  while ( len > 0 ) Pattern [--top] = Pattern [--len] ;
375  }
376  }
377 
378  /* compute numerical values kth row of L (a sparse triangular solve) */
379  D [k] = Y [k] ; /* get D(k,k) and clear Y(k) */
380  Y [k] = 0.0 ;
381 
382  for ( ; top < n ; top++ )
383  {
384  i = Pattern [top] ; /* Pattern [top:n-1] is pattern of L(:,k) */
385  yi = Y [i] ; /* get and clear Y(i) */
386  Y [i] = 0.0 ;
387  p2 = Lp [i] + Lnz [i] ;
388 
389  for ( p = Lp [i] ; p < p2 ; p++ )
390  {
391  Y [Li [p]] -= Lx [p] * yi ;
392  }
393 
394  l_ki = yi / D [i] ; /* the nonzero entry L(k,i) */
395  D [k] -= l_ki * yi ;
396  Li [p] = k ; /* store L(k,i) in column form of L */
397  Lx [p] = l_ki ;
398  Lnz [i]++ ; /* increment count of nonzeros in col i */
399  }
400 
401  if ( D [k] == 0.0 ) return ( k ) ; /* failure, D(k,k) is zero */
402  }
403 
404  return ( n ) ; /* success, diagonal of D is all nonzero */
405 }
406 
407 
408 /* ========================================================================== */
409 /* === ldl::lsolve: solve Lx=b ============================================== */
410 /* ========================================================================== */
411 template<typename T>
412 void ldl<T>::lsolve( ublas::vector<value_type>& X )
413 {
414  size_type j, p, p2 ;
415 
416  for ( j = 0 ; j < n ; j++ )
417  {
418  p2 = Lp [j+1] ;
419 
420  for ( p = Lp [j] ; p < p2 ; p++ )
421  {
422  X [Li [p]] -= Lx [p] * X [j] ;
423  }
424  }
425 }
426 
427 
428 /* ========================================================================== */
429 /* === ldl::dsolve: solve Dx=b ============================================== */
430 /* ========================================================================== */
431 template<typename T>
432 void ldl<T>::dsolve( ublas::vector<value_type>& X )
433 {
434  size_type j ;
435 
436  for ( j = 0 ; j < n ; j++ )
437  {
438  X [j] /= D [j] ;
439  }
440 }
441 
442 
443 /* ========================================================================== */
444 /* === ldl::ltsolve: solve L'x=b ============================================ */
445 /* ========================================================================== */
446 template<typename T>
447 void ldl<T>::ltsolve( ublas::vector<value_type>& X )
448 {
449  size_type p, p2 ;
450 
451  for ( int j = n-1 ; j >= 0 ; j-- )
452  {
453  p2 = Lp [j+1] ;
454 
455  for ( p = Lp [j] ; p < p2 ; p++ )
456  {
457  X [j] -= Lx [p] * X [Li [p]] ;
458  }
459  }
460 }
461 
462 
463 /* ========================================================================== */
464 /* === ldl::perm: permute a vector, x=Pb ===================================== */
465 /* ========================================================================== */
466 template<typename T>
467 void ldl<T>::perm( ublas::vector<value_type> & X,
468  ublas::vector<value_type> const& B )
469 {
470  size_type j ;
471 
472  for ( j = 0 ; j < n ; j++ )
473  {
474  X [j] = B [P [j]] ;
475  }
476 }
477 
478 
479 /* ========================================================================== */
480 /* === ldl::permt: permute a vector, x=P'b =================================== */
481 /* ========================================================================== */
482 template<typename T>
483 void ldl<T>::permt( ublas::vector<value_type> & X,
484  ublas::vector<value_type> const& B )
485 {
486  size_type j ;
487 
488  for ( j = 0 ; j < n ; j++ )
489  {
490  X [P [j]] = B [j] ;
491  }
492 }
493 
494 
495 /* ========================================================================== */
496 /* === ldl::valid_perm: check if a permutation vector is valid =============== */
497 /* ========================================================================== */
498 
499 /* returns 1 if valid, 0 otherwise */
500 template<typename T>
501 bool ldl<T>::valid_perm()
502 {
503  size_type j, k ;
504 
505  if ( Flag.size() == 0 )
506  {
507  return false; /* n must be >= 0, and Flag must be present */
508  }
509 
510  if ( !P )
511  {
512  return true; /* If NULL, P is assumed to be the identity perm. */
513  }
514 
515  for ( j = 0 ; j < n ; j++ )
516  {
517  Flag [j] = 0 ; /* clear the Flag array */
518  }
519 
520  for ( k = 0 ; k < n ; k++ )
521  {
522  j = P [k] ;
523 
524  if ( j >= n || Flag [j] != 0 )
525  {
526  return false; /* P is not valid */
527  }
528 
529  Flag [j] = 1 ;
530  }
531 
532  return true; /* P is valid */
533 }
534 
535 
536 /* ========================================================================== */
537 /* === ldl::valid_matrix: check if a sparse matrix is valid ================== */
538 /* ========================================================================== */
539 
540 /* This routine checks to see if a sparse matrix A is valid for input to
541  * ldl::symbolic and ldl::numeric. It returns 1 if the matrix is valid, 0
542  * otherwise. A is in sparse column form. The numerical values in column j
543  * are stored in Ax [Ap [j] ... Ap [j+1]-1], with row indices in
544  * Ai [Ap [j] ... Ap [j+1]-1]. The Ax array is not checked.
545  */
546 template<typename T>
547 bool ldl<T>::valid_matrix()
548 {
549  size_type j, p ;
550 
551  /*
552  if (n < 0 || !Ap || !Ai || Ap [0] != 0)
553  {
554  return false ; // n must be >= 0, and Ap and Ai must be present
555  }*/
556  for ( j = 0 ; j < n ; j++ )
557  {
558  if ( Ap [j] > Ap [j+1] )
559  {
560  return false; /* Ap must be monotonically nondecreasing */
561  }
562  }
563 
564  for ( p = 0 ; p < Ap [n] ; p++ )
565  {
566  if ( Ai [p] >= n )
567  {
568  return false; /* row indices must be in the range 0 to n-1 */
569  }
570  }
571 
572  return true; /* matrix is valid */
573 }
574 
575 } // namespace glas
576 } // namespace Feel
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319

Generated on Sun Dec 22 2013 13:11:06 for Feel++ by doxygen 1.8.5