amesos_cholmod_l_analyze.c

Go to the documentation of this file.

/* ========================================================================== */
/* === 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

/* This file should make the long int version of CHOLMOD */
#define DLONG 1

#include "amesos_cholmod_internal.h"
#include "amesos_cholmod_cholesky.h"

#ifndef NPARTITION
#include "amesos_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 amesos_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 = amesos_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 ;

  amesos_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