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


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:41:37