1 /* ========================================================================== */
2 /* === Cholesky/cholmod_rowcolcounts ======================================== */
3 /* ========================================================================== */
4
5 /* -----------------------------------------------------------------------------
6 * CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
7 * -------------------------------------------------------------------------- */
8
9 /* Compute the row and column counts of the Cholesky factor L of the matrix
10 * A or A*A'. The etree and its postordering must already be computed (see
11 * cholmod_etree and cholmod_postorder) and given as inputs to this routine.
12 *
13 * For the symmetric case (LL'=A), A is accessed by column. Only the lower
14 * triangular part of A is used. Entries not in this part of the matrix are
15 * ignored. This is the same as storing the upper triangular part of A by
16 * rows, with entries in the lower triangular part being ignored. NOTE: this
17 * representation is the TRANSPOSE of the input to cholmod_etree.
18 *
19 * For the unsymmetric case (LL'=AA'), A is accessed by column. Equivalently,
20 * if A is viewed as a matrix in compressed-row form, this routine computes
21 * the row and column counts for L where LL'=A'A. If the input vector f is
22 * present, then F*F' is analyzed instead, where F = A(:,f).
23 *
24 * The set f is held in fset and fsize.
25 * fset = NULL means ":" in MATLAB. fset is ignored.
26 * fset != NULL means f = fset [0..fset-1].
27 * fset != NULL and fsize = 0 means f is the empty set.
28 * Common->status is set to CHOLMOD_INVALID if fset is invalid.
29 *
30 * In both cases, the columns of A need not be sorted.
31 * A can be packed or unpacked.
32 *
33 * References:
34 * J. Gilbert, E. Ng, B. Peyton, "An efficient algorithm to compute row and
35 * column counts for sparse Cholesky factorization", SIAM J. Matrix Analysis &
36 * Applic., vol 15, 1994, pp. 1075-1091.
37 *
38 * J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
39 * sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
40 *
41 * workspace:
42 * if symmetric: Flag (nrow), Iwork (2*nrow)
43 * if unsymmetric: Flag (nrow), Iwork (2*nrow+ncol), Head (nrow+1)
44 *
45 * Supports any xtype (pattern, real, complex, or zomplex).
46 */
47
48 #ifndef NCHOLESKY
49
50 #include "cholmod_internal.h"
51 #include "cholmod_cholesky.h"
52
53 /* ========================================================================== */
54 /* === initialize_node ====================================================== */
55 /* ========================================================================== */
56
initialize_node(Int k,Int Post[],Int Parent[],Int ColCount[],Int PrevNbr[])57 static Int initialize_node /* initial work for kth node in postordered etree */
58 (
59 Int k, /* at the kth step of the algorithm (and kth node) */
60 Int Post [ ], /* Post [k] = i, the kth node in postordered etree */
61 Int Parent [ ], /* Parent [i] is the parent of i in the etree */
62 Int ColCount [ ], /* ColCount [c] is the current weight of node c */
63 Int PrevNbr [ ] /* PrevNbr [u] = k if u was last considered at step k */
64 )
65 {
66 Int p, parent ;
67 /* determine p, the kth node in the postordered etree */
68 p = Post [k] ;
69 /* adjust the weight if p is not a root of the etree */
70 parent = Parent [p] ;
71 if (parent != EMPTY)
72 {
73 ColCount [parent]-- ;
74 }
75 /* flag node p to exclude self edges (p,p) */
76 PrevNbr [p] = k ;
77 return (p) ;
78 }
79
80
81 /* ========================================================================== */
82 /* === process_edge ========================================================= */
83 /* ========================================================================== */
84
85 /* edge (p,u) is being processed. p < u is a descendant of its ancestor u in
86 * the etree. node p is the kth node in the postordered etree. */
87
process_edge(Int p,Int u,Int k,Int First[],Int PrevNbr[],Int ColCount[],Int PrevLeaf[],Int RowCount[],Int SetParent[],Int Level[])88 static void process_edge
89 (
90 Int p, /* process edge (p,u) of the matrix */
91 Int u,
92 Int k, /* we are at the kth node in the postordered etree */
93 Int First [ ], /* First [i] = k if the postordering of first
94 * descendent of node i is k */
95 Int PrevNbr [ ], /* u was last considered at step k = PrevNbr [u] */
96 Int ColCount [ ], /* ColCount [c] is the current weight of node c */
97 Int PrevLeaf [ ], /* s = PrevLeaf [u] means that s was the last leaf
98 * seen in the subtree rooted at u. */
99 Int RowCount [ ], /* RowCount [i] is # of nonzeros in row i of L,
100 * including the diagonal. Not computed if NULL. */
101 Int SetParent [ ], /* the FIND/UNION data structure, which forms a set
102 * of trees. A root i has i = SetParent [i]. Following
103 * a path from i to the root q of the subtree containing
104 * i means that q is the SetParent representative of i.
105 * All nodes in the tree could have their SetParent
106 * equal to the root q; the tree representation is used
107 * to save time. When a path is traced from i to its
108 * root q, the path is re-traversed to set the SetParent
109 * of the whole path to be the root q. */
110 Int Level [ ] /* Level [i] = length of path from node i to root */
111 )
112 {
113 Int prevleaf, q, s, sparent ;
114 if (First [p] > PrevNbr [u])
115 {
116 /* p is a leaf of the subtree of u */
117 ColCount [p]++ ;
118 prevleaf = PrevLeaf [u] ;
119 if (prevleaf == EMPTY)
120 {
121 /* p is the first leaf of subtree of u; RowCount will be incremented
122 * by the length of the path in the etree from p up to u. */
123 q = u ;
124 }
125 else
126 {
127 /* q = FIND (prevleaf): find the root q of the
128 * SetParent tree containing prevleaf */
129 for (q = prevleaf ; q != SetParent [q] ; q = SetParent [q])
130 {
131 ;
132 }
133 /* the root q has been found; re-traverse the path and
134 * perform path compression */
135 s = prevleaf ;
136 for (s = prevleaf ; s != q ; s = sparent)
137 {
138 sparent = SetParent [s] ;
139 SetParent [s] = q ;
140 }
141 /* adjust the RowCount and ColCount; RowCount will be incremented by
142 * the length of the path from p to the SetParent root q, and
143 * decrement the ColCount of q by one. */
144 ColCount [q]-- ;
145 }
146 if (RowCount != NULL)
147 {
148 /* if RowCount is being computed, increment it by the length of
149 * the path from p to q */
150 RowCount [u] += (Level [p] - Level [q]) ;
151 }
152 /* p is a leaf of the subtree of u, so mark PrevLeaf [u] to be p */
153 PrevLeaf [u] = p ;
154 }
155 /* flag u has having been processed at step k */
156 PrevNbr [u] = k ;
157 }
158
159
160 /* ========================================================================== */
161 /* === finalize_node ======================================================== */
162 /* ========================================================================== */
163
finalize_node(Int p,Int Parent[],Int SetParent[])164 static void finalize_node /* compute UNION (p, Parent [p]) */
165 (
166 Int p,
167 Int Parent [ ], /* Parent [p] is the parent of p in the etree */
168 Int SetParent [ ] /* see process_edge, above */
169 )
170 {
171 /* all nodes in the SetParent tree rooted at p now have as their final
172 * root the node Parent [p]. This computes UNION (p, Parent [p]) */
173 if (Parent [p] != EMPTY)
174 {
175 SetParent [p] = Parent [p] ;
176 }
177 }
178
179
180 /* ========================================================================== */
181 /* === cholmod_rowcolcounts ================================================= */
182 /* ========================================================================== */
183
CHOLMOD(rowcolcounts)184 int CHOLMOD(rowcolcounts)
185 (
186 /* ---- input ---- */
187 cholmod_sparse *A, /* matrix to analyze */
188 Int *fset, /* subset of 0:(A->ncol)-1 */
189 size_t fsize, /* size of fset */
190 Int *Parent, /* size nrow. Parent [i] = p if p is the parent of i */
191 Int *Post, /* size nrow. Post [k] = i if i is the kth node in
192 * the postordered etree. */
193 /* ---- output --- */
194 Int *RowCount, /* size nrow. RowCount [i] = # entries in the ith row of
195 * L, including the diagonal. */
196 Int *ColCount, /* size nrow. ColCount [i] = # entries in the ith
197 * column of L, including the diagonal. */
198 Int *First, /* size nrow. First [i] = k is the least postordering
199 * of any descendant of i. */
200 Int *Level, /* size nrow. Level [i] is the length of the path from
201 * i to the root, with Level [root] = 0. */
202 /* --------------- */
203 cholmod_common *Common
204 )
205 {
206 double fl, ff ;
207 Int *Ap, *Ai, *Anz, *PrevNbr, *SetParent, *Head, *PrevLeaf, *Anext, *Ipost,
208 *Iwork ;
209 Int i, j, r, k, len, s, p, pend, inew, stype, nf, anz, inode, parent,
210 nrow, ncol, packed, use_fset, jj ;
211 size_t w ;
212 int ok = TRUE ;
213
214 /* ---------------------------------------------------------------------- */
215 /* check inputs */
216 /* ---------------------------------------------------------------------- */
217
218 RETURN_IF_NULL_COMMON (FALSE) ;
219 RETURN_IF_NULL (A, FALSE) ;
220 RETURN_IF_NULL (Parent, FALSE) ;
221 RETURN_IF_NULL (Post, FALSE) ;
222 RETURN_IF_NULL (ColCount, FALSE) ;
223 RETURN_IF_NULL (First, FALSE) ;
224 RETURN_IF_NULL (Level, FALSE) ;
225 RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
226 stype = A->stype ;
227 if (stype > 0)
228 {
229 /* symmetric with upper triangular part not supported */
230 ERROR (CHOLMOD_INVALID, "symmetric upper not supported") ;
231 return (FALSE) ;
232 }
233 Common->status = CHOLMOD_OK ;
234
235 /* ---------------------------------------------------------------------- */
236 /* allocate workspace */
237 /* ---------------------------------------------------------------------- */
238
239 nrow = A->nrow ; /* the number of rows of A */
240 ncol = A->ncol ; /* the number of columns of A */
241
242 /* w = 2*nrow + (stype ? 0 : ncol) */
243 w = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
244 w = CHOLMOD(add_size_t) (w, (stype ? 0 : ncol), &ok) ;
245 if (!ok)
246 {
247 ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
248 return (FALSE) ;
249 }
250
251 CHOLMOD(allocate_work) (nrow, w, 0, Common) ;
252 if (Common->status < CHOLMOD_OK)
253 {
254 return (FALSE) ;
255 }
256
257 ASSERT (CHOLMOD(dump_perm) (Post, nrow, nrow, "Post", Common)) ;
258 ASSERT (CHOLMOD(dump_parent) (Parent, nrow, "Parent", Common)) ;
259
260 /* ---------------------------------------------------------------------- */
261 /* get inputs */
262 /* ---------------------------------------------------------------------- */
263
264 Ap = A->p ; /* size ncol+1, column pointers for A */
265 Ai = A->i ; /* the row indices of A, of size nz=Ap[ncol+1] */
266 Anz = A->nz ;
267 packed = A->packed ;
268 ASSERT (IMPLIES (!packed, Anz != NULL)) ;
269
270 /* ---------------------------------------------------------------------- */
271 /* get workspace */
272 /* ---------------------------------------------------------------------- */
273
274 Iwork = Common->Iwork ;
275 SetParent = Iwork ; /* size nrow (i/i/l) */
276 PrevNbr = Iwork + nrow ; /* size nrow (i/i/l) */
277 Anext = Iwork + 2*((size_t) nrow) ; /* size ncol (i/i/l) (unsym only) */
278 PrevLeaf = Common->Flag ; /* size nrow */
279 Head = Common->Head ; /* size nrow+1 (unsym only)*/
280
281 /* ---------------------------------------------------------------------- */
282 /* find the first descendant and level of each node in the tree */
283 /* ---------------------------------------------------------------------- */
284
285 /* First [i] = k if the postordering of first descendent of node i is k */
286 /* Level [i] = length of path from node i to the root (Level [root] = 0) */
287
288 for (i = 0 ; i < nrow ; i++)
289 {
290 First [i] = EMPTY ;
291 }
292
293 /* postorder traversal of the etree */
294 for (k = 0 ; k < nrow ; k++)
295 {
296 /* node i of the etree is the kth node in the postordered etree */
297 i = Post [k] ;
298
299 /* i is a leaf if First [i] is still EMPTY */
300 /* ColCount [i] starts at 1 if i is a leaf, zero otherwise */
301 ColCount [i] = (First [i] == EMPTY) ? 1 : 0 ;
302
303 /* traverse the path from node i to the root, stopping if we find a
304 * node r whose First [r] is already defined. */
305 len = 0 ;
306 for (r = i ; (r != EMPTY) && (First [r] == EMPTY) ; r = Parent [r])
307 {
308 First [r] = k ;
309 len++ ;
310 }
311 if (r == EMPTY)
312 {
313 /* we hit a root node, the level of which is zero */
314 len-- ;
315 }
316 else
317 {
318 /* we stopped at node r, where Level [r] is already defined */
319 len += Level [r] ;
320 }
321 /* re-traverse the path from node i to r; set the level of each node */
322 for (s = i ; s != r ; s = Parent [s])
323 {
324 Level [s] = len-- ;
325 }
326 }
327
328 /* ---------------------------------------------------------------------- */
329 /* AA' case: sort columns of A according to first postordered row index */
330 /* ---------------------------------------------------------------------- */
331
332 fl = 0.0 ;
333 if (stype == 0)
334 {
335 /* [ use PrevNbr [0..nrow-1] as workspace for Ipost */
336 Ipost = PrevNbr ;
337 /* Ipost [i] = k if i is the kth node in the postordered etree. */
338 for (k = 0 ; k < nrow ; k++)
339 {
340 Ipost [Post [k]] = k ;
341 }
342 use_fset = (fset != NULL) ;
343 if (use_fset)
344 {
345 nf = fsize ;
346 /* clear Anext to check fset */
347 for (j = 0 ; j < ncol ; j++)
348 {
349 Anext [j] = -2 ;
350 }
351 /* find the first postordered row in each column of A (post,f)
352 * and place the column in the corresponding link list */
353 for (jj = 0 ; jj < nf ; jj++)
354 {
355 j = fset [jj] ;
356 if (j < 0 || j > ncol || Anext [j] != -2)
357 {
358 /* out-of-range or duplicate entry in fset */
359 ERROR (CHOLMOD_INVALID, "fset invalid") ;
360 return (FALSE) ;
361 }
362 /* flag column j as having been seen */
363 Anext [j] = EMPTY ;
364 }
365 /* fset is now valid */
366 ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
367 }
368 else
369 {
370 nf = ncol ;
371 }
372 for (jj = 0 ; jj < nf ; jj++)
373 {
374 j = (use_fset) ? (fset [jj]) : jj ;
375 /* column j is in the fset; find the smallest row (if any) */
376 p = Ap [j] ;
377 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
378 ff = (double) MAX (0, pend - p) ;
379 fl += ff*ff + ff ;
380 if (pend > p)
381 {
382 k = Ipost [Ai [p]] ;
383 for ( ; p < pend ; p++)
384 {
385 inew = Ipost [Ai [p]] ;
386 k = MIN (k, inew) ;
387 }
388 /* place column j in link list k */
389 ASSERT (k >= 0 && k < nrow) ;
390 Anext [j] = Head [k] ;
391 Head [k] = j ;
392 }
393 }
394 /* Ipost no longer needed for inverse postordering ]
395 * Head [k] contains a link list of all columns whose first
396 * postordered row index is equal to k, for k = 0 to nrow-1. */
397 }
398
399 /* ---------------------------------------------------------------------- */
400 /* compute the row counts and node weights */
401 /* ---------------------------------------------------------------------- */
402
403 if (RowCount != NULL)
404 {
405 for (i = 0 ; i < nrow ; i++)
406 {
407 RowCount [i] = 1 ;
408 }
409 }
410 for (i = 0 ; i < nrow ; i++)
411 {
412 PrevLeaf [i] = EMPTY ;
413 PrevNbr [i] = EMPTY ;
414 SetParent [i] = i ; /* every node is in its own set, by itself */
415 }
416
417 if (stype != 0)
418 {
419
420 /* ------------------------------------------------------------------ */
421 /* symmetric case: LL' = A */
422 /* ------------------------------------------------------------------ */
423
424 /* also determine the number of entries in triu(A) */
425 anz = nrow ;
426 for (k = 0 ; k < nrow ; k++)
427 {
428 /* j is the kth node in the postordered etree */
429 j = initialize_node (k, Post, Parent, ColCount, PrevNbr) ;
430
431 /* for all nonzeros A(i,j) below the diagonal, in column j of A */
432 p = Ap [j] ;
433 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
434 for ( ; p < pend ; p++)
435 {
436 i = Ai [p] ;
437 if (i > j)
438 {
439 /* j is a descendant of i in etree(A) */
440 anz++ ;
441 process_edge (j, i, k, First, PrevNbr, ColCount,
442 PrevLeaf, RowCount, SetParent, Level) ;
443 }
444 }
445 /* update SetParent: UNION (j, Parent [j]) */
446 finalize_node (j, Parent, SetParent) ;
447 }
448 Common->anz = anz ;
449 }
450 else
451 {
452
453 /* ------------------------------------------------------------------ */
454 /* unsymmetric case: LL' = AA' */
455 /* ------------------------------------------------------------------ */
456
457 for (k = 0 ; k < nrow ; k++)
458 {
459 /* inode is the kth node in the postordered etree */
460 inode = initialize_node (k, Post, Parent, ColCount, PrevNbr) ;
461
462 /* for all cols j whose first postordered row is k: */
463 for (j = Head [k] ; j != EMPTY ; j = Anext [j])
464 {
465 /* k is the first postordered row in column j of A */
466 /* for all rows i in column j: */
467 p = Ap [j] ;
468 pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
469 for ( ; p < pend ; p++)
470 {
471 i = Ai [p] ;
472 /* has i already been considered at this step k */
473 if (PrevNbr [i] < k)
474 {
475 /* inode is a descendant of i in etree(AA') */
476 /* process edge (inode,i) and set PrevNbr[i] to k */
477 process_edge (inode, i, k, First, PrevNbr, ColCount,
478 PrevLeaf, RowCount, SetParent, Level) ;
479 }
480 }
481 }
482 /* clear link list k */
483 Head [k] = EMPTY ;
484 /* update SetParent: UNION (inode, Parent [inode]) */
485 finalize_node (inode, Parent, SetParent) ;
486 }
487 }
488
489 /* ---------------------------------------------------------------------- */
490 /* finish computing the column counts */
491 /* ---------------------------------------------------------------------- */
492
493 for (j = 0 ; j < nrow ; j++)
494 {
495 parent = Parent [j] ;
496 if (parent != EMPTY)
497 {
498 /* add the ColCount of j to its parent */
499 ColCount [parent] += ColCount [j] ;
500 }
501 }
502
503 /* ---------------------------------------------------------------------- */
504 /* clear workspace */
505 /* ---------------------------------------------------------------------- */
506
507 Common->mark = EMPTY ;
508 /* CHOLMOD(clear_flag) (Common) ; */
509 CHOLMOD_CLEAR_FLAG (Common) ;
510
511 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
512
513 /* ---------------------------------------------------------------------- */
514 /* flop count and nnz(L) for subsequent LL' numerical factorization */
515 /* ---------------------------------------------------------------------- */
516
517 /* use double to avoid integer overflow. lnz cannot be NaN. */
518 Common->aatfl = fl ;
519 Common->lnz = 0. ;
520 fl = 0 ;
521 for (j = 0 ; j < nrow ; j++)
522 {
523 ff = (double) (ColCount [j]) ;
524 Common->lnz += ff ;
525 fl += ff*ff ;
526 }
527
528 Common->fl = fl ;
529 PRINT1 (("rowcol fl %g lnz %g\n", Common->fl, Common->lnz)) ;
530
531 return (TRUE) ;
532 }
533 #endif
534