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