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