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
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 StorageIndex>
cs_wclear(StorageIndex mark,StorageIndex lemax,StorageIndex * w,StorageIndex n)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>
cs_tdfs(StorageIndex j,StorageIndex k,StorageIndex * head,const StorageIndex * next,StorageIndex * post,StorageIndex * stack)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
84 /** \internal
85 * \ingroup OrderingMethods_Module
86 * Approximate minimum degree ordering algorithm.
87 *
88 * \param[in] C the input selfadjoint matrix stored in compressed column major format.
89 * \param[out] perm the permutation P reducing the fill-in of the input matrix \a C
90 *
91 * Note that the input matrix \a C must be complete, that is both the upper and lower parts have to be stored, as well as the diagonal entries.
92 * On exit the values of C are destroyed */
93 template<typename Scalar, typename StorageIndex>
minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,StorageIndex> & C,PermutationMatrix<Dynamic,Dynamic,StorageIndex> & perm)94 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,StorageIndex>& C, PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm)
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
446