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 internal {
00052   
00053 
00054 #define CS_FLIP(i) (-(i)-2)
00055 #define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
00056 #define CS_MARKED(w,j) (w[j] < 0)
00057 #define CS_MARK(w,j) { w[j] = CS_FLIP (w[j]); }
00058 
00059 /* clear w */
00060 template<typename Index>
00061 static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
00062 {
00063   Index k;
00064   if(mark < 2 || (mark + lemax < 0))
00065   {
00066     for(k = 0; k < n; k++)
00067       if(w[k] != 0)
00068         w[k] = 1;
00069     mark = 2;
00070   }
00071   return (mark);     /* at this point, w[0..n-1] < mark holds */
00072 }
00073 
00074 /* depth-first search and postorder of a tree rooted at node j */
00075 template<typename Index>
00076 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
00077 {
00078   int i, p, top = 0;
00079   if(!head || !next || !post || !stack) return (-1);    /* check inputs */
00080   stack[0] = j;                 /* place j on the stack */
00081   while (top >= 0)                /* while (stack is not empty) */
00082   {
00083     p = stack[top];           /* p = top of stack */
00084     i = head[p];              /* i = youngest child of p */
00085     if(i == -1)
00086     {
00087       top--;                 /* p has no unordered children left */
00088       post[k++] = p;        /* node p is the kth postordered node */
00089     }
00090     else
00091     {
00092       head[p] = next[i];   /* remove i from children of p */
00093       stack[++top] = i;     /* start dfs on child node i */
00094     }
00095   }
00096   return k;
00097 }
00098 
00099 
00105 template<typename Scalar, typename Index>
00106 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic>& perm)
00107 {
00108   typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
00109   
00110   int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
00111       k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
00112       ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
00113   unsigned int h;
00114   
00115   Index n = C.cols();
00116   dense = std::max<Index> (16, 10 * sqrt ((double) n));   /* find dense threshold */
00117   dense = std::min<Index> (n-2, dense);
00118   
00119   Index cnz = C.nonZeros();
00120   perm.resize(n+1);
00121   t = cnz + cnz/5 + 2*n;                 /* add elbow room to C */
00122   C.resizeNonZeros(t);
00123   
00124   Index* W       = new Index[8*(n+1)]; /* get workspace */
00125   Index* len     = W;
00126   Index* nv      = W +   (n+1);
00127   Index* next    = W + 2*(n+1);
00128   Index* head    = W + 3*(n+1);
00129   Index* elen    = W + 4*(n+1);
00130   Index* degree  = W + 5*(n+1);
00131   Index* w       = W + 6*(n+1);
00132   Index* hhead   = W + 7*(n+1);
00133   Index* last    = perm.indices().data();                              /* use P as workspace for last */
00134   
00135   /* --- Initialize quotient graph ---------------------------------------- */
00136   Index* Cp = C._outerIndexPtr();
00137   Index* Ci = C._innerIndexPtr();
00138   for(k = 0; k < n; k++)
00139     len[k] = Cp[k+1] - Cp[k];
00140   len[n] = 0;
00141   nzmax = t;
00142   
00143   for(i = 0; i <= n; i++)
00144   {
00145     head[i]   = -1;                     // degree list i is empty
00146     last[i]   = -1;
00147     next[i]   = -1;
00148     hhead[i]  = -1;                     // hash list i is empty 
00149     nv[i]     = 1;                      // node i is just one node
00150     w[i]      = 1;                      // node i is alive
00151     elen[i]   = 0;                      // Ek of node i is empty
00152     degree[i] = len[i];                 // degree of node i
00153   }
00154   mark = cs_wclear (0, 0, w, n);         /* clear w */
00155   elen[n] = -2;                         /* n is a dead element */
00156   Cp[n] = -1;                           /* n is a root of assembly tree */
00157   w[n] = 0;                             /* n is a dead element */
00158   
00159   /* --- Initialize degree lists ------------------------------------------ */
00160   for(i = 0; i < n; i++)
00161   {
00162     d = degree[i];
00163     if(d == 0)                         /* node i is empty */
00164     {
00165       elen[i] = -2;                 /* element i is dead */
00166       nel++;
00167       Cp[i] = -1;                   /* i is a root of assembly tree */
00168       w[i] = 0;
00169     }
00170     else if(d > dense)                 /* node i is dense */
00171     {
00172       nv[i] = 0;                    /* absorb i into element n */
00173       elen[i] = -1;                 /* node i is dead */
00174       nel++;
00175       Cp[i] = CS_FLIP (n);
00176       nv[n]++;
00177     }
00178     else
00179     {
00180       if(head[d] != -1) last[head[d]] = i;
00181       next[i] = head[d];           /* put node i in degree list d */
00182       head[d] = i;
00183     }
00184   }
00185   
00186   while (nel < n)                         /* while (selecting pivots) do */
00187   {
00188     /* --- Select node of minimum approximate degree -------------------- */
00189     for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
00190     if(next[k] != -1) last[next[k]] = -1;
00191     head[mindeg] = next[k];          /* remove k from degree list */
00192     elenk = elen[k];                  /* elenk = |Ek| */
00193     nvk = nv[k];                      /* # of nodes k represents */
00194     nel += nvk;                        /* nv[k] nodes of A eliminated */
00195     
00196     /* --- Garbage collection ------------------------------------------- */
00197     if(elenk > 0 && cnz + mindeg >= nzmax)
00198     {
00199       for(j = 0; j < n; j++)
00200       {
00201         if((p = Cp[j]) >= 0)      /* j is a live node or element */
00202         {
00203           Cp[j] = Ci[p];          /* save first entry of object */
00204           Ci[p] = CS_FLIP (j);    /* first entry is now CS_FLIP(j) */
00205         }
00206       }
00207       for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
00208       {
00209         if((j = CS_FLIP (Ci[p++])) >= 0)  /* found object j */
00210         {
00211           Ci[q] = Cp[j];       /* restore first entry of object */
00212           Cp[j] = q++;          /* new pointer to object j */
00213           for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
00214         }
00215       }
00216       cnz = q;                       /* Ci[cnz...nzmax-1] now free */
00217     }
00218     
00219     /* --- Construct new element ---------------------------------------- */
00220     dk = 0;
00221     nv[k] = -nvk;                     /* flag k as in Lk */
00222     p = Cp[k];
00223     pk1 = (elenk == 0) ? p : cnz;      /* do in place if elen[k] == 0 */
00224     pk2 = pk1;
00225     for(k1 = 1; k1 <= elenk + 1; k1++)
00226     {
00227       if(k1 > elenk)
00228       {
00229         e = k;                     /* search the nodes in k */
00230         pj = p;                    /* list of nodes starts at Ci[pj]*/
00231         ln = len[k] - elenk;      /* length of list of nodes in k */
00232       }
00233       else
00234       {
00235         e = Ci[p++];              /* search the nodes in e */
00236         pj = Cp[e];
00237         ln = len[e];              /* length of list of nodes in e */
00238       }
00239       for(k2 = 1; k2 <= ln; k2++)
00240       {
00241         i = Ci[pj++];
00242         if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
00243         dk += nvi;                 /* degree[Lk] += size of node i */
00244         nv[i] = -nvi;             /* negate nv[i] to denote i in Lk*/
00245         Ci[pk2++] = i;            /* place i in Lk */
00246         if(next[i] != -1) last[next[i]] = last[i];
00247         if(last[i] != -1)         /* remove i from degree list */
00248         {
00249           next[last[i]] = next[i];
00250         }
00251         else
00252         {
00253           head[degree[i]] = next[i];
00254         }
00255       }
00256       if(e != k)
00257       {
00258         Cp[e] = CS_FLIP (k);      /* absorb e into k */
00259         w[e] = 0;                 /* e is now a dead element */
00260       }
00261     }
00262     if(elenk != 0) cnz = pk2;         /* Ci[cnz...nzmax] is free */
00263     degree[k] = dk;                   /* external degree of k - |Lk\i| */
00264     Cp[k] = pk1;                      /* element k is in Ci[pk1..pk2-1] */
00265     len[k] = pk2 - pk1;
00266     elen[k] = -2;                     /* k is now an element */
00267     
00268     /* --- Find set differences ----------------------------------------- */
00269     mark = cs_wclear (mark, lemax, w, n);  /* clear w if necessary */
00270     for(pk = pk1; pk < pk2; pk++)    /* scan 1: find |Le\Lk| */
00271     {
00272       i = Ci[pk];
00273       if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
00274       nvi = -nv[i];                      /* nv[i] was negated */
00275       wnvi = mark - nvi;
00276       for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)  /* scan Ei */
00277       {
00278         e = Ci[p];
00279         if(w[e] >= mark)
00280         {
00281           w[e] -= nvi;          /* decrement |Le\Lk| */
00282         }
00283         else if(w[e] != 0)        /* ensure e is a live element */
00284         {
00285           w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
00286         }
00287       }
00288     }
00289     
00290     /* --- Degree update ------------------------------------------------ */
00291     for(pk = pk1; pk < pk2; pk++)    /* scan2: degree update */
00292     {
00293       i = Ci[pk];                   /* consider node i in Lk */
00294       p1 = Cp[i];
00295       p2 = p1 + elen[i] - 1;
00296       pn = p1;
00297       for(h = 0, d = 0, p = p1; p <= p2; p++)    /* scan Ei */
00298       {
00299         e = Ci[p];
00300         if(w[e] != 0)             /* e is an unabsorbed element */
00301         {
00302           dext = w[e] - mark;   /* dext = |Le\Lk| */
00303           if(dext > 0)
00304           {
00305             d += dext;         /* sum up the set differences */
00306             Ci[pn++] = e;     /* keep e in Ei */
00307             h += e;            /* compute the hash of node i */
00308           }
00309           else
00310           {
00311             Cp[e] = CS_FLIP (k);  /* aggressive absorb. e->k */
00312             w[e] = 0;             /* e is a dead element */
00313           }
00314         }
00315       }
00316       elen[i] = pn - p1 + 1;        /* elen[i] = |Ei| */
00317       p3 = pn;
00318       p4 = p1 + len[i];
00319       for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
00320       {
00321         j = Ci[p];
00322         if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
00323         d += nvj;                  /* degree(i) += |j| */
00324         Ci[pn++] = j;             /* place j in node list of i */
00325         h += j;                    /* compute hash for node i */
00326       }
00327       if(d == 0)                     /* check for mass elimination */
00328       {
00329         Cp[i] = CS_FLIP (k);      /* absorb i into k */
00330         nvi = -nv[i];
00331         dk -= nvi;                 /* |Lk| -= |i| */
00332         nvk += nvi;                /* |k| += nv[i] */
00333         nel += nvi;
00334         nv[i] = 0;
00335         elen[i] = -1;             /* node i is dead */
00336       }
00337       else
00338       {
00339         degree[i] = std::min<Index> (degree[i], d);   /* update degree(i) */
00340         Ci[pn] = Ci[p3];         /* move first node to end */
00341         Ci[p3] = Ci[p1];         /* move 1st el. to end of Ei */
00342         Ci[p1] = k;               /* add k as 1st element in of Ei */
00343         len[i] = pn - p1 + 1;     /* new len of adj. list of node i */
00344         h %= n;                    /* finalize hash of i */
00345         next[i] = hhead[h];      /* place i in hash bucket */
00346         hhead[h] = i;
00347         last[i] = h;              /* save hash of i in last[i] */
00348       }
00349     }                                   /* scan2 is done */
00350     degree[k] = dk;                   /* finalize |Lk| */
00351     lemax = std::max<Index>(lemax, dk);
00352     mark = cs_wclear (mark+lemax, lemax, w, n);    /* clear w */
00353     
00354     /* --- Supernode detection ------------------------------------------ */
00355     for(pk = pk1; pk < pk2; pk++)
00356     {
00357       i = Ci[pk];
00358       if(nv[i] >= 0) continue;         /* skip if i is dead */
00359       h = last[i];                      /* scan hash bucket of node i */
00360       i = hhead[h];
00361       hhead[h] = -1;                    /* hash bucket will be empty */
00362       for(; i != -1 && next[i] != -1; i = next[i], mark++)
00363       {
00364         ln = len[i];
00365         eln = elen[i];
00366         for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
00367         jlast = i;
00368         for(j = next[i]; j != -1; ) /* compare i with all j */
00369         {
00370           ok = (len[j] == ln) && (elen[j] == eln);
00371           for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
00372           {
00373             if(w[Ci[p]] != mark) ok = 0;    /* compare i and j*/
00374           }
00375           if(ok)                     /* i and j are identical */
00376           {
00377             Cp[j] = CS_FLIP (i);  /* absorb j into i */
00378             nv[i] += nv[j];
00379             nv[j] = 0;
00380             elen[j] = -1;         /* node j is dead */
00381             j = next[j];          /* delete j from hash bucket */
00382             next[jlast] = j;
00383           }
00384           else
00385           {
00386             jlast = j;             /* j and i are different */
00387             j = next[j];
00388           }
00389         }
00390       }
00391     }
00392     
00393     /* --- Finalize new element------------------------------------------ */
00394     for(p = pk1, pk = pk1; pk < pk2; pk++)   /* finalize Lk */
00395     {
00396       i = Ci[pk];
00397       if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
00398       nv[i] = nvi;                      /* restore nv[i] */
00399       d = degree[i] + dk - nvi;         /* compute external degree(i) */
00400       d = std::min<Index> (d, n - nel - nvi);
00401       if(head[d] != -1) last[head[d]] = i;
00402       next[i] = head[d];               /* put i back in degree list */
00403       last[i] = -1;
00404       head[d] = i;
00405       mindeg = std::min<Index> (mindeg, d);       /* find new minimum degree */
00406       degree[i] = d;
00407       Ci[p++] = i;                      /* place i in Lk */
00408     }
00409     nv[k] = nvk;                      /* # nodes absorbed into k */
00410     if((len[k] = p-pk1) == 0)         /* length of adj list of element k*/
00411     {
00412       Cp[k] = -1;                   /* k is a root of the tree */
00413       w[k] = 0;                     /* k is now a dead element */
00414     }
00415     if(elenk != 0) cnz = p;           /* free unused space in Lk */
00416   }
00417   
00418   /* --- Postordering ----------------------------------------------------- */
00419   for(i = 0; i < n; i++) Cp[i] = CS_FLIP (Cp[i]);/* fix assembly tree */
00420   for(j = 0; j <= n; j++) head[j] = -1;
00421   for(j = n; j >= 0; j--)              /* place unordered nodes in lists */
00422   {
00423     if(nv[j] > 0) continue;          /* skip if j is an element */
00424     next[j] = head[Cp[j]];          /* place j in list of its parent */
00425     head[Cp[j]] = j;
00426   }
00427   for(e = n; e >= 0; e--)              /* place elements in lists */
00428   {
00429     if(nv[e] <= 0) continue;         /* skip unless e is an element */
00430     if(Cp[e] != -1)
00431     {
00432       next[e] = head[Cp[e]];      /* place e in list of its parent */
00433       head[Cp[e]] = e;
00434     }
00435   }
00436   for(k = 0, i = 0; i <= n; i++)       /* postorder the assembly tree */
00437   {
00438     if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, perm.indices().data(), w);
00439   }
00440   
00441   perm.indices().conservativeResize(n);
00442 
00443   delete[] W;
00444 }
00445 
00446 } // namespace internal
00447 
00448 #endif // EIGEN_SPARSE_AMD_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:28