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