Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_etree.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_etree =============================================== */
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 elimination tree of A or A'*A
00014  *
00015  * In the symmetric case, the upper triangular part of A is used.  Entries not
00016  * in this part of the matrix are ignored.  Computing the etree of a symmetric
00017  * matrix from just its lower triangular entries is not supported.
00018  *
00019  * In the unsymmetric case, all of A is used, and the etree of A'*A is computed.
00020  *
00021  * References:
00022  *
00023  * J. Liu, "A compact row storage scheme for Cholesky factors", ACM Trans.
00024  * Math. Software, vol 12, 1986, pp. 127-148.
00025  *
00026  * J. Liu, "The role of elimination trees in sparse factorization", SIAM J.
00027  * Matrix Analysis & Applic., vol 11, 1990, pp. 134-172.
00028  *
00029  * J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
00030  * sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
00031  *
00032  * workspace: symmetric: Iwork (nrow), unsymmetric: Iwork (nrow+ncol)
00033  *
00034  * Supports any xtype (pattern, real, complex, or zomplex)
00035  */
00036 
00037 #ifndef NCHOLESKY
00038 
00039 #include "amesos_cholmod_internal.h"
00040 #include "amesos_cholmod_cholesky.h"
00041 
00042 /* ========================================================================== */
00043 /* === update_etree ========================================================= */
00044 /* ========================================================================== */
00045 
00046 static void amesos_update_etree
00047 (
00048     /* inputs, not modified */
00049     Int k,    /* process the edge (k,i) in the input graph */
00050     Int i,
00051     /* inputs, modified on output */
00052     Int Parent [ ], /* Parent [t] = p if p is the parent of t */
00053     Int Ancestor [ ]  /* Ancestor [t] is the ancestor of node t in the
00054          partially-constructed etree */
00055 )
00056 {
00057     Int a ;
00058     for ( ; ; )   /* traverse the path from k to the root of the tree */
00059     {
00060   a = Ancestor [k] ;
00061   if (a == i)
00062   {
00063       /* final ancestor reached; no change to tree */
00064       return ;
00065   }
00066   /* perform path compression */
00067   Ancestor [k] = i ;
00068   if (a == EMPTY)
00069   {
00070       /* final ancestor undefined; this is a new edge in the tree */
00071       Parent [k] = i ;
00072       return ;
00073   }
00074   /* traverse up to the ancestor of k */
00075   k = a ;
00076     }
00077 }
00078 
00079 /* ========================================================================== */
00080 /* === cholmod_etree ======================================================== */
00081 /* ========================================================================== */
00082 
00083 /* Find the elimination tree of A or A'*A */
00084 
00085 int CHOLMOD(etree)
00086 (
00087     /* ---- input ---- */
00088     cholmod_sparse *A,
00089     /* ---- output --- */
00090     Int *Parent,  /* size ncol.  Parent [j] = p if p is the parent of j */
00091     /* --------------- */
00092     cholmod_common *Common
00093 )
00094 {
00095     Int *Ap, *Ai, *Anz, *Ancestor, *Prev, *Iwork ;
00096     Int i, j, jprev, p, pend, nrow, ncol, packed, stype ;
00097     size_t s ;
00098     int ok = TRUE ;
00099 
00100     /* ---------------------------------------------------------------------- */
00101     /* check inputs */
00102     /* ---------------------------------------------------------------------- */
00103 
00104     RETURN_IF_NULL_COMMON (FALSE) ;
00105     RETURN_IF_NULL (A, FALSE) ;
00106     RETURN_IF_NULL (Parent, FALSE) ;
00107     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00108     Common->status = CHOLMOD_OK ;
00109 
00110     /* ---------------------------------------------------------------------- */
00111     /* allocate workspace */
00112     /* ---------------------------------------------------------------------- */
00113 
00114     stype = A->stype ;
00115 
00116     /* s = A->nrow + (stype ? 0 : A->ncol) */
00117     s = CHOLMOD(add_size_t) (A->nrow, (stype ? 0 : A->ncol), &ok) ;
00118     if (!ok)
00119     {
00120   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00121   return (FALSE) ;
00122     }
00123 
00124     CHOLMOD(allocate_work) (0, s, 0, Common) ;
00125     if (Common->status < CHOLMOD_OK)
00126     {
00127   return (FALSE) ;  /* out of memory */
00128     }
00129 
00130     ASSERT (CHOLMOD(dump_sparse) (A, "etree", Common) >= 0) ;
00131     Iwork = Common->Iwork ;
00132 
00133     /* ---------------------------------------------------------------------- */
00134     /* get inputs */
00135     /* ---------------------------------------------------------------------- */
00136 
00137     ncol = A->ncol ;  /* the number of columns of A */
00138     nrow = A->nrow ;  /* the number of rows of A */
00139     Ap = A->p ;   /* size ncol+1, column pointers for A */
00140     Ai = A->i ;   /* the row indices of A */
00141     Anz = A->nz ; /* number of nonzeros in each column of A */
00142     packed = A->packed ;
00143     Ancestor = Iwork ;  /* size ncol (i/i/l) */
00144 
00145     for (j = 0 ; j < ncol ; j++)
00146     {
00147   Parent [j] = EMPTY ;
00148   Ancestor [j] = EMPTY ;
00149     }
00150 
00151     /* ---------------------------------------------------------------------- */
00152     /* compute the etree */
00153     /* ---------------------------------------------------------------------- */
00154 
00155     if (stype > 0)
00156     {
00157 
00158   /* ------------------------------------------------------------------ */
00159   /* symmetric (upper) case: compute etree (A) */
00160   /* ------------------------------------------------------------------ */
00161 
00162   for (j = 0 ; j < ncol ; j++)
00163   {
00164       /* for each row i in column j of triu(A), excluding the diagonal */
00165       p = Ap [j] ;
00166       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00167       for ( ; p < pend ; p++)
00168       {
00169     i = Ai [p] ;
00170     if (i < j)
00171     {
00172         amesos_update_etree (i, j, Parent, Ancestor) ;
00173     }
00174       }
00175   }
00176 
00177     }
00178     else if (stype == 0)
00179     {
00180 
00181   /* ------------------------------------------------------------------ */
00182   /* unsymmetric case: compute etree (A'*A) */
00183   /* ------------------------------------------------------------------ */
00184 
00185   Prev = Iwork + ncol ; /* size nrow (i/i/l) */
00186   for (i = 0 ; i < nrow ; i++)
00187   {
00188       Prev [i] = EMPTY ;
00189   }
00190   for (j = 0 ; j < ncol ; j++)
00191   {
00192       /* for each row i in column j of A */
00193       p = Ap [j] ;
00194       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00195       for ( ; p < pend ; p++)
00196       {
00197     /* a graph is constructed dynamically with one path per row
00198      * of A.  If the ith row of A contains column indices
00199      * (j1,j2,j3,j4) then the new graph has edges (j1,j2), (j2,j3),
00200      * and (j3,j4).  When at node i of this path-graph, all edges
00201      * (jprev,j) are considered, where jprev<j */
00202     i = Ai [p] ;
00203     jprev = Prev [i] ;
00204     if (jprev != EMPTY)
00205     {
00206         amesos_update_etree (jprev, j, Parent, Ancestor) ;
00207     }
00208     Prev [i] = j ;
00209       }
00210   }
00211 
00212     }
00213     else
00214     {
00215 
00216   /* ------------------------------------------------------------------ */
00217   /* symmetric case with lower triangular part not supported */
00218   /* ------------------------------------------------------------------ */
00219 
00220   ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
00221   return (FALSE) ;
00222     }
00223 
00224     ASSERT (CHOLMOD(dump_parent) (Parent, ncol, "Parent", Common)) ;
00225     return (TRUE) ;
00226 }
00227 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines