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 {
109  using std::sqrt;
111 
112  int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
113  k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
114  ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
115  unsigned int h;
116 
117  Index n = C.cols();
118  dense = std::max<Index> (16, Index(10 * sqrt(double(n)))); /* find dense threshold */
119  dense = std::min<Index> (n-2, dense);
120 
121  Index cnz = C.nonZeros();
122  perm.resize(n+1);
123  t = cnz + cnz/5 + 2*n; /* add elbow room to C */
124  C.resizeNonZeros(t);
125 
126  Index* W = new Index[8*(n+1)]; /* get workspace */
127  Index* len = W;
128  Index* nv = W + (n+1);
129  Index* next = W + 2*(n+1);
130  Index* head = W + 3*(n+1);
131  Index* elen = W + 4*(n+1);
132  Index* degree = W + 5*(n+1);
133  Index* w = W + 6*(n+1);
134  Index* hhead = W + 7*(n+1);
135  Index* last = perm.indices().data(); /* use P as workspace for last */
136 
137  /* --- Initialize quotient graph ---------------------------------------- */
138  Index* Cp = C.outerIndexPtr();
139  Index* Ci = C.innerIndexPtr();
140  for(k = 0; k < n; k++)
141  len[k] = Cp[k+1] - Cp[k];
142  len[n] = 0;
143  nzmax = t;
144 
145  for(i = 0; i <= n; i++)
146  {
147  head[i] = -1; // degree list i is empty
148  last[i] = -1;
149  next[i] = -1;
150  hhead[i] = -1; // hash list i is empty
151  nv[i] = 1; // node i is just one node
152  w[i] = 1; // node i is alive
153  elen[i] = 0; // Ek of node i is empty
154  degree[i] = len[i]; // degree of node i
155  }
156  mark = internal::cs_wclear<Index>(0, 0, w, n); /* clear w */
157  elen[n] = -2; /* n is a dead element */
158  Cp[n] = -1; /* n is a root of assembly tree */
159  w[n] = 0; /* n is a dead element */
160 
161  /* --- Initialize degree lists ------------------------------------------ */
162  for(i = 0; i < n; i++)
163  {
164  d = degree[i];
165  if(d == 0) /* node i is empty */
166  {
167  elen[i] = -2; /* element i is dead */
168  nel++;
169  Cp[i] = -1; /* i is a root of assembly tree */
170  w[i] = 0;
171  }
172  else if(d > dense) /* node i is dense */
173  {
174  nv[i] = 0; /* absorb i into element n */
175  elen[i] = -1; /* node i is dead */
176  nel++;
177  Cp[i] = amd_flip (n);
178  nv[n]++;
179  }
180  else
181  {
182  if(head[d] != -1) last[head[d]] = i;
183  next[i] = head[d]; /* put node i in degree list d */
184  head[d] = i;
185  }
186  }
187 
188  while (nel < n) /* while (selecting pivots) do */
189  {
190  /* --- Select node of minimum approximate degree -------------------- */
191  for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
192  if(next[k] != -1) last[next[k]] = -1;
193  head[mindeg] = next[k]; /* remove k from degree list */
194  elenk = elen[k]; /* elenk = |Ek| */
195  nvk = nv[k]; /* # of nodes k represents */
196  nel += nvk; /* nv[k] nodes of A eliminated */
197 
198  /* --- Garbage collection ------------------------------------------- */
199  if(elenk > 0 && cnz + mindeg >= nzmax)
200  {
201  for(j = 0; j < n; j++)
202  {
203  if((p = Cp[j]) >= 0) /* j is a live node or element */
204  {
205  Cp[j] = Ci[p]; /* save first entry of object */
206  Ci[p] = amd_flip (j); /* first entry is now amd_flip(j) */
207  }
208  }
209  for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
210  {
211  if((j = amd_flip (Ci[p++])) >= 0) /* found object j */
212  {
213  Ci[q] = Cp[j]; /* restore first entry of object */
214  Cp[j] = q++; /* new pointer to object j */
215  for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
216  }
217  }
218  cnz = q; /* Ci[cnz...nzmax-1] now free */
219  }
220 
221  /* --- Construct new element ---------------------------------------- */
222  dk = 0;
223  nv[k] = -nvk; /* flag k as in Lk */
224  p = Cp[k];
225  pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */
226  pk2 = pk1;
227  for(k1 = 1; k1 <= elenk + 1; k1++)
228  {
229  if(k1 > elenk)
230  {
231  e = k; /* search the nodes in k */
232  pj = p; /* list of nodes starts at Ci[pj]*/
233  ln = len[k] - elenk; /* length of list of nodes in k */
234  }
235  else
236  {
237  e = Ci[p++]; /* search the nodes in e */
238  pj = Cp[e];
239  ln = len[e]; /* length of list of nodes in e */
240  }
241  for(k2 = 1; k2 <= ln; k2++)
242  {
243  i = Ci[pj++];
244  if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
245  dk += nvi; /* degree[Lk] += size of node i */
246  nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/
247  Ci[pk2++] = i; /* place i in Lk */
248  if(next[i] != -1) last[next[i]] = last[i];
249  if(last[i] != -1) /* remove i from degree list */
250  {
251  next[last[i]] = next[i];
252  }
253  else
254  {
255  head[degree[i]] = next[i];
256  }
257  }
258  if(e != k)
259  {
260  Cp[e] = amd_flip (k); /* absorb e into k */
261  w[e] = 0; /* e is now a dead element */
262  }
263  }
264  if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */
265  degree[k] = dk; /* external degree of k - |Lk\i| */
266  Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */
267  len[k] = pk2 - pk1;
268  elen[k] = -2; /* k is now an element */
269 
270  /* --- Find set differences ----------------------------------------- */
271  mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */
272  for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */
273  {
274  i = Ci[pk];
275  if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
276  nvi = -nv[i]; /* nv[i] was negated */
277  wnvi = mark - nvi;
278  for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */
279  {
280  e = Ci[p];
281  if(w[e] >= mark)
282  {
283  w[e] -= nvi; /* decrement |Le\Lk| */
284  }
285  else if(w[e] != 0) /* ensure e is a live element */
286  {
287  w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
288  }
289  }
290  }
291 
292  /* --- Degree update ------------------------------------------------ */
293  for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */
294  {
295  i = Ci[pk]; /* consider node i in Lk */
296  p1 = Cp[i];
297  p2 = p1 + elen[i] - 1;
298  pn = p1;
299  for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */
300  {
301  e = Ci[p];
302  if(w[e] != 0) /* e is an unabsorbed element */
303  {
304  dext = w[e] - mark; /* dext = |Le\Lk| */
305  if(dext > 0)
306  {
307  d += dext; /* sum up the set differences */
308  Ci[pn++] = e; /* keep e in Ei */
309  h += e; /* compute the hash of node i */
310  }
311  else
312  {
313  Cp[e] = amd_flip (k); /* aggressive absorb. e->k */
314  w[e] = 0; /* e is a dead element */
315  }
316  }
317  }
318  elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */
319  p3 = pn;
320  p4 = p1 + len[i];
321  for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
322  {
323  j = Ci[p];
324  if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
325  d += nvj; /* degree(i) += |j| */
326  Ci[pn++] = j; /* place j in node list of i */
327  h += j; /* compute hash for node i */
328  }
329  if(d == 0) /* check for mass elimination */
330  {
331  Cp[i] = amd_flip (k); /* absorb i into k */
332  nvi = -nv[i];
333  dk -= nvi; /* |Lk| -= |i| */
334  nvk += nvi; /* |k| += nv[i] */
335  nel += nvi;
336  nv[i] = 0;
337  elen[i] = -1; /* node i is dead */
338  }
339  else
340  {
341  degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */
342  Ci[pn] = Ci[p3]; /* move first node to end */
343  Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */
344  Ci[p1] = k; /* add k as 1st element in of Ei */
345  len[i] = pn - p1 + 1; /* new len of adj. list of node i */
346  h %= n; /* finalize hash of i */
347  next[i] = hhead[h]; /* place i in hash bucket */
348  hhead[h] = i;
349  last[i] = h; /* save hash of i in last[i] */
350  }
351  } /* scan2 is done */
352  degree[k] = dk; /* finalize |Lk| */
353  lemax = std::max<Index>(lemax, dk);
354  mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */
355 
356  /* --- Supernode detection ------------------------------------------ */
357  for(pk = pk1; pk < pk2; pk++)
358  {
359  i = Ci[pk];
360  if(nv[i] >= 0) continue; /* skip if i is dead */
361  h = last[i]; /* scan hash bucket of node i */
362  i = hhead[h];
363  hhead[h] = -1; /* hash bucket will be empty */
364  for(; i != -1 && next[i] != -1; i = next[i], mark++)
365  {
366  ln = len[i];
367  eln = elen[i];
368  for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
369  jlast = i;
370  for(j = next[i]; j != -1; ) /* compare i with all j */
371  {
372  ok = (len[j] == ln) && (elen[j] == eln);
373  for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
374  {
375  if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/
376  }
377  if(ok) /* i and j are identical */
378  {
379  Cp[j] = amd_flip (i); /* absorb j into i */
380  nv[i] += nv[j];
381  nv[j] = 0;
382  elen[j] = -1; /* node j is dead */
383  j = next[j]; /* delete j from hash bucket */
384  next[jlast] = j;
385  }
386  else
387  {
388  jlast = j; /* j and i are different */
389  j = next[j];
390  }
391  }
392  }
393  }
394 
395  /* --- Finalize new element------------------------------------------ */
396  for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */
397  {
398  i = Ci[pk];
399  if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
400  nv[i] = nvi; /* restore nv[i] */
401  d = degree[i] + dk - nvi; /* compute external degree(i) */
402  d = std::min<Index> (d, n - nel - nvi);
403  if(head[d] != -1) last[head[d]] = i;
404  next[i] = head[d]; /* put i back in degree list */
405  last[i] = -1;
406  head[d] = i;
407  mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */
408  degree[i] = d;
409  Ci[p++] = i; /* place i in Lk */
410  }
411  nv[k] = nvk; /* # nodes absorbed into k */
412  if((len[k] = p-pk1) == 0) /* length of adj list of element k*/
413  {
414  Cp[k] = -1; /* k is a root of the tree */
415  w[k] = 0; /* k is now a dead element */
416  }
417  if(elenk != 0) cnz = p; /* free unused space in Lk */
418  }
419 
420  /* --- Postordering ----------------------------------------------------- */
421  for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */
422  for(j = 0; j <= n; j++) head[j] = -1;
423  for(j = n; j >= 0; j--) /* place unordered nodes in lists */
424  {
425  if(nv[j] > 0) continue; /* skip if j is an element */
426  next[j] = head[Cp[j]]; /* place j in list of its parent */
427  head[Cp[j]] = j;
428  }
429  for(e = n; e >= 0; e--) /* place elements in lists */
430  {
431  if(nv[e] <= 0) continue; /* skip unless e is an element */
432  if(Cp[e] != -1)
433  {
434  next[e] = head[Cp[e]]; /* place e in list of its parent */
435  head[Cp[e]] = e;
436  }
437  }
438  for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */
439  {
440  if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w);
441  }
442 
443  perm.indices().conservativeResize(n);
444 
445  delete[] W;
446 }
447 
448 } // namespace internal
449 
450 } // end namespace Eigen
451 
452 #endif // EIGEN_SPARSE_AMD_H