The R Project SVN R-packages

Rev

Rev 4428 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed 

/* ========================================================================== */
/* === Cholesky/cholmod_factorize =========================================== */
/* ========================================================================== */

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

/* Computes the numerical factorization of a symmetric matrix.  The primary
 * inputs to this routine are a sparse matrix A and the symbolic factor L from
 * cholmod_analyze or a prior numerical factor L.  If A is symmetric, this
 * routine factorizes A(p,p)+beta*I (beta can be zero), where p is the
 * fill-reducing permutation (L->Perm).  If A is unsymmetric, either
 * A(p,:)*A(p,:)'+beta*I or A(p,f)*A(p,f)'+beta*I is factorized.  The set f and
 * the nonzero pattern of the matrix A must be the same as the matrix passed to
 * cholmod_analyze for the supernodal case.  For the simplicial case, it can
 * be different, but it should be the same for best performance.  beta is real.
 *
 * A simplicial factorization or supernodal factorization is chosen, based on
 * the type of the factor L.  If L->is_super is TRUE, a supernodal LL'
 * factorization is computed.  Otherwise, a simplicial numeric factorization
 * is computed, either LL' or LDL', depending on Common->final_ll.
 *
 * Once the factorization is complete, it can be left as is or optionally
 * converted into any simplicial numeric type, depending on the
 * Common->final_* parameters.  If converted from a supernodal to simplicial
 * type, and the Common->final_resymbol parameter is true, then numerically
 * zero entries in L due to relaxed supernodal amalgamation are removed from
 * the simplicial factor (they are always left in the supernodal form of L).
 * Entries that are numerically zero but present in the simplicial symbolic
 * pattern of L are left in place (that is, the graph of L remains chordal).
 * This is required for the update/downdate/rowadd/rowdel routines to work
 * properly.
 *
 * workspace: Flag (nrow), Head (nrow+1),
 *  if symmetric:   Iwork (2*nrow+2*nsuper)
 *  if unsymmetric: Iwork (2*nrow+MAX(2*nsuper,ncol))
 *      where nsuper is 0 if simplicial, or the # of relaxed supernodes in
 *      L otherwise (nsuper <= nrow).
 *  if simplicial: W (nrow).
 *  Allocates up to two temporary copies of its input matrix (including
 *  both pattern and numerical values).
 *
 * If the matrix is not positive definite the routine returns TRUE, but
 * sets Common->status to CHOLMOD_NOT_POSDEF and L->minor is set to the
 * column at which the failure occurred.  Columns L->minor to L->n-1 are
 * set to zero.
 *
 * Supports any xtype (pattern, real, complex, or zomplex), except that the
 * input matrix A cannot be pattern-only.  If L is simplicial, its numeric
 * xtype matches A on output.  If L is supernodal, its xtype is real if A is
 * real, or complex if A is complex or zomplex.
 */

#ifndef NCHOLESKY

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

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


/* ========================================================================== */
/* === cholmod_factorize ==================================================== */
/* ========================================================================== */

/* Factorizes PAP' (or PAA'P' if A->stype is 0), using a factor obtained
 * from cholmod_analyze.  The analysis can be re-used simply by calling this
 * routine a second time with another matrix.  A must have the same nonzero
 * pattern as that passed to cholmod_analyze. */

int CHOLMOD(factorize)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to factorize */
    /* ---- in/out --- */
    cholmod_factor *L,  /* resulting factorization */
    /* --------------- */
    cholmod_common *Common
)
{
    double zero [2] ;
    zero [0] = 0 ;
    zero [1] = 0 ;
    return (CHOLMOD(factorize_p) (A, zero, NULL, 0, L, Common)) ;
}


/* ========================================================================== */
/* === cholmod_factorize_p ================================================== */
/* ========================================================================== */

/* Same as cholmod_factorize, but with more options. */

int CHOLMOD(factorize_p)
(
    /* ---- input ---- */
    cholmod_sparse *A,  /* matrix to factorize */
    double beta [2],    /* factorize beta*I+A or beta*I+A'*A */
    Int *fset,      /* subset of 0:(A->ncol)-1 */
    size_t fsize,   /* size of fset */
    /* ---- in/out --- */
    cholmod_factor *L,  /* resulting factorization */
    /* --------------- */
    cholmod_common *Common
)
{
    cholmod_sparse *S, *F, *A1, *A2 ;
    Int nrow, ncol, stype, convert, n, nsuper, grow2, status ;
    size_t s, t, uncol ;
    int ok = TRUE ;

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

    RETURN_IF_NULL_COMMON (FALSE) ;
    RETURN_IF_NULL (A, FALSE) ;
    RETURN_IF_NULL (L, FALSE) ;
    RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
    RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
    nrow = A->nrow ;
    ncol = A->ncol ;
    n = L->n ;
    stype = A->stype ;
    if (L->n != A->nrow)
    {
    ERROR (CHOLMOD_INVALID, "A and L dimensions do not match") ;
    return (FALSE) ;
    }
    if (stype != 0 && nrow != ncol)
    {
    ERROR (CHOLMOD_INVALID, "matrix invalid") ;
    return (FALSE) ;
    }
    DEBUG (CHOLMOD(dump_sparse) (A, "A for cholmod_factorize", Common)) ;
    Common->status = CHOLMOD_OK ;

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

    nsuper = (L->is_super ? L->nsuper : 0) ;
    uncol = ((stype != 0) ? 0 : ncol) ;

    /* s = 2*nrow + MAX (uncol, 2*nsuper) */
    s = CHOLMOD(mult_size_t) (nsuper, 2, &ok) ;
    s = MAX (uncol, s) ;
    t = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
    s = CHOLMOD(add_size_t) (s, t, &ok) ;
    if (!ok)
    {
    ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
    return (FALSE) ;
    }

    CHOLMOD(allocate_work) (nrow, s, 0, Common) ;
    if (Common->status < CHOLMOD_OK)
    {
    return (FALSE) ;
    }

    S  = NULL ;
    F  = NULL ;
    A1 = NULL ;
    A2 = NULL ;

    /* convert to another form when done, if requested */
    convert = !(Common->final_asis) ;

    /* ---------------------------------------------------------------------- */
    /* perform supernodal LL' or simplicial LDL' factorization */
    /* ---------------------------------------------------------------------- */

    if (L->is_super)
    {

#ifndef NSUPERNODAL

    /* ------------------------------------------------------------------ */
    /* supernodal factorization */
    /* ------------------------------------------------------------------ */

    if (L->ordering == CHOLMOD_NATURAL)
    {

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

        if (stype > 0)
        {
        /* S = tril (A'), F not needed */
        /* workspace: Iwork (nrow) */
        A1 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
        S = A1 ;
        }
        else if (stype < 0)
        {
        /* This is the fastest option for the natural ordering */
        /* S = A; F not needed */
        S = A ;
        }
        else
        {
        /* F = A(:,f)' */
        /* workspace: Iwork (nrow) */
        /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
        A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
        F = A1 ;
        /* S = A */
        S = A ;
        }

    }
    else
    {

        /* -------------------------------------------------------------- */
        /* permute the input matrix before factorization */
        /* -------------------------------------------------------------- */

        if (stype > 0)
        {
        /* This is the fastest option for factoring a permuted matrix */
        /* S = tril (PAP'); F not needed */
        /* workspace: Iwork (2*nrow) */
        A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
        S = A1 ;
        }
        else if (stype < 0)
        {
        /* A2 = triu (PAP') */
        /* workspace: Iwork (2*nrow) */
        A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
        /* S = tril (A2'); F not needed */
        /* workspace: Iwork (nrow) */
        A1 = CHOLMOD(ptranspose) (A2, 2, NULL, NULL, 0, Common) ;
        S = A1 ;
        CHOLMOD(free_sparse) (&A2, Common) ;
        ASSERT (A2 == NULL) ;
        }
        else
        {
        /* F = A(p,f)' */
        /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
        A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
        F = A1 ;
        /* S = F' */
        /* workspace: Iwork (nrow) */
        A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
        S = A2 ;
        }
    }

    /* ------------------------------------------------------------------ */
    /* supernodal factorization */
    /* ------------------------------------------------------------------ */

    /* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow+2*nsuper) */
    if (Common->status == CHOLMOD_OK)
    {
        CHOLMOD(super_numeric) (S, F, beta, L, Common) ;
    }
    status = Common->status ;
    ASSERT (IMPLIES (status >= CHOLMOD_OK, L->xtype != CHOLMOD_PATTERN)) ;

    /* ------------------------------------------------------------------ */
    /* convert to final form, if requested */
    /* ------------------------------------------------------------------ */

    if (Common->status >= CHOLMOD_OK && convert)
    {
        /* workspace: none */
        ok = CHOLMOD(change_factor) (L->xtype, Common->final_ll,
            Common->final_super, Common->final_pack,
            Common->final_monotonic, L, Common) ;
        if (ok && Common->final_resymbol && !(L->is_super))
        {
        /* workspace: Flag (nrow), Head (nrow+1),
         *  if symmetric:   Iwork (2*nrow)
         *  if unsymmetric: Iwork (2*nrow+ncol) */
        CHOLMOD(resymbol_noperm) (S, fset, fsize, Common->final_pack,
            L, Common) ;
        }
    }

#else

    /* ------------------------------------------------------------------ */
    /* CHOLMOD Supernodal module not installed */
    /* ------------------------------------------------------------------ */

    status = CHOLMOD_NOT_INSTALLED ;
    ERROR (CHOLMOD_NOT_INSTALLED,"Supernodal module not installed") ;

#endif

    }
    else
    {

    /* ------------------------------------------------------------------ */
    /* simplicial LDL' factorization */
    /* ------------------------------------------------------------------ */

    /* Permute the input matrix A if necessary.  cholmod_rowfac requires
     * triu(A) in column form for the symmetric case, and A in column form
     * for the unsymmetric case (the matrix S).  The unsymmetric case
     * requires A in row form, or equivalently A' in column form (the
     * matrix F).
     */

    if (L->ordering == CHOLMOD_NATURAL)
    {

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

        if (stype > 0)
        {
        /* F is not needed, S = A */
        S = A ;
        }
        else if (stype < 0)
        {
        /* F is not needed, S = A' */
        /* workspace: Iwork (nrow) */
        A2 = CHOLMOD(ptranspose) (A, 2, NULL, NULL, 0, Common) ;
        S = A2 ;
        }
        else
        {
        /* F = A (:,f)' */
        /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
        A1 = CHOLMOD(ptranspose) (A, 2, NULL, fset, fsize, Common) ;
        F = A1 ;
        S = A ;
        }

    }
    else
    {

        /* -------------------------------------------------------------- */
        /* permute the input matrix before factorization */
        /* -------------------------------------------------------------- */

        if (stype > 0)
        {
        /* F = tril (A (p,p)') */
        /* workspace: Iwork (2*nrow) */
        A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
        /* A2 = triu (F') */
        /* workspace: Iwork (nrow) */
        A2 = CHOLMOD(ptranspose) (A1, 2, NULL, NULL, 0, Common) ;
        /* the symmetric case does not need F, free it and set to NULL*/
        CHOLMOD(free_sparse) (&A1, Common) ;
        }
        else if (stype < 0)
        {
        /* A2 = triu (A (p,p)'), F not needed.  This is the fastest
         * way to factorize a matrix using the simplicial routine
         * (cholmod_rowfac). */
        /* workspace: Iwork (2*nrow) */
        A2 = CHOLMOD(ptranspose) (A, 2, L->Perm, NULL, 0, Common) ;
        }
        else
        {
        /* F = A (p,f)' */
        /* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
        A1 = CHOLMOD(ptranspose) (A, 2, L->Perm, fset, fsize, Common) ;
        F = A1 ;
        /* A2 = F' */
        /* workspace: Iwork (nrow) */
        A2 = CHOLMOD(ptranspose) (F, 2, NULL, NULL, 0, Common) ;
        }
        S = A2 ;
    }

    /* ------------------------------------------------------------------ */
    /* simplicial LDL' or LL' factorization */
    /* ------------------------------------------------------------------ */

    /* factorize beta*I+S (symmetric) or beta*I+F*F' (unsymmetric) */
    /* workspace: Flag (nrow), W (nrow), Iwork (2*nrow) */
    if (Common->status == CHOLMOD_OK)
    {
        grow2 = Common->grow2 ;
        L->is_ll = BOOLEAN (Common->final_ll) ;
        if (L->xtype == CHOLMOD_PATTERN && Common->final_pack)
        {
        /* allocate a factor with exactly the space required */
        Common->grow2 = 0 ;
        }
        CHOLMOD(rowfac) (S, F, beta, 0, nrow, L, Common) ;
        Common->grow2 = grow2 ;
    }
    status = Common->status ;

    /* ------------------------------------------------------------------ */
    /* convert to final form, if requested */
    /* ------------------------------------------------------------------ */

    if (Common->status >= CHOLMOD_OK && convert)
    {
        /* workspace: none */
        CHOLMOD(change_factor) (L->xtype, L->is_ll, FALSE,
            Common->final_pack, Common->final_monotonic, L, Common) ;
    }
    }

    /* ---------------------------------------------------------------------- */
    /* free A1 and A2 if they exist */
    /* ---------------------------------------------------------------------- */

    CHOLMOD(free_sparse) (&A1, Common) ;
    CHOLMOD(free_sparse) (&A2, Common) ;
    Common->status = MAX (Common->status, status) ;
    return (Common->status >= CHOLMOD_OK) ;
}
#endif