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


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:07