/* ========================================================================== */
/* === Tcov/test_ops ======================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* 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 matrix operators. */
#include "cm.h"
/* ========================================================================== */
/* === nzdiag =============================================================== */
/* ========================================================================== */
/* Count the entries on the diagonal */
Int nzdiag (cholmod_sparse *A)
{
Int *Ap, *Ai, *Anz ;
Int nnzdiag, packed, j, p, i, pend, ncol ;
if (A == NULL)
{
return (EMPTY) ;
}
nnzdiag = 0 ;
ncol = A->ncol ;
Ap = A->p ;
Ai = A->i ;
Anz = A->nz ;
packed = A->packed ;
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i == j)
{
nnzdiag++ ;
}
}
}
return (nnzdiag) ;
}
/* ========================================================================== */
/* === check_partition ====================================================== */
/* ========================================================================== */
/* Check a node separator, and return the # of nodes in the node separator or
* -1 if the separater is invalid. A node j is in the left part if
* Part [j] = 0, in the right part if Part [j] = 1, and in the separator if
* Part [j] = 2.
*/
Int check_partition (cholmod_sparse *A, Int *Part)
{
Int *Ap, *Ai, *Anz ;
Int chek [3], which, i, j, p, n, pend, packed ;
if (A == NULL || Part == NULL || A->nrow != A->ncol)
{
return (EMPTY) ;
}
n = A->nrow ;
Ap = A->p ;
Ai = A->i ;
Anz = A->nz ;
packed = A->packed ;
chek [0] = 0 ;
chek [1] = 0 ;
chek [2] = 0 ;
for (j = 0 ; j < n ; j++)
{
which = Part [j] ;
p = Ap [j] ;
pend = (packed) ? (Ap [j+1]) : (p + Anz [j]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (which == 0)
{
if (Part [i] == 1)
{
return (EMPTY) ;
}
}
else if (which == 1)
{
if (Part [i] == 0)
{
return (EMPTY) ;
}
}
}
if (which < 0 || which > 2)
{
return (EMPTY) ;
}
chek [which]++ ;
}
return (chek [2]) ;
}
/* ========================================================================== */
/* === check_equality ======================================================= */
/* ========================================================================== */
/* Ensure two sparse matrices are identical. */
static void check_equality (cholmod_sparse *E, cholmod_sparse *D, Int xtype)
{
double *Ex, *Ez, *Dx, *Dz ;
Int *Ep, *Ei, *Dp, *Di ;
Int j, nz, p, ncol ;
if (E == NULL || D == NULL || D->xtype != xtype || E->xtype != xtype)
{
return ;
}
Ep = E->p ;
Ei = E->i ;
Ex = E->x ;
Ez = E->z ;
Dp = D->p ;
Di = D->i ;
Dx = D->x ;
Dz = D->z ;
OK (E->ncol == D->ncol) ;
OK (E->nrow == D->nrow) ;
ncol = E->ncol ;
for (j = 0 ; j <= ncol ; j++)
{
OK (Ep [j] == Dp [j]) ;
}
nz = Ep [ncol] ;
for (p = 0 ; p < nz ; p++)
{
OK (Ei [p] == Di [p]) ;
}
if (xtype == CHOLMOD_REAL)
{
for (p = 0 ; p < nz ; p++)
{
OK (Ex [p] == Dx [p]) ;
}
}
else if (xtype == CHOLMOD_COMPLEX)
{
for (p = 0 ; p < 2*nz ; p++)
{
OK (Ex [p] == Dx [p]) ;
}
}
else if (xtype == CHOLMOD_ZOMPLEX)
{
for (p = 0 ; p < nz ; p++)
{
OK (Ex [p] == Dx [p]) ;
OK (Ez [p] == Dz [p]) ;
}
}
}
/* ========================================================================== */
/* === test_ops ============================================================= */
/* ========================================================================== */
/* Test various matrix operations. */
double test_ops (cholmod_sparse *A)
{
double maxerr = 0, r, x, anorm, r1, rinf ;
double *Sx, *Sz ;
Int *Pinv, *P, *Si, *Sj, *Q, *Qinv, *fset, *Partition ;
cholmod_triplet *S ;
cholmod_sparse *C, *D, *E, *F, *G, *H, *AT, *Zs ;
cholmod_dense *X, *Y ;
Int n, kk, k, nrow, ncol, len, nz, ok, i, j, stype, nmin, mode, isreal,
xtype, xtype2, mtype, asym, xmatched, pmatched, nzoffdiag, nz_diag ;
size_t nz1, nz2 ;
void (*save) (int, char *, int, char *) ;
double alpha [2], beta [2], *Xx ;
FILE *f ;
int option, save3 ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "nothing for test_ops") ;
return (1) ;
}
nrow = A->nrow ;
ncol = A->ncol ;
H = NULL ;
n = MAX (nrow, ncol) ;
nmin = MIN (nrow, ncol) ;
xtype = A->xtype ;
isreal = (A->xtype == CHOLMOD_REAL) ;
/* ---------------------------------------------------------------------- */
/* norm */
/* ---------------------------------------------------------------------- */
CHOLMOD(print_sparse) (A, "A for testops", cm) ;
r1 = CHOLMOD(norm_sparse) (A, 1, cm) ;
rinf = CHOLMOD(norm_sparse) (A, 0, cm) ;
anorm = r1 ;
/* E = pattern of A */
E = CHOLMOD(copy) (A, 0, 0, cm) ;
CHOLMOD(print_sparse) (E, "E = pattern of A", cm) ;
r1 = CHOLMOD(norm_sparse) (E, 1, cm) ;
rinf = CHOLMOD(norm_sparse) (E, 0, cm) ;
OK (r1 <= nrow) ;
OK (rinf <= ncol) ;
CHOLMOD(free_sparse) (&E, cm) ;
/* E = pattern of A, but exclude the diagonal */
E = CHOLMOD(copy) (A, 0, -1, cm) ;
CHOLMOD(print_sparse) (E, "E = spones (A), excl diag", cm) ;
r1 = CHOLMOD(norm_sparse) (E, 1, cm) ;
rinf = CHOLMOD(norm_sparse) (E, 0, cm) ;
if (nrow < ncol)
{
OK (r1 <= nrow) ;
OK (rinf < MAX (1,ncol)) ;
}
else
{
OK (r1 < MAX (1,nrow)) ;
OK (rinf <= ncol) ;
}
CHOLMOD(free_sparse) (&E, cm) ;
/* ---------------------------------------------------------------------- */
/* copy */
/* ---------------------------------------------------------------------- */
if (A->stype)
{
/* E = tril (A), no diagonal */
E = CHOLMOD(copy) (A, -1, -1, cm) ;
CHOLMOD(print_sparse) (E, "E, lower and no diagonal", cm) ;
CHOLMOD(band_inplace) (0, n, 0, E, cm) ;
CHOLMOD(print_sparse) (E, "E Empty", cm) ;
nz = CHOLMOD(nnz) (E, cm) ;
if (E != NULL)
{
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&E, cm) ;
}
/* ---------------------------------------------------------------------- */
/* read/write */
/* ---------------------------------------------------------------------- */
/*
i = cm->print ;
cm->print = 4 ;
CHOLMOD(print_sparse) (A, "A for read/write", cm) ;
cm->print = i ;
*/
/* delete the contents of the temp1.mtx and temp2.mtx file */
f = fopen ("temp1.mtx", "w") ;
fprintf (f, "temp1\n") ;
fclose (f) ;
f = fopen ("temp3.mtx", "w") ;
fprintf (f, "temp3\n") ;
fclose (f) ;
CHOLMOD(free_work) (cm) ;
f = fopen ("temp1.mtx", "w") ;
asym = CHOLMOD(write_sparse) (f, A, NULL, "comments.txt", cm) ;
fclose (f) ;
printf ("write_sparse, asym: %d\n", asym) ;
OK (IMPLIES (A != NULL, asym > EMPTY)) ;
f = fopen ("temp1.mtx", "r") ;
C = CHOLMOD(read_sparse) (f, cm) ;
fclose (f) ;
printf ("got_sparse\n") ;
CHOLMOD(free_sparse) (&C, cm) ;
save3 = A->xtype ;
A->xtype = CHOLMOD_PATTERN ;
f = fopen ("temp3.mtx", "w") ;
asym = CHOLMOD(write_sparse) (f, A, NULL, "comments.txt", cm) ;
A->xtype = save3 ;
fclose (f) ;
printf ("write_sparse3, asym: %d\n", asym) ;
f = fopen ("temp3.mtx", "r") ;
C = CHOLMOD(read_sparse) (f, cm) ;
fclose (f) ;
printf ("got_sparse3\n") ;
CHOLMOD(free_sparse) (&C, cm) ;
for (i = 0 ; i <= 1 ; i++)
{
f = fopen ("temp2.mtx", "w") ;
fprintf (f, "temp2\n") ;
fclose (f) ;
X = CHOLMOD(ones) (4, 4, CHOLMOD_REAL, cm) ;
if (X != NULL)
{
Xx = X->x ;
Xx [0] = (i == 0) ? 1.1e308 : -1.1e308 ;
}
f = fopen ("temp2.mtx", "w") ;
ok = CHOLMOD(write_dense) (f, X, "comments.txt", cm) ;
fclose (f) ;
printf ("wrote dense\n") ;
f = fopen ("temp2.mtx", "r") ;
Y = CHOLMOD(read_dense) (f, cm) ;
fclose (f) ;
printf ("got dense\n") ;
CHOLMOD(free_dense) (&X, cm) ;
CHOLMOD(free_dense) (&Y, cm) ;
}
/* ---------------------------------------------------------------------- */
/* symmetry */
/* ---------------------------------------------------------------------- */
CHOLMOD(free_work) (cm) ;
xmatched = 0 ;
pmatched = 0 ;
nzoffdiag = 0 ;
nz_diag = 0 ;
for (option = 0 ; option <= 2 ; option++)
{
asym = CHOLMOD(symmetry) (A, option, &xmatched, &pmatched, &nzoffdiag,
&nz_diag, cm);
printf ("symmetry, asym: %d matched %d %d offdiag %d diag %d\n", asym,
xmatched, pmatched, nzoffdiag, nz_diag) ;
}
/* ---------------------------------------------------------------------- */
/* transpose */
/* ---------------------------------------------------------------------- */
C = CHOLMOD(allocate_sparse) (A->ncol, A->nrow, A->nzmax, TRUE, TRUE,
-(A->stype), A->xtype, cm) ;
D = CHOLMOD(allocate_sparse) (A->nrow, A->ncol, A->nzmax, TRUE, TRUE,
A->stype, A->xtype, cm) ;
CHOLMOD(free_work) (cm) ;
ok = (C != NULL && D != NULL) ;
/* C = A' */
if (ok)
{
if (A->stype)
{
ok = CHOLMOD(transpose_sym) (A, 2, NULL, C, cm) ;
}
else
{
ok = CHOLMOD(transpose_unsym) (A, 2, NULL, NULL, 0, C, cm) ;
}
OK (ok || cm->status == CHOLMOD_OUT_OF_MEMORY) ;
}
/* D = C' */
if (ok)
{
if (A->stype)
{
ok = CHOLMOD(transpose_sym) (C, 2, NULL, D, cm) ;
}
else
{
ok = CHOLMOD(transpose_unsym) (C, 2, NULL, NULL, 0, D, cm) ;
}
OK (ok || cm->status == CHOLMOD_OUT_OF_MEMORY) ;
}
if (ok)
{
ok = CHOLMOD(check_sparse) (D, cm) ;
OK (ok || cm->status == CHOLMOD_OUT_OF_MEMORY) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
/* ---------------------------------------------------------------------- */
/* C = A with jumbled triplets */
/* ---------------------------------------------------------------------- */
S = CHOLMOD(sparse_to_triplet) (A, cm) ; /* [ */
if (S != NULL && nmin > 0)
{
/* double the number of entries in S */
nz1 = S->nzmax ;
nz2 = 2*nz1 ;
ok = CHOLMOD(reallocate_triplet) (nz2, S, cm) ;
if (ok)
{
/* add duplicate entries, but keep the matrix the same */
OK (S->nzmax == nz2) ;
Si = S->i ;
Sj = S->j ;
Sx = S->x ;
Sz = S->z ;
for (k = nz1 ; k < ((Int) nz2) ; k++)
{
kk = nrand (k) ; /* RAND */
Si [k] = Si [kk] ;
Sj [k] = Sj [kk] ;
if (S->xtype == CHOLMOD_REAL)
{
x = Sx [kk] * (xrand (4.) - 2) ; /* RAND */
Sx [k] = x ;
Sx [kk] -= x ;
}
else if (S->xtype == CHOLMOD_COMPLEX)
{
x = Sx [2*kk] * (xrand (4.) - 2) ; /* RAND */
Sx [2*k] = x ;
Sx [2*kk] -= x ;
x = Sx [2*kk+1] * (xrand (4.) - 2) ; /* RAND */
Sx [2*k+1] = x ;
Sx [2*kk+1] -= x ;
}
else
{
x = Sx [kk] * (xrand (4.) - 2) ; /* RAND */
Sx [k] = x ;
Sx [kk] -= x ;
x = Sz [kk] * (xrand (4.) - 2) ; /* RAND */
Sz [k] = x ;
Sz [kk] -= x ;
}
}
/* randomly jumble the entries */
for (k = 0 ; k < ((Int) (nz2-1)) ; k++)
{
kk = k + nrand (nz2-k) ; /* RAND */
i = Si [k] ;
Si [k] = Si [kk] ;
Si [kk] = i ;
j = Sj [k] ;
Sj [k] = Sj [kk] ;
Sj [kk] = j ;
if (S->xtype == CHOLMOD_REAL)
{
x = Sx [k] ;
Sx [k] = Sx [kk] ;
Sx [kk] = x ;
}
else if (S->xtype == CHOLMOD_COMPLEX)
{
x = Sx [2*k] ;
Sx [2*k] = Sx [2*kk] ;
Sx [2*kk] = x ;
x = Sx [2*k+1] ;
Sx [2*k+1] = Sx [2*kk+1] ;
Sx [2*kk+1] = x ;
}
else
{
x = Sx [k] ;
Sx [k] = Sx [kk] ;
Sx [kk] = x ;
x = Sz [k] ;
Sz [k] = Sz [kk] ;
Sz [kk] = x ;
}
}
S->nnz = nz2 ;
}
else
{
OK (S->nzmax == nz1) ;
OK (S->nnz == nz1) ;
}
}
CHOLMOD(print_triplet) (S, "S jumbled", cm) ;
C = CHOLMOD(triplet_to_sparse) (S, 0, cm) ; /* [ */
CHOLMOD(print_sparse) (A, "A", cm) ;
CHOLMOD(print_sparse) (C, "C", cm) ;
Zs = CHOLMOD(spzeros) (nrow, ncol, 1, xtype, cm) ; /* [ */
G = NULL ;
F = NULL ;
if (isreal)
{
/* G=A+0 */
G = CHOLMOD(add) (A, Zs, one, one, TRUE, TRUE, cm) ; /* [ */
/* F = G-C */
F = CHOLMOD(add) (G, C, one, minusone, TRUE, TRUE, cm) ; /* [ */
CHOLMOD(print_sparse) (F, "F", cm) ;
r = CHOLMOD(norm_sparse) (F, 1, cm) ;
MAXERR (maxerr, r, anorm) ;
r = CHOLMOD(norm_sparse) (F, 0, cm) ;
CHOLMOD(drop) (0, F, cm) ;
rinf = CHOLMOD(norm_sparse) (F, 0, cm) ;
if (F != NULL)
{
OK (r == rinf) ;
}
MAXERR (maxerr, r, anorm) ;
MAXERR (maxerr, rinf, anorm) ;
/* E = F, with change of type and dropping small entries */
for (stype = -1 ; stype <= 1 ; stype++)
{
if (stype != 0 && (F != NULL && F->nrow != F->ncol))
{
continue ;
}
for (mode = 0 ; mode <= 1 ; mode++)
{
E = CHOLMOD(copy) (F, stype, mode, cm) ; /* [ */
CHOLMOD(drop) (1e-16, E, cm) ;
r1 = CHOLMOD(norm_sparse) (E, 1, cm) ;
rinf = CHOLMOD(norm_sparse) (E, 0, cm) ;
if (E != NULL)
{
if (mode == 0)
{
/* pattern only */
OK (r1 <= nrow) ;
OK (rinf <= ncol) ;
}
else
{
MAXERR (maxerr, r1, anorm) ;
MAXERR (maxerr, rinf, anorm) ;
}
}
CHOLMOD(free_sparse) (&E, cm) ; /* ] */
}
}
CHOLMOD(free_sparse) (&F, cm) ; /* ] */
CHOLMOD(free_sparse) (&G, cm) ; /* ] */
}
Y = CHOLMOD(ones) (nrow, 1, xtype, cm) ; /* [ */
X = CHOLMOD(ones) (ncol, 1, xtype, cm) ; /* [ */
alpha [0] = 0 ;
alpha [1] = 0 ;
beta [0] = 2 ;
beta [1] = 0 ;
/* Y = 0*A*X + 2*Y */
CHOLMOD(sdmult) (A, FALSE, alpha, beta, X, Y, cm) ;
r = CHOLMOD(norm_dense) (Y, 0, cm) ;
if (Y != NULL && X != NULL && A != NULL)
{
OK ((nrow == 0) ? (r == 0) : (r == 2)) ;
}
alpha [0] = 1 ;
/* Y = 1*(0)*X + 2*Y */
CHOLMOD(sdmult) (Zs, FALSE, alpha, beta, X, Y, cm) ;
r = CHOLMOD(norm_dense) (Y, 0, cm) ;
if (Y != NULL && X != NULL && A != NULL)
{
OK ((nrow == 0) ? (r == 0) : (r == 4)) ;
}
CHOLMOD(free_dense) (&X, cm) ; /* ] */
CHOLMOD(free_dense) (&Y, cm) ; /* ] */
CHOLMOD(free_sparse) (&Zs, cm) ; /* ] */
CHOLMOD(free_sparse) (&C, cm) ; /* ] */
CHOLMOD(free_triplet) (&S, cm) ; /* ] */
/* ---------------------------------------------------------------------- */
/* C = P*A*Q in triplet form */
/* ---------------------------------------------------------------------- */
S = CHOLMOD(sparse_to_triplet) (A, cm) ;
P = prand (nrow) ; /* RAND */
if (A->stype == 0)
{
Q = prand (ncol) ; /* RAND */
}
else
{
/* if A is symmetric, and stored in either upper or lower form, then
* the following code only works if P = Q */
Q = P ;
}
Pinv = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
Qinv = CHOLMOD(malloc) (ncol, sizeof (Int), cm) ;
Partition = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
if (Pinv != NULL && Qinv != NULL && P != NULL && Q != NULL && S != NULL)
{
Si = S->i ;
Sj = S->j ;
nz = S->nnz ;
for (k = 0 ; k < nrow ; k++)
{
Pinv [P [k]] = k ;
}
for (k = 0 ; k < ncol ; k++)
{
Qinv [Q [k]] = k ;
}
for (k = 0 ; k < nz ; k++)
{
Si [k] = Pinv [Si [k]] ;
}
for (k = 0 ; k < nz ; k++)
{
Sj [k] = Qinv [Sj [k]] ;
}
}
C = CHOLMOD(triplet_to_sparse) (S, 0, cm) ;
D = NULL ;
E = NULL ;
F = NULL ;
G = NULL ;
if (isreal)
{
/* ------------------------------------------------------------------ */
/* E = P*A*Q in sparse form */
/* ------------------------------------------------------------------ */
D = CHOLMOD(copy) (A, 0, 1, cm) ;
E = CHOLMOD(submatrix) (D, P, nrow, Q, ncol, TRUE, FALSE, cm) ;
CHOLMOD(sort) (E, cm) ;
/* ------------------------------------------------------------------ */
/* F = E-G */
/* ------------------------------------------------------------------ */
G = CHOLMOD(copy) (C, 0, 1, cm) ;
F = CHOLMOD(add) (E, G, one, minusone, TRUE, TRUE, cm) ;
CHOLMOD(drop) (0, F, cm) ;
nz = CHOLMOD(nnz) (F, cm) ;
if (F != NULL)
{
OK (nz == 0) ;
}
}
CHOLMOD(free_sparse) (&F, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&H, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_triplet) (&S, cm) ;
/* ---------------------------------------------------------------------- */
/* submatrix */
/* ---------------------------------------------------------------------- */
if (A->stype == 0 && isreal)
{
/* E = A(:,:) */
E = CHOLMOD(submatrix) (A, NULL, -1, NULL, -1, TRUE, TRUE, cm) ;
/* C = A-E */
C = CHOLMOD(add) (A, E, one, minusone, TRUE, TRUE, cm) ;
ok = CHOLMOD(drop) (0., C, cm) ;
nz = CHOLMOD(nnz) (C, cm) ;
if (C != NULL)
{
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
}
/* ---------------------------------------------------------------------- */
/* test band and add, unsymmetric */
/* ---------------------------------------------------------------------- */
if (isreal)
{
CHOLMOD(print_sparse) (A, "A for do triplet", cm) ;
/* E = A */
E = CHOLMOD(copy) (A, 0, 1, cm) ;
CHOLMOD(print_sparse) (E, "E=triu(A)", cm) ;
/* E = triu (E) */
CHOLMOD(band_inplace) (0, ncol, 1, E, cm) ;
/* F = A */
F = CHOLMOD(copy) (A, 0, 1, cm) ;
CHOLMOD(print_sparse) (F, "F=tril(A)", cm) ;
/* F = tril(F,-1) */
CHOLMOD(band_inplace) (-nrow, -1, 1, F, cm) ;
CHOLMOD(print_sparse) (F, "Ftril", cm) ;
/* G = E+F */
G = CHOLMOD(add) (E, F, one, one, TRUE, TRUE, cm) ;
CHOLMOD(print_sparse) (G, "G=E+F", cm) ;
/* D = A-G, which should be empty */
D = CHOLMOD(add) (G, A, one, minusone, TRUE, TRUE, cm) ;
CHOLMOD(print_sparse) (D, "D=A-G", cm) ;
CHOLMOD(drop) (0, D, cm) ;
CHOLMOD(print_sparse) (D, "D drop", cm) ;
nz = CHOLMOD(nnz) (D, cm) ;
if (D != NULL)
{
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&F, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
D = CHOLMOD(band) (A, 1, -1, 0, cm) ;
nz = CHOLMOD(nnz) (D, cm) ;
if (D != NULL)
{
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&D, cm) ;
D = CHOLMOD(band) (A, 0, 0, 0, cm) ;
nz = CHOLMOD(nnz) (D, cm) ;
if (D != NULL)
{
OK (nz == nzdiag (D)) ;
}
CHOLMOD(free_sparse) (&D, cm) ;
}
/* ---------------------------------------------------------------------- */
/* test band, add and copy_sparse (symmetric) */
/* ---------------------------------------------------------------------- */
if (A->stype && isreal)
{
/* E = A, in symmetric/upper form */
E = CHOLMOD(copy) (A, 1, 1, cm) ;
CHOLMOD(print_sparse) (E, "E=A in sym/upper form", cm) ;
/* E = -E */
CHOLMOD(scale) (M1, CHOLMOD_SCALAR, E, cm) ;
/* F = A, in symmetric/lower form */
F = CHOLMOD(copy) (A, -1, 1, cm) ;
CHOLMOD(print_sparse) (F, "F=A in sym/lower form", cm) ;
/* C = F (exact copy) */
C = CHOLMOD(copy_sparse) (F, cm) ;
/* G = E+C */
G = CHOLMOD(add) (E, C, one, one, TRUE, FALSE, cm) ;
CHOLMOD(print_sparse) (G, "G=E+F", cm) ;
CHOLMOD(sort) (G, cm) ;
CHOLMOD(drop) (0, G, cm) ;
CHOLMOD(print_sparse) (G, "G drop", cm) ;
nz = CHOLMOD(nnz) (G, cm) ;
if (G != NULL)
{
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&F, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
}
/* ---------------------------------------------------------------------- */
/* try a dense identity matrix */
/* ---------------------------------------------------------------------- */
X = CHOLMOD(eye) (3, 4, CHOLMOD_REAL, cm) ;
CHOLMOD(print_dense) (X, "Dense identity", cm) ;
CHOLMOD(free_dense) (&X, cm) ;
/* ---------------------------------------------------------------------- */
/* bisector and nested_dissection */
/* ---------------------------------------------------------------------- */
#ifndef NPARTITION
if (A != NULL && A->nrow == A->ncol)
{
UF_long nc, nc_new ;
Int cnz, csep, save2 ;
Int *Cnw, *Cew, *Cmember, *CParent, *Perm ;
double save1 ;
/* try CHOLMOD's interface to METIS_NodeComputeSeparator */
cm->metis_memory = 2.0 ;
CHOLMOD(print_sparse) (A, "A for bisect", cm) ;
csep = CHOLMOD(bisect) (A, NULL, 0, TRUE, Partition, cm) ;
if (csep != EMPTY)
{
OK (csep == check_partition (A, Partition)) ;
}
/* try the raw interface to METIS_NodeComputeSeparator */
CHOLMOD(print_sparse) (A, "A for metis bisect", cm) ;
/* C = A+A', remove the diagonal */
AT = CHOLMOD(transpose) (A, 0, cm) ;
E = CHOLMOD(add) (A, AT, one, one, FALSE, TRUE, cm) ;
CHOLMOD(free_sparse) (&AT, cm) ;
C = CHOLMOD(copy) (E, 0, -1, cm) ;
CHOLMOD(print_sparse) (C, "C for metis bisect", cm) ;
cnz = (C != NULL) ? (C->nzmax) : 0 ;
Cew = CHOLMOD(malloc) (cnz, sizeof (Int), cm) ;
Cnw = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
if (Cnw != NULL)
{
for (j = 0 ; j < (Int) (A->nrow) ; j++)
{
Cnw [j] = 1 ;
}
}
if (Cew != NULL)
{
for (j = 0 ; j < cnz ; j++)
{
Cew [j] = 1 ;
}
}
csep = CHOLMOD(metis_bisector) (C, Cnw, Cew, Partition, cm) ;
if (csep != EMPTY)
{
OK (csep == check_partition (C, Partition)) ;
}
CHOLMOD(free) (nrow, sizeof (Int), Cnw, cm) ;
CHOLMOD(free) (cnz, sizeof (Int), Cew, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
Cmember = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
CParent = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
Perm = CHOLMOD(malloc) (nrow, sizeof (Int), cm) ;
save1 = cm->method [cm->current].nd_small ;
save2 = cm->method [cm->current].nd_oksep ;
cm->method [cm->current].nd_small = 1 ;
cm->method [cm->current].nd_oksep = 1.0 ;
nc = CHOLMOD(nested_dissection) (A, NULL, 0, Perm, CParent, Cmember,
cm);
if (nc > 0)
{
OK (CHOLMOD(check_perm) (Perm, n, n, cm)) ;
}
CHOLMOD(free_work) (cm) ;
/* collapse the septree */
if (nc > 0 && n > 0)
{
nc_new = CHOLMOD(collapse_septree) (n, nc, 0.1, 400,
CParent, Cmember, cm) ;
/* error checks */
save = cm->error_handler ;
cm->error_handler = NULL ;
nc_new = CHOLMOD(collapse_septree) (n, nc, 0.1, 400,
CParent, Cmember, NULL) ;
OK (nc_new == EMPTY) ;
nc_new = CHOLMOD(collapse_septree) (n, nc, 0.1, 400,
NULL, Cmember, cm) ;
OK (nc_new == EMPTY) ;
nc_new = CHOLMOD(collapse_septree) (n, nc, 0.1, 400,
CParent, NULL, cm) ;
OK (nc_new == EMPTY) ;
nc_new = CHOLMOD(collapse_septree) (0, 1, 0.1, 400,
CParent, Cmember, cm) ;
OK (nc_new == EMPTY) ;
nc_new = CHOLMOD(collapse_septree) (1, 1, 0.1, 400,
CParent, Cmember, cm) ;
OK (nc_new == 1 || nc_new == EMPTY) ;
nc_new = CHOLMOD(collapse_septree) (Int_max, Int_max,
0.1, 400, CParent, Cmember, cm) ;
OK (nc_new == EMPTY) ;
cm->error_handler = save ;
}
CHOLMOD(free) (nrow, sizeof (Int), Cmember, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), CParent, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), Perm, cm) ;
cm->method [cm->current].nd_small = save1 ;
cm->method [cm->current].nd_oksep = save2 ;
}
#endif
/* ---------------------------------------------------------------------- */
/* dense to/from sparse conversions */
/* ---------------------------------------------------------------------- */
/* convert A to real, remove zero entries, and then convert to pattern */
if (MAX (nrow, ncol) < 1000)
{
C = CHOLMOD(copy_sparse) (A, cm) ;
CHOLMOD(sparse_xtype) (CHOLMOD_REAL, C, cm) ;
CHOLMOD(drop) (0., C, cm) ;
D = CHOLMOD(copy) (C, 0, 0, cm) ;
CHOLMOD(sort) (D, cm) ;
/* X = dense copy of C */
CHOLMOD(sparse_xtype) (CHOLMOD_PATTERN, C, cm) ;
X = CHOLMOD(sparse_to_dense) (C, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
/* change X to sparse pattern and then real/complex/zomplex, it should
* equal D */
for (xtype2 = CHOLMOD_REAL ; xtype2 <= CHOLMOD_ZOMPLEX ; xtype2++)
{
E = CHOLMOD(dense_to_sparse) (X, FALSE, cm) ;
ok = CHOLMOD(sparse_xtype) (xtype2, E, cm) ;
ok = CHOLMOD(sparse_xtype) (xtype2, D, cm) ;
if (xtype2 == CHOLMOD_REAL)
{
F = CHOLMOD(add) (E, D, one, minusone, TRUE, TRUE, cm) ;
r = CHOLMOD(norm_sparse) (F, 0, cm) ;
if (F != NULL)
{
OK (r == 0) ;
}
CHOLMOD(free_sparse) (&F, cm) ;
}
else
{
check_equality (E, D, xtype2) ;
}
CHOLMOD(free_sparse) (&E, cm) ;
}
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_dense) (&X, cm) ;
}
/* ---------------------------------------------------------------------- */
/* unsymmetric transpose */
/* ---------------------------------------------------------------------- */
len = ncol/2 ;
fset = prand (ncol) ; /* RAND */
CHOLMOD(print_perm) (P, nrow, nrow, "P", cm) ;
CHOLMOD(print_subset) (fset, ncol, ncol, "fset", cm) ;
if (isreal)
{
C = CHOLMOD(copy) (A, 0, 1, cm) ;
D = CHOLMOD(ptranspose) (C, 1, P, fset, len, cm) ;
E = CHOLMOD(transpose) (D, 1, cm) ;
F = CHOLMOD(transpose) (E, 1, cm) ;
G = CHOLMOD(add) (D, F, one, minusone, TRUE, FALSE, cm) ;
r = CHOLMOD(norm_sparse) (G, 0, cm) ;
if (G != NULL)
{
OK (r == 0) ;
}
CHOLMOD(drop) (0, G, cm) ;
r = CHOLMOD(norm_sparse) (G, 0, cm) ;
nz = CHOLMOD(nnz) (G, cm) ;
if (G != NULL)
{
OK (r == 0) ;
OK (nz == 0) ;
}
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&F, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
}
/* ---------------------------------------------------------------------- */
/* symmetric array transpose */
/* ---------------------------------------------------------------------- */
if (A->stype != 0)
{
/* C = A(p,p).' */
C = CHOLMOD(ptranspose) (A, 1, P, NULL, 0, cm) ;
/* D = C(pinv,pinv).' */
D = CHOLMOD(ptranspose) (C, 1, Pinv, NULL, 0, cm) ;
CHOLMOD(sort) (D, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
/* C = A, sorted */
C = CHOLMOD(copy_sparse) (A, cm) ;
CHOLMOD(sort) (C, cm) ;
/* C and D should be equal */
check_equality (C, D, xtype) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
/* C = A.' */
C = CHOLMOD(transpose) (A, 1, cm) ;
/* D = C.' */
D = CHOLMOD(transpose) (C, 1, cm) ;
CHOLMOD(sort) (D, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
/* C = A, sorted */
C = CHOLMOD(copy_sparse) (A, cm) ;
CHOLMOD(sort) (C, cm) ;
/* C and D should be equal */
check_equality (C, D, xtype) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&D, cm) ;
}
/* ---------------------------------------------------------------------- */
/* matrix multiply */
/* ---------------------------------------------------------------------- */
if (isreal)
{
/* this fails for a large arrowhead matrix, so turn off error hanlder */
save = cm->error_handler ;
cm->error_handler = NULL ;
AT = CHOLMOD(transpose) (A, 1, cm) ;
D = CHOLMOD(copy) (A, 0, 1, cm) ;
if (n > NLARGE) progress (1, '.') ;
C = CHOLMOD(aat) (D, NULL, 0, 1, cm) ;
if (n > NLARGE) progress (1, '.') ;
CHOLMOD(print_common) ("After A*A'", cm) ;
for (stype = -1 ; stype <= 1 ; stype++)
{
if (n > NLARGE) progress (1, '.') ;
E = CHOLMOD(ssmult) (A, AT, stype, TRUE, TRUE, cm) ;
if (n > NLARGE) progress (1, '.') ;
G = CHOLMOD(add) (C, E, one, minusone, TRUE, FALSE, cm) ;
if (n > NLARGE) progress (1, '.') ;
r = CHOLMOD(norm_sparse) (G, 0, cm) ;
if (G != NULL)
{
MAXERR (maxerr, r, anorm) ;
}
CHOLMOD(drop) (0, G, cm) ;
r = CHOLMOD(norm_sparse) (G, 0, cm) ;
if (G != NULL)
{
MAXERR (maxerr, r, anorm) ;
}
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
}
if (nrow == ncol)
{
/* E = pattern of A */
E = CHOLMOD(copy) (A, 0, 0, cm) ;
/* G = E*E */
if (n > NLARGE) progress (1, '.') ;
G = CHOLMOD(ssmult) (E, E, 0, FALSE, FALSE, cm) ;
if (n > NLARGE) progress (1, '.') ;
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&G, cm) ;
}
cm->error_handler = save ;
CHOLMOD(free_sparse) (&D, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&AT, cm) ;
}
/* ---------------------------------------------------------------------- */
/* free P, Q, and their inverses */
/* ---------------------------------------------------------------------- */
CHOLMOD(free) (ncol, sizeof (Int), fset, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), P, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), Pinv, cm) ;
if (A->stype == 0)
{
CHOLMOD(free) (ncol, sizeof (Int), Q, cm) ;
}
CHOLMOD(free) (ncol, sizeof (Int), Qinv, cm) ;
CHOLMOD(free) (nrow, sizeof (Int), Partition, cm) ;
progress (0, '.') ;
return (maxerr) ;
}
syntax highlighted by Code2HTML, v. 0.9.1