Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_camd.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Partition/cholmod_camd =============================================== */
00003 /* ========================================================================== */
00004 
00005 /* -----------------------------------------------------------------------------
00006  * CHOLMOD/Partition Module.  Copyright (C) 2005-2006, Timothy A. Davis
00007  * The CHOLMOD/Partition 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 /* CHOLMOD interface to the CAMD ordering routine.  Orders A if the matrix is
00014  * symmetric.  On output, Perm [k] = i if row/column i of A is the kth
00015  * row/column of P*A*P'.  This corresponds to A(p,p) in MATLAB notation.
00016  *
00017  * If A is unsymmetric, cholmod_camd orders A*A'.  On output, Perm [k] = i if
00018  * row/column i of A*A' is the kth row/column of P*A*A'*P'.  This corresponds to
00019  * A(p,:)*A(p,:)' in MATLAB notation.  If f is present, A(p,f)*A(p,f)' is
00020  * ordered.
00021  *
00022  * Computes the flop count for a subsequent LL' factorization, the number
00023  * of nonzeros in L, and the number of nonzeros in the matrix ordered (A,
00024  * A*A' or A(:,f)*A(:,f)').
00025  *
00026  * workspace: Iwork (4*nrow). Head (nrow).
00027  *
00028  * Allocates a temporary copy of A+A' or A*A' (with
00029  * both upper and lower triangular parts) as input to CAMD.
00030  * Also allocates 3*(n+1) additional integer workspace (not in Common).
00031  *
00032  * Supports any xtype (pattern, real, complex, or zomplex)
00033  */
00034 
00035 #ifndef NPARTITION
00036 
00037 #include "amesos_cholmod_internal.h"
00038 #include "amesos_camd.h"
00039 #include "amesos_cholmod_partition.h"
00040 
00041 #if (CAMD_VERSION < CAMD_VERSION_CODE (2,0))
00042 #error "CAMD v2.0 or later is required"
00043 #endif
00044 
00045 /* ========================================================================== */
00046 /* === cholmod_camd ========================================================= */
00047 /* ========================================================================== */
00048 
00049 int CHOLMOD(camd)
00050 (
00051     /* ---- input ---- */
00052     cholmod_sparse *A,  /* matrix to order */
00053     Int *fset,    /* subset of 0:(A->ncol)-1 */
00054     size_t fsize, /* size of fset */
00055     Int *Cmember, /* size nrow.  see cholmod_ccolamd.c for description.*/
00056     /* ---- output ---- */
00057     Int *Perm,    /* size A->nrow, output permutation */
00058     /* --------------- */
00059     cholmod_common *Common
00060 )
00061 {
00062     double Info [CAMD_INFO], Control2 [CAMD_CONTROL], *Control ;
00063     Int *Cp, *Len, *Nv, *Head, *Elen, *Degree, *Wi, *Next, *BucketSet,
00064   *Work3n, *p ;
00065     cholmod_sparse *C ;
00066     Int j, n, cnz ;
00067     size_t s ;
00068     int ok = TRUE ;
00069 
00070     /* ---------------------------------------------------------------------- */
00071     /* get inputs */
00072     /* ---------------------------------------------------------------------- */
00073 
00074     RETURN_IF_NULL_COMMON (FALSE) ;
00075     RETURN_IF_NULL (A, FALSE) ;
00076     n = A->nrow ;
00077 
00078     /* s = 4*n */
00079     s = CHOLMOD(mult_size_t) (n, 4, &ok) ;
00080     if (!ok)
00081     {
00082   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00083   return (FALSE) ;
00084     }
00085 
00086     RETURN_IF_NULL (Perm, FALSE) ;
00087     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00088     Common->status = CHOLMOD_OK ;
00089     if (n == 0)
00090     {
00091   /* nothing to do */
00092   Common->fl = 0 ;
00093   Common->lnz = 0 ;
00094   Common->anz = 0 ;
00095   return (TRUE) ;
00096     }
00097 
00098     /* ---------------------------------------------------------------------- */
00099     /* get workspace */
00100     /* ---------------------------------------------------------------------- */
00101 
00102     /* cholmod_analyze has allocated Cmember at Common->Iwork + 5*n+uncol, and
00103      * CParent at Common->Iwork + 4*n+uncol, where uncol is 0 if A is symmetric
00104      * or A->ncol otherwise.  Thus, only the first 4n integers in Common->Iwork
00105      * can be used here. */
00106 
00107     CHOLMOD(allocate_work) (n, s, 0, Common) ;
00108     if (Common->status < CHOLMOD_OK)
00109     {
00110   return (FALSE) ;
00111     }
00112 
00113     p = Common->Iwork ;
00114     Degree = p ; p += n ; /* size n */
00115     Elen   = p ; p += n ; /* size n */
00116     Len    = p ; p += n ; /* size n */
00117     Nv     = p ; p += n ; /* size n */
00118 
00119     Work3n = CHOLMOD(malloc) (n+1, 3*sizeof (Int), Common) ;
00120     if (Common->status < CHOLMOD_OK)
00121     {
00122   return (FALSE) ;
00123     }
00124     p = Work3n ;
00125     Next = p ; p += n ;   /* size n */
00126     Wi   = p ; p += (n+1) ; /* size n+1 */
00127     BucketSet = p ;   /* size n */
00128 
00129     Head = Common->Head ; /* size n+1 */
00130 
00131     /* ---------------------------------------------------------------------- */
00132     /* construct the input matrix for CAMD */
00133     /* ---------------------------------------------------------------------- */
00134 
00135     if (A->stype == 0)
00136     {
00137   /* C = A*A' or A(:,f)*A(:,f)', add extra space of nnz(C)/2+n to C */
00138   C = CHOLMOD(aat) (A, fset, fsize, -2, Common) ;
00139     }
00140     else
00141     {
00142   /* C = A+A', but use only the upper triangular part of A if A->stype = 1
00143    * and only the lower part of A if A->stype = -1.  Add extra space of
00144    * nnz(C)/2+n to C. */
00145   C = CHOLMOD(copy) (A, 0, -2, Common) ;
00146     }
00147 
00148     if (Common->status < CHOLMOD_OK)
00149     {
00150   /* out of memory, fset invalid, or other error */
00151   CHOLMOD(free) (n+1, 3*sizeof (Int), Work3n, Common) ;
00152   return (FALSE) ;
00153     }
00154 
00155     Cp = C->p ;
00156     for (j = 0 ; j < n ; j++)
00157     {
00158   Len [j] = Cp [j+1] - Cp [j] ;
00159     }
00160 
00161     /* C does not include the diagonal, and both upper and lower parts.
00162      * Common->anz includes the diagonal, and just the lower part of C */
00163     cnz = Cp [n] ;
00164     Common->anz = cnz / 2 + n ;
00165 
00166     /* ---------------------------------------------------------------------- */
00167     /* order C using CAMD */
00168     /* ---------------------------------------------------------------------- */
00169 
00170     /* get parameters */
00171     if (Common->current < 0 || Common->current >= CHOLMOD_MAXMETHODS)
00172     {
00173   /* use CAMD defaults */
00174   Control = NULL ;
00175     }
00176     else
00177     {
00178   Control = Control2 ;
00179   Control [CAMD_DENSE] = Common->method [Common->current].prune_dense ;
00180   Control [CAMD_AGGRESSIVE] = Common->method [Common->current].aggressive;
00181     }
00182 
00183     /* CAMD_2 does not use camd_malloc and camd_free, but set these pointers
00184      * just be safe. */
00185     amesos_camd_malloc = Common->malloc_memory ;
00186     amesos_camd_free = Common->free_memory ;
00187     amesos_camd_calloc = Common->calloc_memory ;
00188     amesos_camd_realloc = Common->realloc_memory ;
00189 
00190     /* CAMD_2 doesn't print anything either, but future versions might,
00191      * so set the camd_printf pointer too. */
00192     amesos_camd_printf = Common->print_function ;
00193 
00194 #ifdef LONG
00195     /* DEBUG (camd_l_debug_init ("cholmod_l_camd")) ; */
00196     amesos_camd_l2 (n, C->p,  C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen,
00197       Degree, Wi, Control, Info, Cmember, BucketSet) ;
00198 #else
00199     /* DEBUG (camd_debug_init ("cholmod_camd")) ; */
00200     amesos_camd_2 (n, C->p,  C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen,
00201       Degree, Wi, Control, Info, Cmember, BucketSet) ;
00202 #endif
00203 
00204     /* LL' flop count.  Need to subtract n for LL' flop count.  Note that this
00205      * is a slight upper bound which is often exact (see CAMD/Source/camd_2.c
00206      * for details).  cholmod_analyze computes an exact flop count and
00207      * fill-in. */
00208     Common->fl = Info [CAMD_NDIV] + 2 * Info [CAMD_NMULTSUBS_LDL] + n ;
00209 
00210     /* Info [CAMD_LNZ] excludes the diagonal */
00211     Common->lnz = n + Info [CAMD_LNZ] ;
00212 
00213     /* ---------------------------------------------------------------------- */
00214     /* free the CAMD workspace and clear the persistent workspace in Common */
00215     /* ---------------------------------------------------------------------- */
00216 
00217     ASSERT (IMPLIES (Common->status == CHOLMOD_OK,
00218     CHOLMOD(dump_perm) (Perm, n, n, "CAMD2 perm", Common))) ;
00219     CHOLMOD(free_sparse) (&C, Common) ;
00220     for (j = 0 ; j <= n ; j++)
00221     {
00222   Head [j] = EMPTY ;
00223     }
00224     CHOLMOD(free) (n+1, 3*sizeof (Int), Work3n, Common) ;
00225     return (TRUE) ;
00226 }
00227 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines