• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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