Amd.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 /*
00026 
00027 NOTE: this routine has been adapted from the CSparse library:
00028 
00029 Copyright (c) 2006, Timothy A. Davis.
00030 http://www.cise.ufl.edu/research/sparse/CSparse
00031 
00032 CSparse is free software; you can redistribute it and/or
00033 modify it under the terms of the GNU Lesser General Public
00034 License as published by the Free Software Foundation; either
00035 version 2.1 of the License, or (at your option) any later version.
00036 
00037 CSparse is distributed in the hope that it will be useful,
00038 but WITHOUT ANY WARRANTY; without even the implied warranty of
00039 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00040 Lesser General Public License for more details.
00041 
00042 You should have received a copy of the GNU Lesser General Public
00043 License along with this Module; if not, write to the Free Software
00044 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00045 
00046 */
00047 
00048 #ifndef EIGEN_SPARSE_AMD_H
00049 #define EIGEN_SPARSE_AMD_H
00050 
00051 namespace Eigen { 
00052 
00053 namespace internal {
00054   
00055 template<typename T> inline T amd_flip(const T& i) { return -i-2; }
00056 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; }
00057 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; }
00058 template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); }
00059 
00060 /* clear w */
00061 template<typename Index>
00062 static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
00063 {
00064   Index k;
00065   if(mark < 2 || (mark + lemax < 0))
00066   {
00067     for(k = 0; k < n; k++)
00068       if(w[k] != 0)
00069         w[k] = 1;
00070     mark = 2;
00071   }
00072   return (mark);     /* at this point, w[0..n-1] < mark holds */
00073 }
00074 
00075 /* depth-first search and postorder of a tree rooted at node j */
00076 template<typename Index>
00077 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
00078 {
00079   int i, p, top = 0;
00080   if(!head || !next || !post || !stack) return (-1);    /* check inputs */
00081   stack[0] = j;                 /* place j on the stack */
00082   while (top >= 0)                /* while (stack is not empty) */
00083   {
00084     p = stack[top];           /* p = top of stack */
00085     i = head[p];              /* i = youngest child of p */
00086     if(i == -1)
00087     {
00088       top--;                 /* p has no unordered children left */
00089       post[k++] = p;        /* node p is the kth postordered node */
00090     }
00091     else
00092     {
00093       head[p] = next[i];   /* remove i from children of p */
00094       stack[++top] = i;     /* start dfs on child node i */
00095     }
00096   }
00097   return k;
00098 }
00099 
00100 
00106 template<typename Scalar, typename Index>
00107 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm)
00108 {
00109   typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
00110   
00111   int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
00112       k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
00113       ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
00114   unsigned int h;
00115   
00116   Index n = C.cols();
00117   dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
00118   dense = std::min<Index> (n-2, dense);
00119   
00120   Index cnz = C.nonZeros();
00121   perm.resize(n+1);
00122   t = cnz + cnz/5 + 2*n;                 /* add elbow room to C */
00123   C.resizeNonZeros(t);
00124   
00125   Index* W       = new Index[8*(n+1)]; /* get workspace */
00126   Index* len     = W;
00127   Index* nv      = W +   (n+1);
00128   Index* next    = W + 2*(n+1);
00129   Index* head    = W + 3*(n+1);
00130   Index* elen    = W + 4*(n+1);
00131   Index* degree  = W + 5*(n+1);
00132   Index* w       = W + 6*(n+1);
00133   Index* hhead   = W + 7*(n+1);
00134   Index* last    = perm.indices().data();                              /* use P as workspace for last */
00135   
00136   /* --- Initialize quotient graph ---------------------------------------- */
00137   Index* Cp = C.outerIndexPtr();
00138   Index* Ci = C.innerIndexPtr();
00139   for(k = 0; k < n; k++)
00140     len[k] = Cp[k+1] - Cp[k];
00141   len[n] = 0;
00142   nzmax = t;
00143   
00144   for(i = 0; i <= n; i++)
00145   {
00146     head[i]   = -1;                     // degree list i is empty
00147     last[i]   = -1;
00148     next[i]   = -1;
00149     hhead[i]  = -1;                     // hash list i is empty 
00150     nv[i]     = 1;                      // node i is just one node
00151     w[i]      = 1;                      // node i is alive
00152     elen[i]   = 0;                      // Ek of node i is empty
00153     degree[i] = len[i];                 // degree of node i
00154   }
00155   mark = internal::cs_wclear<Index>(0, 0, w, n);         /* clear w */
00156   elen[n] = -2;                         /* n is a dead element */
00157   Cp[n] = -1;                           /* n is a root of assembly tree */
00158   w[n] = 0;                             /* n is a dead element */
00159   
00160   /* --- Initialize degree lists ------------------------------------------ */
00161   for(i = 0; i < n; i++)
00162   {
00163     d = degree[i];
00164     if(d == 0)                         /* node i is empty */
00165     {
00166       elen[i] = -2;                 /* element i is dead */
00167       nel++;
00168       Cp[i] = -1;                   /* i is a root of assembly tree */
00169       w[i] = 0;
00170     }
00171     else if(d > dense)                 /* node i is dense */
00172     {
00173       nv[i] = 0;                    /* absorb i into element n */
00174       elen[i] = -1;                 /* node i is dead */
00175       nel++;
00176       Cp[i] = amd_flip (n);
00177       nv[n]++;
00178     }
00179     else
00180     {
00181       if(head[d] != -1) last[head[d]] = i;
00182       next[i] = head[d];           /* put node i in degree list d */
00183       head[d] = i;
00184     }
00185   }
00186   
00187   while (nel < n)                         /* while (selecting pivots) do */
00188   {
00189     /* --- Select node of minimum approximate degree -------------------- */
00190     for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
00191     if(next[k] != -1) last[next[k]] = -1;
00192     head[mindeg] = next[k];          /* remove k from degree list */
00193     elenk = elen[k];                  /* elenk = |Ek| */
00194     nvk = nv[k];                      /* # of nodes k represents */
00195     nel += nvk;                        /* nv[k] nodes of A eliminated */
00196     
00197     /* --- Garbage collection ------------------------------------------- */
00198     if(elenk > 0 && cnz + mindeg >= nzmax)
00199     {
00200       for(j = 0; j < n; j++)
00201       {
00202         if((p = Cp[j]) >= 0)      /* j is a live node or element */
00203         {
00204           Cp[j] = Ci[p];          /* save first entry of object */
00205           Ci[p] = amd_flip (j);    /* first entry is now amd_flip(j) */
00206         }
00207       }
00208       for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
00209       {
00210         if((j = amd_flip (Ci[p++])) >= 0)  /* found object j */
00211         {
00212           Ci[q] = Cp[j];       /* restore first entry of object */
00213           Cp[j] = q++;          /* new pointer to object j */
00214           for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
00215         }
00216       }
00217       cnz = q;                       /* Ci[cnz...nzmax-1] now free */
00218     }
00219     
00220     /* --- Construct new element ---------------------------------------- */
00221     dk = 0;
00222     nv[k] = -nvk;                     /* flag k as in Lk */
00223     p = Cp[k];
00224     pk1 = (elenk == 0) ? p : cnz;      /* do in place if elen[k] == 0 */
00225     pk2 = pk1;
00226     for(k1 = 1; k1 <= elenk + 1; k1++)
00227     {
00228       if(k1 > elenk)
00229       {
00230         e = k;                     /* search the nodes in k */
00231         pj = p;                    /* list of nodes starts at Ci[pj]*/
00232         ln = len[k] - elenk;      /* length of list of nodes in k */
00233       }
00234       else
00235       {
00236         e = Ci[p++];              /* search the nodes in e */
00237         pj = Cp[e];
00238         ln = len[e];              /* length of list of nodes in e */
00239       }
00240       for(k2 = 1; k2 <= ln; k2++)
00241       {
00242         i = Ci[pj++];
00243         if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
00244         dk += nvi;                 /* degree[Lk] += size of node i */
00245         nv[i] = -nvi;             /* negate nv[i] to denote i in Lk*/
00246         Ci[pk2++] = i;            /* place i in Lk */
00247         if(next[i] != -1) last[next[i]] = last[i];
00248         if(last[i] != -1)         /* remove i from degree list */
00249         {
00250           next[last[i]] = next[i];
00251         }
00252         else
00253         {
00254           head[degree[i]] = next[i];
00255         }
00256       }
00257       if(e != k)
00258       {
00259         Cp[e] = amd_flip (k);      /* absorb e into k */
00260         w[e] = 0;                 /* e is now a dead element */
00261       }
00262     }
00263     if(elenk != 0) cnz = pk2;         /* Ci[cnz...nzmax] is free */
00264     degree[k] = dk;                   /* external degree of k - |Lk\i| */
00265     Cp[k] = pk1;                      /* element k is in Ci[pk1..pk2-1] */
00266     len[k] = pk2 - pk1;
00267     elen[k] = -2;                     /* k is now an element */
00268     
00269     /* --- Find set differences ----------------------------------------- */
00270     mark = internal::cs_wclear<Index>(mark, lemax, w, n);  /* clear w if necessary */
00271     for(pk = pk1; pk < pk2; pk++)    /* scan 1: find |Le\Lk| */
00272     {
00273       i = Ci[pk];
00274       if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
00275       nvi = -nv[i];                      /* nv[i] was negated */
00276       wnvi = mark - nvi;
00277       for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)  /* scan Ei */
00278       {
00279         e = Ci[p];
00280         if(w[e] >= mark)
00281         {
00282           w[e] -= nvi;          /* decrement |Le\Lk| */
00283         }
00284         else if(w[e] != 0)        /* ensure e is a live element */
00285         {
00286           w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
00287         }
00288       }
00289     }
00290     
00291     /* --- Degree update ------------------------------------------------ */
00292     for(pk = pk1; pk < pk2; pk++)    /* scan2: degree update */
00293     {
00294       i = Ci[pk];                   /* consider node i in Lk */
00295       p1 = Cp[i];
00296       p2 = p1 + elen[i] - 1;
00297       pn = p1;
00298       for(h = 0, d = 0, p = p1; p <= p2; p++)    /* scan Ei */
00299       {
00300         e = Ci[p];
00301         if(w[e] != 0)             /* e is an unabsorbed element */
00302         {
00303           dext = w[e] - mark;   /* dext = |Le\Lk| */
00304           if(dext > 0)
00305           {
00306             d += dext;         /* sum up the set differences */
00307             Ci[pn++] = e;     /* keep e in Ei */
00308             h += e;            /* compute the hash of node i */
00309           }
00310           else
00311           {
00312             Cp[e] = amd_flip (k);  /* aggressive absorb. e->k */
00313             w[e] = 0;             /* e is a dead element */
00314           }
00315         }
00316       }
00317       elen[i] = pn - p1 + 1;        /* elen[i] = |Ei| */
00318       p3 = pn;
00319       p4 = p1 + len[i];
00320       for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
00321       {
00322         j = Ci[p];
00323         if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
00324         d += nvj;                  /* degree(i) += |j| */
00325         Ci[pn++] = j;             /* place j in node list of i */
00326         h += j;                    /* compute hash for node i */
00327       }
00328       if(d == 0)                     /* check for mass elimination */
00329       {
00330         Cp[i] = amd_flip (k);      /* absorb i into k */
00331         nvi = -nv[i];
00332         dk -= nvi;                 /* |Lk| -= |i| */
00333         nvk += nvi;                /* |k| += nv[i] */
00334         nel += nvi;
00335         nv[i] = 0;
00336         elen[i] = -1;             /* node i is dead */
00337       }
00338       else
00339       {
00340         degree[i] = std::min<Index> (degree[i], d);   /* update degree(i) */
00341         Ci[pn] = Ci[p3];         /* move first node to end */
00342         Ci[p3] = Ci[p1];         /* move 1st el. to end of Ei */
00343         Ci[p1] = k;               /* add k as 1st element in of Ei */
00344         len[i] = pn - p1 + 1;     /* new len of adj. list of node i */
00345         h %= n;                    /* finalize hash of i */
00346         next[i] = hhead[h];      /* place i in hash bucket */
00347         hhead[h] = i;
00348         last[i] = h;              /* save hash of i in last[i] */
00349       }
00350     }                                   /* scan2 is done */
00351     degree[k] = dk;                   /* finalize |Lk| */
00352     lemax = std::max<Index>(lemax, dk);
00353     mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n);    /* clear w */
00354     
00355     /* --- Supernode detection ------------------------------------------ */
00356     for(pk = pk1; pk < pk2; pk++)
00357     {
00358       i = Ci[pk];
00359       if(nv[i] >= 0) continue;         /* skip if i is dead */
00360       h = last[i];                      /* scan hash bucket of node i */
00361       i = hhead[h];
00362       hhead[h] = -1;                    /* hash bucket will be empty */
00363       for(; i != -1 && next[i] != -1; i = next[i], mark++)
00364       {
00365         ln = len[i];
00366         eln = elen[i];
00367         for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
00368         jlast = i;
00369         for(j = next[i]; j != -1; ) /* compare i with all j */
00370         {
00371           ok = (len[j] == ln) && (elen[j] == eln);
00372           for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
00373           {
00374             if(w[Ci[p]] != mark) ok = 0;    /* compare i and j*/
00375           }
00376           if(ok)                     /* i and j are identical */
00377           {
00378             Cp[j] = amd_flip (i);  /* absorb j into i */
00379             nv[i] += nv[j];
00380             nv[j] = 0;
00381             elen[j] = -1;         /* node j is dead */
00382             j = next[j];          /* delete j from hash bucket */
00383             next[jlast] = j;
00384           }
00385           else
00386           {
00387             jlast = j;             /* j and i are different */
00388             j = next[j];
00389           }
00390         }
00391       }
00392     }
00393     
00394     /* --- Finalize new element------------------------------------------ */
00395     for(p = pk1, pk = pk1; pk < pk2; pk++)   /* finalize Lk */
00396     {
00397       i = Ci[pk];
00398       if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
00399       nv[i] = nvi;                      /* restore nv[i] */
00400       d = degree[i] + dk - nvi;         /* compute external degree(i) */
00401       d = std::min<Index> (d, n - nel - nvi);
00402       if(head[d] != -1) last[head[d]] = i;
00403       next[i] = head[d];               /* put i back in degree list */
00404       last[i] = -1;
00405       head[d] = i;
00406       mindeg = std::min<Index> (mindeg, d);       /* find new minimum degree */
00407       degree[i] = d;
00408       Ci[p++] = i;                      /* place i in Lk */
00409     }
00410     nv[k] = nvk;                      /* # nodes absorbed into k */
00411     if((len[k] = p-pk1) == 0)         /* length of adj list of element k*/
00412     {
00413       Cp[k] = -1;                   /* k is a root of the tree */
00414       w[k] = 0;                     /* k is now a dead element */
00415     }
00416     if(elenk != 0) cnz = p;           /* free unused space in Lk */
00417   }
00418   
00419   /* --- Postordering ----------------------------------------------------- */
00420   for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */
00421   for(j = 0; j <= n; j++) head[j] = -1;
00422   for(j = n; j >= 0; j--)              /* place unordered nodes in lists */
00423   {
00424     if(nv[j] > 0) continue;          /* skip if j is an element */
00425     next[j] = head[Cp[j]];          /* place j in list of its parent */
00426     head[Cp[j]] = j;
00427   }
00428   for(e = n; e >= 0; e--)              /* place elements in lists */
00429   {
00430     if(nv[e] <= 0) continue;         /* skip unless e is an element */
00431     if(Cp[e] != -1)
00432     {
00433       next[e] = head[Cp[e]];      /* place e in list of its parent */
00434       head[Cp[e]] = e;
00435     }
00436   }
00437   for(k = 0, i = 0; i <= n; i++)       /* postorder the assembly tree */
00438   {
00439     if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w);
00440   }
00441   
00442   perm.indices().conservativeResize(n);
00443 
00444   delete[] W;
00445 }
00446 
00447 } // namespace internal
00448 
00449 } // end namespace Eigen
00450 
00451 #endif // EIGEN_SPARSE_AMD_H