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