Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_rowcolcounts.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_rowcolcounts ======================================== */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
00007  * The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
00008  * Lesser General Public License.  See lesser.txt for a text of the license.
00009  * CHOLMOD is also available under other licenses; contact authors for details.
00010  * http://www.cise.ufl.edu/research/sparse
00011  * -------------------------------------------------------------------------- */
00012 
00013 /* Compute the row and column counts of the Cholesky factor L of the matrix
00014  * A or A*A'.  The etree and its postordering must already be computed (see
00015  * cholmod_etree and cholmod_postorder) and given as inputs to this routine.
00016  *
00017  * For the symmetric case (LL'=A), A is accessed by column.  Only the lower
00018  * triangular part of A is used.  Entries not in this part of the matrix are
00019  * ignored.  This is the same as storing the upper triangular part of A by
00020  * rows, with entries in the lower triangular part being ignored.  NOTE: this
00021  * representation is the TRANSPOSE of the input to cholmod_etree.
00022  *
00023  * For the unsymmetric case (LL'=AA'), A is accessed by column.  Equivalently,
00024  * if A is viewed as a matrix in compressed-row form, this routine computes
00025  * the row and column counts for L where LL'=A'A.  If the input vector f is
00026  * present, then F*F' is analyzed instead, where F = A(:,f).
00027  *
00028  * The set f is held in fset and fsize.
00029  *  fset = NULL means ":" in MATLAB. fset is ignored.
00030  *  fset != NULL means f = fset [0..fset-1].
00031  *  fset != NULL and fsize = 0 means f is the empty set.
00032  *  Common->status is set to CHOLMOD_INVALID if fset is invalid.
00033  *
00034  * In both cases, the columns of A need not be sorted.
00035  * A can be packed or unpacked.
00036  *
00037  * References:
00038  * J. Gilbert, E. Ng, B. Peyton, "An efficient algorithm to compute row and
00039  * column counts for sparse Cholesky factorization", SIAM J. Matrix Analysis &
00040  * Applic., vol 15, 1994, pp. 1075-1091.
00041  *
00042  * J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
00043  * sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
00044  *
00045  * workspace:
00046  *  if symmetric:   Flag (nrow), Iwork (2*nrow)
00047  *  if unsymmetric: Flag (nrow), Iwork (2*nrow+ncol), Head (nrow+1)
00048  *
00049  * Supports any xtype (pattern, real, complex, or zomplex).
00050  */
00051 
00052 #ifndef NCHOLESKY
00053 
00054 #include "amesos_cholmod_internal.h"
00055 #include "amesos_cholmod_cholesky.h"
00056 
00057 /* ========================================================================== */
00058 /* === initialize_node ====================================================== */
00059 /* ========================================================================== */
00060 
00061 static int amesos_initialize_node  /* initial work for kth node in postordered etree */
00062 (
00063     Int k,    /* at the kth step of the algorithm (and kth node) */
00064     Int Post [ ], /* Post [k] = i, the kth node in postordered etree */
00065     Int Parent [ ], /* Parent [i] is the parent of i in the etree */
00066     Int ColCount [ ], /* ColCount [c] is the current weight of node c */
00067     Int PrevNbr [ ] /* PrevNbr [u] = k if u was last considered at step k */
00068 )
00069 {
00070     Int p, parent ;
00071     /* determine p, the kth node in the postordered etree */
00072     p = Post [k] ;
00073     /* adjust the weight if p is not a root of the etree */
00074     parent = Parent [p] ;
00075     if (parent != EMPTY)
00076     {
00077   ColCount [parent]-- ;
00078     }
00079     /* flag node p to exclude self edges (p,p) */
00080     PrevNbr [p] = k ;
00081     return (p) ;
00082 }
00083 
00084 
00085 /* ========================================================================== */
00086 /* === process_edge ========================================================= */
00087 /* ========================================================================== */
00088 
00089 /* edge (p,u) is being processed.  p < u is a descendant of its ancestor u in
00090  * the etree.  node p is the kth node in the postordered etree.  */
00091 
00092 static void amesos_process_edge
00093 (
00094     Int p,    /* process edge (p,u) of the matrix */
00095     Int u,
00096     Int k,    /* we are at the kth node in the postordered etree */
00097     Int First [ ],  /* First [i] = k if the postordering of first
00098        * descendent of node i is k */
00099     Int PrevNbr [ ],  /* u was last considered at step k = PrevNbr [u] */
00100     Int ColCount [ ], /* ColCount [c] is the current weight of node c */
00101     Int PrevLeaf [ ], /* s = PrevLeaf [u] means that s was the last leaf
00102        * seen in the subtree rooted at u.  */
00103     Int RowCount [ ], /* RowCount [i] is # of nonzeros in row i of L,
00104        * including the diagonal.  Not computed if NULL. */
00105     Int SetParent [ ],  /* the FIND/UNION data structure, which forms a set
00106        * of trees.  A root i has i = SetParent [i].  Following
00107        * a path from i to the root q of the subtree containing
00108        * i means that q is the SetParent representative of i.
00109        * All nodes in the tree could have their SetParent
00110        * equal to the root q; the tree representation is used
00111        * to save time.  When a path is traced from i to its
00112        * root q, the path is re-traversed to set the SetParent
00113        * of the whole path to be the root q. */
00114     Int Level [ ]  /* Level [i] = length of path from node i to root */
00115 )
00116 {
00117     Int prevleaf, q, s, sparent ;
00118     if (First [p] > PrevNbr [u])
00119     {
00120   /* p is a leaf of the subtree of u */
00121   ColCount [p]++ ;
00122   prevleaf = PrevLeaf [u] ;
00123   if (prevleaf == EMPTY)
00124   {
00125       /* p is the first leaf of subtree of u; RowCount will be incremented
00126        * by the length of the path in the etree from p up to u. */
00127       q = u ;
00128   }
00129   else
00130   {
00131       /* q = FIND (prevleaf): find the root q of the
00132        * SetParent tree containing prevleaf */
00133       for (q = prevleaf ; q != SetParent [q] ; q = SetParent [q])
00134       {
00135     ;
00136       }
00137       /* the root q has been found; re-traverse the path and
00138        * perform path compression */
00139       s = prevleaf ;
00140       for (s = prevleaf ; s != q ; s = sparent)
00141       {
00142     sparent = SetParent [s] ;
00143     SetParent [s] = q ;
00144       }
00145       /* adjust the RowCount and ColCount; RowCount will be incremented by
00146        * the length of the path from p to the SetParent root q, and
00147        * decrement the ColCount of q by one. */
00148       ColCount [q]-- ;
00149   }
00150   if (RowCount != NULL)
00151   {
00152       /* if RowCount is being computed, increment it by the length of
00153        * the path from p to q */
00154       RowCount [u] += (Level [p] - Level [q]) ;
00155   }
00156   /* p is a leaf of the subtree of u, so mark PrevLeaf [u] to be p */
00157   PrevLeaf [u] = p ;
00158     }
00159     /* flag u has having been processed at step k */
00160     PrevNbr [u] = k ;
00161 }
00162 
00163 
00164 /* ========================================================================== */
00165 /* === finalize_node ======================================================== */
00166 /* ========================================================================== */
00167 
00168 static void amesos_finalize_node    /* compute UNION (p, Parent [p]) */
00169 (
00170     Int p,
00171     Int Parent [ ], /* Parent [p] is the parent of p in the etree */
00172     Int SetParent [ ] /* see process_edge, above */
00173 )
00174 {
00175     /* all nodes in the SetParent tree rooted at p now have as their final
00176      * root the node Parent [p].  This computes UNION (p, Parent [p]) */
00177     if (Parent [p] != EMPTY)
00178     {
00179   SetParent [p] = Parent [p] ;
00180     }
00181 }
00182 
00183 
00184 /* ========================================================================== */
00185 /* === cholmod_rowcolcounts ================================================= */
00186 /* ========================================================================== */
00187 
00188 int CHOLMOD(rowcolcounts)
00189 (
00190     /* ---- input ---- */
00191     cholmod_sparse *A,  /* matrix to analyze */
00192     Int *fset,    /* subset of 0:(A->ncol)-1 */
00193     size_t fsize, /* size of fset */
00194     Int *Parent,  /* size nrow.  Parent [i] = p if p is the parent of i */
00195     Int *Post,    /* size nrow.  Post [k] = i if i is the kth node in
00196        * the postordered etree. */
00197     /* ---- output --- */
00198     Int *RowCount,  /* size nrow. RowCount [i] = # entries in the ith row of
00199        * L, including the diagonal. */
00200     Int *ColCount,  /* size nrow. ColCount [i] = # entries in the ith
00201        * column of L, including the diagonal. */
00202     Int *First,   /* size nrow.  First [i] = k is the least postordering
00203        * of any descendant of i. */
00204     Int *Level,   /* size nrow.  Level [i] is the length of the path from
00205        * i to the root, with Level [root] = 0. */
00206     /* --------------- */
00207     cholmod_common *Common
00208 )
00209 {
00210     double fl, ff ;
00211     Int *Ap, *Ai, *Anz, *PrevNbr, *SetParent, *Head, *PrevLeaf, *Anext, *Ipost,
00212   *Iwork ;
00213     Int i, j, r, k, len, s, p, pend, inew, stype, nf, anz, inode, parent,
00214   nrow, ncol, packed, use_fset, jj ;
00215     size_t w ;
00216     int ok = TRUE ;
00217 
00218     /* ---------------------------------------------------------------------- */
00219     /* check inputs */
00220     /* ---------------------------------------------------------------------- */
00221 
00222     RETURN_IF_NULL_COMMON (FALSE) ;
00223     RETURN_IF_NULL (A, FALSE) ;
00224     RETURN_IF_NULL (Parent, FALSE) ;
00225     RETURN_IF_NULL (Post, FALSE) ;
00226     RETURN_IF_NULL (ColCount, FALSE) ;
00227     RETURN_IF_NULL (First, FALSE) ;
00228     RETURN_IF_NULL (Level, FALSE) ;
00229     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00230     stype = A->stype ;
00231     if (stype > 0)
00232     {
00233   /* symmetric with upper triangular part not supported */
00234   ERROR (CHOLMOD_INVALID, "symmetric upper not supported") ;
00235   return (FALSE) ;
00236     }
00237     Common->status = CHOLMOD_OK ;
00238 
00239     /* ---------------------------------------------------------------------- */
00240     /* allocate workspace */
00241     /* ---------------------------------------------------------------------- */
00242 
00243     nrow = A->nrow ;  /* the number of rows of A */
00244     ncol = A->ncol ;  /* the number of columns of A */
00245 
00246     /* w = 2*nrow + (stype ? 0 : ncol) */
00247     w = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
00248     w = CHOLMOD(add_size_t) (w, (stype ? 0 : ncol), &ok) ;
00249     if (!ok)
00250     {
00251   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00252   return (FALSE) ;
00253     }
00254 
00255     CHOLMOD(allocate_work) (nrow, w, 0, Common) ;
00256     if (Common->status < CHOLMOD_OK)
00257     {
00258   return (FALSE) ;
00259     }
00260 
00261     ASSERT (CHOLMOD(dump_perm) (Post, nrow, nrow, "Post", Common)) ;
00262     ASSERT (CHOLMOD(dump_parent) (Parent, nrow, "Parent", Common)) ;
00263 
00264     /* ---------------------------------------------------------------------- */
00265     /* get inputs */
00266     /* ---------------------------------------------------------------------- */
00267 
00268     Ap = A->p ; /* size ncol+1, column pointers for A */
00269     Ai = A->i ; /* the row indices of A, of size nz=Ap[ncol+1] */
00270     Anz = A->nz ;
00271     packed = A->packed ;
00272     ASSERT (IMPLIES (!packed, Anz != NULL)) ;
00273 
00274     /* ---------------------------------------------------------------------- */
00275     /* get workspace */
00276     /* ---------------------------------------------------------------------- */
00277 
00278     Iwork = Common->Iwork ;
00279     SetParent = Iwork ;       /* size nrow (i/i/l) */
00280     PrevNbr   = Iwork + nrow ;      /* size nrow (i/i/l) */
00281     Anext     = Iwork + 2*((size_t) nrow) ;    /* size ncol (i/i/l) (unsym only) */
00282     PrevLeaf  = Common->Flag ;      /* size nrow */
00283     Head      = Common->Head ;      /* size nrow+1 (unsym only)*/
00284 
00285     /* ---------------------------------------------------------------------- */
00286     /* find the first descendant and level of each node in the tree */
00287     /* ---------------------------------------------------------------------- */
00288 
00289     /* First [i] = k if the postordering of first descendent of node i is k */
00290     /* Level [i] = length of path from node i to the root (Level [root] = 0) */
00291 
00292     for (i = 0 ; i < nrow ; i++)
00293     {
00294   First [i] = EMPTY ;
00295     }
00296 
00297     /* postorder traversal of the etree */
00298     for (k = 0 ; k < nrow ; k++)
00299     {
00300   /* node i of the etree is the kth node in the postordered etree */
00301   i = Post [k] ;
00302 
00303   /* i is a leaf if First [i] is still EMPTY */
00304   /* ColCount [i] starts at 1 if i is a leaf, zero otherwise */
00305   ColCount [i] = (First [i] == EMPTY) ? 1 : 0 ;
00306 
00307   /* traverse the path from node i to the root, stopping if we find a
00308    * node r whose First [r] is already defined. */
00309   len = 0 ;
00310   for (r = i ; (r != EMPTY) && (First [r] == EMPTY) ; r = Parent [r])
00311   {
00312       First [r] = k ;
00313       len++ ;
00314   }
00315   if (r == EMPTY)
00316   {
00317       /* we hit a root node, the level of which is zero */
00318       len-- ;
00319   }
00320   else
00321   {
00322       /* we stopped at node r, where Level [r] is already defined */
00323       len += Level [r] ;
00324   }
00325   /* re-traverse the path from node i to r; set the level of each node */
00326   for (s = i ; s != r ; s = Parent [s])
00327   {
00328       Level [s] = len-- ;
00329   }
00330     }
00331 
00332     /* ---------------------------------------------------------------------- */
00333     /* AA' case: sort columns of A according to first postordered row index */
00334     /* ---------------------------------------------------------------------- */
00335 
00336     fl = 0.0 ;
00337     if (stype == 0)
00338     {
00339   /* [ use PrevNbr [0..nrow-1] as workspace for Ipost */
00340   Ipost = PrevNbr ;
00341   /* Ipost [i] = k if i is the kth node in the postordered etree. */
00342   for (k = 0 ; k < nrow ; k++)
00343   {
00344       Ipost [Post [k]] = k ;
00345   }
00346   use_fset = (fset != NULL) ;
00347   if (use_fset)
00348   {
00349       nf = fsize ;
00350       /* clear Anext to check fset */
00351       for (j = 0 ; j < ncol ; j++)
00352       {
00353     Anext [j] = -2 ;
00354       }
00355       /* find the first postordered row in each column of A (post,f)
00356        * and place the column in the corresponding link list */
00357       for (jj = 0 ; jj < nf ; jj++)
00358       {
00359     j = fset [jj] ;
00360     if (j < 0 || j > ncol || Anext [j] != -2)
00361     {
00362         /* out-of-range or duplicate entry in fset */
00363         ERROR (CHOLMOD_INVALID, "fset invalid") ;
00364         return (FALSE) ;
00365     }
00366     /* flag column j as having been seen */
00367     Anext [j] = EMPTY ;
00368       }
00369       /* fset is now valid */
00370       ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
00371   }
00372   else
00373   {
00374       nf = ncol ;
00375   }
00376   for (jj = 0 ; jj < nf ; jj++)
00377   {
00378       j = (use_fset) ? (fset [jj]) : jj ;
00379       /* column j is in the fset; find the smallest row (if any) */
00380       p = Ap [j] ;
00381       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00382       ff = (double) MAX (0, pend - p) ;
00383       fl += ff*ff + ff ;
00384       if (pend > p)
00385       {
00386     k = Ipost [Ai [p]] ;
00387     for ( ; p < pend ; p++)
00388     {
00389         inew = Ipost [Ai [p]] ;
00390         k = MIN (k, inew) ;
00391     }
00392     /* place column j in link list k */
00393     ASSERT (k >= 0 && k < nrow) ;
00394     Anext [j] = Head [k] ;
00395     Head [k] = j ;
00396       }
00397   }
00398   /* Ipost no longer needed for inverse postordering ]
00399    * Head [k] contains a link list of all columns whose first
00400    * postordered row index is equal to k, for k = 0 to nrow-1. */
00401     }
00402 
00403     /* ---------------------------------------------------------------------- */
00404     /* compute the row counts and node weights */
00405     /* ---------------------------------------------------------------------- */
00406 
00407     if (RowCount != NULL)
00408     {
00409   for (i = 0 ; i < nrow ; i++)
00410   {
00411       RowCount [i] = 1 ;
00412   }
00413     }
00414     for (i = 0 ; i < nrow ; i++)
00415     {
00416   PrevLeaf [i] = EMPTY ;
00417   PrevNbr [i] = EMPTY ;
00418   SetParent [i] = i ; /* every node is in its own set, by itself */
00419     }
00420 
00421     if (stype != 0)
00422     {
00423 
00424   /* ------------------------------------------------------------------ */
00425   /* symmetric case: LL' = A */
00426   /* ------------------------------------------------------------------ */
00427 
00428   /* also determine the number of entries in triu(A) */
00429   anz = nrow ;
00430   for (k = 0 ; k < nrow ; k++)
00431   {
00432       /* j is the kth node in the postordered etree */
00433       j = amesos_initialize_node (k, Post, Parent, ColCount, PrevNbr) ;
00434 
00435       /* for all nonzeros A(i,j) below the diagonal, in column j of A */
00436       p = Ap [j] ;
00437       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00438       for ( ; p < pend ; p++)
00439       {
00440     i = Ai [p] ;
00441     if (i > j)
00442     {
00443         /* j is a descendant of i in etree(A) */
00444         anz++ ;
00445         amesos_process_edge (j, i, k, First, PrevNbr, ColCount,
00446           PrevLeaf, RowCount, SetParent, Level) ;
00447     }
00448       }
00449       /* update SetParent: UNION (j, Parent [j]) */
00450       amesos_finalize_node (j, Parent, SetParent) ;
00451   }
00452   Common->anz = anz ;
00453     }
00454     else
00455     {
00456 
00457   /* ------------------------------------------------------------------ */
00458   /* unsymmetric case: LL' = AA' */
00459   /* ------------------------------------------------------------------ */
00460 
00461   for (k = 0 ; k < nrow ; k++)
00462   {
00463       /* inode is the kth node in the postordered etree */
00464       inode = amesos_initialize_node (k, Post, Parent, ColCount, PrevNbr) ;
00465 
00466       /* for all cols j whose first postordered row is k: */
00467       for (j = Head [k] ; j != EMPTY ; j = Anext [j])
00468       {
00469     /* k is the first postordered row in column j of A */
00470     /* for all rows i in column j: */
00471     p = Ap [j] ;
00472     pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00473     for ( ; p < pend ; p++)
00474     {
00475         i = Ai [p] ;
00476         /* has i already been considered at this step k */
00477         if (PrevNbr [i] < k)
00478         {
00479       /* inode is a descendant of i in etree(AA') */
00480       /* process edge (inode,i) and set PrevNbr[i] to k */
00481       amesos_process_edge (inode, i, k, First, PrevNbr, ColCount,
00482         PrevLeaf, RowCount, SetParent, Level) ;
00483         }
00484     }
00485       }
00486       /* clear link list k */
00487       Head [k] = EMPTY ;
00488       /* update SetParent: UNION (inode, Parent [inode]) */
00489       amesos_finalize_node (inode, Parent, SetParent) ;
00490   }
00491     }
00492 
00493     /* ---------------------------------------------------------------------- */
00494     /* finish computing the column counts */
00495     /* ---------------------------------------------------------------------- */
00496 
00497     for (j = 0 ; j < nrow ; j++)
00498     {
00499   parent = Parent [j] ;
00500   if (parent != EMPTY)
00501   {
00502       /* add the ColCount of j to its parent */
00503       ColCount [parent] += ColCount [j] ;
00504   }
00505     }
00506 
00507     /* ---------------------------------------------------------------------- */
00508     /* clear workspace */
00509     /* ---------------------------------------------------------------------- */
00510 
00511     Common->mark = EMPTY ;
00512     CHOLMOD(clear_flag) (Common) ;
00513     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00514 
00515     /* ---------------------------------------------------------------------- */
00516     /* flop count and nnz(L) for subsequent LL' numerical factorization */
00517     /* ---------------------------------------------------------------------- */
00518 
00519     /* use double to avoid integer overflow.  lnz cannot be NaN. */
00520     Common->aatfl = fl ;
00521     Common->lnz = 0. ;
00522     fl = 0 ;
00523     for (j = 0 ; j < nrow ; j++)
00524     {
00525   ff = (double) (ColCount [j]) ;
00526   Common->lnz += ff ;
00527   fl += ff*ff ;
00528     }
00529 
00530     Common->fl = fl ;
00531     PRINT1 (("rowcol fl %g lnz %g\n", Common->fl, Common->lnz)) ;
00532 
00533     return (TRUE) ;
00534 }
00535 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines