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


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