1 /* Compute the row counts of the Cholesky factor L of the matrix A. Uses
2 * the lower triangular part of A. */
3
4 #include "cs_mex.h"
5
6 static
firstdesc(CS_INT n,CS_INT * parent,CS_INT * post,CS_INT * first,CS_INT * level)7 void firstdesc (CS_INT n, CS_INT *parent, CS_INT *post, CS_INT *first,
8 CS_INT *level)
9 {
10 CS_INT len, i, k, r, s ;
11 for (i = 0 ; i < n ; i++) first [i] = -1 ;
12 for (k = 0 ; k < n ; k++)
13 {
14 i = post [k] ; /* node i of etree is kth postordered node */
15 len = 0 ; /* traverse from i towards the root */
16 for (r = i ; r != -1 && first [r] == -1 ; r = parent [r], len++)
17 first [r] = k ;
18 len += (r == -1) ? (-1) : level [r] ; /* root node or end of path */
19 for (s = i ; s != r ; s = parent [s]) level [s] = len-- ;
20 }
21 }
22
23 /* return rowcount [0..n-1] */
24 static
rowcnt(cs_dl * A,CS_INT * parent,CS_INT * post)25 CS_INT *rowcnt (cs_dl *A, CS_INT *parent, CS_INT *post)
26 {
27 CS_INT i, j, k, len, s, p, jprev, q, n, sparent, jleaf, *Ap, *Ai, *maxfirst,
28 *ancestor, *prevleaf, *w, *first, *level, *rowcount ;
29 n = A->n ; Ap = A->p ; Ai = A->i ; /* get A */
30 w = cs_dl_malloc (5*n, sizeof (CS_INT)) ; /* get workspace */
31 ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ; first = w+3*n ;
32 level = w+4*n ;
33 rowcount = cs_dl_malloc (n, sizeof (CS_INT)) ; /* allocate result */
34 firstdesc (n, parent, post, first, level) ; /* find first and level */
35 for (i = 0 ; i < n ; i++)
36 {
37 rowcount [i] = 1 ; /* count the diagonal of L */
38 prevleaf [i] = -1 ; /* no previous leaf of the ith row subtree */
39 maxfirst [i] = -1 ; /* max first[j] for node j in ith subtree */
40 ancestor [i] = i ; /* every node is in its own set, by itself */
41 }
42 for (k = 0 ; k < n ; k++)
43 {
44 j = post [k] ; /* j is the kth node in the postordered etree */
45 for (p = Ap [j] ; p < Ap [j+1] ; p++)
46 {
47 i = Ai [p] ;
48 q = cs_dl_leaf (i, j, first, maxfirst, prevleaf, ancestor, &jleaf) ;
49 if (jleaf) rowcount [i] += (level [j] - level [q]) ;
50 }
51 if (parent [j] != -1) ancestor [j] = parent [j] ;
52 }
53 cs_dl_free (w) ;
54 return (rowcount) ;
55 }
56
mexFunction(int nargout,mxArray * pargout[],int nargin,const mxArray * pargin[])57 void mexFunction
58 (
59 int nargout,
60 mxArray *pargout [ ],
61 int nargin,
62 const mxArray *pargin [ ]
63 )
64 {
65 cs_dl *A, Amatrix ;
66 double *x ;
67 CS_INT i, m, n, *parent, *post, *rowcount ;
68
69 if (nargout > 1 || nargin != 3)
70 {
71 mexErrMsgTxt ("Usage: r = cs_rowcnt(A,parent,post)") ;
72 }
73
74 /* get inputs */
75 A = cs_dl_mex_get_sparse (&Amatrix, 1, 0, pargin [0]) ;
76 n = A->n ;
77
78 parent = cs_dl_mex_get_int (n, pargin [1], &i, 0) ;
79 post = cs_dl_mex_get_int (n, pargin [2], &i, 1) ;
80
81 rowcount = rowcnt (A, parent, post) ;
82
83 pargout [0] = mxCreateDoubleMatrix (1, n, mxREAL) ;
84 x = mxGetPr (pargout [0]) ;
85 for (i = 0 ; i < n ; i++) x [i] = rowcount [i] ;
86
87 cs_dl_free (rowcount) ;
88 }
89