EVSL  1.1.0 EigenValues Slicing Library
cs_counts.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* column counts of LL'=A or LL'=A'A, given parent & post ordering */
4 #define NEXT(J) (ata ? next [J] : -1)
5 static void init_ata (cs *AT, const CS_INT *post, CS_INT *w, CS_INT **head, CS_INT **next)
6 {
7  CS_INT i, k, p, m = AT->n, n = AT->m, *ATp = AT->p, *ATi = AT->i ;
8  *head = w+4*n, *next = w+5*n+1 ;
9  for (k = 0 ; k < n ; k++) w [post [k]] = k ; /* invert post */
10  for (i = 0 ; i < m ; i++)
11  {
12  for (k = n, p = ATp[i] ; p < ATp[i+1] ; p++) k = CS_MIN (k, w [ATi[p]]);
13  (*next) [i] = (*head) [k] ; /* place row i in linked list k */
14  (*head) [k] = i ;
15  }
16 }
17 CS_INT *cs_counts (const cs *A, const CS_INT *parent, const CS_INT *post, CS_INT ata)
18 {
19  CS_INT i, j, k, n, m, J, s, p, q, jleaf, *ATp, *ATi, *maxfirst, *prevleaf,
20  *ancestor, *head = NULL, *next = NULL, *colcount, *w, *first, *delta ;
21  cs *AT ;
22  if (!CS_CSC (A) || !parent || !post) return (NULL) ; /* check inputs */
23  m = A->m ; n = A->n ;
24  s = 4*n + (ata ? (n+m+1) : 0) ;
25  delta = colcount = cs_malloc (n, sizeof (CS_INT)) ; /* allocate result */
26  w = cs_malloc (s, sizeof (CS_INT)) ; /* get workspace */
27  AT = cs_transpose (A, 0) ; /* AT = A' */
28  if (!AT || !colcount || !w) return (cs_idone (colcount, AT, w, 0)) ;
29  ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ; first = w+3*n ;
30  for (k = 0 ; k < s ; k++) w [k] = -1 ; /* clear workspace w [0..s-1] */
31  for (k = 0 ; k < n ; k++) /* find first [j] */
32  {
33  j = post [k] ;
34  delta [j] = (first [j] == -1) ? 1 : 0 ; /* delta[j]=1 if j is a leaf */
35  for ( ; j != -1 && first [j] == -1 ; j = parent [j]) first [j] = k ;
36  }
37  ATp = AT->p ; ATi = AT->i ;
38  if (ata) init_ata (AT, post, w, &head, &next) ;
39  for (i = 0 ; i < n ; i++) ancestor [i] = i ; /* each node in its own set */
40  for (k = 0 ; k < n ; k++)
41  {
42  j = post [k] ; /* j is the kth node in postordered etree */
43  if (parent [j] != -1) delta [parent [j]]-- ; /* j is not a root */
44  for (J = HEAD (k,j) ; J != -1 ; J = NEXT (J)) /* J=j for LL'=A case */
45  {
46  for (p = ATp [J] ; p < ATp [J+1] ; p++)
47  {
48  i = ATi [p] ;
49  q = cs_leaf (i, j, first, maxfirst, prevleaf, ancestor, &jleaf);
50  if (jleaf >= 1) delta [j]++ ; /* A(i,j) is in skeleton */
51  if (jleaf == 2) delta [q]-- ; /* account for overlap in q */
52  }
53  }
54  if (parent [j] != -1) ancestor [j] = parent [j] ;
55  }
56  for (j = 0 ; j < n ; j++) /* sum up delta's of each child */
57  {
58  if (parent [j] != -1) colcount [parent [j]] += colcount [j] ;
59  }
60  return (cs_idone (colcount, AT, w, 1)) ; /* success: free workspace */
61 }
CS_INT cs_leaf(CS_INT i, CS_INT j, const CS_INT *first, CS_INT *maxfirst, CS_INT *prevleaf, CS_INT *ancestor, CS_INT *jleaf)
Definition: cs_leaf.c:3
#define cs
Definition: cs.h:637
Definition: cs_counts.c:3
#define CS_CSC(A)
Definition: cs.h:659
cs * cs_transpose(const cs *A, CS_INT values)
Definition: cs_transpose.c:3
#define NEXT(J)
Definition: cs_counts.c:4
#define CS_INT
Definition: cs.h:627
CS_INT * cs_counts(const cs *A, const CS_INT *parent, const CS_INT *post, CS_INT ata)
Definition: cs_counts.c:17
#define CS_MIN(a, b)
Definition: cs.h:654
void * cs_malloc(CS_INT n, size_t size)
Definition: cs_malloc.c:10
CS_INT * cs_idone(CS_INT *p, cs *C, void *w, CS_INT ok)
Definition: cs_util.c:98