Amd.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 /*
26 
27 NOTE: this routine has been adapted from the CSparse library:
28 
29 Copyright (c) 2006, Timothy A. Davis.
30 http://www.cise.ufl.edu/research/sparse/CSparse
31 
32 CSparse is free software; you can redistribute it and/or
33 modify it under the terms of the GNU Lesser General Public
34 License as published by the Free Software Foundation; either
35 version 2.1 of the License, or (at your option) any later version.
36 
37 CSparse is distributed in the hope that it will be useful,
38 but WITHOUT ANY WARRANTY; without even the implied warranty of
39 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
40 Lesser General Public License for more details.
41 
42 You should have received a copy of the GNU Lesser General Public
43 License along with this Module; if not, write to the Free Software
44 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
45 
46 */
47 
48 #ifndef EIGEN_SPARSE_AMD_H
49 #define EIGEN_SPARSE_AMD_H
50 
51 namespace Eigen {
52 
53 namespace internal {
54 
55 template<typename T> inline T amd_flip(const T& i) { return -i-2; }
56 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; }
57 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; }
58 template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); }
59 
60 /* clear w */
61 template<typename Index>
62 static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
63 {
64  Index k;
65  if(mark < 2 || (mark + lemax < 0))
66  {
67  for(k = 0; k < n; k++)
68  if(w[k] != 0)
69  w[k] = 1;
70  mark = 2;
71  }
72  return (mark); /* at this point, w[0..n-1] < mark holds */
73 }
74 
75 /* depth-first search and postorder of a tree rooted at node j */
76 template<typename Index>
77 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
78 {
79  int i, p, top = 0;
80  if(!head || !next || !post || !stack) return (-1); /* check inputs */
81  stack[0] = j; /* place j on the stack */
82  while (top >= 0) /* while (stack is not empty) */
83  {
84  p = stack[top]; /* p = top of stack */
85  i = head[p]; /* i = youngest child of p */
86  if(i == -1)
87  {
88  top--; /* p has no unordered children left */
89  post[k++] = p; /* node p is the kth postordered node */
90  }
91  else
92  {
93  head[p] = next[i]; /* remove i from children of p */
94  stack[++top] = i; /* start dfs on child node i */
95  }
96  }
97  return k;
98 }
99 
100 
106 template<typename Scalar, typename Index>
108 {
110 
111  int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
112  k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
113  ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
114  unsigned int h;
115 
116  Index n = C.cols();
117  dense = std::max<Index> (16, 10 * sqrt ((double) n)); /* find dense threshold */
118  dense = std::min<Index> (n-2, dense);
119 
120  Index cnz = C.nonZeros();
121  perm.resize(n+1);
122  t = cnz + cnz/5 + 2*n; /* add elbow room to C */
123  C.resizeNonZeros(t);
124 
125  Index* W = new Index[8*(n+1)]; /* get workspace */
126  Index* len = W;
127  Index* nv = W + (n+1);
128  Index* next = W + 2*(n+1);
129  Index* head = W + 3*(n+1);
130  Index* elen = W + 4*(n+1);
131  Index* degree = W + 5*(n+1);
132  Index* w = W + 6*(n+1);
133  Index* hhead = W + 7*(n+1);
134  Index* last = perm.indices().data(); /* use P as workspace for last */
135 
136  /* --- Initialize quotient graph ---------------------------------------- */
137  Index* Cp = C.outerIndexPtr();
138  Index* Ci = C.innerIndexPtr();
139  for(k = 0; k < n; k++)
140  len[k] = Cp[k+1] - Cp[k];
141  len[n] = 0;
142  nzmax = t;
143 
144  for(i = 0; i <= n; i++)
145  {
146  head[i] = -1; // degree list i is empty
147  last[i] = -1;
148  next[i] = -1;
149  hhead[i] = -1; // hash list i is empty
150  nv[i] = 1; // node i is just one node
151  w[i] = 1; // node i is alive
152  elen[i] = 0; // Ek of node i is empty
153  degree[i] = len[i]; // degree of node i
154  }
155  mark = internal::cs_wclear<Index>(0, 0, w, n); /* clear w */
156  elen[n] = -2; /* n is a dead element */
157  Cp[n] = -1; /* n is a root of assembly tree */
158  w[n] = 0; /* n is a dead element */
159 
160  /* --- Initialize degree lists ------------------------------------------ */
161  for(i = 0; i < n; i++)
162  {
163  d = degree[i];
164  if(d == 0) /* node i is empty */
165  {
166  elen[i] = -2; /* element i is dead */
167  nel++;
168  Cp[i] = -1; /* i is a root of assembly tree */
169  w[i] = 0;
170  }
171  else if(d > dense) /* node i is dense */
172  {
173  nv[i] = 0; /* absorb i into element n */
174  elen[i] = -1; /* node i is dead */
175  nel++;
176  Cp[i] = amd_flip (n);
177  nv[n]++;
178  }
179  else
180  {
181  if(head[d] != -1) last[head[d]] = i;
182  next[i] = head[d]; /* put node i in degree list d */
183  head[d] = i;
184  }
185  }
186 
187  while (nel < n) /* while (selecting pivots) do */
188  {
189  /* --- Select node of minimum approximate degree -------------------- */
190  for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
191  if(next[k] != -1) last[next[k]] = -1;
192  head[mindeg] = next[k]; /* remove k from degree list */
193  elenk = elen[k]; /* elenk = |Ek| */
194  nvk = nv[k]; /* # of nodes k represents */
195  nel += nvk; /* nv[k] nodes of A eliminated */
196 
197  /* --- Garbage collection ------------------------------------------- */
198  if(elenk > 0 && cnz + mindeg >= nzmax)
199  {
200  for(j = 0; j < n; j++)
201  {
202  if((p = Cp[j]) >= 0) /* j is a live node or element */
203  {
204  Cp[j] = Ci[p]; /* save first entry of object */
205  Ci[p] = amd_flip (j); /* first entry is now amd_flip(j) */
206  }
207  }
208  for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
209  {
210  if((j = amd_flip (Ci[p++])) >= 0) /* found object j */
211  {
212  Ci[q] = Cp[j]; /* restore first entry of object */
213  Cp[j] = q++; /* new pointer to object j */
214  for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
215  }
216  }
217  cnz = q; /* Ci[cnz...nzmax-1] now free */
218  }
219 
220  /* --- Construct new element ---------------------------------------- */
221  dk = 0;
222  nv[k] = -nvk; /* flag k as in Lk */
223  p = Cp[k];
224  pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */
225  pk2 = pk1;
226  for(k1 = 1; k1 <= elenk + 1; k1++)
227  {
228  if(k1 > elenk)
229  {
230  e = k; /* search the nodes in k */
231  pj = p; /* list of nodes starts at Ci[pj]*/
232  ln = len[k] - elenk; /* length of list of nodes in k */
233  }
234  else
235  {
236  e = Ci[p++]; /* search the nodes in e */
237  pj = Cp[e];
238  ln = len[e]; /* length of list of nodes in e */
239  }
240  for(k2 = 1; k2 <= ln; k2++)
241  {
242  i = Ci[pj++];
243  if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
244  dk += nvi; /* degree[Lk] += size of node i */
245  nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/
246  Ci[pk2++] = i; /* place i in Lk */
247  if(next[i] != -1) last[next[i]] = last[i];
248  if(last[i] != -1) /* remove i from degree list */
249  {
250  next[last[i]] = next[i];
251  }
252  else
253  {
254  head[degree[i]] = next[i];
255  }
256  }
257  if(e != k)
258  {
259  Cp[e] = amd_flip (k); /* absorb e into k */
260  w[e] = 0; /* e is now a dead element */
261  }
262  }
263  if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */
264  degree[k] = dk; /* external degree of k - |Lk\i| */
265  Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */
266  len[k] = pk2 - pk1;
267  elen[k] = -2; /* k is now an element */
268 
269  /* --- Find set differences ----------------------------------------- */
270  mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */
271  for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */
272  {
273  i = Ci[pk];
274  if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
275  nvi = -nv[i]; /* nv[i] was negated */
276  wnvi = mark - nvi;
277  for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */
278  {
279  e = Ci[p];
280  if(w[e] >= mark)
281  {
282  w[e] -= nvi; /* decrement |Le\Lk| */
283  }
284  else if(w[e] != 0) /* ensure e is a live element */
285  {
286  w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
287  }
288  }
289  }
290 
291  /* --- Degree update ------------------------------------------------ */
292  for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */
293  {
294  i = Ci[pk]; /* consider node i in Lk */
295  p1 = Cp[i];
296  p2 = p1 + elen[i] - 1;
297  pn = p1;
298  for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */
299  {
300  e = Ci[p];
301  if(w[e] != 0) /* e is an unabsorbed element */
302  {
303  dext = w[e] - mark; /* dext = |Le\Lk| */
304  if(dext > 0)
305  {
306  d += dext; /* sum up the set differences */
307  Ci[pn++] = e; /* keep e in Ei */
308  h += e; /* compute the hash of node i */
309  }
310  else
311  {
312  Cp[e] = amd_flip (k); /* aggressive absorb. e->k */
313  w[e] = 0; /* e is a dead element */
314  }
315  }
316  }
317  elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */
318  p3 = pn;
319  p4 = p1 + len[i];
320  for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
321  {
322  j = Ci[p];
323  if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
324  d += nvj; /* degree(i) += |j| */
325  Ci[pn++] = j; /* place j in node list of i */
326  h += j; /* compute hash for node i */
327  }
328  if(d == 0) /* check for mass elimination */
329  {
330  Cp[i] = amd_flip (k); /* absorb i into k */
331  nvi = -nv[i];
332  dk -= nvi; /* |Lk| -= |i| */
333  nvk += nvi; /* |k| += nv[i] */
334  nel += nvi;
335  nv[i] = 0;
336  elen[i] = -1; /* node i is dead */
337  }
338  else
339  {
340  degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */
341  Ci[pn] = Ci[p3]; /* move first node to end */
342  Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */
343  Ci[p1] = k; /* add k as 1st element in of Ei */
344  len[i] = pn - p1 + 1; /* new len of adj. list of node i */
345  h %= n; /* finalize hash of i */
346  next[i] = hhead[h]; /* place i in hash bucket */
347  hhead[h] = i;
348  last[i] = h; /* save hash of i in last[i] */
349  }
350  } /* scan2 is done */
351  degree[k] = dk; /* finalize |Lk| */
352  lemax = std::max<Index>(lemax, dk);
353  mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */
354 
355  /* --- Supernode detection ------------------------------------------ */
356  for(pk = pk1; pk < pk2; pk++)
357  {
358  i = Ci[pk];
359  if(nv[i] >= 0) continue; /* skip if i is dead */
360  h = last[i]; /* scan hash bucket of node i */
361  i = hhead[h];
362  hhead[h] = -1; /* hash bucket will be empty */
363  for(; i != -1 && next[i] != -1; i = next[i], mark++)
364  {
365  ln = len[i];
366  eln = elen[i];
367  for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
368  jlast = i;
369  for(j = next[i]; j != -1; ) /* compare i with all j */
370  {
371  ok = (len[j] == ln) && (elen[j] == eln);
372  for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
373  {
374  if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/
375  }
376  if(ok) /* i and j are identical */
377  {
378  Cp[j] = amd_flip (i); /* absorb j into i */
379  nv[i] += nv[j];
380  nv[j] = 0;
381  elen[j] = -1; /* node j is dead */
382  j = next[j]; /* delete j from hash bucket */
383  next[jlast] = j;
384  }
385  else
386  {
387  jlast = j; /* j and i are different */
388  j = next[j];
389  }
390  }
391  }
392  }
393 
394  /* --- Finalize new element------------------------------------------ */
395  for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */
396  {
397  i = Ci[pk];
398  if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
399  nv[i] = nvi; /* restore nv[i] */
400  d = degree[i] + dk - nvi; /* compute external degree(i) */
401  d = std::min<Index> (d, n - nel - nvi);
402  if(head[d] != -1) last[head[d]] = i;
403  next[i] = head[d]; /* put i back in degree list */
404  last[i] = -1;
405  head[d] = i;
406  mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */
407  degree[i] = d;
408  Ci[p++] = i; /* place i in Lk */
409  }
410  nv[k] = nvk; /* # nodes absorbed into k */
411  if((len[k] = p-pk1) == 0) /* length of adj list of element k*/
412  {
413  Cp[k] = -1; /* k is a root of the tree */
414  w[k] = 0; /* k is now a dead element */
415  }
416  if(elenk != 0) cnz = p; /* free unused space in Lk */
417  }
418 
419  /* --- Postordering ----------------------------------------------------- */
420  for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */
421  for(j = 0; j <= n; j++) head[j] = -1;
422  for(j = n; j >= 0; j--) /* place unordered nodes in lists */
423  {
424  if(nv[j] > 0) continue; /* skip if j is an element */
425  next[j] = head[Cp[j]]; /* place j in list of its parent */
426  head[Cp[j]] = j;
427  }
428  for(e = n; e >= 0; e--) /* place elements in lists */
429  {
430  if(nv[e] <= 0) continue; /* skip unless e is an element */
431  if(Cp[e] != -1)
432  {
433  next[e] = head[Cp[e]]; /* place e in list of its parent */
434  head[Cp[e]] = e;
435  }
436  }
437  for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */
438  {
439  if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w);
440  }
441 
442  perm.indices().conservativeResize(n);
443 
444  delete[] W;
445 }
446 
447 } // namespace internal
448 
449 } // end namespace Eigen
450 
451 #endif // EIGEN_SPARSE_AMD_H