Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_colamd.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_colamd ============================================== */
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 /* CHOLMOD interface to the COLAMD ordering routine (version 2.4 or later).
00014  * Finds a permutation p such that the Cholesky factorization of PAA'P' is
00015  * sparser than AA' using colamd.  If the postorder input parameter is TRUE,
00016  * the column etree is found and postordered, and the colamd ordering is then
00017  * combined with its postordering.  A must be unsymmetric.
00018  *
00019  * There can be no duplicate entries in f.
00020  * f can be length 0 to n if A is m-by-n.
00021  *
00022  * workspace: Iwork (4*nrow+ncol), Head (nrow+1), Flag (nrow)
00023  *  Allocates a copy of its input matrix, which
00024  *  is then used as CCOLAMD's workspace.
00025  *
00026  * Supports any xtype (pattern, real, complex, or zomplex)
00027  */
00028 
00029 #ifndef NCHOLESKY
00030 
00031 /* This file should make the long int version of CHOLMOD */
00032 #define DLONG 1
00033 
00034 #include "amesos_cholmod_internal.h"
00035 #include "amesos_colamd.h"
00036 #include "amesos_cholmod_cholesky.h"
00037 
00038 #if (!defined (COLAMD_VERSION) || (COLAMD_VERSION < COLAMD_VERSION_CODE (2,5)))
00039 #error "COLAMD v2.5 or later is required"
00040 #endif
00041 
00042 /* ========================================================================== */
00043 /* === cholmod_colamd ======================================================= */
00044 /* ========================================================================== */
00045 
00046 int CHOLMOD(colamd)
00047 (
00048     /* ---- input ---- */
00049     cholmod_sparse *A,  /* matrix to order */
00050     Int *fset,    /* subset of 0:(A->ncol)-1 */
00051     size_t fsize, /* size of fset */
00052     int postorder,  /* if TRUE, follow with a coletree postorder */
00053     /* ---- output --- */
00054     Int *Perm,    /* size A->nrow, output permutation */
00055     /* --------------- */
00056     cholmod_common *Common
00057 )
00058 {
00059     double knobs [COLAMD_KNOBS] ;
00060     cholmod_sparse *C ;
00061     Int *NewPerm, *Parent, *Post, *Work2n ;
00062     Int k, nrow, ncol ;
00063     size_t s, alen ;
00064     int ok = TRUE ;
00065 
00066     /* ---------------------------------------------------------------------- */
00067     /* check inputs */
00068     /* ---------------------------------------------------------------------- */
00069 
00070     RETURN_IF_NULL_COMMON (FALSE) ;
00071     RETURN_IF_NULL (A, FALSE) ;
00072     RETURN_IF_NULL (Perm, FALSE) ;
00073     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
00074     if (A->stype != 0)
00075     {
00076   ERROR (CHOLMOD_INVALID, "matrix must be unsymmetric") ;
00077   return (FALSE) ;
00078     }
00079     Common->status = CHOLMOD_OK ;
00080 
00081     /* ---------------------------------------------------------------------- */
00082     /* get inputs */
00083     /* ---------------------------------------------------------------------- */
00084 
00085     nrow = A->nrow ;
00086     ncol = A->ncol ;
00087 
00088     /* ---------------------------------------------------------------------- */
00089     /* allocate workspace */
00090     /* ---------------------------------------------------------------------- */
00091 
00092     /* Note: this is less than the space used in cholmod_analyze, so if
00093      * cholmod_colamd is being called by that routine, no space will be
00094      * allocated.
00095      */
00096 
00097     /* s = 4*nrow + ncol */
00098     s = CHOLMOD(mult_size_t) (nrow, 4, &ok) ;
00099     s = CHOLMOD(add_size_t) (s, ncol, &ok) ;
00100 
00101 #ifdef LONG
00102     alen = amesos_colamd_l_recommended (A->nzmax, ncol, nrow) ;
00103     amesos_colamd_l_set_defaults (knobs) ;
00104 #else
00105     alen = amesos_colamd_recommended (A->nzmax, ncol, nrow) ;
00106     amesos_colamd_set_defaults (knobs) ;
00107 #endif
00108 
00109     if (!ok || alen == 0)
00110     {
00111   ERROR (CHOLMOD_TOO_LARGE, "matrix invalid or too large") ;
00112   return (FALSE) ;
00113     }
00114 
00115     CHOLMOD(allocate_work) (0, s, 0, Common) ;
00116     if (Common->status < CHOLMOD_OK)
00117     {
00118   return (FALSE) ;
00119     }
00120 
00121     /* ---------------------------------------------------------------------- */
00122     /* allocate COLAMD workspace */
00123     /* ---------------------------------------------------------------------- */
00124 
00125     /* colamd_printf is only available in colamd v2.4 or later */
00126     amesos_colamd_printf = Common->print_function ;
00127 
00128     C = CHOLMOD(allocate_sparse) (ncol, nrow, alen, TRUE, TRUE, 0,
00129       CHOLMOD_PATTERN, Common) ;
00130 
00131     /* ---------------------------------------------------------------------- */
00132     /* copy (and transpose) the input matrix A into the colamd workspace */
00133     /* ---------------------------------------------------------------------- */
00134 
00135     /* C = A (:,f)', which also packs A if needed. */
00136     /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset) */
00137     ok = CHOLMOD(transpose_unsym) (A, 0, NULL, fset, fsize, C, Common) ;
00138 
00139     /* ---------------------------------------------------------------------- */
00140     /* order the matrix (destroys the contents of C->i and C->p) */
00141     /* ---------------------------------------------------------------------- */
00142 
00143     /* get parameters */
00144     if (Common->current < 0 || Common->current >= CHOLMOD_MAXMETHODS)
00145     {
00146   /* this is the CHOLMOD default, not the COLAMD default */
00147   knobs [COLAMD_DENSE_ROW] = -1 ;
00148     }
00149     else
00150     {
00151   /* get the knobs from the Common parameters */
00152   knobs [COLAMD_DENSE_COL] = Common->method[Common->current].prune_dense ;
00153   knobs [COLAMD_DENSE_ROW] = Common->method[Common->current].prune_dense2;
00154   knobs [COLAMD_AGGRESSIVE] = Common->method[Common->current].aggressive ;
00155     }
00156 
00157     if (ok)
00158     {
00159   Int *Cp ;
00160   Int stats [COLAMD_STATS] ;
00161   Cp = C->p ;
00162 
00163 #ifdef LONG
00164   amesos_colamd_l (ncol, nrow, alen, C->i, Cp, knobs, stats) ;
00165 #else
00166   amesos_colamd (ncol, nrow, alen, C->i, Cp, knobs, stats) ;
00167 #endif
00168 
00169   ok = stats [COLAMD_STATUS] ;
00170   ok = (ok == COLAMD_OK || ok == COLAMD_OK_BUT_JUMBLED) ;
00171   /* permutation returned in C->p, if the ordering succeeded */
00172   for (k = 0 ; k < nrow ; k++)
00173   {
00174       Perm [k] = Cp [k] ;
00175   }
00176     }
00177 
00178     CHOLMOD(free_sparse) (&C, Common) ;
00179 
00180     /* ---------------------------------------------------------------------- */
00181     /* column etree postordering */
00182     /* ---------------------------------------------------------------------- */
00183 
00184     if (postorder)
00185     {
00186   /* use the last 2*n space in Iwork for Parent and Post */
00187   Work2n = Common->Iwork ;
00188   Work2n += 2*((size_t) nrow) + ncol ;
00189   Parent = Work2n ;   /* size nrow (i/i/l) */
00190   Post   = Work2n + nrow ;  /* size nrow (i/i/l) */
00191 
00192   /* workspace: Iwork (2*nrow+ncol), Flag (nrow), Head (nrow+1) */
00193   ok = ok && CHOLMOD(analyze_ordering) (A, CHOLMOD_COLAMD, Perm, fset,
00194     fsize, Parent, Post, NULL, NULL, NULL, Common) ;
00195 
00196   /* combine the colamd permutation with its postordering */
00197   if (ok)
00198   {
00199       NewPerm = Common->Iwork ;   /* size nrow (i/i/l) */
00200       for (k = 0 ; k < nrow ; k++)
00201       {
00202     NewPerm [k] = Perm [Post [k]] ;
00203       }
00204       for (k = 0 ; k < nrow ; k++)
00205       {
00206     Perm [k] = NewPerm [k] ;
00207       }
00208   }
00209     }
00210 
00211     return (ok) ;
00212 }
00213 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines