/* ========================================================================== */
/* === Tcov/solve =========================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Tcov Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Tcov 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
* -------------------------------------------------------------------------- */
/* Test CHOLMOD for solving various systems of linear equations. */
#include "cm.h"
#define NFTYPES 17
Int ll_types [NFTYPES] = { 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0 } ;
Int pk_types [NFTYPES] = { 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0 } ;
Int mn_types [NFTYPES] = { 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0 } ;
Int co_types [NFTYPES] = { 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1 } ;
#define NRHS 9
#ifndef NDEBUG
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN int cholmod_dump, cholmod_l_dump ;
#endif
/* ========================================================================== */
/* === test_solver ========================================================== */
/* ========================================================================== */
/* Test solve(A) with various control parameters */
double test_solver (cholmod_sparse *A)
{
double err, maxerr = 0 ;
int save ;
for (cm->postorder = 0 ; cm->postorder <= 1 ; cm->postorder++)
{
my_srand (42) ; /* RAND reset */
/* simplicial, no extra memory */
printf ("test_solver: simplicial, no extra memory\n") ;
cm->supernodal = CHOLMOD_SIMPLICIAL ;
cm->grow2 = 0 ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* simplicial, extra space in columns of L */
printf ("test_solver: simplicial, extra space in columns of L\n") ;
cm->grow2 = 5 ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* supernodal */
printf ("test_solver: supernodal\n") ;
cm->supernodal = CHOLMOD_SUPERNODAL ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* supernodal, without final resymbol */
printf ("test_solver: supernodal, without final resymbol\n") ;
cm->final_resymbol = FALSE ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* supernodal, with resymbol, final_super false */
printf ("test_solver: supernodal, with resymbol\n") ;
cm->supernodal = CHOLMOD_SUPERNODAL ;
cm->final_asis = FALSE ;
cm->final_resymbol = TRUE ;
cm->final_super = FALSE ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* supernodal, with resymbol, final_super tree */
printf ("test_solver: supernodal, with resymbol\n") ;
cm->supernodal = CHOLMOD_SUPERNODAL ;
cm->final_asis = FALSE ;
cm->final_resymbol = TRUE ;
cm->final_super = TRUE ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* simplicial LL' */
printf ("test_solver: simplicial LL', try NESDIS instead of METIS\n") ;
cm->supernodal = CHOLMOD_SIMPLICIAL ;
cm->final_ll = TRUE ;
cm->default_nesdis = TRUE ;
err = solve (A) ;
cm->default_nesdis = FALSE ;
MAXERR (maxerr, err, 1) ;
cm->final_ll = FALSE ;
printf ("test_solver err: %6.2e\n", err) ;
/* supernodal, without final resymbol, and no relaxed supernodes */
printf (
"test_solver: supernodal, without final resymbol, and no relaxed supernodes\n") ;
cm->supernodal = CHOLMOD_SUPERNODAL ;
cm->final_asis = TRUE ;
cm->nrelax [0] = 0 ;
cm->nrelax [1] = 0 ;
cm->nrelax [2] = 0 ;
cm->zrelax [0] = 0 ;
cm->zrelax [1] = 0 ;
cm->zrelax [2] = 0 ;
cm->grow0 = 1 ;
cm->grow1 = 1 ;
cm->grow2 = 0 ;
err = solve (A) ;
MAXERR (maxerr, err, 1) ;
printf ("test_solver err: %6.2e\n", err) ;
/* ------------------------------------------------------------------ */
/* restore defaults */
/* ------------------------------------------------------------------ */
cm->dbound = 0.0 ;
cm->grow0 = 1.2 ;
cm->grow1 = 1.2 ;
cm->grow2 = 5 ;
cm->final_asis = TRUE ;
cm->final_super = TRUE ;
cm->final_ll = FALSE ;
cm->final_pack = TRUE ;
cm->final_monotonic = TRUE ;
cm->final_resymbol = FALSE ;
cm->supernodal = CHOLMOD_AUTO ;
cm->nrelax [0] = 4 ;
cm->nrelax [1] = 16 ;
cm->nrelax [2] = 48 ;
cm->zrelax [0] = 0.8 ;
cm->zrelax [1] = 0.1 ;
cm->zrelax [2] = 0.05 ;
/* do not restore these defaults: */
/*
cm->maxrank = ...
cm->metis_memory = 2.0
cm->metis_nswitch = 3000
cm->metis_dswitch = 0.66
cm->print = 3
cm->precise = FALSE
*/
}
progress (1, '.') ;
return (maxerr) ;
}
/* ========================================================================== */
/* === solve ================================================================ */
/* ========================================================================== */
/* solve Ax=b or AA'x=b, systems involving just L, D, or L', and update/downdate
* the system. Returns the worst-case residual. This routine keeps going if
* it runs out of memory (unless the error handler terminates it), because
* it is used both normally and in the memory tests.
*/
double solve (cholmod_sparse *A)
{
double r, enorm, snorm, maxerr = 0, gnorm, anorm, bnorm, xnorm, norm, rcond;
cholmod_factor *L, *Lcopy, *L2 ;
cholmod_sparse *Lmat, *Lo, *D, *Up, *S, *LD, *LDL, *E, *I, *C, *CC, *Ct,
*Ssym, *Cperm, *C2, *S2, *H, *F, *Lo1, *Aboth ;
cholmod_dense *X, *Cdense, *B, *Bcomplex, *Bzomplex, *Breal, *W, *R,
*A3, *C3, *E3 ;
cholmod_sparse *AFt, *AF, *G, *RowK, *Bsparse, *Xsparse ;
double *Cx ;
Int *P, *cset, *fset, *Parent, *Post, *RowCount, *ColCount,
*First, *Level, *rcount, *ccount, *Lp, *Li ;
Int p, i, j, k, n, nrhs, save, save2, csize, rank, nrow, ncol, is_ll,
xtype, isreal, prefer_zomplex, Lxtype, xtype2, save3 ;
int blas_ok = TRUE ;
if (cm->print > 1)
{
printf ("============================================== in solve:\n") ;
}
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "nothing to solve") ;
return (1) ;
}
/* ---------------------------------------------------------------------- */
/* construct right-hand-side (Ax=b if symmetric, AA'x=b otherwise) */
/* ---------------------------------------------------------------------- */
n = A->nrow ;
nrow = A->nrow ;
ncol = A->ncol ;
xtype = A->xtype ;
isreal = (xtype == CHOLMOD_REAL) ;
B = rhs (A, NRHS, n) ;
anorm = CHOLMOD(norm_sparse) (A, 1, cm) ;
save = cm->final_asis ;
cm->final_asis = TRUE ;
/* contents of these will be revised later */
Bzomplex = CHOLMOD(copy_dense) (B, cm) ;
Bcomplex = CHOLMOD(copy_dense) (Bzomplex, cm) ;
Breal = CHOLMOD(copy_dense) (Bzomplex, cm) ;
/* ---------------------------------------------------------------------- */
/* analyze */
/* ---------------------------------------------------------------------- */
if (n < 100 && cm->nmethods == 1 && cm->method[0].ordering == CHOLMOD_GIVEN)
{
Int *UserPerm = prand (nrow) ; /* RAND */
L = CHOLMOD(analyze_p) (A, UserPerm, NULL, 0, cm) ;
OK (CHOLMOD(print_common) ("with UserPerm", cm)) ;
CHOLMOD(free) (nrow, sizeof (Int), UserPerm, cm) ;
}
else
{
L = CHOLMOD(analyze) (A, cm) ;
}
/* test rowadd on a symbolic factor */
if (isreal)
{
RowK = CHOLMOD(spzeros) (n, 1, 0, CHOLMOD_REAL, cm) ;
Lcopy = CHOLMOD(copy_factor) (L, cm) ;
if (n > 0)
{
CHOLMOD(rowadd) (0, RowK, Lcopy, cm) ;
CHOLMOD(check_factor) (Lcopy, cm) ;
CHOLMOD(print_factor) (Lcopy, "Lcopy, now numeric", cm) ;
}
CHOLMOD(free_sparse) (&RowK, cm) ;
CHOLMOD(free_factor) (&Lcopy, cm) ;
}
/* ---------------------------------------------------------------------- */
/* factorize */
/* ---------------------------------------------------------------------- */
CHOLMOD(factorize) (A, L, cm) ;
/* ---------------------------------------------------------------------- */
/* various solves */
/* ---------------------------------------------------------------------- */
if (B != NULL)
{
/* B->ncol = 1 ; */
prefer_zomplex = (B->xtype == CHOLMOD_ZOMPLEX) ;
}
else
{
prefer_zomplex = FALSE ;
}
for (k = 0 ; k <= 1 ; k++)
{
cm->prefer_zomplex = k ;
/* compute the residual, X complex if L or B not real */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
r = resid (A, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* zomplex right-hand-side */
CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, Bzomplex, cm) ;
if (Bzomplex != NULL && B != NULL && B->xtype == CHOLMOD_REAL
&& Bzomplex->xtype == CHOLMOD_ZOMPLEX)
{
/* add an arbitrary imaginary part */
double *Bz = Bzomplex->z ;
for (j = 0 ; j < NRHS ; j++)
{
for (i = 0 ; i < n ; i++)
{
Bz [i+j*n] = (double) (i+j*n) ;
}
}
}
X = CHOLMOD(solve) (CHOLMOD_A, L, Bzomplex, cm) ;
r = resid (A, X, Bzomplex) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* complex right-hand-side */
CHOLMOD(dense_xtype) (CHOLMOD_COMPLEX, Bcomplex, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, Bcomplex, cm) ;
r = resid (A, X, Bcomplex) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* real right-hand-side */
CHOLMOD(dense_xtype) (CHOLMOD_REAL, Breal, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, Breal, cm) ;
r = resid (A, X, Breal) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* sparse solve of Ax=b, b real */
Bsparse = CHOLMOD(dense_to_sparse) (Breal, TRUE, cm) ;
Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
r = resid_sparse (A, Xsparse, Bsparse) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_sparse) (&Bsparse, cm) ;
CHOLMOD(free_sparse) (&Xsparse, cm) ;
/* sparse solve of Ax=b, b complex */
Bsparse = CHOLMOD(dense_to_sparse) (Bcomplex, TRUE, cm) ;
Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
r = resid_sparse (A, Xsparse, Bsparse) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_sparse) (&Bsparse, cm) ;
CHOLMOD(free_sparse) (&Xsparse, cm) ;
/* sparse solve of Ax=b, b zomplex */
Bsparse = CHOLMOD(dense_to_sparse) (Bzomplex, TRUE, cm) ;
Xsparse = CHOLMOD(spsolve) (CHOLMOD_A, L, Bsparse, cm) ;
r = resid_sparse (A, Xsparse, Bsparse) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_sparse) (&Bsparse, cm) ;
CHOLMOD(free_sparse) (&Xsparse, cm) ;
}
cm->prefer_zomplex = FALSE ;
/* ---------------------------------------------------------------------- */
/* sparse solve to compute inv(A) */
/* ---------------------------------------------------------------------- */
CHOLMOD(print_sparse) (A, "A", cm) ;
CHOLMOD(print_factor) (L, "L", cm) ;
rcond = CHOLMOD(rcond) (L, cm) ;
if (cm->print > 1)
{
printf ("rcond: %g\n", rcond) ;
}
if (n < 100 && A->stype != 0)
{
/* solve A*C=I, so C should equal A inverse */
I = CHOLMOD(speye) (n, n, CHOLMOD_REAL, cm) ;
C = CHOLMOD(spsolve) (CHOLMOD_A, L, I, cm) ;
/* compute norm of A*C-I */
if (isreal && n > 10)
{
/* A and C are large and real */
E = CHOLMOD(ssmult) (A, C, 0, TRUE, FALSE, cm) ;
F = CHOLMOD(add) (E, I, minusone, one, TRUE, FALSE, cm) ;
r = CHOLMOD(norm_sparse) (F, 1, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&F, cm) ;
}
else
{
/* There is no complex ssmult or add, so use the BLAS.
* Also test sparse_to_dense for small symmetric matrices. */
A3 = CHOLMOD(sparse_to_dense) (A, cm) ;
C3 = CHOLMOD(sparse_to_dense) (C, cm) ;
xtype2 = isreal ? CHOLMOD_REAL : CHOLMOD_COMPLEX ;
CHOLMOD(dense_xtype) (xtype2, A3, cm) ;
CHOLMOD(dense_xtype) (xtype2, C3, cm) ;
E3 = CHOLMOD(eye) (n, n, xtype2, cm) ;
if (A3 != NULL && C3 != NULL && E3 != NULL)
{
/* E3 = A3*C3-I */
if (isreal)
{
BLAS_dgemm ("N", "N", n, n, n, one, A3->x, n, C3->x, n,
minusone, E3->x, n) ;
}
else
{
BLAS_zgemm ("N", "N", n, n, n, one, A3->x, n, C3->x, n,
minusone, E3->x, n) ;
}
OK (blas_ok) ;
}
r = CHOLMOD(norm_dense) (E3, 1, cm) ;
CHOLMOD(free_dense) (&A3, cm) ;
CHOLMOD(free_dense) (&C3, cm) ;
CHOLMOD(free_dense) (&E3, cm) ;
}
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_sparse) (&I, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
}
/* ---------------------------------------------------------------------- */
/* change complexity of L and solve again; test copy/change routines */
/* ---------------------------------------------------------------------- */
/* change to complex, otherwise leave as */
Lcopy = CHOLMOD(copy_factor) (L, cm) ;
CHOLMOD(factor_xtype) (CHOLMOD_COMPLEX, Lcopy, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, Lcopy, B, cm) ;
r = resid (A, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_factor) (&Lcopy, cm) ;
/* change to zomplex LDL' */
Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
Lcopy = CHOLMOD(copy_factor) (L, cm) ;
CHOLMOD(check_factor) (L, cm) ;
CHOLMOD(print_factor) (Lcopy, "Lcopy", cm) ;
CHOLMOD(change_factor) (Lxtype, FALSE, FALSE, TRUE, TRUE, Lcopy, cm) ;
CHOLMOD(factor_xtype) (CHOLMOD_ZOMPLEX, Lcopy, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, Lcopy, B, cm) ;
r = resid (A, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_factor) (&Lcopy, cm) ;
Lcopy = CHOLMOD(copy_factor) (L, cm) ;
CHOLMOD(change_factor) (Lxtype, TRUE, FALSE, FALSE, FALSE, Lcopy, cm) ;
CHOLMOD(check_factor) (L, cm) ;
CHOLMOD(print_factor) (Lcopy, "Lcopy LL unpacked", cm) ;
CHOLMOD(free_factor) (&Lcopy, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
cm->final_asis = save ;
/* ---------------------------------------------------------------------- */
/* solve again, but use cm->final_asis as given */
/* ---------------------------------------------------------------------- */
if (n < 100 && cm->nmethods == 1 && cm->method[0].ordering == CHOLMOD_GIVEN)
{
Int *UserPerm = prand (nrow) ; /* RAND */
L = CHOLMOD(analyze_p) (A, UserPerm, NULL, 0, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), UserPerm, cm) ;
}
else
{
L = CHOLMOD(analyze) (A, cm) ;
}
CHOLMOD(print_factor) (L, "Lsymbolic", cm) ;
CHOLMOD(factorize) (A, L, cm) ;
/* turn off memory tests [ */
save3 = my_tries ;
my_tries = -1 ;
CHOLMOD(print_factor) (L, "Lnumeric for solver tests", cm) ;
CHOLMOD(print_dense) (B, "B for solver tests", cm) ;
CHOLMOD(print_dense) (Breal, "Breal for solver tests", cm) ;
CHOLMOD(print_dense) (Bcomplex, "Bcomplex for solver tests", cm) ;
CHOLMOD(print_dense) (Bzomplex, "Bzomplex for solver tests", cm) ;
if (B != NULL && Breal != NULL && Bcomplex != NULL && Bzomplex != NULL)
{
for (nrhs = 1 ; nrhs <= NRHS ; nrhs++)
{
for (cm->prefer_zomplex = 0 ; cm->prefer_zomplex <= 1 ;
cm->prefer_zomplex++)
{
B->ncol = nrhs ;
Breal->ncol = nrhs ;
Bcomplex->ncol = nrhs ;
Bzomplex->ncol = nrhs ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
r = resid (A, X, B) ;
CHOLMOD(free_dense) (&X, cm) ;
MAXERR (maxerr, r, 1) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, Breal, cm) ;
r = resid (A, X, Breal) ;
CHOLMOD(free_dense) (&X, cm) ;
MAXERR (maxerr, r, 1) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, Bcomplex, cm) ;
r = resid (A, X, Bcomplex) ;
CHOLMOD(free_dense) (&X, cm) ;
MAXERR (maxerr, r, 1) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, Bzomplex, cm) ;
r = resid (A, X, Bzomplex) ;
CHOLMOD(free_dense) (&X, cm) ;
MAXERR (maxerr, r, 1) ;
}
}
}
cm->prefer_zomplex = FALSE ;
/* turn memory tests back on, where we left off ] */
my_tries = save3 ;
Lcopy = CHOLMOD(copy_factor) (L, cm) ;
/* ---------------------------------------------------------------------- */
/* convert L to LDL' packed or LL packed */
/* ---------------------------------------------------------------------- */
is_ll = (L == NULL) ? FALSE : (L->is_ll) ;
Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
CHOLMOD(change_factor) (Lxtype, is_ll, FALSE, TRUE, TRUE, Lcopy, cm) ;
/* ---------------------------------------------------------------------- */
/* extract L, D, and L' as matrices from Lcopy */
/* ---------------------------------------------------------------------- */
CHOLMOD(resymbol) (A, NULL, 0, TRUE, Lcopy, cm) ;
Lmat = CHOLMOD(factor_to_sparse) (Lcopy, cm) ;
CHOLMOD(check_sparse) (Lmat, cm) ;
I = CHOLMOD(speye) (n, n, CHOLMOD_REAL, cm) ;
Lo = NULL ;
D = NULL ;
Lxtype = (Lmat == NULL) ? CHOLMOD_REAL : (Lmat->xtype) ;
if (isreal)
{
/* use band and add */
if (!is_ll)
{
/* factorization is LDL' = Lo*D*Up */
Lo1 = CHOLMOD(band) (Lmat, -n, -1, 1, cm) ;
Lo = CHOLMOD(add) (Lo1, I, one, one, TRUE, TRUE, cm) ;
CHOLMOD(free_sparse) (&Lo1, cm) ;
D = CHOLMOD(band) (Lmat, 0, 0, 1, cm) ;
}
else
{
/* factorization is LL' = Lo*D*Up */
Lo = CHOLMOD(band) (Lmat, -n, 0, 1, cm) ;
D = CHOLMOD(speye) (n, n, Lxtype, cm) ;
}
}
else
{
/* band and add do not work for c/zomplex matrices*/
D = CHOLMOD(speye) (n, n, Lxtype, cm) ;
Lo = CHOLMOD(copy_sparse) (Lmat, cm) ;
if (!is_ll && D != NULL && Lo != NULL)
{
/* factorization is LDL' = Lo*D*Up */
double *Dx = D->x ;
double *Lx = Lo->x ;
Lp = Lo->p ;
for (k = 0 ; k < n ; k++)
{
p = Lp [k] ;
if (Lxtype == CHOLMOD_COMPLEX)
{
Dx [2*k] = Lx [2*p] ;
Lx [2*p] = 1 ;
}
else
{
Dx [k] = Lx [p] ;
Lx [p] = 1 ;
}
}
}
}
Up = CHOLMOD(transpose) (Lo, 2, cm) ;
/* ---------------------------------------------------------------------- */
/* compute 1-norm of (Lo*D*Up - PAP') or (Lo*D*Up - PAA'P') */
/* ---------------------------------------------------------------------- */
P = (L != NULL) ? (L->Perm) : NULL ;
S = NULL ;
G = NULL ;
if (isreal)
{
if (A->stype == 0)
{
/* G = A*A', try with fset = prand (ncol) */
fset = prand (ncol) ; /* RAND */
AFt = CHOLMOD(ptranspose) (A, 1, NULL, fset, ncol, cm) ;
AF = CHOLMOD(transpose) (AFt, 1, cm) ;
CHOLMOD(free) (ncol, sizeof (Int), fset, cm) ;
G = CHOLMOD(ssmult) (AF, AFt, 0, TRUE, TRUE, cm) ;
/* also try aat */
H = CHOLMOD(aat) (AF, NULL, 0, 1, cm) ;
E = CHOLMOD(add) (G, H, one, minusone, TRUE, FALSE, cm) ;
enorm = CHOLMOD(norm_sparse) (E, 0, cm) ;
gnorm = CHOLMOD(norm_sparse) (G, 0, cm) ;
MAXERR (maxerr, enorm, gnorm) ;
if (cm->print > 1)
{
printf ("enorm %g gnorm %g hnorm %g\n", enorm, gnorm,
CHOLMOD(norm_sparse) (H, 0, cm)) ;
}
if (gnorm > 0)
{
enorm /= gnorm ;
}
OK (enorm < 1e-8) ;
CHOLMOD(free_sparse) (&AFt, cm) ;
CHOLMOD(free_sparse) (&AF, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&H, cm) ;
}
else
{
/* G = A */
G = CHOLMOD(copy) (A, 0, 1, cm) ;
}
if (A->stype == 0)
{
/* S = PAA'P' */
S = CHOLMOD(submatrix) (G, P, n, P, n, TRUE, FALSE, cm) ;
}
else
{
/* S = PAP' */
Aboth = CHOLMOD(copy) (A, 0, 1, cm) ;
S = CHOLMOD(submatrix) (Aboth, P, n, P, n, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&Aboth, cm) ;
}
if (n < NSMALL)
{
/* only do this for small test matrices, since L*D*L' can have many
* nonzero entries */
/* E = L*D*L' - S */
LD = CHOLMOD(ssmult) (Lo, D, 0, TRUE, FALSE, cm) ;
LDL = CHOLMOD(ssmult) (LD, Up, 0, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&LD, cm) ;
E = CHOLMOD(add) (LDL, S, one, minusone, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&LDL, cm) ;
/* e = norm (E) / norm (S) */
enorm = CHOLMOD(norm_sparse) (E, 1, cm) ;
snorm = CHOLMOD(norm_sparse) (S, 0, cm) ;
MAXERR (maxerr, enorm, snorm) ;
CHOLMOD(free_sparse) (&E, cm) ;
}
/* check the row/col counts */
RowCount = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
ColCount = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
Parent = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
Post = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
First = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
Level = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
rcount = CHOLMOD(calloc) (n, sizeof (Int), cm) ;
ccount = CHOLMOD(calloc) (n, sizeof (Int), cm) ;
if (S != NULL && Lmat != NULL && RowCount != NULL && ColCount != NULL &&
Parent != NULL && Post != NULL && First != NULL && Level != NULL &&
rcount != NULL && ccount != NULL)
{
S->stype = 1 ;
CHOLMOD(etree) (S, Parent, cm) ;
CHOLMOD(print_parent) (Parent, n, "Parent", cm) ;
CHOLMOD(postorder) (Parent, n, NULL, Post, cm) ;
CHOLMOD(print_perm) (Post, n, n, "Post", cm) ;
S->stype = -1 ;
CHOLMOD(rowcolcounts) (S, NULL, 0, Parent, Post, RowCount, ColCount,
First, Level, cm) ;
Lp = Lmat->p ;
Li = Lmat->i ;
OK (Lmat->packed) ;
for (j = 0 ; j < n ; j++)
{
for (p = Lp [j] ; p < Lp [j+1] ; p++)
{
i = Li [p] ;
rcount [i]++ ;
ccount [j]++ ;
}
}
/* a singular matrix will only be partially factorized */
if (L->minor == (size_t) n)
{
for (j = 0 ; j < n ; j++)
{
OK (ccount [j] == ColCount [j]) ;
}
}
for (i = 0 ; i < (Int) (L->minor) ; i++)
{
OK (rcount [i] == RowCount [i]) ;
}
}
CHOLMOD(free) (n, sizeof (Int), RowCount, cm) ;
CHOLMOD(free) (n, sizeof (Int), ColCount, cm) ;
CHOLMOD(free) (n, sizeof (Int), Parent, cm) ;
CHOLMOD(free) (n, sizeof (Int), Post, cm) ;
CHOLMOD(free) (n, sizeof (Int), First, cm) ;
CHOLMOD(free) (n, sizeof (Int), Level, cm) ;
CHOLMOD(free) (n, sizeof (Int), rcount, cm) ;
CHOLMOD(free) (n, sizeof (Int), ccount, cm) ;
CHOLMOD(free_sparse) (&S, cm) ;
}
CHOLMOD(free_factor) (&Lcopy, cm) ;
/* ---------------------------------------------------------------------- */
/* solve other systems */
/* ---------------------------------------------------------------------- */
/* turn off memory tests [ */
save3 = my_tries ;
my_tries = -1 ;
for (nrhs = 1 ; nrhs <= 4 ; nrhs++) /* reduced here (6 to 4) */
{
if (B == NULL)
{
break ;
}
B->ncol = nrhs ;
/* solve LDL'x=b */
X = CHOLMOD(solve) (CHOLMOD_LDLt, L, B, cm) ;
r = resid3 (Lo, D, Up, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* solve LDx=b */
X = CHOLMOD(solve) (CHOLMOD_LD, L, B, cm) ;
r = resid3 (Lo, D, NULL, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* solve DL'x=b */
X = CHOLMOD(solve) (CHOLMOD_DLt, L, B, cm) ;
r = resid3 (D, Up, NULL, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* solve Lx=b */
X = CHOLMOD(solve) (CHOLMOD_L, L, B, cm) ;
r = resid3 (Lo, NULL, NULL, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* solve L'x=b */
X = CHOLMOD(solve) (CHOLMOD_Lt, L, B, cm) ;
r = resid3 (Up, NULL, NULL, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* solve Dx=b */
X = CHOLMOD(solve) (CHOLMOD_D, L, B, cm) ;
r = resid3 (D, NULL, NULL, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
save2 = cm->prefer_zomplex ;
for (k = 0 ; k <= 1 ; k++)
{
cm->prefer_zomplex = k ;
/* x=Pb */
X = CHOLMOD(solve) (CHOLMOD_P, L, B, cm) ;
r = pnorm (X, P, B, FALSE) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* x=P'b */
X = CHOLMOD(solve) (CHOLMOD_Pt, L, B, cm) ;
r = pnorm (X, P, B, TRUE) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
}
cm->prefer_zomplex = save2 ;
}
/* turn memory tests back on, where we left off ] */
my_tries = save3 ;
CHOLMOD(free_dense) (&B, cm) ;
CHOLMOD(free_sparse) (&I, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_sparse) (&Lo, cm) ;
CHOLMOD(free_sparse) (&Up, cm) ;
CHOLMOD(free_sparse) (&Lmat, cm) ;
/* ---------------------------------------------------------------------- */
/* update the factorization */
/* ---------------------------------------------------------------------- */
/* turn off memory tests [ */
save3 = my_tries ;
my_tries = -1 ;
B = rhs (A, 1, n) ;
for (rank = 1 ; isreal && rank <= ((n < 100) ? 33 : 2) ; rank++)
{
/* pick a random C */
Cdense = CHOLMOD(zeros) (n, rank, CHOLMOD_REAL, cm) ;
if (Cdense != NULL)
{
Cx = Cdense->x ;
for (k = 0 ; k < 10*rank ; k++)
{
i = nrand (n) ; /* RAND */
j = nrand (rank) ; /* RAND */
Cx [i+j*n] += xrand (1.) ; /* RAND */
}
}
C = CHOLMOD(dense_to_sparse) (Cdense, TRUE, cm) ;
CHOLMOD(free_dense) (&Cdense, cm) ;
/* permute the rows according to L->Perm */
Cperm = CHOLMOD(submatrix) (C, P, n, NULL, -1, TRUE, TRUE, cm) ;
/* update */
CHOLMOD(updown) (TRUE, Cperm, L, cm) ;
CHOLMOD(free_sparse) (&Cperm, cm) ;
/* solve (G+C*C')x=b */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
/* an alternative method would be to compute (G*x + C*(C'*x) - b) */
/* compute S = G+C*C', with no sort */
Ct = CHOLMOD(transpose) (C, 1, cm) ;
CC = CHOLMOD(ssmult) (C, Ct, 0, TRUE, FALSE, cm) ;
S = CHOLMOD(add) (G, CC, one, one, TRUE, TRUE, cm) ;
Ssym = CHOLMOD(copy) (S, 1, 1, cm) ;
CHOLMOD(free_sparse) (&CC, cm) ;
CHOLMOD(free_sparse) (&Ct, cm) ;
/* compute norm (S*x-b) */
r = resid (Ssym, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* ------------------------------------------------------------------ */
/* factor A+CC' from scratch, using same permutation */
/* ------------------------------------------------------------------ */
if (rank == 1)
{
save = cm->nmethods ;
save2 = cm->method [0].ordering ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_GIVEN ;
L2 = CHOLMOD(analyze_p) (Ssym, P, NULL, 0, cm) ;
cm->nmethods = save ;
cm->method [0].ordering = save2 ;
CHOLMOD(factorize) (Ssym, L2, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, L2, B, cm) ;
r = resid (Ssym, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_factor) (&L2, cm) ;
}
CHOLMOD(free_sparse) (&Ssym, cm) ;
/* ------------------------------------------------------------------ */
/* downdate, with just the first half of C */
/* ------------------------------------------------------------------ */
csize = MAX (1, rank / 2) ;
cset = CHOLMOD(malloc) (csize, sizeof (Int), cm) ;
if (cset != NULL)
{
for (i = 0 ; i < csize ; i++)
{
cset [i] = i ;
}
}
C2 = CHOLMOD(submatrix) (C, NULL, -1, cset, csize, TRUE, TRUE, cm) ;
Cperm = CHOLMOD(submatrix) (C2, P, n, NULL, -1, TRUE, TRUE, cm) ;
CHOLMOD(free) (csize, sizeof (Int), cset, cm) ;
CHOLMOD(updown) (FALSE, Cperm, L, cm) ;
CHOLMOD(free_sparse) (&Cperm, cm) ;
/* solve (G+C*C'-C2*C2')x=b */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
/* This is an expensive way to compute the residual. A better
* method would be to compute (G*x + C*(C'*x) - C2*(C2'*x) - b),
* using sdmult. It is done just to test the ssmult
* routine. */
/* compute S2 = G+C*C'-C2*C2', with no sort */
Ct = CHOLMOD(transpose) (C2, 1, cm) ;
CC = CHOLMOD(ssmult) (C2, Ct, 0, TRUE, FALSE, cm) ;
S2 = CHOLMOD(add) (S, CC, one, minusone, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&CC, cm) ;
CHOLMOD(free_sparse) (&Ct, cm) ;
CHOLMOD(free_sparse) (&C2, cm) ;
/* Ssym is a symmetric/upper copy of S2 */
Ssym = CHOLMOD(copy) (S2, 1, 1, cm) ;
CHOLMOD(free_sparse) (&S2, cm) ;
/* compute norm (S2*x-b) */
r = resid (Ssym, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* ------------------------------------------------------------------ */
/* factor S2 scratch, using same permutation */
/* ------------------------------------------------------------------ */
if (rank == 1)
{
save = cm->nmethods ;
save2 = cm->method [0].ordering ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_GIVEN ;
L2 = CHOLMOD(analyze_p) (Ssym, P, NULL, 0, cm) ;
cm->nmethods = save ;
cm->method [0].ordering = save2 ;
CHOLMOD(factorize) (Ssym, L2, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, L2, B, cm) ;
r = resid (Ssym, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_factor) (&L2, cm) ;
}
/* ------------------------------------------------------------------ */
/* re-do the symbolic factorization on L */
/* ------------------------------------------------------------------ */
CHOLMOD(resymbol) (Ssym, NULL, 0, TRUE, L, cm) ;
/* solve (G+C*C'-C2*C2')x=b again */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
/* compute norm (S2*x-b) */
r = resid (Ssym, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_sparse) (&Ssym, cm) ;
/* ------------------------------------------------------------------ */
/* downdate, with the remaining part of C, to get original L */
/* ------------------------------------------------------------------ */
if (rank > csize)
{
cset = CHOLMOD(malloc) (rank-csize, sizeof (Int), cm) ;
if (cset != NULL)
{
for (i = csize ; i < rank ; i++)
{
cset [i-csize] = i ;
}
}
if (rank - csize == 4)
{
/* test the subset print/check routine */
CHOLMOD(print_subset) (cset, rank-csize, rank, "cset", cm) ;
}
C2 = CHOLMOD(submatrix) (C, NULL, -1, cset, rank-csize, TRUE, TRUE,
cm) ;
Cperm = CHOLMOD(submatrix) (C2, P, n, NULL, -1, TRUE, TRUE, cm) ;
CHOLMOD(free) (rank-csize, sizeof (Int), cset, cm) ;
CHOLMOD(updown) (FALSE, Cperm, L, cm) ;
CHOLMOD(free_sparse) (&Cperm, cm) ;
CHOLMOD(free_sparse) (&C2, cm) ;
/* solve the original system */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
/* compute the residual */
r = resid (A, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&S, cm) ;
}
/* turn memory tests back on, where we left off ] */
my_tries = save3 ;
CHOLMOD(free_dense) (&B, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
/* ---------------------------------------------------------------------- */
/* factor again, then change the factor type many times */
/* ---------------------------------------------------------------------- */
C = CHOLMOD(copy_sparse) (A, cm) ;
if (C != NULL)
{
C->sorted = FALSE ;
}
L = CHOLMOD(analyze) (C, cm) ;
CHOLMOD(factorize) (C, L, cm) ;
if (L != NULL && !(L->is_super))
{
CHOLMOD(resymbol) (C, NULL, 0, TRUE, L, cm) ;
}
B = rhs (C, 1, n) ;
cm->prefer_zomplex = prefer_zomplex ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
cm->prefer_zomplex = FALSE ;
r = resid (C, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
for (k = 0 ; k < NFTYPES ; k++)
{
if (co_types [k] && n > 1)
{
/* reallocate column zero of L, to make it non-monotonic */
CHOLMOD(reallocate_column) (0, n, L, cm) ;
}
Lxtype = (L == NULL) ? CHOLMOD_REAL : (L->xtype) ;
CHOLMOD(change_factor) (Lxtype, ll_types [k], FALSE, pk_types [k],
mn_types [k], L, cm) ;
cm->prefer_zomplex = prefer_zomplex ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
cm->prefer_zomplex = FALSE ;
r = resid (C, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
}
/* reallocate a column and solve again */
if (n > 3)
{
CHOLMOD(change_factor) (CHOLMOD_REAL, FALSE, FALSE, TRUE, TRUE, L, cm) ;
CHOLMOD(reallocate_column) (0, n, L, cm) ;
cm->prefer_zomplex = prefer_zomplex ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
cm->prefer_zomplex = FALSE ;
r = resid (C, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
CHOLMOD(free_dense) (&B, cm) ;
CHOLMOD(free_dense) (&Breal, cm) ;
CHOLMOD(free_dense) (&Bcomplex, cm) ;
CHOLMOD(free_dense) (&Bzomplex, cm) ;
/* ---------------------------------------------------------------------- */
/* factorize and solve (A*A'+beta*I)x=b or A'x=b */
/* ---------------------------------------------------------------------- */
if (A->stype == 0)
{
double *Rx, *Rz, *Xx, *Xz ;
double beta [2] ;
beta [0] = 3.14159 ;
beta [1] = 0 ;
L = CHOLMOD(analyze) (A, cm) ;
CHOLMOD(factorize_p) (A, beta, NULL, 0, L, cm) ;
B = rhs (A, 1, n) ;
cm->prefer_zomplex = prefer_zomplex ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
cm->prefer_zomplex = FALSE ;
/* compute the residual */
/* W = A'*X */
W = CHOLMOD(zeros) (ncol, 1, xtype, cm) ;
CHOLMOD(sdmult) (A, 2, one, zero, X, W, cm) ;
/* R = B */
R = CHOLMOD(copy_dense) (B, cm) ;
/* R = A*W - R */
CHOLMOD(sdmult) (A, 0, one, minusone, W, R, cm) ;
/* R = R + beta*X */
if (R != NULL && X != NULL)
{
Rx = R->x ;
Rz = R->z ;
Xx = X->x ;
Xz = X->z ;
for (i = 0 ; i < nrow ; i++)
{
switch (xtype)
{
case CHOLMOD_REAL:
Rx [i] += beta [0] * Xx [i] ;
break ;
case CHOLMOD_COMPLEX:
Rx [2*i ] += beta [0] * Xx [2*i ] ;
Rx [2*i+1] += beta [0] * Xx [2*i+1] ;
break ;
case CHOLMOD_ZOMPLEX:
Rx [i] += beta [0] * Xx [i] ;
Rz [i] += beta [0] * Xz [i] ;
break ;
}
}
}
r = CHOLMOD(norm_dense) (R, 2, cm) ;
bnorm = CHOLMOD(norm_dense) (B, 2, cm) ;
xnorm = CHOLMOD(norm_dense) (X, 2, cm) ;
norm = MAX (r, xnorm) ;
if (norm > 0)
{
r /= norm ;
}
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_dense) (&R, cm) ;
CHOLMOD(free_dense) (&W, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
CHOLMOD(free_dense) (&B, cm) ;
}
/* ---------------------------------------------------------------------- */
/* test rowdel and updown */
/* ---------------------------------------------------------------------- */
if (isreal && A->stype == 1 && n > 0 && n < NLARGE)
{
Int save4, save5, save6 ;
save4 = cm->nmethods ;
save5 = cm->method [0].ordering ;
save6 = cm->supernodal ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_NATURAL ;
cm->supernodal = CHOLMOD_SUPERNODAL ;
B = rhs (A, 1, n) ;
L = CHOLMOD(analyze) (A, cm) ;
CHOLMOD(factorize) (A, L, cm) ;
/* solve Ax=b */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
r = resid (A, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* determine the new system with row/column k missing */
k = n/2 ;
S = CHOLMOD(copy) (A, 0, 1, cm) ;
RowK = CHOLMOD(submatrix) (S, NULL, -1, &k, 1, TRUE, TRUE, cm) ;
CHOLMOD(print_sparse) (S, "S", cm) ;
CHOLMOD(print_sparse) (RowK, "RowK of S", cm) ;
CHOLMOD(free_sparse) (&RowK, cm) ;
prune_row (S, k) ;
if (S != NULL)
{
S->stype = 1 ;
}
/* delete row k of L (converts to LDL') */
/* printf ("rowdel here:\n") ; */
CHOLMOD(rowdel) (k, NULL, L, cm) ;
CHOLMOD(resymbol) (S, NULL, 0, TRUE, L, cm) ;
/* solve with row k missing */
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
r = resid (S, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_sparse) (&S, cm) ;
/* factorize again */
CHOLMOD(free_factor) (&L, cm) ;
L = CHOLMOD(analyze) (A, cm) ;
CHOLMOD(factorize) (A, L, cm) ;
/* rank-3 update (converts to LDL') and solve */
C = CHOLMOD(speye) (n, 3, CHOLMOD_REAL, cm) ;
CC = CHOLMOD(aat) (C, NULL, 0, 1, cm) ;
S = CHOLMOD(add) (A, CC, one, one, TRUE, TRUE, cm) ;
if (S != NULL)
{
S->stype = 1 ;
}
CHOLMOD(updown) (TRUE, C, L, cm) ;
X = CHOLMOD(solve) (CHOLMOD_A, L, B, cm) ;
r = resid (S, X, B) ;
MAXERR (maxerr, r, 1) ;
CHOLMOD(free_dense) (&X, cm) ;
/* free everything */
CHOLMOD(free_sparse) (&S, cm) ;
CHOLMOD(free_sparse) (&CC, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&S, cm) ;
CHOLMOD(free_factor) (&L, cm) ;
CHOLMOD(free_dense) (&B, cm) ;
cm->nmethods = save4 ;
cm->method [0].ordering = save5 ;
cm->supernodal = save6 ;
}
/* ---------------------------------------------------------------------- */
/* free remaining workspace */
/* ---------------------------------------------------------------------- */
OK (CHOLMOD(print_common) ("cm", cm)) ;
progress (0, '.') ;
return (maxerr) ;
}
syntax highlighted by Code2HTML, v. 0.9.1