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 */
3 #define HEAD(k,j) (ata ? head [k] : j)
4 #define NEXT(J) (ata ? next [J] : -1)
5 static void init_ata (cs *AT, const int *post, int *w, int **head, int **next)
6 {
7  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 int *cs_counts (const cs *A, const int *parent, const int *post, int ata)
18 {
19  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 = (int*) cs_malloc (n, sizeof (int)) ; /* allocate result */
26  w = (int*) cs_malloc (s, sizeof (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 }
#define HEAD(k, j)
Definition: cs_counts.c:3
int n
Definition: cs.h:20
cs * cs_transpose(const cs *A, int values)
Definition: cs_transpose.c:3
int * cs_counts(const cs *A, const int *parent, const int *post, int ata)
Definition: cs_counts.c:17
int cs_leaf(int i, int j, const int *first, int *maxfirst, int *prevleaf, int *ancestor, int *jleaf)
Definition: cs_leaf.c:3
#define CS_CSC(A)
Definition: cs.h:140
int * cs_idone(int *p, cs *C, void *w, int ok)
Definition: cs_util.c:97
int * p
Definition: cs.h:21
Definition: cs.h:16
int * i
Definition: cs.h:22
#define NEXT(J)
Definition: cs_counts.c:4
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
SegmentReturnType head(Index vecSize)
Definition: BlockMethods.h:781
Expression next(const Expression &arg)
static void init_ata(cs *AT, const int *post, int *w, int **head, int **next)
Definition: cs_counts.c:5
#define CS_MIN(a, b)
Definition: cs.h:135
int m
Definition: cs.h:19


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:31