Amesos Package Browser (Single Doxygen Collection) Development
amesos_cholmod_l_analyze.c
Go to the documentation of this file.
00001 /* ========================================================================== */
00002 /* === Cholesky/cholmod_analyze ============================================= */
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 /* Order and analyze a matrix (either simplicial or supernodal), in prepartion
00014  * for numerical factorization via cholmod_factorize or via the "expert"
00015  * routines cholmod_rowfac and cholmod_super_numeric.
00016  *
00017  * symmetric case:    A or A(p,p)
00018  * unsymmetric case:  AA', A(p,:)*A(p,:)', A(:,f)*A(:,f)', or A(p,f)*A(p,f)'
00019  *
00020  * For the symmetric case, only the upper or lower triangular part of A is
00021  * accessed (depending on the type of A).  LL'=A (or permuted A) is analzed.
00022  * For the unsymmetric case (LL'=AA' or permuted A).
00023  *
00024  * There can be no duplicate entries in p or f.  p is of length m if A is
00025  * m-by-n.  f can be length 0 to n.
00026  *
00027  * In both cases, the columns of A need not be sorted.  A can be in packed
00028  * or unpacked form.
00029  *
00030  * Ordering options include:
00031  *
00032  *  natural:    A is not permuted to reduce fill-in
00033  *  given:      a permutation can be provided to this routine (UserPerm)
00034  *  AMD:      approximate minumum degree (AMD for the symmetric case,
00035  *        COLAMD for the AA' case).
00036  *  METIS:      nested dissection with METIS_NodeND
00037  *  NESDIS:     nested dissection using METIS_NodeComputeSeparator,
00038  *        typically followed by a constrained minimum degree
00039  *        (CAMD for the symmetric case, CCOLAMD for the AA' case).
00040  *
00041  * Multiple ordering options can be tried (up to 9 of them), and the best one
00042  * is selected (the one that gives the smallest number of nonzeros in the
00043  * simplicial factor L).  If one method fails, cholmod_analyze keeps going, and
00044  * picks the best among the methods that succeeded.  This routine fails (and
00045  * returns NULL) if either initial memory allocation fails, all ordering methods
00046  * fail, or the supernodal analysis (if requested) fails.  By default, the 9
00047  * methods available are:
00048  *
00049  *  1) given permutation (skipped if UserPerm is NULL)
00050  *  2) AMD (symmetric case) or COLAMD (unsymmetric case)
00051  *  3) METIS with default parameters
00052  *  4) NESDIS with default parameters (stopping the partitioning when
00053  *      the graph is of size nd_small = 200 or less, remove nodes with
00054  *      more than max (16, prune_dense * sqrt (n)) nodes where
00055  *      prune_dense = 10, and follow partitioning with CCOLAMD, a
00056  *      constrained minimum degree ordering).
00057  *  5) natural
00058  *  6) NESDIS, nd_small = 20000, prune_dense = 10
00059  *  7) NESDIS, nd_small =     4, prune_dense = 10, no min degree
00060  *  8) NESDIS, nd_small =   200, prune_dense = 0
00061  *  9) COLAMD for A*A' or AMD for A
00062  *
00063  * By default, the first two are tried, and METIS is tried if AMD reports a high
00064  * flop count and fill-in.  Let fl denote the flop count for the AMD, ordering,
00065  * nnz(L) the # of nonzeros in L, and nnz(tril(A)) (or A*A').  If
00066  * fl/nnz(L) >= 500 and nnz(L)/nnz(tril(A)) >= 5, then METIS is attempted.  The
00067  * best ordering is used (UserPerm if given, AMD, and METIS if attempted).  If
00068  * you do not have METIS, only the first two will be tried (user permutation,
00069  * if provided, and AMD/COLAMD).  This default behavior is obtained when
00070  * Common->nmethods is zero.  In this case, methods 0, 1, and 2 in
00071  * Common->method [..] are reset to User-provided, AMD, and METIS (or NESDIS
00072  * if Common->default_nesdis is set to the non-default value of TRUE),
00073  * respectively.
00074  *
00075  * You can modify these 9 methods and the number of methods tried by changing
00076  * parameters in the Common argument.  If you know the best ordering for your
00077  * matrix, set Common->nmethods to 1 and set Common->method[0].ordering to the
00078  * requested ordering method.  Parameters for each method can also be modified
00079  * (refer to cholmod.h for details).
00080  *
00081  * Note that it is possible for METIS to terminate your program if it runs out
00082  * of memory.  This is not the case for any CHOLMOD or minimum degree ordering
00083  * routine (AMD, COLAMD, CAMD, CCOLAMD, or CSYMAMD).  Since NESDIS relies on
00084  * METIS, it too can terminate your program.
00085  *
00086  * The factor L is returned as simplicial symbolic (L->is_super FALSE) if
00087  * Common->supernodal <= CHOLMOD_SIMPLICIAL (0) or as supernodal symbolic if
00088  * Common->supernodal >= CHOLMOD_SUPERNODAL (2).  If Common->supernodal is
00089  * equal to CHOLMOD_AUTO (1), then L is simplicial if the flop count per
00090  * nonzero in L is less than Common->supernodal_switch (default: 40), and
00091  * is returned as a supernodal factor otherwise.
00092  *
00093  * In both cases, L->xtype is CHOLMOD_PATTERN.
00094  * A subsequent call to cholmod_factorize will perform a
00095  * simplicial or supernodal factorization, depending on the type of L.
00096  *
00097  * For the simplicial case, L contains the fill-reducing permutation (L->Perm)
00098  * and the counts of nonzeros in each column of L (L->ColCount).  For the
00099  * supernodal case, L also contains the nonzero pattern of each supernode.
00100  *
00101  * workspace: Flag (nrow), Head (nrow+1)
00102  *  if symmetric:   Iwork (6*nrow)
00103  *  if unsymmetric: Iwork (6*nrow+ncol).
00104  *  calls various ordering routines, which typically allocate O(nnz(A))
00105  *  temporary workspace ((2 to 3)*nnz(A) * sizeof (Int) is typical, but it
00106  *  can be much higher if A*A' must be explicitly formed for METIS).  Also
00107  *  allocates up to 2 temporary (permuted/transpose) copies of the nonzero
00108  *  pattern of A, and up to 3*n*sizeof(Int) additional workspace.
00109  *
00110  * Supports any xtype (pattern, real, complex, or zomplex)
00111  */
00112 
00113 #ifndef NCHOLESKY
00114 
00115 /* This file should make the long int version of CHOLMOD */
00116 #define DLONG 1
00117 
00118 #include "amesos_cholmod_internal.h"
00119 #include "amesos_cholmod_cholesky.h"
00120 
00121 #ifndef NPARTITION
00122 #include "amesos_cholmod_partition.h"
00123 #endif
00124 
00125 
00126 /* ========================================================================== */
00127 /* === cholmod_analyze ====================================================== */
00128 /* ========================================================================== */
00129 
00130 /* Orders and analyzes A, AA', PAP', or PAA'P' and returns a symbolic factor
00131  * that can later be passed to cholmod_factorize. */
00132 
00133 cholmod_factor *CHOLMOD(analyze)
00134 (
00135     /* ---- input ---- */
00136     cholmod_sparse *A,  /* matrix to order and analyze */
00137     /* --------------- */
00138     cholmod_common *Common
00139 )
00140 {
00141     return (CHOLMOD(analyze_p) (A, NULL, NULL, 0, Common)) ;
00142 }
00143 
00144 
00145 /* ========================================================================== */
00146 /* === permute_matrices ===================================================== */
00147 /* ========================================================================== */
00148 
00149 /* Permute and transpose a matrix.  Allocates the A1 and A2 matrices, if needed,
00150  * or returns them as NULL if not needed.
00151  */
00152 
00153 static int amesos_permute_matrices
00154 (
00155     /* ---- input ---- */
00156     cholmod_sparse *A,  /* matrix to permute */
00157     Int ordering, /* ordering method used */
00158     Int *Perm,    /* fill-reducing permutation */
00159     Int *fset,    /* subset of 0:(A->ncol)-1 */
00160     size_t fsize, /* size of fset */
00161     Int do_rowcolcounts,/* if TRUE, compute both S and F.  If FALSE, only
00162        * S is needed for the symmetric case, and only F for
00163        * the unsymmetric case */
00164     /* ---- output --- */
00165     cholmod_sparse **A1_handle,     /* see comments below for A1, A2, S, F */
00166     cholmod_sparse **A2_handle,
00167     cholmod_sparse **S_handle,
00168     cholmod_sparse **F_handle,
00169     /* --------------- */
00170     cholmod_common *Common
00171 )
00172 {
00173     cholmod_sparse *A1, *A2, *S, *F ;
00174 
00175     *A1_handle = NULL ;
00176     *A2_handle = NULL ;
00177     *S_handle = NULL ;
00178     *F_handle = NULL ;
00179     A1 = NULL ;
00180     A2 = NULL ;
00181 
00182     if (ordering == CHOLMOD_NATURAL)
00183     {
00184 
00185   /* ------------------------------------------------------------------ */
00186   /* natural ordering of A */
00187   /* ------------------------------------------------------------------ */
00188 
00189   if (A->stype < 0)
00190   {
00191       /* symmetric lower case: A already in lower form, so S=A' */
00192       /* workspace: Iwork (nrow) */
00193       A2 = CHOLMOD(ptranspose) (A, 0, NULL, NULL, 0, Common) ;
00194       F = A ;
00195       S = A2 ;
00196   }
00197   else if (A->stype > 0)
00198   {
00199       /* symmetric upper case: F = pattern of triu (A)', S = A */
00200       /* workspace: Iwork (nrow) */
00201       if (do_rowcolcounts)
00202       {
00203     /* F not needed for symmetric case if do_rowcolcounts FALSE */
00204     A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
00205       }
00206       F = A1 ;
00207       S = A ;
00208   }
00209   else
00210   {
00211       /* unsymmetric case: F = pattern of A (:,f)',  S = A */
00212       /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
00213       A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
00214       F = A1 ;
00215       S = A ;
00216   }
00217 
00218     }
00219     else
00220     {
00221 
00222   /* ------------------------------------------------------------------ */
00223   /* A is permuted */
00224   /* ------------------------------------------------------------------ */
00225 
00226   if (A->stype < 0)
00227   {
00228       /* symmetric lower case: S = tril (A (p,p))' and F = S' */
00229       /* workspace: Iwork (2*nrow) */
00230       A2 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
00231       S = A2 ;
00232       /* workspace: Iwork (nrow) */
00233       if (do_rowcolcounts)
00234       {
00235     /* F not needed for symmetric case if do_rowcolcounts FALSE */
00236     A1 = CHOLMOD(ptranspose) (A2, 0, NULL, NULL, 0, Common) ;
00237       }
00238       F = A1 ;
00239   }
00240   else if (A->stype > 0)
00241   {
00242       /* symmetric upper case: F = triu (A (p,p))' and S = F' */
00243       /* workspace: Iwork (2*nrow) */
00244       A1 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
00245       F = A1 ;
00246       /* workspace: Iwork (nrow) */
00247       A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
00248       S = A2 ;
00249   }
00250   else
00251   {
00252       /* unsymmetric case:     F = A (p,f)'         and S = F' */
00253       /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
00254       A1 = CHOLMOD(ptranspose) (A, 0, Perm, fset, fsize, Common) ;
00255       F = A1 ;
00256       if (do_rowcolcounts)
00257       {
00258     /* S not needed for unsymmetric case if do_rowcolcounts FALSE */
00259     /* workspace: Iwork (nrow) */
00260     A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
00261       }
00262       S = A2 ;
00263   }
00264     }
00265 
00266     /* If any cholmod_*transpose fails, one or more matrices will be NULL */
00267     *A1_handle = A1 ;
00268     *A2_handle = A2 ;
00269     *S_handle = S ;
00270     *F_handle = F ;
00271     return (Common->status == CHOLMOD_OK) ;
00272 }
00273 
00274 
00275 /* ========================================================================== */
00276 /* === cholmod_analyze_ordering ============================================= */
00277 /* ========================================================================== */
00278 
00279 /* Given a matrix A and its fill-reducing permutation, compute the elimination
00280  * tree, its (non-weighted) postordering, and the number of nonzeros in each
00281  * column of L.  Also computes the flop count, the total nonzeros in L, and
00282  * the nonzeros in A (Common->fl, Common->lnz, and Common->anz).
00283  *
00284  * The column counts of L, flop count, and other statistics from
00285  * cholmod_rowcolcounts are not computed if ColCount is NULL.
00286  *
00287  * workspace: Iwork (2*nrow if symmetric, 2*nrow+ncol if unsymmetric),
00288  *  Flag (nrow), Head (nrow+1)
00289  */
00290 
00291 int CHOLMOD(analyze_ordering)
00292 (
00293     /* ---- input ---- */
00294     cholmod_sparse *A,  /* matrix to analyze */
00295     int ordering, /* ordering method used */
00296     Int *Perm,    /* size n, fill-reducing permutation to analyze */
00297     Int *fset,    /* subset of 0:(A->ncol)-1 */
00298     size_t fsize, /* size of fset */
00299     /* ---- output --- */
00300     Int *Parent,  /* size n, elimination tree */
00301     Int *Post,    /* size n, postordering of elimination tree */
00302     Int *ColCount,  /* size n, nnz in each column of L */
00303     /* ---- workspace  */
00304     Int *First,   /* size nworkspace for cholmod_postorder */
00305     Int *Level,   /* size n workspace for cholmod_postorder */
00306     /* --------------- */
00307     cholmod_common *Common
00308 )
00309 {
00310     cholmod_sparse *A1, *A2, *S, *F ;
00311     Int n, ok, do_rowcolcounts ;
00312 
00313     /* check inputs */
00314     RETURN_IF_NULL_COMMON (FALSE) ;
00315     RETURN_IF_NULL (A, FALSE) ;
00316 
00317     n = A->nrow ;
00318 
00319     do_rowcolcounts = (ColCount != NULL) ;
00320 
00321     /* permute A according to Perm and fset */
00322     ok = amesos_permute_matrices (A, ordering, Perm, fset, fsize, do_rowcolcounts,
00323       &A1, &A2, &S, &F, Common) ;
00324 
00325     /* find etree of S (symmetric upper/lower case) or F (unsym case) */
00326     /* workspace: symmmetric: Iwork (nrow), unsym: Iwork (nrow+ncol) */
00327     ok = ok && CHOLMOD(etree) (A->stype ? S:F, Parent, Common) ;
00328 
00329     /* postorder the etree (required by cholmod_rowcolcounts) */
00330     /* workspace: Iwork (2*nrow) */
00331     ok = ok && (CHOLMOD(postorder) (Parent, n, NULL, Post, Common) == n) ;
00332 
00333     /* cholmod_postorder doesn't set Common->status if it returns < n */
00334     Common->status = (!ok && Common->status == CHOLMOD_OK) ?
00335   CHOLMOD_INVALID : Common->status ;
00336 
00337     /* analyze LL'=S or SS' or S(:,f)*S(:,f)' */
00338     /* workspace:
00339      *  if symmetric:   Flag (nrow), Iwork (2*nrow)
00340      *  if unsymmetric: Flag (nrow), Iwork (2*nrow+ncol), Head (nrow+1)
00341      */
00342     if (do_rowcolcounts)
00343     {
00344   ok = ok && CHOLMOD(rowcolcounts) (A->stype ? F:S, fset, fsize, Parent,
00345       Post, NULL, ColCount, First, Level, Common) ;
00346     }
00347 
00348     /* free temporary matrices and return result */
00349     CHOLMOD(free_sparse) (&A1, Common) ;
00350     CHOLMOD(free_sparse) (&A2, Common) ;
00351     return (ok) ;
00352 }
00353 
00354 
00355 /* ========================================================================== */
00356 /* === Free workspace and return L ========================================== */
00357 /* ========================================================================== */
00358 
00359 #define FREE_WORKSPACE_AND_RETURN \
00360 { \
00361     Common->no_workspace_reallocate = FALSE ; \
00362     CHOLMOD(free) (n, sizeof (Int), Lparent,  Common) ; \
00363     CHOLMOD(free) (n, sizeof (Int), Perm,     Common) ; \
00364     CHOLMOD(free) (n, sizeof (Int), ColCount, Common) ; \
00365     if (Common->status < CHOLMOD_OK) \
00366     { \
00367   CHOLMOD(free_factor) (&L, Common) ; \
00368     } \
00369     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ; \
00370     return (L) ; \
00371 }
00372 
00373 
00374 /* ========================================================================== */
00375 /* === cholmod_analyze_p ==================================================== */
00376 /* ========================================================================== */
00377 
00378 /* Orders and analyzes A, AA', PAP', PAA'P', FF', or PFF'P and returns a
00379  * symbolic factor that can later be passed to cholmod_factorize, where
00380  * F = A(:,fset) if fset is not NULL and A->stype is zero.
00381  * UserPerm is tried if non-NULL.  */
00382 
00383 cholmod_factor *CHOLMOD(analyze_p)
00384 (
00385     /* ---- input ---- */
00386     cholmod_sparse *A,  /* matrix to order and analyze */
00387     Int *UserPerm,  /* user-provided permutation, size A->nrow */
00388     Int *fset,    /* subset of 0:(A->ncol)-1 */
00389     size_t fsize, /* size of fset */
00390     /* --------------- */
00391     cholmod_common *Common
00392 )
00393 {
00394     double lnz_best ;
00395     Int *First, *Level, *Work4n, *Cmember, *CParent, *ColCount, *Lperm, *Parent,
00396   *Post, *Perm, *Lparent, *Lcolcount ;
00397     cholmod_factor *L ;
00398     Int k, n, ordering, method, nmethods, status, default_strategy, ncol, uncol,
00399   skip_analysis, skip_best ;
00400     size_t s ;
00401     int ok = TRUE ;
00402 
00403     /* ---------------------------------------------------------------------- */
00404     /* check inputs */
00405     /* ---------------------------------------------------------------------- */
00406 
00407     RETURN_IF_NULL_COMMON (NULL) ;
00408     RETURN_IF_NULL (A, NULL) ;
00409     RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
00410     Common->status = CHOLMOD_OK ;
00411     status = CHOLMOD_OK ;
00412     Common->selected = EMPTY ;
00413     Common->called_nd = FALSE ;
00414 
00415     /* ---------------------------------------------------------------------- */
00416     /* get inputs */
00417     /* ---------------------------------------------------------------------- */
00418 
00419     n = A->nrow ;
00420     ncol = A->ncol ;
00421     uncol = (A->stype == 0) ? (A->ncol) : 0 ;
00422 
00423     /* ---------------------------------------------------------------------- */
00424     /* set the default strategy */
00425     /* ---------------------------------------------------------------------- */
00426 
00427     lnz_best = (double) EMPTY ;
00428     skip_best = FALSE ;
00429     nmethods = MIN (Common->nmethods, CHOLMOD_MAXMETHODS) ;
00430     nmethods = MAX (0, nmethods) ;
00431     PRINT1 (("nmethods "ID"\n", nmethods)) ;
00432 
00433     default_strategy = (nmethods == 0) ;
00434     if (default_strategy)
00435     {
00436   /* default strategy: try UserPerm, if given.  Try AMD for A, or COLAMD
00437    * to order A*A'.  Try METIS for the symmetric case only if AMD reports
00438    * a high degree of fill-in and flop count.  Always try METIS for the
00439    * unsymmetric case.  METIS is not tried if the Partition Module
00440    * isn't installed.   If Common->default_nesdis is TRUE, then NESDIS
00441    * is used as the 3rd ordering instead. */
00442   Common->method [0].ordering = CHOLMOD_GIVEN ;/* skip if UserPerm NULL */
00443   Common->method [1].ordering = CHOLMOD_AMD ;
00444   Common->method [2].ordering = 
00445       (Common->default_nesdis ? CHOLMOD_NESDIS : CHOLMOD_METIS) ;
00446 #ifndef NPARTITION
00447   nmethods = 3 ;
00448 #else
00449   nmethods = 2 ;
00450 #endif
00451     }
00452 
00453 #ifdef NSUPERNODAL
00454     /* CHOLMOD Supernodal module not installed, just do simplicial analysis */
00455     Common->supernodal = CHOLMOD_SIMPLICIAL ;
00456 #endif
00457 
00458     /* ---------------------------------------------------------------------- */
00459     /* allocate workspace */
00460     /* ---------------------------------------------------------------------- */
00461 
00462     /* Note: enough space needs to be allocated here so that routines called by
00463      * cholmod_analyze do not reallocate the space.
00464      */
00465 
00466     /* s = 6*n + uncol */
00467     s = CHOLMOD(mult_size_t) (n, 6, &ok) ;
00468     s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
00469     if (!ok)
00470     {
00471   ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
00472   return (NULL) ;
00473     }
00474 
00475     CHOLMOD(allocate_work) (n, s, 0, Common) ;
00476     if (Common->status < CHOLMOD_OK)
00477     {
00478   return (NULL) ;     /* out of memory */
00479     }
00480     ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00481 
00482     /* ensure that subsequent routines, called by cholmod_analyze, do not
00483      * reallocate any workspace.  This is set back to FALSE in the
00484      * FREE_WORKSPACE_AND_RETURN macro, which is the only way this function
00485      * returns to its caller. */
00486     Common->no_workspace_reallocate = TRUE ;
00487 
00488     /* Use the last 4*n Int's in Iwork for Parent, First, Level, and Post, since
00489      * other CHOLMOD routines will use the first 2n+uncol space.  The ordering
00490      * routines (cholmod_amd, cholmod_colamd, cholmod_ccolamd, cholmod_metis)
00491      * are an exception.  They can use all 6n + ncol space, since the contents
00492      * of Parent, First, Level, and Post are not needed across calls to those
00493      * routines. */
00494     Work4n = Common->Iwork ;
00495     Work4n += 2*((size_t) n) + uncol ;
00496     Parent = Work4n ;
00497     First  = Work4n + n ;
00498     Level  = Work4n + 2*((size_t) n) ;
00499     Post   = Work4n + 3*((size_t) n) ;
00500 
00501     /* note that this assignment means that cholmod_nested_dissection,
00502      * cholmod_ccolamd, and cholmod_camd can use only the first 4n+uncol
00503      * space in Common->Iwork */
00504     Cmember = Post ;
00505     CParent = Level ;
00506 
00507     /* ---------------------------------------------------------------------- */
00508     /* allocate more workspace, and an empty simplicial symbolic factor */
00509     /* ---------------------------------------------------------------------- */
00510 
00511     L = CHOLMOD(allocate_factor) (n, Common) ;
00512     Lparent  = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
00513     Perm     = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
00514     ColCount = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
00515     if (Common->status < CHOLMOD_OK)
00516     {
00517   /* out of memory */
00518   FREE_WORKSPACE_AND_RETURN ;
00519     }
00520     Lperm = L->Perm ;
00521     Lcolcount = L->ColCount ;
00522     Common->anz = EMPTY ;
00523 
00524     /* ---------------------------------------------------------------------- */
00525     /* try all the requested ordering options and backup to AMD if needed */
00526     /* ---------------------------------------------------------------------- */
00527 
00528     /* turn off error handling [ */
00529     Common->try_catch = TRUE ;
00530 
00531     for (method = 0 ; method <= nmethods ; method++)
00532     {
00533 
00534   /* ------------------------------------------------------------------ */
00535   /* determine the method to try */
00536   /* ------------------------------------------------------------------ */
00537 
00538   Common->fl = EMPTY ;
00539   Common->lnz = EMPTY ;
00540   skip_analysis = FALSE ;
00541 
00542   if (method == nmethods)
00543   {
00544       /* All methods failed: backup to AMD */
00545       if (Common->selected == EMPTY && !default_strategy)
00546       {
00547     PRINT1 (("All methods requested failed: backup to AMD\n")) ;
00548     ordering = CHOLMOD_AMD ;
00549       }
00550       else
00551       {
00552     break ;
00553       }
00554   }
00555   else
00556   {
00557       ordering = Common->method [method].ordering ;
00558   }
00559   Common->current = method ;
00560   PRINT1 (("method "ID": Try method: "ID"\n", method, ordering)) ;
00561 
00562   /* ------------------------------------------------------------------ */
00563   /* find the fill-reducing permutation */
00564   /* ------------------------------------------------------------------ */
00565 
00566   if (ordering == CHOLMOD_NATURAL)
00567   {
00568 
00569       /* -------------------------------------------------------------- */
00570       /* natural ordering */
00571       /* -------------------------------------------------------------- */
00572 
00573       for (k = 0 ; k < n ; k++)
00574       {
00575     Perm [k] = k ;
00576       }
00577 
00578   }
00579   else if (ordering == CHOLMOD_GIVEN)
00580   {
00581 
00582       /* -------------------------------------------------------------- */
00583       /* use given ordering of A, if provided */
00584       /* -------------------------------------------------------------- */
00585 
00586       if (UserPerm == NULL)
00587       {
00588     /* this is not an error condition */
00589     PRINT1 (("skip, no user perm given\n")) ;
00590     continue ;
00591       }
00592       for (k = 0 ; k < n ; k++)
00593       {
00594     /* UserPerm is checked in cholmod_ptranspose */
00595     Perm [k] = UserPerm [k] ;
00596       }
00597 
00598   }
00599   else if (ordering == CHOLMOD_AMD)
00600   {
00601 
00602       /* -------------------------------------------------------------- */
00603       /* AMD ordering of A, A*A', or A(:,f)*A(:,f)' */
00604       /* -------------------------------------------------------------- */
00605 
00606       CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
00607       skip_analysis = TRUE ;
00608 
00609   }
00610   else if (ordering == CHOLMOD_COLAMD)
00611   {
00612 
00613       /* -------------------------------------------------------------- */
00614       /* AMD for symmetric case, COLAMD for A*A' or A(:,f)*A(:,f)' */
00615       /* -------------------------------------------------------------- */
00616 
00617       if (A->stype)
00618       {
00619     CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
00620     skip_analysis = TRUE ;
00621       }
00622       else
00623       {
00624     /* Alternative:
00625     CHOLMOD(ccolamd) (A, fset, fsize, NULL, Perm, Common) ;
00626     */
00627     /* do not postorder, it is done later, below */
00628     /* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1)*/
00629     CHOLMOD(colamd) (A, fset, fsize, FALSE, Perm, Common) ;
00630       }
00631 
00632   }
00633   else if (ordering == CHOLMOD_METIS)
00634   {
00635 
00636       /* -------------------------------------------------------------- */
00637       /* use METIS_NodeND directly (via a CHOLMOD wrapper) */
00638       /* -------------------------------------------------------------- */
00639 
00640 #ifndef NPARTITION
00641       /* postorder parameter is false, because it will be later, below */
00642       /* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1) */
00643       Common->called_nd = TRUE ;
00644       CHOLMOD(metis) (A, fset, fsize, FALSE, Perm, Common) ;
00645 #else
00646       Common->status = CHOLMOD_NOT_INSTALLED ;
00647 #endif
00648 
00649   }
00650   else if (ordering == CHOLMOD_NESDIS)
00651   {
00652 
00653       /* -------------------------------------------------------------- */
00654       /* use CHOLMOD's nested dissection */
00655       /* -------------------------------------------------------------- */
00656 
00657       /* this method is based on METIS' node bissection routine
00658        * (METIS_NodeComputeSeparator).  In contrast to METIS_NodeND,
00659        * it calls CAMD or CCOLAMD on the whole graph, instead of MMD
00660        * on just the leaves. */
00661 #ifndef NPARTITION
00662       /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow) */
00663       Common->called_nd = TRUE ;
00664       CHOLMOD(nested_dissection) (A, fset, fsize, Perm, CParent, Cmember,
00665         Common) ;
00666 #else
00667       Common->status = CHOLMOD_NOT_INSTALLED ;
00668 #endif
00669 
00670   }
00671   else
00672   {
00673 
00674       /* -------------------------------------------------------------- */
00675       /* invalid ordering method */
00676       /* -------------------------------------------------------------- */
00677 
00678       Common->status = CHOLMOD_INVALID ;
00679       PRINT1 (("No such ordering: "ID"\n", ordering)) ;
00680   }
00681 
00682   ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
00683 
00684   if (Common->status < CHOLMOD_OK)
00685   {
00686       /* out of memory, or method failed */
00687       status = MIN (status, Common->status) ;
00688       Common->status = CHOLMOD_OK ;
00689       continue ;
00690   }
00691 
00692   /* ------------------------------------------------------------------ */
00693   /* analyze the ordering */
00694   /* ------------------------------------------------------------------ */
00695 
00696   if (!skip_analysis)
00697   {
00698       if (!CHOLMOD(analyze_ordering) (A, ordering, Perm, fset, fsize,
00699         Parent, Post, ColCount, First, Level, Common))
00700       {
00701     /* ordering method failed; clear status and try next method */
00702     status = MIN (status, Common->status) ;
00703     Common->status = CHOLMOD_OK ;
00704     continue ;
00705       }
00706   }
00707 
00708   ASSERT (Common->fl >= 0 && Common->lnz >= 0) ;
00709   Common->method [method].fl  = Common->fl ;
00710   Common->method [method].lnz = Common->lnz ;
00711   PRINT1 (("lnz %g fl %g\n", Common->lnz, Common->fl)) ;
00712 
00713   /* ------------------------------------------------------------------ */
00714   /* pick the best method */
00715   /* ------------------------------------------------------------------ */
00716 
00717   /* fl.pt. compare, but lnz can never be NaN */
00718   if (Common->selected == EMPTY || Common->lnz < lnz_best)
00719   {
00720       Common->selected = method ;
00721       PRINT1 (("this is best so far, method "ID"\n", method)) ;
00722       L->ordering = ordering ;
00723       lnz_best = Common->lnz ;
00724       for (k = 0 ; k < n ; k++)
00725       {
00726     Lperm [k] = Perm [k] ;
00727       }
00728       /* save the results of cholmod_analyze_ordering, if it was called */
00729       skip_best = skip_analysis ;
00730       if (!skip_analysis)
00731       {
00732     /* save the column count; becomes permanent part of L */
00733     for (k = 0 ; k < n ; k++)
00734     {
00735         Lcolcount [k] = ColCount [k] ;
00736     }
00737     /* Parent is needed for weighted postordering and for supernodal
00738      * analysis.  Does not become a permanent part of L */
00739     for (k = 0 ; k < n ; k++)
00740     {
00741         Lparent [k] = Parent [k] ;
00742     }
00743       }
00744   }
00745 
00746   /* ------------------------------------------------------------------ */
00747   /* determine if METIS is to be skipped */
00748   /* ------------------------------------------------------------------ */
00749 
00750   if (default_strategy && ordering == CHOLMOD_AMD)
00751   {
00752       if ((Common->fl < 500 * Common->lnz) ||
00753     (Common->lnz < 5 * Common->anz))
00754       {
00755     /* AMD found an ordering with less than 500 flops per nonzero in
00756      * L, or one with a fill-in ratio (nnz(L)/nnz(A)) of less than
00757      * 5.  This is pretty good, and it's unlikely that METIS will do
00758      * better (this heuristic is based on tests on all symmetric
00759      * positive definite matrices in the UF sparse matrix
00760      * collection, and it works well across a wide range of
00761      * problems).  METIS can take much more time and AMD. */
00762     break ;
00763       }
00764   }
00765     }
00766 
00767     /* turn error printing back on ] */
00768     Common->try_catch = FALSE ;
00769 
00770     /* ---------------------------------------------------------------------- */
00771     /* return if no ordering method succeeded */
00772     /* ---------------------------------------------------------------------- */
00773 
00774     if (Common->selected == EMPTY)
00775     {
00776   /* All methods failed.  
00777    * If two or more methods failed, they may have failed for different
00778    * reasons.  Both would clear Common->status and skip to the next
00779    * method.  Common->status needs to be restored here to the worst error
00780    * obtained in any of the methods.  CHOLMOD_INVALID is worse
00781    * than CHOLMOD_OUT_OF_MEMORY, since the former implies something may
00782    * be wrong with the user's input.  CHOLMOD_OUT_OF_MEMORY is simply an
00783    * indication of lack of resources. */
00784   ASSERT (status < CHOLMOD_OK) ;
00785   ERROR (status, "all methods failed") ;
00786   FREE_WORKSPACE_AND_RETURN ;
00787     }
00788 
00789     /* ---------------------------------------------------------------------- */
00790     /* do the analysis for AMD, if skipped */
00791     /* ---------------------------------------------------------------------- */
00792 
00793     Common->fl  = Common->method [Common->selected].fl  ;
00794     Common->lnz = Common->method [Common->selected].lnz ;
00795     ASSERT (Common->lnz >= 0) ;
00796 
00797     if (skip_best)
00798     {
00799   if (!CHOLMOD(analyze_ordering) (A, L->ordering, Lperm, fset, fsize,
00800     Lparent, Post, Lcolcount, First, Level, Common))
00801   {
00802       /* out of memory, or method failed */
00803       FREE_WORKSPACE_AND_RETURN ;
00804   }
00805     }
00806 
00807     /* ---------------------------------------------------------------------- */
00808     /* postorder the etree, weighted by the column counts */
00809     /* ---------------------------------------------------------------------- */
00810 
00811     if (Common->postorder)
00812     {
00813   /* combine the fill-reducing ordering with the weighted postorder */
00814   /* workspace: Iwork (2*nrow) */
00815   if (CHOLMOD(postorder) (Lparent, n, Lcolcount, Post, Common) == n)
00816   {
00817       /* use First and Level as workspace [ */
00818       Int *Wi = First, *InvPost = Level ;
00819       Int newchild, oldchild, newparent, oldparent ;
00820 
00821       for (k = 0 ; k < n ; k++)
00822       {
00823     Wi [k] = Lperm [Post [k]] ;
00824       }
00825       for (k = 0 ; k < n ; k++)
00826       {
00827     Lperm [k] = Wi [k] ;
00828       }
00829 
00830       for (k = 0 ; k < n ; k++)
00831       {
00832     Wi [k] = Lcolcount [Post [k]] ;
00833       }
00834       for (k = 0 ; k < n ; k++)
00835       {
00836     Lcolcount [k] = Wi [k] ;
00837       }
00838       for (k = 0 ; k < n ; k++)
00839       {
00840     InvPost [Post [k]] = k ;
00841       }
00842 
00843       /* updated Lparent needed only for supernodal case */
00844       for (newchild = 0 ; newchild < n ; newchild++)
00845       {
00846     oldchild = Post [newchild] ;
00847     oldparent = Lparent [oldchild] ;
00848     newparent = (oldparent == EMPTY) ? EMPTY : InvPost [oldparent] ;
00849     Wi [newchild] = newparent ;
00850       }
00851       for (k = 0 ; k < n ; k++)
00852       {
00853     Lparent [k] = Wi [k] ;
00854       }
00855       /* done using Iwork as workspace ] */
00856 
00857       /* L is now postordered, no longer in natural ordering */
00858       if (L->ordering == CHOLMOD_NATURAL)
00859       {
00860     L->ordering = CHOLMOD_POSTORDERED ;
00861       }
00862   }
00863     }
00864 
00865     /* ---------------------------------------------------------------------- */
00866     /* supernodal analysis, if requested or if selected automatically */
00867     /* ---------------------------------------------------------------------- */
00868 
00869 #ifndef NSUPERNODAL
00870     if (Common->supernodal > CHOLMOD_AUTO
00871     || (Common->supernodal == CHOLMOD_AUTO &&
00872   Common->lnz > 0 &&
00873   (Common->fl / Common->lnz) >= Common->supernodal_switch))
00874     {
00875   cholmod_sparse *S, *F, *A2, *A1 ;
00876 
00877   amesos_permute_matrices (A, L->ordering, Lperm, fset, fsize, TRUE,
00878     &A1, &A2, &S, &F, Common) ;
00879 
00880   /* workspace: Flag (nrow), Head (nrow), Iwork (5*nrow) */
00881   CHOLMOD(super_symbolic) (S, F, Lparent, L, Common) ;
00882   PRINT1 (("status %d\n", Common->status)) ;
00883 
00884   CHOLMOD(free_sparse) (&A1, Common) ;
00885   CHOLMOD(free_sparse) (&A2, Common) ;
00886     }
00887 #endif
00888 
00889     /* ---------------------------------------------------------------------- */
00890     /* free temporary matrices and workspace, and return result L */
00891     /* ---------------------------------------------------------------------- */
00892 
00893     FREE_WORKSPACE_AND_RETURN ;
00894 }
00895 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines