cs_amd.c
Go to the documentation of this file.
00001 #include "cs.h"
00002 /* clear w */
00003 static int cs_wclear (int mark, int lemax, int *w, int n)
00004 {
00005     int k ;
00006     if (mark < 2 || (mark + lemax < 0))
00007     {
00008         for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ;
00009         mark = 2 ;
00010     }
00011     return (mark) ;     /* at this point, w [0..n-1] < mark holds */
00012 }
00013 
00014 /* keep off-diagonal entries; drop diagonal entries */
00015 static int cs_diag (int i, int j, double aij, void *other) { return (i != j) ; }
00016 
00017 /* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
00018 int *cs_amd (int order, const cs *A)  /* order 0:natural, 1:Chol, 2:LU, 3:QR */
00019 {
00020     cs *C, *A2, *AT ;
00021     int *Cp, *Ci, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w,
00022         *hhead, *ATp, *ATi, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
00023         k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
00024         ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ;
00025     unsigned int h ;
00026     /* --- Construct matrix C ----------------------------------------------- */
00027     if (!CS_CSC (A) || order <= 0 || order > 3) return (NULL) ; /* check */
00028     AT = cs_transpose (A, 0) ;              /* compute A' */
00029     if (!AT) return (NULL) ;
00030     m = A->m ; n = A->n ;
00031     dense = CS_MAX (16, 10 * (int)sqrt ((double) n)) ;   /* find dense threshold */
00032     dense = CS_MIN (n-2, dense) ;
00033     if (order == 1 && n == m)
00034     {
00035         C = cs_add (A, AT, 0, 0) ;          /* C = A+A' */
00036     }
00037     else if (order == 2)
00038     {
00039         ATp = AT->p ;                       /* drop dense columns from AT */
00040         ATi = AT->i ;
00041         for (p2 = 0, j = 0 ; j < m ; j++)
00042         {
00043             p = ATp [j] ;                   /* column j of AT starts here */
00044             ATp [j] = p2 ;                  /* new column j starts here */
00045             if (ATp [j+1] - p > dense) continue ;   /* skip dense col j */
00046             for ( ; p < ATp [j+1] ; p++) ATi [p2++] = ATi [p] ;
00047         }
00048         ATp [m] = p2 ;                      /* finalize AT */
00049         A2 = cs_transpose (AT, 0) ;         /* A2 = AT' */
00050         C = A2 ? cs_multiply (AT, A2) : NULL ;  /* C=A'*A with no dense rows */
00051         cs_spfree (A2) ;
00052     }
00053     else
00054     {
00055         C = cs_multiply (AT, A) ;           /* C=A'*A */
00056     }
00057     cs_spfree (AT) ;
00058     if (!C) return (NULL) ;
00059     cs_fkeep (C, &cs_diag, NULL) ;          /* drop diagonal entries */
00060     Cp = C->p ;
00061     cnz = Cp [n] ;
00062     P = (int*) cs_malloc (n+1, sizeof (int)) ;     /* allocate result */
00063     W = (int*) cs_malloc (8*(n+1), sizeof (int)) ; /* get workspace */
00064     t = cnz + cnz/5 + 2*n ;                 /* add elbow room to C */
00065     if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ;
00066     len  = W           ; nv     = W +   (n+1) ; next   = W + 2*(n+1) ;
00067     head = W + 3*(n+1) ; elen   = W + 4*(n+1) ; degree = W + 5*(n+1) ;
00068     w    = W + 6*(n+1) ; hhead  = W + 7*(n+1) ;
00069     last = P ;                              /* use P as workspace for last */
00070     /* --- Initialize quotient graph ---------------------------------------- */
00071     for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
00072     len [n] = 0 ;
00073     nzmax = C->nzmax ;
00074     Ci = C->i ;
00075     for (i = 0 ; i <= n ; i++)
00076     {
00077         head [i] = -1 ;                     /* degree list i is empty */
00078         last [i] = -1 ;
00079         next [i] = -1 ;
00080         hhead [i] = -1 ;                    /* hash list i is empty */
00081         nv [i] = 1 ;                        /* node i is just one node */
00082         w [i] = 1 ;                         /* node i is alive */
00083         elen [i] = 0 ;                      /* Ek of node i is empty */
00084         degree [i] = len [i] ;              /* degree of node i */
00085     }
00086     mark = cs_wclear (0, 0, w, n) ;         /* clear w */
00087     elen [n] = -2 ;                         /* n is a dead element */
00088     Cp [n] = -1 ;                           /* n is a root of assembly tree */
00089     w [n] = 0 ;                             /* n is a dead element */
00090     /* --- Initialize degree lists ------------------------------------------ */
00091     for (i = 0 ; i < n ; i++)
00092     {
00093         d = degree [i] ;
00094         if (d == 0)                         /* node i is empty */
00095         {
00096             elen [i] = -2 ;                 /* element i is dead */
00097             nel++ ;
00098             Cp [i] = -1 ;                   /* i is a root of assembly tree */
00099             w [i] = 0 ;
00100         }
00101         else if (d > dense)                 /* node i is dense */
00102         {
00103             nv [i] = 0 ;                    /* absorb i into element n */
00104             elen [i] = -1 ;                 /* node i is dead */
00105             nel++ ;
00106             Cp [i] = CS_FLIP (n) ;
00107             nv [n]++ ;
00108         }
00109         else
00110         {
00111             if (head [d] != -1) last [head [d]] = i ;
00112             next [i] = head [d] ;           /* put node i in degree list d */
00113             head [d] = i ;
00114         }
00115     }
00116     while (nel < n)                         /* while (selecting pivots) do */
00117     {
00118         /* --- Select node of minimum approximate degree -------------------- */
00119         for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
00120         if (next [k] != -1) last [next [k]] = -1 ;
00121         head [mindeg] = next [k] ;          /* remove k from degree list */
00122         elenk = elen [k] ;                  /* elenk = |Ek| */
00123         nvk = nv [k] ;                      /* # of nodes k represents */
00124         nel += nvk ;                        /* nv[k] nodes of A eliminated */
00125         /* --- Garbage collection ------------------------------------------- */
00126         if (elenk > 0 && cnz + mindeg >= nzmax)
00127         {
00128             for (j = 0 ; j < n ; j++)
00129             {
00130                 if ((p = Cp [j]) >= 0)      /* j is a live node or element */
00131                 {
00132                     Cp [j] = Ci [p] ;       /* save first entry of object */
00133                     Ci [p] = CS_FLIP (j) ;  /* first entry is now CS_FLIP(j) */
00134                 }
00135             }
00136             for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */
00137             {
00138                 if ((j = CS_FLIP (Ci [p++])) >= 0)  /* found object j */
00139                 {
00140                     Ci [q] = Cp [j] ;       /* restore first entry of object */
00141                     Cp [j] = q++ ;          /* new pointer to object j */
00142                     for (k3 = 0 ; k3 < len [j]-1 ; k3++) Ci [q++] = Ci [p++] ;
00143                 }
00144             }
00145             cnz = q ;                       /* Ci [cnz...nzmax-1] now free */
00146         }
00147         /* --- Construct new element ---------------------------------------- */
00148         dk = 0 ;
00149         nv [k] = -nvk ;                     /* flag k as in Lk */
00150         p = Cp [k] ;
00151         pk1 = (elenk == 0) ? p : cnz ;      /* do in place if elen[k] == 0 */
00152         pk2 = pk1 ;
00153         for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
00154         {
00155             if (k1 > elenk)
00156             {
00157                 e = k ;                     /* search the nodes in k */
00158                 pj = p ;                    /* list of nodes starts at Ci[pj]*/
00159                 ln = len [k] - elenk ;      /* length of list of nodes in k */
00160             }
00161             else
00162             {
00163                 e = Ci [p++] ;              /* search the nodes in e */
00164                 pj = Cp [e] ;
00165                 ln = len [e] ;              /* length of list of nodes in e */
00166             }
00167             for (k2 = 1 ; k2 <= ln ; k2++)
00168             {
00169                 i = Ci [pj++] ;
00170                 if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */
00171                 dk += nvi ;                 /* degree[Lk] += size of node i */
00172                 nv [i] = -nvi ;             /* negate nv[i] to denote i in Lk*/
00173                 Ci [pk2++] = i ;            /* place i in Lk */
00174                 if (next [i] != -1) last [next [i]] = last [i] ;
00175                 if (last [i] != -1)         /* remove i from degree list */
00176                 {
00177                     next [last [i]] = next [i] ;
00178                 }
00179                 else
00180                 {
00181                     head [degree [i]] = next [i] ;
00182                 }
00183             }
00184             if (e != k)
00185             {
00186                 Cp [e] = CS_FLIP (k) ;      /* absorb e into k */
00187                 w [e] = 0 ;                 /* e is now a dead element */
00188             }
00189         }
00190         if (elenk != 0) cnz = pk2 ;         /* Ci [cnz...nzmax] is free */
00191         degree [k] = dk ;                   /* external degree of k - |Lk\i| */
00192         Cp [k] = pk1 ;                      /* element k is in Ci[pk1..pk2-1] */
00193         len [k] = pk2 - pk1 ;
00194         elen [k] = -2 ;                     /* k is now an element */
00195         /* --- Find set differences ----------------------------------------- */
00196         mark = cs_wclear (mark, lemax, w, n) ;  /* clear w if necessary */
00197         for (pk = pk1 ; pk < pk2 ; pk++)    /* scan 1: find |Le\Lk| */
00198         {
00199             i = Ci [pk] ;
00200             if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */
00201             nvi = -nv [i] ;                      /* nv [i] was negated */
00202             wnvi = mark - nvi ;
00203             for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++)  /* scan Ei */
00204             {
00205                 e = Ci [p] ;
00206                 if (w [e] >= mark)
00207                 {
00208                     w [e] -= nvi ;          /* decrement |Le\Lk| */
00209                 }
00210                 else if (w [e] != 0)        /* ensure e is a live element */
00211                 {
00212                     w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */
00213                 }
00214             }
00215         }
00216         /* --- Degree update ------------------------------------------------ */
00217         for (pk = pk1 ; pk < pk2 ; pk++)    /* scan2: degree update */
00218         {
00219             i = Ci [pk] ;                   /* consider node i in Lk */
00220             p1 = Cp [i] ;
00221             p2 = p1 + elen [i] - 1 ;
00222             pn = p1 ;
00223             for (h = 0, d = 0, p = p1 ; p <= p2 ; p++)    /* scan Ei */
00224             {
00225                 e = Ci [p] ;
00226                 if (w [e] != 0)             /* e is an unabsorbed element */
00227                 {
00228                     dext = w [e] - mark ;   /* dext = |Le\Lk| */
00229                     if (dext > 0)
00230                     {
00231                         d += dext ;         /* sum up the set differences */
00232                         Ci [pn++] = e ;     /* keep e in Ei */
00233                         h += e ;            /* compute the hash of node i */
00234                     }
00235                     else
00236                     {
00237                         Cp [e] = CS_FLIP (k) ;  /* aggressive absorb. e->k */
00238                         w [e] = 0 ;             /* e is a dead element */
00239                     }
00240                 }
00241             }
00242             elen [i] = pn - p1 + 1 ;        /* elen[i] = |Ei| */
00243             p3 = pn ;
00244             p4 = p1 + len [i] ;
00245             for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */
00246             {
00247                 j = Ci [p] ;
00248                 if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
00249                 d += nvj ;                  /* degree(i) += |j| */
00250                 Ci [pn++] = j ;             /* place j in node list of i */
00251                 h += j ;                    /* compute hash for node i */
00252             }
00253             if (d == 0)                     /* check for mass elimination */
00254             {
00255                 Cp [i] = CS_FLIP (k) ;      /* absorb i into k */
00256                 nvi = -nv [i] ;
00257                 dk -= nvi ;                 /* |Lk| -= |i| */
00258                 nvk += nvi ;                /* |k| += nv[i] */
00259                 nel += nvi ;
00260                 nv [i] = 0 ;
00261                 elen [i] = -1 ;             /* node i is dead */
00262             }
00263             else
00264             {
00265                 degree [i] = CS_MIN (degree [i], d) ;   /* update degree(i) */
00266                 Ci [pn] = Ci [p3] ;         /* move first node to end */
00267                 Ci [p3] = Ci [p1] ;         /* move 1st el. to end of Ei */
00268                 Ci [p1] = k ;               /* add k as 1st element in of Ei */
00269                 len [i] = pn - p1 + 1 ;     /* new len of adj. list of node i */
00270                 h %= n ;                    /* finalize hash of i */
00271                 next [i] = hhead [h] ;      /* place i in hash bucket */
00272                 hhead [h] = i ;
00273                 last [i] = h ;              /* save hash of i in last[i] */
00274             }
00275         }                                   /* scan2 is done */
00276         degree [k] = dk ;                   /* finalize |Lk| */
00277         lemax = CS_MAX (lemax, dk) ;
00278         mark = cs_wclear (mark+lemax, lemax, w, n) ;    /* clear w */
00279         /* --- Supernode detection ------------------------------------------ */
00280         for (pk = pk1 ; pk < pk2 ; pk++)
00281         {
00282             i = Ci [pk] ;
00283             if (nv [i] >= 0) continue ;         /* skip if i is dead */
00284             h = last [i] ;                      /* scan hash bucket of node i */
00285             i = hhead [h] ;
00286             hhead [h] = -1 ;                    /* hash bucket will be empty */
00287             for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
00288             {
00289                 ln = len [i] ;
00290                 eln = elen [i] ;
00291                 for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark;
00292                 jlast = i ;
00293                 for (j = next [i] ; j != -1 ; ) /* compare i with all j */
00294                 {
00295                     ok = (len [j] == ln) && (elen [j] == eln) ;
00296                     for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
00297                     {
00298                         if (w [Ci [p]] != mark) ok = 0 ;    /* compare i and j*/
00299                     }
00300                     if (ok)                     /* i and j are identical */
00301                     {
00302                         Cp [j] = CS_FLIP (i) ;  /* absorb j into i */
00303                         nv [i] += nv [j] ;
00304                         nv [j] = 0 ;
00305                         elen [j] = -1 ;         /* node j is dead */
00306                         j = next [j] ;          /* delete j from hash bucket */
00307                         next [jlast] = j ;
00308                     }
00309                     else
00310                     {
00311                         jlast = j ;             /* j and i are different */
00312                         j = next [j] ;
00313                     }
00314                 }
00315             }
00316         }
00317         /* --- Finalize new element------------------------------------------ */
00318         for (p = pk1, pk = pk1 ; pk < pk2 ; pk++)   /* finalize Lk */
00319         {
00320             i = Ci [pk] ;
00321             if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
00322             nv [i] = nvi ;                      /* restore nv[i] */
00323             d = degree [i] + dk - nvi ;         /* compute external degree(i) */
00324             d = CS_MIN (d, n - nel - nvi) ;
00325             if (head [d] != -1) last [head [d]] = i ;
00326             next [i] = head [d] ;               /* put i back in degree list */
00327             last [i] = -1 ;
00328             head [d] = i ;
00329             mindeg = CS_MIN (mindeg, d) ;       /* find new minimum degree */
00330             degree [i] = d ;
00331             Ci [p++] = i ;                      /* place i in Lk */
00332         }
00333         nv [k] = nvk ;                      /* # nodes absorbed into k */
00334         if ((len [k] = p-pk1) == 0)         /* length of adj list of element k*/
00335         {
00336             Cp [k] = -1 ;                   /* k is a root of the tree */
00337             w [k] = 0 ;                     /* k is now a dead element */
00338         }
00339         if (elenk != 0) cnz = p ;           /* free unused space in Lk */
00340     }
00341     /* --- Postordering ----------------------------------------------------- */
00342     for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
00343     for (j = 0 ; j <= n ; j++) head [j] = -1 ;
00344     for (j = n ; j >= 0 ; j--)              /* place unordered nodes in lists */
00345     {
00346         if (nv [j] > 0) continue ;          /* skip if j is an element */
00347         next [j] = head [Cp [j]] ;          /* place j in list of its parent */
00348         head [Cp [j]] = j ;
00349     }
00350     for (e = n ; e >= 0 ; e--)              /* place elements in lists */
00351     {
00352         if (nv [e] <= 0) continue ;         /* skip unless e is an element */
00353         if (Cp [e] != -1)
00354         {
00355             next [e] = head [Cp [e]] ;      /* place e in list of its parent */
00356             head [Cp [e]] = e ;
00357         }
00358     }
00359     for (k = 0, i = 0 ; i <= n ; i++)       /* postorder the assembly tree */
00360     {
00361         if (Cp [i] == -1) k = cs_tdfs (i, k, head, next, P, w) ;
00362     }
00363     return (cs_idone (P, C, W, 1)) ;
00364 }


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:56