Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_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 /* This file should make the long int version of CHOLMOD */
00040 #define DLONG 1
00041 
00042 #include "amesos_cholmod_internal.h"
00043 #include "amesos_cholmod_cholesky.h"
00044 
00045 /* ========================================================================== */
00046 /* === update_etree ========================================================= */
00047 /* ========================================================================== */
00048 
00049 static void amesos_update_etree
00050 (
00051     /* inputs, not modified */
00052     Int k,    /* process the edge (k,i) in the input graph */
00053     Int i,
00054     /* inputs, modified on output */
00055     Int Parent [ ], /* Parent [t] = p if p is the parent of t */
00056     Int Ancestor [ ]  /* Ancestor [t] is the ancestor of node t in the
00057          partially-constructed etree */
00058 )
00059 {
00060     Int a ;
00061     for ( ; ; )   /* traverse the path from k to the root of the tree */
00062     {
00063   a = Ancestor [k] ;
00064   if (a == i)
00065   {
00066       /* final ancestor reached; no change to tree */
00067       return ;
00068   }
00069   /* perform path compression */
00070   Ancestor [k] = i ;
00071   if (a == EMPTY)
00072   {
00073       /* final ancestor undefined; this is a new edge in the tree */
00074       Parent [k] = i ;
00075       return ;
00076   }
00077   /* traverse up to the ancestor of k */
00078   k = a ;
00079     }
00080 }
00081 
00082 /* ========================================================================== */
00083 /* === cholmod_etree ======================================================== */
00084 /* ========================================================================== */
00085 
00086 /* Find the elimination tree of A or A'*A */
00087 
00088 int CHOLMOD(etree)
00089 (
00090     /* ---- input ---- */
00091     cholmod_sparse *A,
00092     /* ---- output --- */
00093     Int *Parent,  /* size ncol.  Parent [j] = p if p is the parent of j */
00094     /* --------------- */
00095     cholmod_common *Common
00096 )
00097 {
00098     Int *Ap, *Ai, *Anz, *Ancestor, *Prev, *Iwork ;
00099     Int i, j, jprev, p, pend, nrow, ncol, packed, stype ;
00100     size_t s ;
00101     int ok = TRUE ;
00102 
00103     /* ---------------------------------------------------------------------- */
00104     /* check inputs */
00105     /* ---------------------------------------------------------------------- */
00106 
00107     RETURN_IF_NULL_COMMON (FALSE) ;
00108     RETURN_IF_NULL (A, FALSE) ;
00109     RETURN_IF_NULL (Parent, FALSE) ;
00110     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00111     Common->status = CHOLMOD_OK ;
00112 
00113     /* ---------------------------------------------------------------------- */
00114     /* allocate workspace */
00115     /* ---------------------------------------------------------------------- */
00116 
00117     stype = A->stype ;
00118 
00119     /* s = A->nrow + (stype ? 0 : A->ncol) */
00120     s = CHOLMOD(add_size_t) (A->nrow, (stype ? 0 : A->ncol), &ok) ;
00121     if (!ok)
00122     {
00123   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00124   return (FALSE) ;
00125     }
00126 
00127     CHOLMOD(allocate_work) (0, s, 0, Common) ;
00128     if (Common->status < CHOLMOD_OK)
00129     {
00130   return (FALSE) ;  /* out of memory */
00131     }
00132 
00133     ASSERT (CHOLMOD(dump_sparse) (A, "etree", Common) >= 0) ;
00134     Iwork = Common->Iwork ;
00135 
00136     /* ---------------------------------------------------------------------- */
00137     /* get inputs */
00138     /* ---------------------------------------------------------------------- */
00139 
00140     ncol = A->ncol ;  /* the number of columns of A */
00141     nrow = A->nrow ;  /* the number of rows of A */
00142     Ap = A->p ;   /* size ncol+1, column pointers for A */
00143     Ai = A->i ;   /* the row indices of A */
00144     Anz = A->nz ; /* number of nonzeros in each column of A */
00145     packed = A->packed ;
00146     Ancestor = Iwork ;  /* size ncol (i/i/l) */
00147 
00148     for (j = 0 ; j < ncol ; j++)
00149     {
00150   Parent [j] = EMPTY ;
00151   Ancestor [j] = EMPTY ;
00152     }
00153 
00154     /* ---------------------------------------------------------------------- */
00155     /* compute the etree */
00156     /* ---------------------------------------------------------------------- */
00157 
00158     if (stype > 0)
00159     {
00160 
00161   /* ------------------------------------------------------------------ */
00162   /* symmetric (upper) case: compute etree (A) */
00163   /* ------------------------------------------------------------------ */
00164 
00165   for (j = 0 ; j < ncol ; j++)
00166   {
00167       /* for each row i in column j of triu(A), excluding the diagonal */
00168       p = Ap [j] ;
00169       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00170       for ( ; p < pend ; p++)
00171       {
00172     i = Ai [p] ;
00173     if (i < j)
00174     {
00175         amesos_update_etree (i, j, Parent, Ancestor) ;
00176     }
00177       }
00178   }
00179 
00180     }
00181     else if (stype == 0)
00182     {
00183 
00184   /* ------------------------------------------------------------------ */
00185   /* unsymmetric case: compute etree (A'*A) */
00186   /* ------------------------------------------------------------------ */
00187 
00188   Prev = Iwork + ncol ; /* size nrow (i/i/l) */
00189   for (i = 0 ; i < nrow ; i++)
00190   {
00191       Prev [i] = EMPTY ;
00192   }
00193   for (j = 0 ; j < ncol ; j++)
00194   {
00195       /* for each row i in column j of A */
00196       p = Ap [j] ;
00197       pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
00198       for ( ; p < pend ; p++)
00199       {
00200     /* a graph is constructed dynamically with one path per row
00201      * of A.  If the ith row of A contains column indices
00202      * (j1,j2,j3,j4) then the new graph has edges (j1,j2), (j2,j3),
00203      * and (j3,j4).  When at node i of this path-graph, all edges
00204      * (jprev,j) are considered, where jprev<j */
00205     i = Ai [p] ;
00206     jprev = Prev [i] ;
00207     if (jprev != EMPTY)
00208     {
00209         amesos_update_etree (jprev, j, Parent, Ancestor) ;
00210     }
00211     Prev [i] = j ;
00212       }
00213   }
00214 
00215     }
00216     else
00217     {
00218 
00219   /* ------------------------------------------------------------------ */
00220   /* symmetric case with lower triangular part not supported */
00221   /* ------------------------------------------------------------------ */
00222 
00223   ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
00224   return (FALSE) ;
00225     }
00226 
00227     ASSERT (CHOLMOD(dump_parent) (Parent, ncol, "Parent", Common)) ;
00228     return (TRUE) ;
00229 }
00230 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines