/* ========================================================================== */
/* === Supernodal/cholmod_super_numeric ===================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Supernodal Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Supernodal Module is licensed under Version 2.0 of the GNU
* General Public License. See gpl.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 Cholesky factorization of A+beta*I or A*F+beta*I. Only the
* the lower triangular part of A+beta*I or A*F+beta*I is accessed. The
* matrices A and F must already be permuted according to the fill-reduction
* permutation L->Perm. cholmod_factorize is an "easy" wrapper for this code
* which applies that permutation. beta is real.
*
* Symmetric case: A is a symmetric (lower) matrix. F is not accessed.
* With a fill-reducing permutation, A(p,p) should be passed instead, where is
* p is L->Perm.
*
* Unsymmetric case: A is unsymmetric, and F must be present. Normally, F=A'.
* With a fill-reducing permutation, A(p,f) and A(p,f)' should be passed as A
* and F, respectively, where f is a list of the subset of the columns of A.
*
* The input factorization L must be supernodal (L->is_super is TRUE). It can
* either be symbolic or numeric. In the first case, L has been analyzed by
* cholmod_analyze or cholmod_super_symbolic, but the matrix has not yet been
* numerically factorized. The numerical values are allocated here and the
* factorization is computed. In the second case, a prior matrix has been
* analyzed and numerically factorized, and a new matrix is being factorized.
* The numerical values of L are replaced with the new numerical factorization.
*
* L->is_ll is ignored, and set to TRUE. This routine always computes an LL'
* factorization. Supernodal LDL' factorization is not (yet) supported.
* FUTURE WORK: perform a supernodal LDL' factorization if L->is_ll is FALSE.
*
* Uses BLAS routines dsyrk, dgemm, dtrsm, and the LAPACK routine dpotrf.
* The supernodal solver uses BLAS routines dtrsv, dgemv, dtrsm, and dgemm.
*
* 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. The supernode containing the non-positive
* diagonal entry is set to zero (this includes columns to the left of L->minor
* in the same supernode), as are all subsequent supernodes.
*
* workspace: Flag (nrow), Head (nrow+1), Iwork (2*nrow + 4*nsuper).
* Allocates temporary space of size L->maxcsize * sizeof(double)
* (twice that for the complex/zomplex case).
*
* If L is supernodal symbolic on input, it is converted to a supernodal numeric
* factor on output, with an xtype of real if A is real, or complex if A is
* complex or zomplex. If L is supernodal numeric on input, its xtype must
* match A (except that L can be complex and A zomplex). The xtype of A and F
* must match.
*/
#ifndef NSUPERNODAL
#include "cholmod_internal.h"
#include "cholmod_supernodal.h"
/* ========================================================================== */
/* === TEMPLATE ============================================================= */
/* ========================================================================== */
#define REAL
#include "t_cholmod_super_numeric.c"
#define COMPLEX
#include "t_cholmod_super_numeric.c"
#define ZOMPLEX
#include "t_cholmod_super_numeric.c"
/* ========================================================================== */
/* === cholmod_super_numeric ================================================ */
/* ========================================================================== */
/* Returns TRUE if successful, or if the matrix is not positive definite.
* Returns FALSE if out of memory, inputs are invalid, or other fatal error
* occurs.
*/
int CHOLMOD(super_numeric)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to factorize */
cholmod_sparse *F, /* F = A' or A(:,f)' */
double beta [2], /* beta*I is added to diagonal of matrix to factorize */
/* ---- in/out --- */
cholmod_factor *L, /* factorization */
/* --------------- */
cholmod_common *Common
)
{
cholmod_dense *C ;
Int *Super, *Map, *SuperMap ;
size_t maxcsize ;
Int nsuper, n, i, k, s, stype, nrow ;
int ok = TRUE, symbolic ;
size_t t, w ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (FALSE) ;
RETURN_IF_NULL (L, FALSE) ;
RETURN_IF_NULL (A, FALSE) ;
RETURN_IF_XTYPE_INVALID (A, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
RETURN_IF_XTYPE_INVALID (L, CHOLMOD_PATTERN, CHOLMOD_COMPLEX, FALSE) ;
stype = A->stype ;
if (stype < 0)
{
if (A->nrow != A->ncol || A->nrow != L->n)
{
ERROR (CHOLMOD_INVALID, "invalid dimensions") ;
return (FALSE) ;
}
}
else if (stype == 0)
{
if (A->nrow != L->n)
{
ERROR (CHOLMOD_INVALID, "invalid dimensions") ;
return (FALSE) ;
}
RETURN_IF_NULL (F, FALSE) ;
RETURN_IF_XTYPE_INVALID (F, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
if (A->nrow != F->ncol || A->ncol != F->nrow || F->stype != 0)
{
ERROR (CHOLMOD_INVALID, "F invalid") ;
return (FALSE) ;
}
if (A->xtype != F->xtype)
{
ERROR (CHOLMOD_INVALID, "A and F must have same xtype") ;
return (FALSE) ;
}
}
else
{
/* symmetric upper case not suppored */
ERROR (CHOLMOD_INVALID, "symmetric upper case not supported") ;
return (FALSE) ;
}
if (!(L->is_super))
{
ERROR (CHOLMOD_INVALID, "L not supernodal") ;
return (FALSE) ;
}
if (L->xtype != CHOLMOD_PATTERN)
{
if (! ((A->xtype == CHOLMOD_REAL && L->xtype == CHOLMOD_REAL)
|| (A->xtype == CHOLMOD_COMPLEX && L->xtype == CHOLMOD_COMPLEX)
|| (A->xtype == CHOLMOD_ZOMPLEX && L->xtype == CHOLMOD_COMPLEX)))
{
ERROR (CHOLMOD_INVALID, "complex type mismatch") ;
return (FALSE) ;
}
}
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* allocate workspace in Common */
/* ---------------------------------------------------------------------- */
nsuper = L->nsuper ;
maxcsize = L->maxcsize ;
nrow = A->nrow ;
n = nrow ;
PRINT1 (("nsuper "ID" maxcsize %g\n", nsuper, (double) maxcsize)) ;
ASSERT (nsuper >= 0 && maxcsize > 0) ;
/* w = 2*n + 4*nsuper */
w = CHOLMOD(mult_size_t) (n, 2, &ok) ;
t = CHOLMOD(mult_size_t) (nsuper, 4, &ok) ;
w = CHOLMOD(add_size_t) (w, t, &ok) ;
if (!ok)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (FALSE) ;
}
CHOLMOD(allocate_work) (n, w, 0, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (FALSE) ;
}
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
/* ---------------------------------------------------------------------- */
/* get the current factor L and allocate numerical part, if needed */
/* ---------------------------------------------------------------------- */
Super = L->super ;
symbolic = (L->xtype == CHOLMOD_PATTERN) ;
if (symbolic)
{
/* convert to supernodal numeric by allocating L->x */
CHOLMOD(change_factor) (
(A->xtype == CHOLMOD_REAL) ? CHOLMOD_REAL : CHOLMOD_COMPLEX,
TRUE, TRUE, TRUE, TRUE, L, Common) ;
if (Common->status < CHOLMOD_OK)
{
/* the factor L remains in symbolic supernodal form */
return (FALSE) ;
}
}
ASSERT (L->dtype == DTYPE) ;
ASSERT (L->xtype == CHOLMOD_REAL || L->xtype == CHOLMOD_COMPLEX) ;
/* supernodal LDL' is not supported */
L->is_ll = TRUE ;
/* ---------------------------------------------------------------------- */
/* get more workspace */
/* ---------------------------------------------------------------------- */
C = CHOLMOD(allocate_dense) (maxcsize, 1, maxcsize, L->xtype, Common) ;
if (Common->status < CHOLMOD_OK)
{
int status = Common->status ;
if (symbolic)
{
/* Change L back to symbolic, since the numeric values are not
* initialized. This cannot fail. */
CHOLMOD(change_factor) (CHOLMOD_PATTERN, TRUE, TRUE, TRUE, TRUE,
L, Common) ;
}
/* the factor L is now back to the form it had on input */
Common->status = status ;
return (FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* get workspace */
/* ---------------------------------------------------------------------- */
SuperMap = Common->Iwork ; /* size n (i/i/l) */
Map = Common->Flag ; /* size n, use Flag as workspace for Map array */
for (i = 0 ; i < n ; i++)
{
Map [i] = EMPTY ;
}
/* ---------------------------------------------------------------------- */
/* find the mapping of nodes to relaxed supernodes */
/* ---------------------------------------------------------------------- */
/* SuperMap [k] = s if column k is contained in supernode s */
for (s = 0 ; s < nsuper ; s++)
{
PRINT1 (("Super ["ID"] "ID" ncols "ID"\n",
s, Super[s], Super[s+1]-Super[s]));
for (k = Super [s] ; k < Super [s+1] ; k++)
{
SuperMap [k] = s ;
PRINT2 (("relaxed SuperMap ["ID"] = "ID"\n", k, SuperMap [k])) ;
}
}
/* ---------------------------------------------------------------------- */
/* supernodal numerical factorization, using template routine */
/* ---------------------------------------------------------------------- */
switch (A->xtype)
{
case CHOLMOD_REAL:
ok = r_cholmod_super_numeric (A, F, beta, L, C, Common) ;
break ;
case CHOLMOD_COMPLEX:
ok = c_cholmod_super_numeric (A, F, beta, L, C, Common) ;
break ;
case CHOLMOD_ZOMPLEX:
/* This operates on complex L, not zomplex */
ok = z_cholmod_super_numeric (A, F, beta, L, C, Common) ;
break ;
}
/* ---------------------------------------------------------------------- */
/* clear Common workspace, free temp workspace C, and return */
/* ---------------------------------------------------------------------- */
/* Flag array was used as workspace, clear it */
Common->mark = EMPTY ;
CHOLMOD(clear_flag) (Common) ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
CHOLMOD(free_dense) (&C, Common) ;
return (ok) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1