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