src/OrderingMethods/Amd.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 /*
0011 NOTE: this routine has been adapted from the CSparse library:
0012
0013 Copyright (c) 2006, Timothy A. Davis.
0014 http://www.suitesparse.com
0015
0016 The author of CSparse, Timothy A. Davis., has executed a license with Google LLC
0017 to permit distribution of this code and derivative works as part of Eigen under
0018 the Mozilla Public License v. 2.0, as stated at the top of this file.
0019 */
0020
0021 #ifndef EIGEN_SPARSE_AMD_H
0022 #define EIGEN_SPARSE_AMD_H
0023
0024 namespace Eigen {
0025
0026 namespace internal {
0027
0028 template<typename T> inline T amd_flip(const T& i) { return -i-2; }
0029 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; }
0030 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; }
0031 template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); }
0032
0033 /* clear w */
0034 template<typename StorageIndex>
0035 static StorageIndex cs_wclear (StorageIndex mark, StorageIndex lemax, StorageIndex *w, StorageIndex n)
0036 {
0037   StorageIndex k;
0038   if(mark < 2 || (mark + lemax < 0))
0039   {
0040     for(k = 0; k < n; k++)
0041       if(w[k] != 0)
0042         w[k] = 1;
0043     mark = 2;
0044   }
0045   return (mark);     /* at this point, w[0..n-1] < mark holds */
0046 }
0047
0048 /* depth-first search and postorder of a tree rooted at node j */
0049 template<typename StorageIndex>
0050 StorageIndex cs_tdfs(StorageIndex j, StorageIndex k, StorageIndex *head, const StorageIndex *next, StorageIndex *post, StorageIndex *stack)
0051 {
0052   StorageIndex i, p, top = 0;
0053   if(!head || !next || !post || !stack) return (-1);    /* check inputs */
0054   stack[0] = j;                 /* place j on the stack */
0055   while (top >= 0)                /* while (stack is not empty) */
0056   {
0057     p = stack[top];           /* p = top of stack */
0058     i = head[p];              /* i = youngest child of p */
0059     if(i == -1)
0060     {
0061       top--;                 /* p has no unordered children left */
0062       post[k++] = p;        /* node p is the kth postordered node */
0063     }
0064     else
0065     {
0066       head[p] = next[i];   /* remove i from children of p */
0067       stack[++top] = i;     /* start dfs on child node i */
0068     }
0069   }
0070   return k;
0071 }
0072
0073
0074 /** \internal
0075   * \ingroup OrderingMethods_Module
0076   * Approximate minimum degree ordering algorithm.
0077   *
0078   * \param[in] C the input selfadjoint matrix stored in compressed column major format.
0079   * \param[out] perm the permutation P reducing the fill-in of the input matrix \a C
0080   *
0081   * 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.
0082   * On exit the values of C are destroyed */
0083 template<typename Scalar, typename StorageIndex>
0084 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,StorageIndex>& C, PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm)
0085 {
0086   using std::sqrt;
0087
0088   StorageIndex d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
0089                 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
0090                 ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t, h;
0091
0092   StorageIndex n = StorageIndex(C.cols());
0093   dense = std::max<StorageIndex> (16, StorageIndex(10 * sqrt(double(n))));   /* find dense threshold */
0094   dense = (std::min)(n-2, dense);
0095
0096   StorageIndex cnz = StorageIndex(C.nonZeros());
0097   perm.resize(n+1);
0098   t = cnz + cnz/5 + 2*n;                 /* add elbow room to C */
0099   C.resizeNonZeros(t);
0100
0101   // get workspace
0102   ei_declare_aligned_stack_constructed_variable(StorageIndex,W,8*(n+1),0);
0103   StorageIndex* len     = W;
0104   StorageIndex* nv      = W +   (n+1);
0105   StorageIndex* next    = W + 2*(n+1);
0106   StorageIndex* head    = W + 3*(n+1);
0107   StorageIndex* elen    = W + 4*(n+1);
0108   StorageIndex* degree  = W + 5*(n+1);
0109   StorageIndex* w       = W + 6*(n+1);
0110   StorageIndex* hhead   = W + 7*(n+1);
0111   StorageIndex* last    = perm.indices().data();                              /* use P as workspace for last */
0112
0113   /* --- Initialize quotient graph ---------------------------------------- */
0114   StorageIndex* Cp = C.outerIndexPtr();
0115   StorageIndex* Ci = C.innerIndexPtr();
0116   for(k = 0; k < n; k++)
0117     len[k] = Cp[k+1] - Cp[k];
0118   len[n] = 0;
0119   nzmax = t;
0120
0121   for(i = 0; i <= n; i++)
0122   {
0123     head[i]   = -1;                     // degree list i is empty
0124     last[i]   = -1;
0125     next[i]   = -1;
0126     hhead[i]  = -1;                     // hash list i is empty
0127     nv[i]     = 1;                      // node i is just one node
0128     w[i]      = 1;                      // node i is alive
0129     elen[i]   = 0;                      // Ek of node i is empty
0130     degree[i] = len[i];                 // degree of node i
0131   }
0132   mark = internal::cs_wclear<StorageIndex>(0, 0, w, n);         /* clear w */
0133
0134   /* --- Initialize degree lists ------------------------------------------ */
0135   for(i = 0; i < n; i++)
0136   {
0137     bool has_diag = false;
0138     for(p = Cp[i]; p<Cp[i+1]; ++p)
0139       if(Ci[p]==i)
0140       {
0141         has_diag = true;
0142         break;
0143       }
0144
0145     d = degree[i];
0146     if(d == 1 && has_diag)           /* node i is empty */
0147     {
0148       elen[i] = -2;                 /* element i is dead */
0149       nel++;
0150       Cp[i] = -1;                   /* i is a root of assembly tree */
0151       w[i] = 0;
0152     }
0153     else if(d > dense || !has_diag)  /* node i is dense or has no structural diagonal element */
0154     {
0155       nv[i] = 0;                    /* absorb i into element n */
0156       elen[i] = -1;                 /* node i is dead */
0157       nel++;
0158       Cp[i] = amd_flip (n);
0159       nv[n]++;
0160     }
0161     else
0162     {
0163       if(head[d] != -1) last[head[d]] = i;
0164       next[i] = head[d];           /* put node i in degree list d */
0165       head[d] = i;
0166     }
0167   }
0168
0169   elen[n] = -2;                         /* n is a dead element */
0170   Cp[n] = -1;                           /* n is a root of assembly tree */
0171   w[n] = 0;                             /* n is a dead element */
0172
0173   while (nel < n)                         /* while (selecting pivots) do */
0174   {
0175     /* --- Select node of minimum approximate degree -------------------- */
0176     for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
0177     if(next[k] != -1) last[next[k]] = -1;
0178     head[mindeg] = next[k];          /* remove k from degree list */
0179     elenk = elen[k];                  /* elenk = |Ek| */
0180     nvk = nv[k];                      /* # of nodes k represents */
0181     nel += nvk;                        /* nv[k] nodes of A eliminated */
0182
0183     /* --- Garbage collection ------------------------------------------- */
0184     if(elenk > 0 && cnz + mindeg >= nzmax)
0185     {
0186       for(j = 0; j < n; j++)
0187       {
0188         if((p = Cp[j]) >= 0)      /* j is a live node or element */
0189         {
0190           Cp[j] = Ci[p];          /* save first entry of object */
0191           Ci[p] = amd_flip (j);    /* first entry is now amd_flip(j) */
0192         }
0193       }
0194       for(q = 0, p = 0; p < cnz; ) /* scan all of memory */
0195       {
0196         if((j = amd_flip (Ci[p++])) >= 0)  /* found object j */
0197         {
0198           Ci[q] = Cp[j];       /* restore first entry of object */
0199           Cp[j] = q++;          /* new pointer to object j */
0200           for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
0201         }
0202       }
0203       cnz = q;                       /* Ci[cnz...nzmax-1] now free */
0204     }
0205
0206     /* --- Construct new element ---------------------------------------- */
0207     dk = 0;
0208     nv[k] = -nvk;                     /* flag k as in Lk */
0209     p = Cp[k];
0210     pk1 = (elenk == 0) ? p : cnz;      /* do in place if elen[k] == 0 */
0211     pk2 = pk1;
0212     for(k1 = 1; k1 <= elenk + 1; k1++)
0213     {
0214       if(k1 > elenk)
0215       {
0216         e = k;                     /* search the nodes in k */
0217         pj = p;                    /* list of nodes starts at Ci[pj]*/
0218         ln = len[k] - elenk;      /* length of list of nodes in k */
0219       }
0220       else
0221       {
0222         e = Ci[p++];              /* search the nodes in e */
0223         pj = Cp[e];
0224         ln = len[e];              /* length of list of nodes in e */
0225       }
0226       for(k2 = 1; k2 <= ln; k2++)
0227       {
0228         i = Ci[pj++];
0229         if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */
0230         dk += nvi;                 /* degree[Lk] += size of node i */
0231         nv[i] = -nvi;             /* negate nv[i] to denote i in Lk*/
0232         Ci[pk2++] = i;            /* place i in Lk */
0233         if(next[i] != -1) last[next[i]] = last[i];
0234         if(last[i] != -1)         /* remove i from degree list */
0235         {
0236           next[last[i]] = next[i];
0237         }
0238         else
0239         {
0240           head[degree[i]] = next[i];
0241         }
0242       }
0243       if(e != k)
0244       {
0245         Cp[e] = amd_flip (k);      /* absorb e into k */
0246         w[e] = 0;                 /* e is now a dead element */
0247       }
0248     }
0249     if(elenk != 0) cnz = pk2;         /* Ci[cnz...nzmax] is free */
0250     degree[k] = dk;                   /* external degree of k - |Lk\i| */
0251     Cp[k] = pk1;                      /* element k is in Ci[pk1..pk2-1] */
0252     len[k] = pk2 - pk1;
0253     elen[k] = -2;                     /* k is now an element */
0254
0255     /* --- Find set differences ----------------------------------------- */
0256     mark = internal::cs_wclear<StorageIndex>(mark, lemax, w, n);  /* clear w if necessary */
0257     for(pk = pk1; pk < pk2; pk++)    /* scan 1: find |Le\Lk| */
0258     {
0259       i = Ci[pk];
0260       if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */
0261       nvi = -nv[i];                      /* nv[i] was negated */
0262       wnvi = mark - nvi;
0263       for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)  /* scan Ei */
0264       {
0265         e = Ci[p];
0266         if(w[e] >= mark)
0267         {
0268           w[e] -= nvi;          /* decrement |Le\Lk| */
0269         }
0270         else if(w[e] != 0)        /* ensure e is a live element */
0271         {
0272           w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */
0273         }
0274       }
0275     }
0276
0277     /* --- Degree update ------------------------------------------------ */
0278     for(pk = pk1; pk < pk2; pk++)    /* scan2: degree update */
0279     {
0280       i = Ci[pk];                   /* consider node i in Lk */
0281       p1 = Cp[i];
0282       p2 = p1 + elen[i] - 1;
0283       pn = p1;
0284       for(h = 0, d = 0, p = p1; p <= p2; p++)    /* scan Ei */
0285       {
0286         e = Ci[p];
0287         if(w[e] != 0)             /* e is an unabsorbed element */
0288         {
0289           dext = w[e] - mark;   /* dext = |Le\Lk| */
0290           if(dext > 0)
0291           {
0292             d += dext;         /* sum up the set differences */
0293             Ci[pn++] = e;     /* keep e in Ei */
0294             h += e;            /* compute the hash of node i */
0295           }
0296           else
0297           {
0298             Cp[e] = amd_flip (k);  /* aggressive absorb. e->k */
0299             w[e] = 0;             /* e is a dead element */
0300           }
0301         }
0302       }
0303       elen[i] = pn - p1 + 1;        /* elen[i] = |Ei| */
0304       p3 = pn;
0305       p4 = p1 + len[i];
0306       for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */
0307       {
0308         j = Ci[p];
0309         if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */
0310         d += nvj;                  /* degree(i) += |j| */
0311         Ci[pn++] = j;             /* place j in node list of i */
0312         h += j;                    /* compute hash for node i */
0313       }
0314       if(d == 0)                     /* check for mass elimination */
0315       {
0316         Cp[i] = amd_flip (k);      /* absorb i into k */
0317         nvi = -nv[i];
0318         dk -= nvi;                 /* |Lk| -= |i| */
0319         nvk += nvi;                /* |k| += nv[i] */
0320         nel += nvi;
0321         nv[i] = 0;
0322         elen[i] = -1;             /* node i is dead */
0323       }
0324       else
0325       {
0326         degree[i] = std::min<StorageIndex> (degree[i], d);   /* update degree(i) */
0327         Ci[pn] = Ci[p3];         /* move first node to end */
0328         Ci[p3] = Ci[p1];         /* move 1st el. to end of Ei */
0329         Ci[p1] = k;               /* add k as 1st element in of Ei */
0330         len[i] = pn - p1 + 1;     /* new len of adj. list of node i */
0331         h %= n;                    /* finalize hash of i */
0332         next[i] = hhead[h];      /* place i in hash bucket */
0333         hhead[h] = i;
0334         last[i] = h;      /* save hash of i in last[i] */
0335       }
0336     }                                   /* scan2 is done */
0337     degree[k] = dk;                   /* finalize |Lk| */
0338     lemax = std::max<StorageIndex>(lemax, dk);
0339     mark = internal::cs_wclear<StorageIndex>(mark+lemax, lemax, w, n);    /* clear w */
0340
0341     /* --- Supernode detection ------------------------------------------ */
0342     for(pk = pk1; pk < pk2; pk++)
0343     {
0344       i = Ci[pk];
0345       if(nv[i] >= 0) continue;         /* skip if i is dead */
0346       h = last[i];                      /* scan hash bucket of node i */
0347       i = hhead[h];
0348       hhead[h] = -1;                    /* hash bucket will be empty */
0349       for(; i != -1 && next[i] != -1; i = next[i], mark++)
0350       {
0351         ln = len[i];
0352         eln = elen[i];
0353         for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
0354         jlast = i;
0355         for(j = next[i]; j != -1; ) /* compare i with all j */
0356         {
0357           ok = (len[j] == ln) && (elen[j] == eln);
0358           for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
0359           {
0360             if(w[Ci[p]] != mark) ok = 0;    /* compare i and j*/
0361           }
0362           if(ok)                     /* i and j are identical */
0363           {
0364             Cp[j] = amd_flip (i);  /* absorb j into i */
0365             nv[i] += nv[j];
0366             nv[j] = 0;
0367             elen[j] = -1;         /* node j is dead */
0368             j = next[j];          /* delete j from hash bucket */
0369             next[jlast] = j;
0370           }
0371           else
0372           {
0373             jlast = j;             /* j and i are different */
0374             j = next[j];
0375           }
0376         }
0377       }
0378     }
0379
0380     /* --- Finalize new element------------------------------------------ */
0381     for(p = pk1, pk = pk1; pk < pk2; pk++)   /* finalize Lk */
0382     {
0383       i = Ci[pk];
0384       if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */
0385       nv[i] = nvi;                      /* restore nv[i] */
0386       d = degree[i] + dk - nvi;         /* compute external degree(i) */
0387       d = std::min<StorageIndex> (d, n - nel - nvi);
0388       if(head[d] != -1) last[head[d]] = i;
0389       next[i] = head[d];               /* put i back in degree list */
0390       last[i] = -1;
0391       head[d] = i;
0392       mindeg = std::min<StorageIndex> (mindeg, d);       /* find new minimum degree */
0393       degree[i] = d;
0394       Ci[p++] = i;                      /* place i in Lk */
0395     }
0396     nv[k] = nvk;                      /* # nodes absorbed into k */
0397     if((len[k] = p-pk1) == 0)         /* length of adj list of element k*/
0398     {
0399       Cp[k] = -1;                   /* k is a root of the tree */
0400       w[k] = 0;                     /* k is now a dead element */
0401     }
0402     if(elenk != 0) cnz = p;           /* free unused space in Lk */
0403   }
0404
0405   /* --- Postordering ----------------------------------------------------- */
0406   for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */
0407   for(j = 0; j <= n; j++) head[j] = -1;
0408   for(j = n; j >= 0; j--)              /* place unordered nodes in lists */
0409   {
0410     if(nv[j] > 0) continue;          /* skip if j is an element */
0411     next[j] = head[Cp[j]];          /* place j in list of its parent */
0412     head[Cp[j]] = j;
0413   }
0414   for(e = n; e >= 0; e--)              /* place elements in lists */
0415   {
0416     if(nv[e] <= 0) continue;         /* skip unless e is an element */
0417     if(Cp[e] != -1)
0418     {
0419       next[e] = head[Cp[e]];      /* place e in list of its parent */
0420       head[Cp[e]] = e;
0421     }
0422   }
0423   for(k = 0, i = 0; i <= n; i++)       /* postorder the assembly tree */
0424   {
0425     if(Cp[i] == -1) k = internal::cs_tdfs<StorageIndex>(i, k, head, next, perm.indices().data(), w);
0426   }
0427
0428   perm.indices().conservativeResize(n);
0429 }
0430
0431 } // namespace internal
0432
0433 } // end namespace Eigen
0434
0435 #endif // EIGEN_SPARSE_AMD_H