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
amd_flip(const T & i)38 template<typename T> inline T amd_flip(const T& i) { return -i-2; }
amd_unflip(const T & i)39 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; }
amd_marked(const T0 * w,const T1 & j)40 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; }
amd_mark(const T0 * w,const T1 & j)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>
cs_wclear(Index mark,Index lemax,Index * w,Index n)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>
cs_tdfs(Index j,Index k,Index * head,const Index * next,Index * post,Index * stack)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
84 /** \internal
85 * \ingroup OrderingMethods_Module
86 * Approximate minimum degree ordering algorithm.
87 * \returns the permutation P reducing the fill-in of the input matrix \a C
88 * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.
89 * On exit the values of C are destroyed */
90 template<typename Scalar, typename Index>
minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index> & C,PermutationMatrix<Dynamic,Dynamic,Index> & perm)91 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm)
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
436