/* ========================================================================== */
/* === 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
syntax highlighted by Code2HTML, v. 0.9.1