/* ========================================================================== */
/* === Cholesky/cholmod_analyze ============================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/Cholesky Module.  Copyright (C) 2005-2006, Timothy A. Davis
 * The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
 * Lesser General Public License.  See lesser.txt for a text of the license.
 * CHOLMOD is also available under other licenses; contact authors for details.
 * http://www.cise.ufl.edu/research/sparse
 * -------------------------------------------------------------------------- */

/* Order and analyze a matrix (either simplicial or supernodal), in prepartion
 * for numerical factorization via cholmod_factorize or via the "expert"
 * routines cholmod_rowfac and cholmod_super_numeric.
 *
 * symmetric case:    A or A(p,p)
 * unsymmetric case:  AA', A(p,:)*A(p,:)', A(:,f)*A(:,f)', or A(p,f)*A(p,f)'
 *
 * For the symmetric case, only the upper or lower triangular part of A is
 * accessed (depending on the type of A).  LL'=A (or permuted A) is analzed.
 * For the unsymmetric case (LL'=AA' or permuted A).
 *
 * There can be no duplicate entries in p or f.  p is of length m if A is
 * m-by-n.  f can be length 0 to n.
 *
 * In both cases, the columns of A need not be sorted.  A can be in packed
 * or unpacked form.
 *
 * Ordering options include:
 *
 *	natural:    A is not permuted to reduce fill-in
 *	given:	    a permutation can be provided to this routine (UserPerm)
 *	AMD:	    approximate minumum degree (AMD for the symmetric case,
 *		    COLAMD for the AA' case).
 *	METIS:	    nested dissection with METIS_NodeND
 *	NESDIS:	    nested dissection using METIS_NodeComputeSeparator,
 *		    typically followed by a constrained minimum degree
 *		    (CAMD for the symmetric case, CCOLAMD for the AA' case).
 *
 * Multiple ordering options can be tried (up to 9 of them), and the best one
 * is selected (the one that gives the smallest number of nonzeros in the
 * simplicial factor L).  If one method fails, cholmod_analyze keeps going, and
 * picks the best among the methods that succeeded.  This routine fails (and
 * returns NULL) if either initial memory allocation fails, all ordering methods
 * fail, or the supernodal analysis (if requested) fails.  By default, the 9
 * methods available are:
 *
 *	1) given permutation (skipped if UserPerm is NULL)
 *	2) AMD (symmetric case) or COLAMD (unsymmetric case)
 *	3) METIS with default parameters
 *	4) NESDIS with default parameters (stopping the partitioning when
 *	    the graph is of size nd_small = 200 or less, remove nodes with
 *	    more than max (16, prune_dense * sqrt (n)) nodes where
 *	    prune_dense = 10, and follow partitioning with CCOLAMD, a
 *	    constrained minimum degree ordering).
 *	5) natural
 *	6) NESDIS, nd_small = 20000, prune_dense = 10
 *	7) NESDIS, nd_small =     4, prune_dense = 10, no min degree
 *	8) NESDIS, nd_small =   200, prune_dense = 0
 *	9) COLAMD for A*A' or AMD for A
 *
 * By default, the first two are tried, and METIS is tried if AMD reports a high
 * flop count and fill-in.  Let fl denote the flop count for the AMD, ordering,
 * nnz(L) the # of nonzeros in L, and nnz(tril(A)) (or A*A').  If
 * fl/nnz(L) >= 500 and nnz(L)/nnz(tril(A)) >= 5, then METIS is attempted.  The
 * best ordering is used (UserPerm if given, AMD, and METIS if attempted).  If
 * you do not have METIS, only the first two will be tried (user permutation,
 * if provided, and AMD/COLAMD).  This default behavior is obtained when
 * Common->nmethods is zero.  In this case, methods 0, 1, and 2 in
 * Common->method [..] are reset to User-provided, AMD, and METIS (or NESDIS
 * if Common->default_nesdis is set to the non-default value of TRUE),
 * respectively.
 *
 * You can modify these 9 methods and the number of methods tried by changing
 * parameters in the Common argument.  If you know the best ordering for your
 * matrix, set Common->nmethods to 1 and set Common->method[0].ordering to the
 * requested ordering method.  Parameters for each method can also be modified
 * (refer to cholmod.h for details).
 *
 * Note that it is possible for METIS to terminate your program if it runs out
 * of memory.  This is not the case for any CHOLMOD or minimum degree ordering
 * routine (AMD, COLAMD, CAMD, CCOLAMD, or CSYMAMD).  Since NESDIS relies on
 * METIS, it too can terminate your program.
 *
 * The factor L is returned as simplicial symbolic (L->is_super FALSE) if
 * Common->supernodal <= CHOLMOD_SIMPLICIAL (0) or as supernodal symbolic if
 * Common->supernodal >= CHOLMOD_SUPERNODAL (2).  If Common->supernodal is
 * equal to CHOLMOD_AUTO (1), then L is simplicial if the flop count per
 * nonzero in L is less than Common->supernodal_switch (default: 40), and
 * is returned as a supernodal factor otherwise.
 *
 * In both cases, L->xtype is CHOLMOD_PATTERN.
 * A subsequent call to cholmod_factorize will perform a
 * simplicial or supernodal factorization, depending on the type of L.
 *
 * For the simplicial case, L contains the fill-reducing permutation (L->Perm)
 * and the counts of nonzeros in each column of L (L->ColCount).  For the
 * supernodal case, L also contains the nonzero pattern of each supernode.
 *
 * workspace: Flag (nrow), Head (nrow+1)
 *	if symmetric:   Iwork (6*nrow)
 *	if unsymmetric: Iwork (6*nrow+ncol).
 *	calls various ordering routines, which typically allocate O(nnz(A))
 *	temporary workspace ((2 to 3)*nnz(A) * sizeof (Int) is typical, but it
 *	can be much higher if A*A' must be explicitly formed for METIS).  Also
 *	allocates up to 2 temporary (permuted/transpose) copies of the nonzero
 *	pattern of A, and up to 3*n*sizeof(Int) additional workspace.
 *
 * Supports any xtype (pattern, real, complex, or zomplex)
 */

#ifndef NCHOLESKY

#include "cholmod_internal.h"
#include "cholmod_cholesky.h"

#ifndef NSUPERNODAL
#include "cholmod_supernodal.h"
#endif

#ifndef NPARTITION
#include "cholmod_partition.h"
#endif


/* ========================================================================== */
/* === cholmod_analyze ====================================================== */
/* ========================================================================== */

/* Orders and analyzes A, AA', PAP', or PAA'P' and returns a symbolic factor
 * that can later be passed to cholmod_factorize. */

cholmod_factor *CHOLMOD(analyze)
(
    /* ---- input ---- */
    cholmod_sparse *A,	/* matrix to order and analyze */
    /* --------------- */
    cholmod_common *Common
)
{
    return (CHOLMOD(analyze_p) (A, NULL, NULL, 0, Common)) ;
}


/* ========================================================================== */
/* === permute_matrices ===================================================== */
/* ========================================================================== */

/* Permute and transpose a matrix.  Allocates the A1 and A2 matrices, if needed,
 * or returns them as NULL if not needed.
 */

static int permute_matrices
(
    /* ---- input ---- */
    cholmod_sparse *A,	/* matrix to permute */
    Int ordering,	/* ordering method used */
    Int *Perm,		/* fill-reducing permutation */
    Int *fset,		/* subset of 0:(A->ncol)-1 */
    size_t fsize,	/* size of fset */
    Int do_rowcolcounts,/* if TRUE, compute both S and F.  If FALSE, only
			 * S is needed for the symmetric case, and only F for
			 * the unsymmetric case */
    /* ---- output --- */
    cholmod_sparse **A1_handle,	    /* see comments below for A1, A2, S, F */
    cholmod_sparse **A2_handle,
    cholmod_sparse **S_handle,
    cholmod_sparse **F_handle,
    /* --------------- */
    cholmod_common *Common
)
{
    cholmod_sparse *A1, *A2, *S, *F ;

    *A1_handle = NULL ;
    *A2_handle = NULL ;
    *S_handle = NULL ;
    *F_handle = NULL ;
    A1 = NULL ;
    A2 = NULL ;

    if (ordering == CHOLMOD_NATURAL)
    {

	/* ------------------------------------------------------------------ */
	/* natural ordering of A */
	/* ------------------------------------------------------------------ */

	if (A->stype < 0)
	{
	    /* symmetric lower case: A already in lower form, so S=A' */
	    /* workspace: Iwork (nrow) */
	    A2 = CHOLMOD(ptranspose) (A, 0, NULL, NULL, 0, Common) ;
	    F = A ;
	    S = A2 ;
	}
	else if (A->stype > 0)
	{
	    /* symmetric upper case: F = pattern of triu (A)', S = A */
	    /* workspace: Iwork (nrow) */
	    if (do_rowcolcounts)
	    {
		/* F not needed for symmetric case if do_rowcolcounts FALSE */
		A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
	    }
	    F = A1 ;
	    S = A ;
	}
	else
	{
	    /* unsymmetric case: F = pattern of A (:,f)',  S = A */
	    /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
	    A1 = CHOLMOD(ptranspose) (A, 0, NULL, fset, fsize, Common) ;
	    F = A1 ;
	    S = A ;
	}

    }
    else
    {

	/* ------------------------------------------------------------------ */
	/* A is permuted */
	/* ------------------------------------------------------------------ */

	if (A->stype < 0)
	{
	    /* symmetric lower case: S = tril (A (p,p))' and F = S' */
	    /* workspace: Iwork (2*nrow) */
	    A2 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
	    S = A2 ;
	    /* workspace: Iwork (nrow) */
	    if (do_rowcolcounts)
	    {
		/* F not needed for symmetric case if do_rowcolcounts FALSE */
		A1 = CHOLMOD(ptranspose) (A2, 0, NULL, NULL, 0, Common) ;
	    }
	    F = A1 ;
	}
	else if (A->stype > 0)
	{
	    /* symmetric upper case: F = triu (A (p,p))' and S = F' */
	    /* workspace: Iwork (2*nrow) */
	    A1 = CHOLMOD(ptranspose) (A, 0, Perm, NULL, 0, Common) ;
	    F = A1 ;
	    /* workspace: Iwork (nrow) */
	    A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
	    S = A2 ;
	}
	else
	{
	    /* unsymmetric case:     F = A (p,f)'         and S = F' */
	    /* workspace: Iwork (nrow if no fset, MAX(nrow,ncol) if fset) */
	    A1 = CHOLMOD(ptranspose) (A, 0, Perm, fset, fsize, Common) ;
	    F = A1 ;
	    if (do_rowcolcounts)
	    {
		/* S not needed for unsymmetric case if do_rowcolcounts FALSE */
		/* workspace: Iwork (nrow) */
		A2 = CHOLMOD(ptranspose) (A1, 0, NULL, NULL, 0, Common) ;
	    }
	    S = A2 ;
	}
    }

    /* If any cholmod_*transpose fails, one or more matrices will be NULL */
    *A1_handle = A1 ;
    *A2_handle = A2 ;
    *S_handle = S ;
    *F_handle = F ;
    return (Common->status == CHOLMOD_OK) ;
}


/* ========================================================================== */
/* === cholmod_analyze_ordering ============================================= */
/* ========================================================================== */

/* Given a matrix A and its fill-reducing permutation, compute the elimination
 * tree, its (non-weighted) postordering, and the number of nonzeros in each
 * column of L.  Also computes the flop count, the total nonzeros in L, and
 * the nonzeros in A (Common->fl, Common->lnz, and Common->anz).
 *
 * The column counts of L, flop count, and other statistics from
 * cholmod_rowcolcounts are not computed if ColCount is NULL.
 *
 * workspace: Iwork (2*nrow if symmetric, 2*nrow+ncol if unsymmetric),
 *	Flag (nrow), Head (nrow+1)
 */

int CHOLMOD(analyze_ordering)
(
    /* ---- input ---- */
    cholmod_sparse *A,	/* matrix to analyze */
    int ordering,	/* ordering method used */
    Int *Perm,		/* size n, fill-reducing permutation to analyze */
    Int *fset,		/* subset of 0:(A->ncol)-1 */
    size_t fsize,	/* size of fset */
    /* ---- output --- */
    Int *Parent,	/* size n, elimination tree */
    Int *Post,		/* size n, postordering of elimination tree */
    Int *ColCount,	/* size n, nnz in each column of L */
    /* ---- workspace  */
    Int *First,		/* size nworkspace for cholmod_postorder */
    Int *Level,		/* size n workspace for cholmod_postorder */
    /* --------------- */
    cholmod_common *Common
)
{
    cholmod_sparse *A1, *A2, *S, *F ;
    Int n, ok, do_rowcolcounts ;

    /* check inputs */
    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;

    n = A->nrow ;

    do_rowcolcounts = (ColCount != NULL) ;

    /* permute A according to Perm and fset */
    ok = permute_matrices (A, ordering, Perm, fset, fsize, do_rowcolcounts,
	    &A1, &A2, &S, &F, Common) ;

    /* find etree of S (symmetric upper/lower case) or F (unsym case) */
    /* workspace: symmmetric: Iwork (nrow), unsym: Iwork (nrow+ncol) */
    ok = ok && CHOLMOD(etree) (A->stype ? S:F, Parent, Common) ;

    /* postorder the etree (required by cholmod_rowcolcounts) */
    /* workspace: Iwork (2*nrow) */
    ok = ok && (CHOLMOD(postorder) (Parent, n, NULL, Post, Common) == n) ;

    /* cholmod_postorder doesn't set Common->status if it returns < n */
    Common->status = (!ok && Common->status == CHOLMOD_OK) ?
	CHOLMOD_INVALID : Common->status ;

    /* analyze LL'=S or SS' or S(:,f)*S(:,f)' */
    /* workspace:
     *	if symmetric:   Flag (nrow), Iwork (2*nrow)
     *	if unsymmetric: Flag (nrow), Iwork (2*nrow+ncol), Head (nrow+1)
     */
    if (do_rowcolcounts)
    {
	ok = ok && CHOLMOD(rowcolcounts) (A->stype ? F:S, fset, fsize, Parent,
	    Post, NULL, ColCount, First, Level, Common) ;
    }

    /* free temporary matrices and return result */
    CHOLMOD(free_sparse) (&A1, Common) ;
    CHOLMOD(free_sparse) (&A2, Common) ;
    return (ok) ;
}


/* ========================================================================== */
/* === Free workspace and return L ========================================== */
/* ========================================================================== */

#define FREE_WORKSPACE_AND_RETURN \
{ \
    Common->no_workspace_reallocate = FALSE ; \
    CHOLMOD(free) (n, sizeof (Int), Lparent,  Common) ; \
    CHOLMOD(free) (n, sizeof (Int), Perm,     Common) ; \
    CHOLMOD(free) (n, sizeof (Int), ColCount, Common) ; \
    if (Common->status < CHOLMOD_OK) \
    { \
	CHOLMOD(free_factor) (&L, Common) ; \
    } \
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ; \
    return (L) ; \
}


/* ========================================================================== */
/* === cholmod_analyze_p ==================================================== */
/* ========================================================================== */

/* Orders and analyzes A, AA', PAP', PAA'P', FF', or PFF'P and returns a
 * symbolic factor that can later be passed to cholmod_factorize, where
 * F = A(:,fset) if fset is not NULL and A->stype is zero.
 * UserPerm is tried if non-NULL.  */

cholmod_factor *CHOLMOD(analyze_p)
(
    /* ---- input ---- */
    cholmod_sparse *A,	/* matrix to order and analyze */
    Int *UserPerm,	/* user-provided permutation, size A->nrow */
    Int *fset,		/* subset of 0:(A->ncol)-1 */
    size_t fsize,	/* size of fset */
    /* --------------- */
    cholmod_common *Common
)
{
    double lnz_best ;
    Int *First, *Level, *Work4n, *Cmember, *CParent, *ColCount, *Lperm, *Parent,
	*Post, *Perm, *Lparent, *Lcolcount ;
    cholmod_factor *L ;
    Int k, n, ordering, method, nmethods, status, default_strategy, ncol, uncol,
	skip_analysis, skip_best ;
    size_t s ;
    int ok = TRUE ;

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    RETURN_IF_NULL_COMMON (NULL) ;
    RETURN_IF_NULL (A, NULL) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, NULL) ;
    Common->status = CHOLMOD_OK ;
    status = CHOLMOD_OK ;
    Common->selected = EMPTY ;
    Common->called_nd = FALSE ;

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    n = A->nrow ;
    ncol = A->ncol ;
    uncol = (A->stype == 0) ? (A->ncol) : 0 ;

    /* ---------------------------------------------------------------------- */
    /* set the default strategy */
    /* ---------------------------------------------------------------------- */

    lnz_best = (double) EMPTY ;
    skip_best = FALSE ;
    nmethods = MIN (Common->nmethods, CHOLMOD_MAXMETHODS) ;
    nmethods = MAX (0, nmethods) ;
    PRINT1 (("nmethods "ID"\n", nmethods)) ;

    default_strategy = (nmethods == 0) ;
    if (default_strategy)
    {
	/* default strategy: try UserPerm, if given.  Try AMD for A, or COLAMD
	 * to order A*A'.  Try METIS for the symmetric case only if AMD reports
	 * a high degree of fill-in and flop count.  Always try METIS for the
	 * unsymmetric case.  METIS is not tried if the Partition Module
	 * isn't installed.   If Common->default_nesdis is TRUE, then NESDIS
	 * is used as the 3rd ordering instead. */
	Common->method [0].ordering = CHOLMOD_GIVEN ;/* skip if UserPerm NULL */
	Common->method [1].ordering = CHOLMOD_AMD ;
	Common->method [2].ordering = 
	    (Common->default_nesdis ? CHOLMOD_NESDIS : CHOLMOD_METIS) ;
#ifndef NPARTITION
	nmethods = 3 ;
#else
	nmethods = 2 ;
#endif
    }

#ifdef NSUPERNODAL
    /* CHOLMOD Supernodal module not installed, just do simplicial analysis */
    Common->supernodal = CHOLMOD_SIMPLICIAL ;
#endif

    /* ---------------------------------------------------------------------- */
    /* allocate workspace */
    /* ---------------------------------------------------------------------- */

    /* Note: enough space needs to be allocated here so that routines called by
     * cholmod_analyze do not reallocate the space.
     */

    /* s = 6*n + uncol */
    s = CHOLMOD(mult_size_t) (n, 6, &ok) ;
    s = CHOLMOD(add_size_t) (s, uncol, &ok) ;
    if (!ok)
    {
	ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
	return (NULL) ;
    }

    CHOLMOD(allocate_work) (n, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
	return (NULL) ;	    /* out of memory */
    }
    ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

    /* ensure that subsequent routines, called by cholmod_analyze, do not
     * reallocate any workspace.  This is set back to FALSE in the
     * FREE_WORKSPACE_AND_RETURN macro, which is the only way this function
     * returns to its caller. */
    Common->no_workspace_reallocate = TRUE ;

    /* Use the last 4*n Int's in Iwork for Parent, First, Level, and Post, since
     * other CHOLMOD routines will use the first 2n+uncol space.  The ordering
     * routines (cholmod_amd, cholmod_colamd, cholmod_ccolamd, cholmod_metis)
     * are an exception.  They can use all 6n + ncol space, since the contents
     * of Parent, First, Level, and Post are not needed across calls to those
     * routines. */
    Work4n = Common->Iwork ;
    Work4n += 2*((size_t) n) + uncol ;
    Parent = Work4n ;
    First  = Work4n + n ;
    Level  = Work4n + 2*((size_t) n) ;
    Post   = Work4n + 3*((size_t) n) ;

    /* note that this assignment means that cholmod_nested_dissection,
     * cholmod_ccolamd, and cholmod_camd can use only the first 4n+uncol
     * space in Common->Iwork */
    Cmember = Post ;
    CParent = Level ;

    /* ---------------------------------------------------------------------- */
    /* allocate more workspace, and an empty simplicial symbolic factor */
    /* ---------------------------------------------------------------------- */

    L = CHOLMOD(allocate_factor) (n, Common) ;
    Lparent  = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
    Perm     = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
    ColCount = CHOLMOD(malloc) (n, sizeof (Int), Common) ;
    if (Common->status < CHOLMOD_OK)
    {
	/* out of memory */
	FREE_WORKSPACE_AND_RETURN ;
    }
    Lperm = L->Perm ;
    Lcolcount = L->ColCount ;
    Common->anz = EMPTY ;

    /* ---------------------------------------------------------------------- */
    /* try all the requested ordering options and backup to AMD if needed */
    /* ---------------------------------------------------------------------- */

    /* turn off error handling [ */
    Common->try_catch = TRUE ;

    for (method = 0 ; method <= nmethods ; method++)
    {

	/* ------------------------------------------------------------------ */
	/* determine the method to try */
	/* ------------------------------------------------------------------ */

	Common->fl = EMPTY ;
	Common->lnz = EMPTY ;
	skip_analysis = FALSE ;

	if (method == nmethods)
	{
	    /* All methods failed: backup to AMD */
	    if (Common->selected == EMPTY && !default_strategy)
	    {
		PRINT1 (("All methods requested failed: backup to AMD\n")) ;
		ordering = CHOLMOD_AMD ;
	    }
	    else
	    {
		break ;
	    }
	}
	else
	{
	    ordering = Common->method [method].ordering ;
	}
	Common->current = method ;
	PRINT1 (("method "ID": Try method: "ID"\n", method, ordering)) ;

	/* ------------------------------------------------------------------ */
	/* find the fill-reducing permutation */
	/* ------------------------------------------------------------------ */

	if (ordering == CHOLMOD_NATURAL)
	{

	    /* -------------------------------------------------------------- */
	    /* natural ordering */
	    /* -------------------------------------------------------------- */

	    for (k = 0 ; k < n ; k++)
	    {
		Perm [k] = k ;
	    }

	}
	else if (ordering == CHOLMOD_GIVEN)
	{

	    /* -------------------------------------------------------------- */
	    /* use given ordering of A, if provided */
	    /* -------------------------------------------------------------- */

	    if (UserPerm == NULL)
	    {
		/* this is not an error condition */
		PRINT1 (("skip, no user perm given\n")) ;
		continue ;
	    }
	    for (k = 0 ; k < n ; k++)
	    {
		/* UserPerm is checked in cholmod_ptranspose */
		Perm [k] = UserPerm [k] ;
	    }

	}
	else if (ordering == CHOLMOD_AMD)
	{

	    /* -------------------------------------------------------------- */
	    /* AMD ordering of A, A*A', or A(:,f)*A(:,f)' */
	    /* -------------------------------------------------------------- */

	    CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
	    skip_analysis = TRUE ;

	}
	else if (ordering == CHOLMOD_COLAMD)
	{

	    /* -------------------------------------------------------------- */
	    /* AMD for symmetric case, COLAMD for A*A' or A(:,f)*A(:,f)' */
	    /* -------------------------------------------------------------- */

	    if (A->stype)
	    {
		CHOLMOD(amd) (A, fset, fsize, Perm, Common) ;
		skip_analysis = TRUE ;
	    }
	    else
	    {
		/* Alternative:
		CHOLMOD(ccolamd) (A, fset, fsize, NULL, Perm, Common) ;
		*/
		/* do not postorder, it is done later, below */
		/* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1)*/
		CHOLMOD(colamd) (A, fset, fsize, FALSE, Perm, Common) ;
	    }

	}
	else if (ordering == CHOLMOD_METIS)
	{

	    /* -------------------------------------------------------------- */
	    /* use METIS_NodeND directly (via a CHOLMOD wrapper) */
	    /* -------------------------------------------------------------- */

#ifndef NPARTITION
	    /* postorder parameter is false, because it will be later, below */
	    /* workspace: Iwork (4*nrow+uncol), Flag (nrow), Head (nrow+1) */
	    Common->called_nd = TRUE ;
	    CHOLMOD(metis) (A, fset, fsize, FALSE, Perm, Common) ;
#else
	    Common->status = CHOLMOD_NOT_INSTALLED ;
#endif

	}
	else if (ordering == CHOLMOD_NESDIS)
	{

	    /* -------------------------------------------------------------- */
	    /* use CHOLMOD's nested dissection */
	    /* -------------------------------------------------------------- */

	    /* this method is based on METIS' node bissection routine
	     * (METIS_NodeComputeSeparator).  In contrast to METIS_NodeND,
	     * it calls CAMD or CCOLAMD on the whole graph, instead of MMD
	     * on just the leaves. */
#ifndef NPARTITION
	    /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow) */
	    Common->called_nd = TRUE ;
	    CHOLMOD(nested_dissection) (A, fset, fsize, Perm, CParent, Cmember,
		    Common) ;
#else
	    Common->status = CHOLMOD_NOT_INSTALLED ;
#endif

	}
	else
	{

	    /* -------------------------------------------------------------- */
	    /* invalid ordering method */
	    /* -------------------------------------------------------------- */

	    Common->status = CHOLMOD_INVALID ;
	    PRINT1 (("No such ordering: "ID"\n", ordering)) ;
	}

	ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;

	if (Common->status < CHOLMOD_OK)
	{
	    /* out of memory, or method failed */
	    status = MIN (status, Common->status) ;
	    Common->status = CHOLMOD_OK ;
	    continue ;
	}

	/* ------------------------------------------------------------------ */
	/* analyze the ordering */
	/* ------------------------------------------------------------------ */

	if (!skip_analysis)
	{
	    if (!CHOLMOD(analyze_ordering) (A, ordering, Perm, fset, fsize,
		    Parent, Post, ColCount, First, Level, Common))
	    {
		/* ordering method failed; clear status and try next method */
		status = MIN (status, Common->status) ;
		Common->status = CHOLMOD_OK ;
		continue ;
	    }
	}

	ASSERT (Common->fl >= 0 && Common->lnz >= 0) ;
	Common->method [method].fl  = Common->fl ;
	Common->method [method].lnz = Common->lnz ;
	PRINT1 (("lnz %g fl %g\n", Common->lnz, Common->fl)) ;

	/* ------------------------------------------------------------------ */
	/* pick the best method */
	/* ------------------------------------------------------------------ */

	/* fl.pt. compare, but lnz can never be NaN */
	if (Common->selected == EMPTY || Common->lnz < lnz_best)
	{
	    Common->selected = method ;
	    PRINT1 (("this is best so far, method "ID"\n", method)) ;
	    L->ordering = ordering ;
	    lnz_best = Common->lnz ;
	    for (k = 0 ; k < n ; k++)
	    {
		Lperm [k] = Perm [k] ;
	    }
	    /* save the results of cholmod_analyze_ordering, if it was called */
	    skip_best = skip_analysis ;
	    if (!skip_analysis)
	    {
		/* save the column count; becomes permanent part of L */
		for (k = 0 ; k < n ; k++)
		{
		    Lcolcount [k] = ColCount [k] ;
		}
		/* Parent is needed for weighted postordering and for supernodal
		 * analysis.  Does not become a permanent part of L */
		for (k = 0 ; k < n ; k++)
		{
		    Lparent [k] = Parent [k] ;
		}
	    }
	}

	/* ------------------------------------------------------------------ */
	/* determine if METIS is to be skipped */
	/* ------------------------------------------------------------------ */

	if (default_strategy && ordering == CHOLMOD_AMD)
	{
	    if ((Common->fl < 500 * Common->lnz) ||
		(Common->lnz < 5 * Common->anz))
	    {
		/* AMD found an ordering with less than 500 flops per nonzero in
		 * L, or one with a fill-in ratio (nnz(L)/nnz(A)) of less than
		 * 5.  This is pretty good, and it's unlikely that METIS will do
		 * better (this heuristic is based on tests on all symmetric
		 * positive definite matrices in the UF sparse matrix
		 * collection, and it works well across a wide range of
		 * problems).  METIS can take much more time and AMD. */
		break ;
	    }
	}
    }

    /* turn error printing back on ] */
    Common->try_catch = FALSE ;

    /* ---------------------------------------------------------------------- */
    /* return if no ordering method succeeded */
    /* ---------------------------------------------------------------------- */

    if (Common->selected == EMPTY)
    {
	/* All methods failed.  
	 * If two or more methods failed, they may have failed for different
	 * reasons.  Both would clear Common->status and skip to the next
	 * method.  Common->status needs to be restored here to the worst error
	 * obtained in any of the methods.  CHOLMOD_INVALID is worse
	 * than CHOLMOD_OUT_OF_MEMORY, since the former implies something may
	 * be wrong with the user's input.  CHOLMOD_OUT_OF_MEMORY is simply an
	 * indication of lack of resources. */
	ASSERT (status < CHOLMOD_OK) ;
	ERROR (status, "all methods failed") ;
	FREE_WORKSPACE_AND_RETURN ;
    }

    /* ---------------------------------------------------------------------- */
    /* do the analysis for AMD, if skipped */
    /* ---------------------------------------------------------------------- */

    Common->fl  = Common->method [Common->selected].fl  ;
    Common->lnz = Common->method [Common->selected].lnz ;
    ASSERT (Common->lnz >= 0) ;

    if (skip_best)
    {
	if (!CHOLMOD(analyze_ordering) (A, L->ordering, Lperm, fset, fsize,
		Lparent, Post, Lcolcount, First, Level, Common))
	{
	    /* out of memory, or method failed */
	    FREE_WORKSPACE_AND_RETURN ;
	}
    }

    /* ---------------------------------------------------------------------- */
    /* postorder the etree, weighted by the column counts */
    /* ---------------------------------------------------------------------- */

    if (Common->postorder)
    {
	/* combine the fill-reducing ordering with the weighted postorder */
	/* workspace: Iwork (2*nrow) */
	if (CHOLMOD(postorder) (Lparent, n, Lcolcount, Post, Common) == n)
	{
	    /* use First and Level as workspace [ */
	    Int *Wi = First, *InvPost = Level ;
	    Int newchild, oldchild, newparent, oldparent ;

	    for (k = 0 ; k < n ; k++)
	    {
		Wi [k] = Lperm [Post [k]] ;
	    }
	    for (k = 0 ; k < n ; k++)
	    {
		Lperm [k] = Wi [k] ;
	    }

	    for (k = 0 ; k < n ; k++)
	    {
		Wi [k] = Lcolcount [Post [k]] ;
	    }
	    for (k = 0 ; k < n ; k++)
	    {
		Lcolcount [k] = Wi [k] ;
	    }
	    for (k = 0 ; k < n ; k++)
	    {
		InvPost [Post [k]] = k ;
	    }

	    /* updated Lparent needed only for supernodal case */
	    for (newchild = 0 ; newchild < n ; newchild++)
	    {
		oldchild = Post [newchild] ;
		oldparent = Lparent [oldchild] ;
		newparent = (oldparent == EMPTY) ? EMPTY : InvPost [oldparent] ;
		Wi [newchild] = newparent ;
	    }
	    for (k = 0 ; k < n ; k++)
	    {
		Lparent [k] = Wi [k] ;
	    }
	    /* done using Iwork as workspace ] */

	    /* L is now postordered, no longer in natural ordering */
	    if (L->ordering == CHOLMOD_NATURAL)
	    {
		L->ordering = CHOLMOD_POSTORDERED ;
	    }
	}
    }

    /* ---------------------------------------------------------------------- */
    /* supernodal analysis, if requested or if selected automatically */
    /* ---------------------------------------------------------------------- */

#ifndef NSUPERNODAL
    if (Common->supernodal > CHOLMOD_AUTO
    || (Common->supernodal == CHOLMOD_AUTO &&
	Common->lnz > 0 &&
	(Common->fl / Common->lnz) >= Common->supernodal_switch))
    {
	cholmod_sparse *S, *F, *A2, *A1 ;

	permute_matrices (A, L->ordering, Lperm, fset, fsize, TRUE,
		&A1, &A2, &S, &F, Common) ;

	/* workspace: Flag (nrow), Head (nrow), Iwork (5*nrow) */
	CHOLMOD(super_symbolic) (S, F, Lparent, L, Common) ;
	PRINT1 (("status %d\n", Common->status)) ;

	CHOLMOD(free_sparse) (&A1, Common) ;
	CHOLMOD(free_sparse) (&A2, Common) ;
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* free temporary matrices and workspace, and return result L */
    /* ---------------------------------------------------------------------- */

    FREE_WORKSPACE_AND_RETURN ;
}
#endif