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