Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_amd.c
Go to the documentation of this file.
```00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_amd ================================================= */
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 /* CHOLMOD interface to the AMD 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_amd 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 (6*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 AMD.
00030  *
00031  * Supports any xtype (pattern, real, complex, or zomplex)
00032  */
00033
00034 #ifndef NCHOLESKY
00035
00036 /* This file should make the long int version of CHOLMOD */
00037 #define DLONG 1
00038
00039 #include "amesos_cholmod_internal.h"
00040 #include "amesos_amd.h"
00041 #include "amesos_cholmod_cholesky.h"
00042
00043 #if (!defined (AMD_VERSION) || (AMD_VERSION < AMD_VERSION_CODE (2,0)))
00044 #error "AMD v2.0 or later is required"
00045 #endif
00046
00047 /* ========================================================================== */
00048 /* === cholmod_amd ========================================================== */
00049 /* ========================================================================== */
00050
00051 int CHOLMOD(amd)
00052 (
00053     /* ---- input ---- */
00054     cholmod_sparse *A,  /* matrix to order */
00055     Int *fset,    /* subset of 0:(A->ncol)-1 */
00056     size_t fsize, /* size of fset */
00057     /* ---- output --- */
00058     Int *Perm,    /* size A->nrow, output permutation */
00059     /* --------------- */
00060     cholmod_common *Common
00061 )
00062 {
00063     double Info [AMD_INFO], Control2 [AMD_CONTROL], *Control ;
00064     Int *Cp, *Len, *Nv, *Head, *Elen, *Degree, *Wi, *Iwork, *Next ;
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     RETURN_IF_NULL (Perm, FALSE) ;
00079     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00080     Common->status = CHOLMOD_OK ;
00081     if (n == 0)
00082     {
00083   /* nothing to do */
00084   Common->fl = 0 ;
00085   Common->lnz = 0 ;
00086   Common->anz = 0 ;
00087   return (TRUE) ;
00088     }
00089
00090     /* ---------------------------------------------------------------------- */
00091     /* get workspace */
00092     /* ---------------------------------------------------------------------- */
00093
00094     /* Note: this is less than the space used in cholmod_analyze, so if
00095      * cholmod_amd is being called by that routine, no space will be
00096      * allocated.
00097      */
00098
00099     /* s = MAX (6*n, A->ncol) */
00100     s = CHOLMOD(mult_size_t) (n, 6, &ok) ;
00101     if (!ok)
00102     {
00103   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00104   return (FALSE) ;
00105     }
00106     s = MAX (s, A->ncol) ;
00107
00108     CHOLMOD(allocate_work) (n, s, 0, Common) ;
00109     if (Common->status < CHOLMOD_OK)
00110     {
00111   return (FALSE) ;
00112     }
00113
00114     Iwork  = Common->Iwork ;
00115     Degree = Iwork ;      /* size n */
00116     Wi     = Iwork + n ;    /* size n */
00117     Len    = Iwork + 2*((size_t) n) ; /* size n */
00118     Nv     = Iwork + 3*((size_t) n) ;   /* size n */
00119     Next   = Iwork + 4*((size_t) n) ;   /* size n */
00120     Elen   = Iwork + 5*((size_t) n) ;   /* size n */
00121
00122     Head = Common->Head ;   /* size n+1, but only n is used */
00123
00124     /* ---------------------------------------------------------------------- */
00125     /* construct the input matrix for AMD */
00126     /* ---------------------------------------------------------------------- */
00127
00128     if (A->stype == 0)
00129     {
00130   /* C = A*A' or A(:,f)*A(:,f)', add extra space of nnz(C)/2+n to C */
00131   C = CHOLMOD(aat) (A, fset, fsize, -2, Common) ;
00132     }
00133     else
00134     {
00135   /* C = A+A', but use only the upper triangular part of A if A->stype = 1
00136    * and only the lower part of A if A->stype = -1.  Add extra space of
00137    * nnz(C)/2+n to C. */
00138   C = CHOLMOD(copy) (A, 0, -2, Common) ;
00139     }
00140
00141     if (Common->status < CHOLMOD_OK)
00142     {
00143   /* out of memory, fset invalid, or other error */
00144   return (FALSE) ;
00145     }
00146
00147     Cp = C->p ;
00148     for (j = 0 ; j < n ; j++)
00149     {
00150   Len [j] = Cp [j+1] - Cp [j] ;
00151     }
00152
00153     /* C does not include the diagonal, and both upper and lower parts.
00154      * Common->anz includes the diagonal, and just the lower part of C */
00155     cnz = Cp [n] ;
00156     Common->anz = cnz / 2 + n ;
00157
00158     /* ---------------------------------------------------------------------- */
00159     /* order C using AMD */
00160     /* ---------------------------------------------------------------------- */
00161
00162     /* get parameters */
00163     if (Common->current < 0 || Common->current >= CHOLMOD_MAXMETHODS)
00164     {
00165   /* use AMD defaults */
00166   Control = NULL ;
00167     }
00168     else
00169     {
00170   Control = Control2 ;
00171   Control [AMD_DENSE] = Common->method [Common->current].prune_dense ;
00172   Control [AMD_AGGRESSIVE] = Common->method [Common->current].aggressive ;
00173     }
00174
00175     /* AMD_2 does not use amd_malloc and amd_free, but set these pointers just
00176      * be safe. */
00177     amesos_amd_malloc = Common->malloc_memory ;
00178     amesos_amd_free = Common->free_memory ;
00179     amesos_amd_calloc = Common->calloc_memory ;
00180     amesos_amd_realloc = Common->realloc_memory ;
00181
00182     /* AMD_2 doesn't print anything either, but future versions might,
00183      * so set the amd_printf pointer too. */
00184     amesos_amd_printf = Common->print_function ;
00185
00186 #ifdef LONG
00187     amesos_amd_l2 (n, C->p,  C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen,
00188       Degree, Wi, Control, Info) ;
00189 #else
00190     amesos_amd_2 (n, C->p,  C->i, Len, C->nzmax, cnz, Nv, Next, Perm, Head, Elen,
00191       Degree, Wi, Control, Info) ;
00192 #endif
00193
00194     /* LL' flop count.  Need to subtract n for LL' flop count.  Note that this
00195      * is a slight upper bound which is often exact (see AMD/Source/amd_2.c for
00196      * details).  cholmod_analyze computes an exact flop count and fill-in. */
00197     Common->fl = Info [AMD_NDIV] + 2 * Info [AMD_NMULTSUBS_LDL] + n ;
00198
00199     /* Info [AMD_LNZ] excludes the diagonal */
00200     Common->lnz = n + Info [AMD_LNZ] ;
00201
00202     /* ---------------------------------------------------------------------- */
00203     /* free the AMD workspace and clear the persistent workspace in Common */
00204     /* ---------------------------------------------------------------------- */
00205
00206     ASSERT (IMPLIES (Common->status == CHOLMOD_OK,
00207     CHOLMOD(dump_perm) (Perm, n, n, "AMD2 perm", Common))) ;
00208     CHOLMOD(free_sparse) (&C, Common) ;
00209     for (j = 0 ; j <= n ; j++)
00210     {
00211   Head [j] = EMPTY ;
00212     }
00213     return (TRUE) ;
00214 }
00215 #endif
```