amesos_cholmod_etree.c

Go to the documentation of this file.

/* ========================================================================== */
/* === Cholesky/cholmod_etree =============================================== */
/* ========================================================================== */

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

/* Compute the elimination tree of A or A'*A
 *
 * In the symmetric case, the upper triangular part of A is used.  Entries not
 * in this part of the matrix are ignored.  Computing the etree of a symmetric
 * matrix from just its lower triangular entries is not supported.
 *
 * In the unsymmetric case, all of A is used, and the etree of A'*A is computed.
 *
 * References:
 *
 * J. Liu, "A compact row storage scheme for Cholesky factors", ACM Trans.
 * Math. Software, vol 12, 1986, pp. 127-148.
 *
 * J. Liu, "The role of elimination trees in sparse factorization", SIAM J.
 * Matrix Analysis & Applic., vol 11, 1990, pp. 134-172.
 *
 * J. Gilbert, X. Li, E. Ng, B. Peyton, "Computing row and column counts for
 * sparse QR and LU factorization", BIT, vol 41, 2001, pp. 693-710.
 *
 * workspace: symmetric: Iwork (nrow), unsymmetric: Iwork (nrow+ncol)
 *
 * Supports any xtype (pattern, real, complex, or zomplex)
 */

#ifndef NCHOLESKY

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

/* ========================================================================== */
/* === update_etree ========================================================= */
/* ========================================================================== */

static void amesos_update_etree
(
    /* inputs, not modified */
    Int k,    /* process the edge (k,i) in the input graph */
    Int i,
    /* inputs, modified on output */
    Int Parent [ ], /* Parent [t] = p if p is the parent of t */
    Int Ancestor [ ]  /* Ancestor [t] is the ancestor of node t in the
         partially-constructed etree */
)
{
    Int a ;
    for ( ; ; )   /* traverse the path from k to the root of the tree */
    {
  a = Ancestor [k] ;
  if (a == i)
  {
      /* final ancestor reached; no change to tree */
      return ;
  }
  /* perform path compression */
  Ancestor [k] = i ;
  if (a == EMPTY)
  {
      /* final ancestor undefined; this is a new edge in the tree */
      Parent [k] = i ;
      return ;
  }
  /* traverse up to the ancestor of k */
  k = a ;
    }
}

/* ========================================================================== */
/* === cholmod_etree ======================================================== */
/* ========================================================================== */

/* Find the elimination tree of A or A'*A */

int CHOLMOD(etree)
(
    /* ---- input ---- */
    cholmod_sparse *A,
    /* ---- output --- */
    Int *Parent,  /* size ncol.  Parent [j] = p if p is the parent of j */
    /* --------------- */
    cholmod_common *Common
)
{
    Int *Ap, *Ai, *Anz, *Ancestor, *Prev, *Iwork ;
    Int i, j, jprev, p, pend, nrow, ncol, packed, stype ;
    size_t s ;
    int ok = TRUE ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_NULL (Parent, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    Common->status = CHOLMOD_OK ;

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

    stype = A->stype ;

    /* s = A->nrow + (stype ? 0 : A->ncol) */
    s = CHOLMOD(add_size_t) (A->nrow, (stype ? 0 : A->ncol), &ok) ;
    if (!ok)
    {
  ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
  return (FALSE) ;
    }

    CHOLMOD(allocate_work) (0, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
  return (FALSE) ;  /* out of memory */
    }

    ASSERT (CHOLMOD(dump_sparse) (A, "etree", Common) >= 0) ;
    Iwork = Common->Iwork ;

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

    ncol = A->ncol ;  /* the number of columns of A */
    nrow = A->nrow ;  /* the number of rows of A */
    Ap = A->p ;   /* size ncol+1, column pointers for A */
    Ai = A->i ;   /* the row indices of A */
    Anz = A->nz ; /* number of nonzeros in each column of A */
    packed = A->packed ;
    Ancestor = Iwork ;  /* size ncol (i/i/l) */

    for (j = 0 ; j < ncol ; j++)
    {
  Parent [j] = EMPTY ;
  Ancestor [j] = EMPTY ;
    }

    /* ---------------------------------------------------------------------- */
    /* compute the etree */
    /* ---------------------------------------------------------------------- */

    if (stype > 0)
    {

  /* ------------------------------------------------------------------ */
  /* symmetric (upper) case: compute etree (A) */
  /* ------------------------------------------------------------------ */

  for (j = 0 ; j < ncol ; j++)
  {
      /* for each row i in column j of triu(A), excluding the diagonal */
      p = Ap [j] ;
      pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
      for ( ; p < pend ; p++)
      {
    i = Ai [p] ;
    if (i < j)
    {
        amesos_update_etree (i, j, Parent, Ancestor) ;
    }
      }
  }

    }
    else if (stype == 0)
    {

  /* ------------------------------------------------------------------ */
  /* unsymmetric case: compute etree (A'*A) */
  /* ------------------------------------------------------------------ */

  Prev = Iwork + ncol ; /* size nrow (i/i/l) */
  for (i = 0 ; i < nrow ; i++)
  {
      Prev [i] = EMPTY ;
  }
  for (j = 0 ; j < ncol ; j++)
  {
      /* for each row i in column j of A */
      p = Ap [j] ;
      pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
      for ( ; p < pend ; p++)
      {
    /* a graph is constructed dynamically with one path per row
     * of A.  If the ith row of A contains column indices
     * (j1,j2,j3,j4) then the new graph has edges (j1,j2), (j2,j3),
     * and (j3,j4).  When at node i of this path-graph, all edges
     * (jprev,j) are considered, where jprev<j */
    i = Ai [p] ;
    jprev = Prev [i] ;
    if (jprev != EMPTY)
    {
        amesos_update_etree (jprev, j, Parent, Ancestor) ;
    }
    Prev [i] = j ;
      }
  }

    }
    else
    {

  /* ------------------------------------------------------------------ */
  /* symmetric case with lower triangular part not supported */
  /* ------------------------------------------------------------------ */

  ERROR (CHOLMOD_INVALID, "symmetric lower not supported") ;
  return (FALSE) ;
    }

    ASSERT (CHOLMOD(dump_parent) (Parent, ncol, "Parent", Common)) ;
    return (TRUE) ;
}
#endif