cs_amd.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* clear w */
3 static int cs_wclear (int mark, int lemax, int *w, int n)
4 {
5  int k ;
6  if (mark < 2 || (mark + lemax < 0))
7  {
8  for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ;
9  mark = 2 ;
10  }
11  return (mark) ; /* at this point, w [0..n-1] < mark holds */
12 }
13 
14 /* keep off-diagonal entries; drop diagonal entries */
15 static int cs_diag (int i, int j, double aij, void *other) { return (i != j) ; }
16 
17 /* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
18 int *cs_amd (int order, const cs *A) /* order 0:natural, 1:Chol, 2:LU, 3:QR */
19 {
20  cs *C, *A2, *AT ;
21  int *Cp, *Ci, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w,
22  *hhead, *ATp, *ATi, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
23  k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
24  ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ;
25  unsigned int h ;
26  /* --- Construct matrix C ----------------------------------------------- */
27  if (!CS_CSC (A) || order <= 0 || order > 3) return (NULL) ; /* check */
28  AT = cs_transpose (A, 0) ; /* compute A' */
29  if (!AT) return (NULL) ;
30  m = A->m ; n = A->n ;
31  dense = CS_MAX (16, 10 * (int)sqrt ((double) n)) ; /* find dense threshold */
32  dense = CS_MIN (n-2, dense) ;
33  if (order == 1 && n == m)
34  {
35  C = cs_add (A, AT, 0, 0) ; /* C = A+A' */
36  }
37  else if (order == 2)
38  {
39  ATp = AT->p ; /* drop dense columns from AT */
40  ATi = AT->i ;
41  for (p2 = 0, j = 0 ; j < m ; j++)
42  {
43  p = ATp [j] ; /* column j of AT starts here */
44  ATp [j] = p2 ; /* new column j starts here */
45  if (ATp [j+1] - p > dense) continue ; /* skip dense col j */
46  for ( ; p < ATp [j+1] ; p++) ATi [p2++] = ATi [p] ;
47  }
48  ATp [m] = p2 ; /* finalize AT */
49  A2 = cs_transpose (AT, 0) ; /* A2 = AT' */
50  C = A2 ? cs_multiply (AT, A2) : NULL ; /* C=A'*A with no dense rows */
51  cs_spfree (A2) ;
52  }
53  else
54  {
55  C = cs_multiply (AT, A) ; /* C=A'*A */
56  }
57  cs_spfree (AT) ;
58  if (!C) return (NULL) ;
59  cs_fkeep (C, &cs_diag, NULL) ; /* drop diagonal entries */
60  Cp = C->p ;
61  cnz = Cp [n] ;
62  P = (int*) cs_malloc (n+1, sizeof (int)) ; /* allocate result */
63  W = (int*) cs_malloc (8*(n+1), sizeof (int)) ; /* get workspace */
64  t = cnz + cnz/5 + 2*n ; /* add elbow room to C */
65  if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ;
66  len = W ; nv = W + (n+1) ; next = W + 2*(n+1) ;
67  head = W + 3*(n+1) ; elen = W + 4*(n+1) ; degree = W + 5*(n+1) ;
68  w = W + 6*(n+1) ; hhead = W + 7*(n+1) ;
69  last = P ; /* use P as workspace for last */
70  /* --- Initialize quotient graph ---------------------------------------- */
71  for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
72  len [n] = 0 ;
73  nzmax = C->nzmax ;
74  Ci = C->i ;
75  for (i = 0 ; i <= n ; i++)
76  {
77  head [i] = -1 ; /* degree list i is empty */
78  last [i] = -1 ;
79  next [i] = -1 ;
80  hhead [i] = -1 ; /* hash list i is empty */
81  nv [i] = 1 ; /* node i is just one node */
82  w [i] = 1 ; /* node i is alive */
83  elen [i] = 0 ; /* Ek of node i is empty */
84  degree [i] = len [i] ; /* degree of node i */
85  }
86  mark = cs_wclear (0, 0, w, n) ; /* clear w */
87  elen [n] = -2 ; /* n is a dead element */
88  Cp [n] = -1 ; /* n is a root of assembly tree */
89  w [n] = 0 ; /* n is a dead element */
90  /* --- Initialize degree lists ------------------------------------------ */
91  for (i = 0 ; i < n ; i++)
92  {
93  d = degree [i] ;
94  if (d == 0) /* node i is empty */
95  {
96  elen [i] = -2 ; /* element i is dead */
97  nel++ ;
98  Cp [i] = -1 ; /* i is a root of assembly tree */
99  w [i] = 0 ;
100  }
101  else if (d > dense) /* node i is dense */
102  {
103  nv [i] = 0 ; /* absorb i into element n */
104  elen [i] = -1 ; /* node i is dead */
105  nel++ ;
106  Cp [i] = CS_FLIP (n) ;
107  nv [n]++ ;
108  }
109  else
110  {
111  if (head [d] != -1) last [head [d]] = i ;
112  next [i] = head [d] ; /* put node i in degree list d */
113  head [d] = i ;
114  }
115  }
116  while (nel < n) /* while (selecting pivots) do */
117  {
118  /* --- Select node of minimum approximate degree -------------------- */
119  for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
120  if (next [k] != -1) last [next [k]] = -1 ;
121  head [mindeg] = next [k] ; /* remove k from degree list */
122  elenk = elen [k] ; /* elenk = |Ek| */
123  nvk = nv [k] ; /* # of nodes k represents */
124  nel += nvk ; /* nv[k] nodes of A eliminated */
125  /* --- Garbage collection ------------------------------------------- */
126  if (elenk > 0 && cnz + mindeg >= nzmax)
127  {
128  for (j = 0 ; j < n ; j++)
129  {
130  if ((p = Cp [j]) >= 0) /* j is a live node or element */
131  {
132  Cp [j] = Ci [p] ; /* save first entry of object */
133  Ci [p] = CS_FLIP (j) ; /* first entry is now CS_FLIP(j) */
134  }
135  }
136  for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */
137  {
138  if ((j = CS_FLIP (Ci [p++])) >= 0) /* found object j */
139  {
140  Ci [q] = Cp [j] ; /* restore first entry of object */
141  Cp [j] = q++ ; /* new pointer to object j */
142  for (k3 = 0 ; k3 < len [j]-1 ; k3++) Ci [q++] = Ci [p++] ;
143  }
144  }
145  cnz = q ; /* Ci [cnz...nzmax-1] now free */
146  }
147  /* --- Construct new element ---------------------------------------- */
148  dk = 0 ;
149  nv [k] = -nvk ; /* flag k as in Lk */
150  p = Cp [k] ;
151  pk1 = (elenk == 0) ? p : cnz ; /* do in place if elen[k] == 0 */
152  pk2 = pk1 ;
153  for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
154  {
155  if (k1 > elenk)
156  {
157  e = k ; /* search the nodes in k */
158  pj = p ; /* list of nodes starts at Ci[pj]*/
159  ln = len [k] - elenk ; /* length of list of nodes in k */
160  }
161  else
162  {
163  e = Ci [p++] ; /* search the nodes in e */
164  pj = Cp [e] ;
165  ln = len [e] ; /* length of list of nodes in e */
166  }
167  for (k2 = 1 ; k2 <= ln ; k2++)
168  {
169  i = Ci [pj++] ;
170  if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */
171  dk += nvi ; /* degree[Lk] += size of node i */
172  nv [i] = -nvi ; /* negate nv[i] to denote i in Lk*/
173  Ci [pk2++] = i ; /* place i in Lk */
174  if (next [i] != -1) last [next [i]] = last [i] ;
175  if (last [i] != -1) /* remove i from degree list */
176  {
177  next [last [i]] = next [i] ;
178  }
179  else
180  {
181  head [degree [i]] = next [i] ;
182  }
183  }
184  if (e != k)
185  {
186  Cp [e] = CS_FLIP (k) ; /* absorb e into k */
187  w [e] = 0 ; /* e is now a dead element */
188  }
189  }
190  if (elenk != 0) cnz = pk2 ; /* Ci [cnz...nzmax] is free */
191  degree [k] = dk ; /* external degree of k - |Lk\i| */
192  Cp [k] = pk1 ; /* element k is in Ci[pk1..pk2-1] */
193  len [k] = pk2 - pk1 ;
194  elen [k] = -2 ; /* k is now an element */
195  /* --- Find set differences ----------------------------------------- */
196  mark = cs_wclear (mark, lemax, w, n) ; /* clear w if necessary */
197  for (pk = pk1 ; pk < pk2 ; pk++) /* scan 1: find |Le\Lk| */
198  {
199  i = Ci [pk] ;
200  if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */
201  nvi = -nv [i] ; /* nv [i] was negated */
202  wnvi = mark - nvi ;
203  for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++) /* scan Ei */
204  {
205  e = Ci [p] ;
206  if (w [e] >= mark)
207  {
208  w [e] -= nvi ; /* decrement |Le\Lk| */
209  }
210  else if (w [e] != 0) /* ensure e is a live element */
211  {
212  w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */
213  }
214  }
215  }
216  /* --- Degree update ------------------------------------------------ */
217  for (pk = pk1 ; pk < pk2 ; pk++) /* scan2: degree update */
218  {
219  i = Ci [pk] ; /* consider node i in Lk */
220  p1 = Cp [i] ;
221  p2 = p1 + elen [i] - 1 ;
222  pn = p1 ;
223  for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) /* scan Ei */
224  {
225  e = Ci [p] ;
226  if (w [e] != 0) /* e is an unabsorbed element */
227  {
228  dext = w [e] - mark ; /* dext = |Le\Lk| */
229  if (dext > 0)
230  {
231  d += dext ; /* sum up the set differences */
232  Ci [pn++] = e ; /* keep e in Ei */
233  h += e ; /* compute the hash of node i */
234  }
235  else
236  {
237  Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */
238  w [e] = 0 ; /* e is a dead element */
239  }
240  }
241  }
242  elen [i] = pn - p1 + 1 ; /* elen[i] = |Ei| */
243  p3 = pn ;
244  p4 = p1 + len [i] ;
245  for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */
246  {
247  j = Ci [p] ;
248  if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
249  d += nvj ; /* degree(i) += |j| */
250  Ci [pn++] = j ; /* place j in node list of i */
251  h += j ; /* compute hash for node i */
252  }
253  if (d == 0) /* check for mass elimination */
254  {
255  Cp [i] = CS_FLIP (k) ; /* absorb i into k */
256  nvi = -nv [i] ;
257  dk -= nvi ; /* |Lk| -= |i| */
258  nvk += nvi ; /* |k| += nv[i] */
259  nel += nvi ;
260  nv [i] = 0 ;
261  elen [i] = -1 ; /* node i is dead */
262  }
263  else
264  {
265  degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */
266  Ci [pn] = Ci [p3] ; /* move first node to end */
267  Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */
268  Ci [p1] = k ; /* add k as 1st element in of Ei */
269  len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */
270  h %= n ; /* finalize hash of i */
271  next [i] = hhead [h] ; /* place i in hash bucket */
272  hhead [h] = i ;
273  last [i] = h ; /* save hash of i in last[i] */
274  }
275  } /* scan2 is done */
276  degree [k] = dk ; /* finalize |Lk| */
277  lemax = CS_MAX (lemax, dk) ;
278  mark = cs_wclear (mark+lemax, lemax, w, n) ; /* clear w */
279  /* --- Supernode detection ------------------------------------------ */
280  for (pk = pk1 ; pk < pk2 ; pk++)
281  {
282  i = Ci [pk] ;
283  if (nv [i] >= 0) continue ; /* skip if i is dead */
284  h = last [i] ; /* scan hash bucket of node i */
285  i = hhead [h] ;
286  hhead [h] = -1 ; /* hash bucket will be empty */
287  for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
288  {
289  ln = len [i] ;
290  eln = elen [i] ;
291  for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark;
292  jlast = i ;
293  for (j = next [i] ; j != -1 ; ) /* compare i with all j */
294  {
295  ok = (len [j] == ln) && (elen [j] == eln) ;
296  for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
297  {
298  if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/
299  }
300  if (ok) /* i and j are identical */
301  {
302  Cp [j] = CS_FLIP (i) ; /* absorb j into i */
303  nv [i] += nv [j] ;
304  nv [j] = 0 ;
305  elen [j] = -1 ; /* node j is dead */
306  j = next [j] ; /* delete j from hash bucket */
307  next [jlast] = j ;
308  }
309  else
310  {
311  jlast = j ; /* j and i are different */
312  j = next [j] ;
313  }
314  }
315  }
316  }
317  /* --- Finalize new element------------------------------------------ */
318  for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */
319  {
320  i = Ci [pk] ;
321  if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
322  nv [i] = nvi ; /* restore nv[i] */
323  d = degree [i] + dk - nvi ; /* compute external degree(i) */
324  d = CS_MIN (d, n - nel - nvi) ;
325  if (head [d] != -1) last [head [d]] = i ;
326  next [i] = head [d] ; /* put i back in degree list */
327  last [i] = -1 ;
328  head [d] = i ;
329  mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */
330  degree [i] = d ;
331  Ci [p++] = i ; /* place i in Lk */
332  }
333  nv [k] = nvk ; /* # nodes absorbed into k */
334  if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/
335  {
336  Cp [k] = -1 ; /* k is a root of the tree */
337  w [k] = 0 ; /* k is now a dead element */
338  }
339  if (elenk != 0) cnz = p ; /* free unused space in Lk */
340  }
341  /* --- Postordering ----------------------------------------------------- */
342  for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
343  for (j = 0 ; j <= n ; j++) head [j] = -1 ;
344  for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */
345  {
346  if (nv [j] > 0) continue ; /* skip if j is an element */
347  next [j] = head [Cp [j]] ; /* place j in list of its parent */
348  head [Cp [j]] = j ;
349  }
350  for (e = n ; e >= 0 ; e--) /* place elements in lists */
351  {
352  if (nv [e] <= 0) continue ; /* skip unless e is an element */
353  if (Cp [e] != -1)
354  {
355  next [e] = head [Cp [e]] ; /* place e in list of its parent */
356  head [Cp [e]] = e ;
357  }
358  }
359  for (k = 0, i = 0 ; i <= n ; i++) /* postorder the assembly tree */
360  {
361  if (Cp [i] == -1) k = cs_tdfs (i, k, head, next, P, w) ;
362  }
363  return (cs_idone (P, C, W, 1)) ;
364 }
IntermediateState sqrt(const Expression &arg)
cs * cs_spfree(cs *A)
Definition: cs_util.c:32
static int cs_wclear(int mark, int lemax, int *w, int n)
Definition: cs_amd.c:3
int n
Definition: cs.h:20
int cs_fkeep(cs *A, int(*fkeep)(int, int, double, void *), void *other)
Definition: cs_fkeep.c:3
int * cs_amd(int order, const cs *A)
Definition: cs_amd.c:18
#define CS_FLIP(i)
Definition: cs.h:136
cs * cs_transpose(const cs *A, int values)
Definition: cs_transpose.c:3
int nzmax
Definition: cs.h:18
#define CS_MAX(a, b)
Definition: cs.h:134
int cs_sprealloc(cs *A, int nzmax)
Definition: cs_util.c:18
#define CS_CSC(A)
Definition: cs.h:140
#define A2
int * cs_idone(int *p, cs *C, void *w, int ok)
Definition: cs_util.c:97
cs * cs_add(const cs *A, const cs *B, double alpha, double beta)
Definition: cs_add.c:3
int * p
Definition: cs.h:21
Definition: cs.h:16
int * i
Definition: cs.h:22
cs * cs_multiply(const cs *A, const cs *B)
Definition: cs_multiply.c:3
int cs_tdfs(int j, int k, int *head, const int *next, int *post, int *stack)
Definition: cs_tdfs.c:3
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
SegmentReturnType head(Index vecSize)
Definition: BlockMethods.h:781
IntermediateState ln(const Expression &arg)
Expression next(const Expression &arg)
#define CS_MIN(a, b)
Definition: cs.h:135
static int cs_diag(int i, int j, double aij, void *other)
Definition: cs_amd.c:15
int m
Definition: cs.h:19


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