/* ========================================================================== */
/* === Tcov/cm ============================================================== */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* 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
* -------------------------------------------------------------------------- */
/* A program for exhaustive statement-coverage for CHOLMOD, AMD, COLAMD, and
* CCOLAMD. It tests every line of code in all three packages.
*
* For a complete test, all CHOLMOD modules, AMD, COLAMD, CCOLAMD, METIS,
* LAPACK, and the BLAS are required. A partial test can be performed without
* the Supernodal and/or Partition modules. METIS is not required if the
* Partition module is not installed. LAPACK and the BLAS are not required
* if the Supernodal module is not installed.
*
* Usage:
*
* cm < input > output
*
* where "input" contains a sparse matrix in triplet form. The first line of
* the file contains four or five integers:
*
* nrow ncol nnz stype complex
*
* where the matrix is nrow-by-ncol. nnz is the number of (i,j,aij) triplets
* in the rest of the file, one triplet per line. stype is -1 (symmetric
* with entries in lower triangular part provided), 0 (unsymmetric), or 1
* (symmetric upper). If the 5th entry is missing, or 0, then the matrix is
* real; if 1 the matrix is complex, if 2 the matrix is zomplex. Each
* subsequent line contains
*
* i j aij
*
* for the row index, column index, and value of A(i,j). Duplicate entries
* are summed. If stype is 2 or 3, the rest of the file is ignored, and a
* special test matrix is constructed (2: arrowhead, 3: tridiagonal plus a
* dense row). Test matrices are located in the Matrix/ subdirectory.
*
* For complex matrices, each line consists of
*
* i j xij zij
*
* where xij is the real part of A(i,j) and zij is the imaginary part.
*
* cm takes one optional parameter. If present (it does not matter what the
* argument is, actually) then extension memory-failure tests are performed.
*/
#include "cm.h"
/* ========================================================================== */
/* === global variables ===================================================== */
/* ========================================================================== */
double zero [2], one [2], minusone [2] ;
cholmod_common Common, *cm ;
cholmod_dense *M1 ;
Int dot = 0 ;
double Zero [2] ;
/* ========================================================================== */
/* === my_rand ============================================================== */
/* ========================================================================== */
/* The POSIX example of rand, duplicated here so that the same sequence will
* be generated on different machines. */
static unsigned long next = 1 ;
/* RAND_MAX assumed to be 32767 */
Int my_rand (void)
{
next = next * 1103515245 + 12345 ;
return ((unsigned)(next/65536) % /* 32768 */ (MY_RAND_MAX + 1)) ;
}
void my_srand (unsigned seed)
{
next = seed ;
}
unsigned long my_seed (void)
{
return (next) ;
}
/* ========================================================================== */
/* === progress ============================================================= */
/* ========================================================================== */
/* print a "." on stderr to indicate progress */
#define LINEWIDTH 70
#define COUNT 100
void progress (Int force, char s)
{
dot++ ;
if (force)
{
dot += (COUNT - (dot % COUNT)) ;
}
if (dot % COUNT == 0)
{
fprintf (stderr, "%c", s) ;
}
if (dot % (COUNT*LINEWIDTH) == 0)
{
fprintf (stderr, "\n") ;
}
fflush (stdout) ;
fflush (stderr) ;
}
/* ========================================================================== */
/* === my_handler =========================================================== */
/* ========================================================================== */
/* An error occurred that should not have occurred */
void my_handler (int status, char *file, int line, char *msg)
{
printf ("Error handler: file %s line %d status %d: %s\n",
file, line, status, msg) ;
if (status < CHOLMOD_OK || status > CHOLMOD_DSMALL)
{
fprintf (stderr, "\n\n************************************************"
"********************************\n"
"*** Test failure: file: %s line: %d\n"
"*** status: %d message: %s\n"
"***********************************************************"
"*********************\n", file, line, status, msg);
fflush (stderr) ;
fflush (stdout) ;
abort ( ) ;
}
}
/* ========================================================================== */
/* === Assert =============================================================== */
/* ========================================================================== */
void Assert (int truth, char *file, int line)
{
if (!truth)
{
my_handler (-1, file, line, "") ;
}
}
/* ========================================================================== */
/* === nrand ================================================================ */
/* ========================================================================== */
/* return a random Int between 0 and n-1 */
Int nrand (Int n)
{
return ((n <= 0) ? (0) : (rand ( ) % n)) ;
}
/* ========================================================================== */
/* === xrand ================================================================ */
/* ========================================================================== */
/* return a random double between 0 and x */
double xrand (double range)
{
return ((range * (double) (my_rand ( ))) / MY_RAND_MAX) ;
}
/* ========================================================================== */
/* === prand ================================================================ */
/* ========================================================================== */
/* allocate and construct a random permutation of 0:n-1 */
Int *prand (Int n)
{
Int *P ;
Int t, j, k ;
P = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
if (P == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot create random perm") ;
return (NULL) ;
}
for (k = 0 ; k < n ; k++)
{
P [k] = k ;
}
for (k = 0 ; k < n-1 ; k++)
{
j = k + nrand (n-k) ;
t = P [j] ;
P [j] = P [k] ;
P [k] = t ;
}
CHOLMOD(print_perm) (P, n, n, "Prandom", cm) ;
return (P) ;
}
/* ========================================================================== */
/* === rand_set ============================================================= */
/* ========================================================================== */
/* allocate and construct a random set of 0:n-1, possibly with duplicates */
Int *rand_set (Int len, Int n)
{
Int *cset ;
Int k ;
cset = CHOLMOD(malloc) (len, sizeof (Int), cm) ;
if (cset == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot create random set") ;
return (NULL) ;
}
for (k = 0 ; k < len ; k++)
{
cset [k] = nrand (n) ;
}
CHOLMOD(print_subset) (cset, len, n, "random cset", cm) ;
return (cset) ;
}
/* ========================================================================== */
/* === read_triplet ========================================================= */
/* ========================================================================== */
/* Read a triplet matrix from a file. */
#define MAXLINE 1024
cholmod_triplet *read_triplet
(
FILE *f
)
{
cholmod_triplet *T ;
double *Tx, *Tz ;
long long x1, x2, x3, x4, x5 ;
Int *Ti, *Tj ;
Int n, j, k, nrow, ncol, nz, stype, arrowhead, tridiag_plus_denserow,
xtype, is_complex ;
char s [MAXLINE] ;
/* ---------------------------------------------------------------------- */
/* read in a triplet matrix from a file */
/* ---------------------------------------------------------------------- */
dot = 0 ;
xtype = 0 ;
if (fgets (s, MAXLINE, f) == NULL)
{
return (NULL) ;
}
x1 = 0 ;
x2 = 0 ;
x3 = 0 ;
x4 = 0 ;
x5 = 0 ;
k = sscanf (s, "%lld %lld %lld %lld %lld\n", &x1, &x2, &x3, &x4, &x5) ;
nrow = x1 ;
ncol = x2 ;
nz = x3 ;
stype = x4 ;
xtype = x5 ;
xtype++ ;
is_complex = (xtype != CHOLMOD_REAL) ;
printf ("read_triplet: nrow "ID" ncol "ID" nz "ID" stype "ID" xtype "ID"\n",
nrow, ncol, nz, stype, xtype) ;
arrowhead = FALSE ;
tridiag_plus_denserow = FALSE ;
n = MAX (nrow, ncol) ;
if (stype == 2)
{
/* ignore nz and the rest of the file, and create an arrowhead matrix */
arrowhead = TRUE ;
nz = nrow + ncol + 1 ;
stype = (nrow == ncol) ? (1) : (0) ;
}
else if (stype == 3)
{
tridiag_plus_denserow = TRUE ;
nrow = n ;
ncol = n ;
nz = 4*n + 4 ;
stype = 0 ;
}
T = CHOLMOD(allocate_triplet) (nrow, ncol, nz, stype,
is_complex ? CHOLMOD_ZOMPLEX : CHOLMOD_REAL, cm) ;
if (T == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot create triplet matrix") ;
return (NULL) ;
}
Ti = T->i ;
Tj = T->j ;
Tx = T->x ;
Tz = T->z ;
if (arrowhead)
{
for (k = 0 ; k < MIN (nrow,ncol) ; k++)
{
Ti [k] = k ;
Tj [k] = k ;
Tx [k] = nrow + xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = nrow + xrand (1) ; /* RAND */
}
}
for (j = 0 ; j < ncol ; j++)
{
Ti [k] = 0 ;
Tj [k] = j ;
Tx [k] = - xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = - xrand (1) ; /* RAND */
}
k++ ;
}
T->nnz = k ;
}
else if (tridiag_plus_denserow)
{
/* dense row, except for the last column */
for (k = 0 ; k < n-1 ; k++)
{
Ti [k] = 0 ;
Tj [k] = k ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
}
/* diagonal */
for (j = 0 ; j < n ; j++)
{
Ti [k] = j ;
Tj [k] = j ;
Tx [k] = nrow + xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = nrow + xrand (1) ; /* RAND */
}
k++ ;
}
/* superdiagonal */
for (j = 1 ; j < n ; j++)
{
Ti [k] = j-1 ;
Tj [k] = j ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
k++ ;
}
/* subdiagonal */
for (j = 0 ; j < n-1 ; j++)
{
Ti [k] = j+1 ;
Tj [k] = j ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
k++ ;
}
/* a few extra terms in the last column */
Ti [k] = MAX (0, n-3) ;
Tj [k] = n-1 ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
k++ ;
Ti [k] = MAX (0, n-4) ;
Tj [k] = n-1 ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
k++ ;
Ti [k] = MAX (0, n-5) ;
Tj [k] = n-1 ;
Tx [k] = xrand (1) ; /* RAND */
if (is_complex)
{
Tz [k] = xrand (1) ; /* RAND */
}
k++ ;
T->nnz = k ;
}
else
{
if (is_complex)
{
for (k = 0 ; k < nz ; k++)
{
if (fscanf (f,""ID" "ID" %lg %lg\n", Ti+k, Tj+k, Tx+k, Tz+k)
== EOF)
{
ERROR (CHOLMOD_INVALID, "Error reading triplet matrix\n") ;
}
}
}
else
{
for (k = 0 ; k < nz ; k++)
{
if (fscanf (f, ""ID" "ID" %lg\n", Ti+k, Tj+k, Tx+k) == EOF)
{
ERROR (CHOLMOD_INVALID, "Error reading triplet matrix\n") ;
}
}
}
T->nnz = nz ;
}
CHOLMOD(triplet_xtype) (xtype, T, cm) ;
/* ---------------------------------------------------------------------- */
/* print the triplet matrix */
/* ---------------------------------------------------------------------- */
cm->print = 4 ;
CHOLMOD(print_triplet) (T, "T input", cm) ;
cm->print = 1 ;
fprintf (stderr, "Test matrix: "ID"-by-"ID" with "ID" entries, stype: "ID
"\n",
(Int) T->nrow, (Int) T->ncol, (Int) T->nnz, (Int) T->stype) ;
printf ("\n\n======================================================\n"
"Test matrix: "ID"-by-"ID" with "ID" entries, stype: "ID"\n",
(Int) T->nrow, (Int) T->ncol, (Int) T->nnz, (Int) T->stype) ;
if (MAX (nrow, ncol) > NLARGE)
{
fprintf (stderr, "Please wait, this will take a while ...") ;
dot = 39*LINEWIDTH ;
}
return (T) ;
}
/* ========================================================================== */
/* === zeros ================================================================ */
/* ========================================================================== */
/* Same as cholmod_zeros or cholmod_l_zeros, except it allows for a leading
* dimension that is different than nrow */
cholmod_dense *zeros (Int nrow, Int ncol, Int d, Int xtype)
{
cholmod_dense *X ;
double *Xx, *Xz ;
Int i, nz ;
X = CHOLMOD(allocate_dense) (nrow, ncol, d, xtype, cm) ;
if (X == NULL)
{
return (NULL) ;
}
Xx = X->x ;
Xz = X->z ;
nz = MAX (1, X->nzmax) ;
switch (X->xtype)
{
case CHOLMOD_REAL:
for (i = 0 ; i < nz ; i++)
{
Xx [i] = 0 ;
}
break ;
case CHOLMOD_COMPLEX:
for (i = 0 ; i < 2*nz ; i++)
{
Xx [i] = 0 ;
}
break ;
case CHOLMOD_ZOMPLEX:
for (i = 0 ; i < nz ; i++)
{
Xx [i] = 0 ;
}
for (i = 0 ; i < nz ; i++)
{
Xz [i] = 0 ;
}
break ;
}
return (X) ;
}
/* ========================================================================== */
/* === xtrue ================================================================ */
/* ========================================================================== */
/* Allocate and construct a dense matrix, X(i,j) = i+j*d+1 */
cholmod_dense *xtrue (Int nrow, Int ncol, Int d, Int xtype)
{
double *x, *z ;
cholmod_dense *X ;
Int j, i ;
X = zeros (nrow, ncol, d, xtype) ;
if (X == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot create dense matrix") ;
return (NULL) ;
}
x = X->x ;
z = X->z ;
if (xtype == CHOLMOD_REAL)
{
for (j = 0 ; j < ncol ; j++)
{
for (i = 0 ; i < nrow ; i++)
{
x [i+j*d] = i+j*d + 1 ;
}
}
}
else if (xtype == CHOLMOD_COMPLEX)
{
for (j = 0 ; j < ncol ; j++)
{
for (i = 0 ; i < nrow ; i++)
{
x [2*(i+j*d) ] = i+j*d + 1 ;
x [2*(i+j*d)+1] = ((double) (j+i*d + 1))/10 ;
}
}
}
else
{
for (j = 0 ; j < ncol ; j++)
{
for (i = 0 ; i < nrow ; i++)
{
x [i+j*d] = i+j*d + 1 ;
z [i+j*d] = ((double) (j+i*d + 1))/10 ;
}
}
}
return (X) ;
}
/* ========================================================================== */
/* === rhs ================================================================== */
/* ========================================================================== */
/* Create a right-hand-side, b = A*x, where x is a known solution */
cholmod_dense *rhs (cholmod_sparse *A, Int nrhs, Int d)
{
Int n ;
cholmod_dense *W, *Z, *B ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute rhs") ;
return (NULL) ;
}
n = A->nrow ;
/* B = zeros (n,rhs) but with leading dimension d */
B = zeros (n, nrhs, d, A->xtype) ;
/* ---------------------------------------------------------------------- */
/* create a known solution */
/* ---------------------------------------------------------------------- */
Z = xtrue (n, nrhs, d, A->xtype) ;
/* ---------------------------------------------------------------------- */
/* compute B = A*Z or A*A'*Z */
/* ---------------------------------------------------------------------- */
if (A->stype == 0)
{
/* W = A'*Z */
W = CHOLMOD(zeros) (A->ncol, nrhs, A->xtype, cm) ;
CHOLMOD(sdmult) (A, TRUE, one, zero, Z, W, cm) ;
/* B = A*W */
CHOLMOD(sdmult) (A, FALSE, one, zero, W, B, cm) ;
CHOLMOD(free_dense) (&W, cm) ;
}
else
{
/* B = A*Z */
CHOLMOD(sdmult) (A, FALSE, one, zero, Z, B, cm) ;
}
CHOLMOD(free_dense) (&Z, cm) ;
return (B) ;
}
/* ========================================================================== */
/* === resid ================================================================ */
/* ========================================================================== */
/* compute r = norm (A*x-b)/norm(b) or r = norm (A*A'*x-b)/norm(b) */
double resid (cholmod_sparse *A, cholmod_dense *X, cholmod_dense *B)
{
double r, bnorm ;
cholmod_dense *R, *X2, *B2 ;
cholmod_sparse *C, *A2 ;
Int d, n, nrhs, xtype ;
if (A == NULL || X == NULL || B == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
return (1) ;
}
cm->status = CHOLMOD_OK ;
n = B->nrow ;
/* ---------------------------------------------------------------------- */
/* convert all inputs to an identical xtype */
/* ---------------------------------------------------------------------- */
xtype = MAX (A->xtype, X->xtype) ;
xtype = MAX (xtype, B->xtype) ;
A2 = NULL ;
B2 = NULL ;
X2 = NULL ;
if (A->xtype != xtype)
{
A2 = CHOLMOD(copy_sparse) (A, cm) ;
CHOLMOD(sparse_xtype) (xtype, A2, cm) ;
A = A2 ;
}
if (X->xtype != xtype)
{
X2 = CHOLMOD(copy_dense) (X, cm) ;
CHOLMOD(dense_xtype) (xtype, X2, cm) ;
X = X2 ;
}
if (B->xtype != xtype)
{
B2 = CHOLMOD(copy_dense) (B, cm) ;
CHOLMOD(dense_xtype) (xtype, B2, cm) ;
B = B2 ;
}
if (cm->status < CHOLMOD_OK)
{
ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
CHOLMOD(free_sparse) (&A2, cm) ;
CHOLMOD(free_dense) (&B2, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
return (1) ;
}
/* ---------------------------------------------------------------------- */
/* get the right-hand-side, B, and its norm */
/* ---------------------------------------------------------------------- */
nrhs = B->ncol ;
d = B->d ;
if (nrhs == 1)
{
/* inf-norm, 1-norm, or 2-norm (random choice) */
bnorm = CHOLMOD(norm_dense) (B, nrand (2), cm) ;
}
else
{
/* inf-norm or 1-norm (random choice) */
bnorm = CHOLMOD(norm_dense) (B, nrand (1), cm) ;
}
/* ---------------------------------------------------------------------- */
/* compute the residual */
/* ---------------------------------------------------------------------- */
if (A->stype == 0)
{
if (n < 10 && A->xtype == CHOLMOD_REAL)
{
/* test cholmod_aat, C = A*A' */
C = CHOLMOD(aat) (A, NULL, 0, 1, cm) ;
/* R = B */
R = CHOLMOD(copy_dense) (B, cm) ;
/* R = C*X - R */
CHOLMOD(sdmult) (C, FALSE, one, minusone, X, R, cm) ;
CHOLMOD(free_sparse) (&C, cm) ;
}
else
{
/* W = A'*X */
cholmod_dense *W ;
W = CHOLMOD(zeros) (A->ncol, nrhs, A->xtype, cm) ;
CHOLMOD(sdmult) (A, TRUE, one, zero, X, W, cm) ;
/* R = B */
R = CHOLMOD(copy_dense) (B, cm) ;
/* R = A*W - R */
CHOLMOD(sdmult) (A, FALSE, one, minusone, W, R, cm) ;
CHOLMOD(free_dense) (&W, cm) ;
}
}
else
{
/* R = B */
R = CHOLMOD(copy_dense) (B, cm) ;
/* R = A*X - R */
CHOLMOD(sdmult) (A, FALSE, one, minusone, X, R, cm) ;
}
/* ---------------------------------------------------------------------- */
/* r = norm (R) / norm (B) */
/* ---------------------------------------------------------------------- */
r = CHOLMOD(norm_dense) (R, 1, cm) ;
CHOLMOD(free_dense) (&R, cm) ;
CHOLMOD(free_sparse) (&A2, cm) ;
CHOLMOD(free_dense) (&B2, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
if (bnorm > 0)
{
r /= bnorm ;
}
return (r) ;
}
/* ========================================================================== */
/* === resid_sparse ========================================================= */
/* ========================================================================== */
/* compute r = norm (A*x-b)/norm(b) or r = norm (A*A'*x-b)/norm(b) */
double resid_sparse (cholmod_sparse *A, cholmod_sparse *X, cholmod_sparse *B)
{
double r, bnorm ;
cholmod_sparse *R, *W, *AT, *W2 ;
cholmod_dense *X2, *B2 ;
Int n, nrhs, xtype ;
if (A == NULL || X == NULL || B == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute resid") ;
return (1) ;
}
/* ---------------------------------------------------------------------- */
/* compute the residual */
/* ---------------------------------------------------------------------- */
xtype = MAX (A->xtype, X->xtype) ;
xtype = MAX (xtype, B->xtype) ;
if (xtype > CHOLMOD_REAL)
{
/* ------------------------------------------------------------------ */
/* convert X and B to dense if any is complex or zomplex */
/* ------------------------------------------------------------------ */
X2 = CHOLMOD(sparse_to_dense) (X, cm) ;
B2 = CHOLMOD(sparse_to_dense) (B, cm) ;
r = resid (A, X2, B2) ;
CHOLMOD(free_dense) (&X2, cm) ;
CHOLMOD(free_dense) (&B2, cm) ;
}
else
{
/* ------------------------------------------------------------------ */
/* all inputs are real */
/* ------------------------------------------------------------------ */
n = B->nrow ;
nrhs = B->ncol ;
/* inf-norm or 1-norm (random choice) */
bnorm = CHOLMOD(norm_sparse) (B, nrand (1), cm) ;
if (A->stype == 0)
{
/* W = A'*X */
AT = CHOLMOD(transpose) (A, 1, cm) ;
W = CHOLMOD(ssmult) (AT, X, 0, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&AT, cm) ;
/* W2 = A*W */
W2 = CHOLMOD(ssmult) (A, W, 0, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&W, cm) ;
/* R = W2 - B */
R = CHOLMOD(add) (W2, B, one, minusone, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&W2, cm) ;
}
else
{
/* W = A*X */
W = CHOLMOD(ssmult) (A, X, 0, TRUE, FALSE, cm) ;
/* R = W - B */
R = CHOLMOD(add) (W, B, one, minusone, TRUE, FALSE, cm) ;
CHOLMOD(free_sparse) (&W, cm) ;
}
r = CHOLMOD(norm_sparse) (R, 1, cm) ;
CHOLMOD(free_sparse) (&R, cm) ;
if (bnorm > 0)
{
r /= bnorm ;
}
}
return (r) ;
}
/* ========================================================================== */
/* === resid3 =============================================================== */
/* ========================================================================== */
/* r = norm (A1*A2*A3*x - b) / norm (b) */
double resid3 (cholmod_sparse *A1, cholmod_sparse *A2, cholmod_sparse *A3,
cholmod_dense *X, cholmod_dense *B)
{
double r, bnorm ;
cholmod_dense *R, *W1, *W2, *X2, *B2 ;
cholmod_sparse *C1, *C2, *C3 ;
Int n, nrhs, d, xtype ;
if (A1 == NULL || X == NULL || B == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute resid3") ;
return (1) ;
}
cm->status = CHOLMOD_OK ;
n = B->nrow ;
/* ---------------------------------------------------------------------- */
/* convert all inputs to an identical xtype */
/* ---------------------------------------------------------------------- */
xtype = MAX (A1->xtype, X->xtype) ;
xtype = MAX (xtype, B->xtype) ;
if (A2 != NULL)
{
xtype = MAX (xtype, A2->xtype) ;
}
if (A3 != NULL)
{
xtype = MAX (xtype, A3->xtype) ;
}
C1 = NULL ;
C2 = NULL ;
C3 = NULL ;
B2 = NULL ;
X2 = NULL ;
if (A1->xtype != xtype)
{
C1 = CHOLMOD(copy_sparse) (A1, cm) ;
CHOLMOD(sparse_xtype) (xtype, C1, cm) ;
A1 = C1 ;
}
if (A2 != NULL && A2->xtype != xtype)
{
C2 = CHOLMOD(copy_sparse) (A2, cm) ;
CHOLMOD(sparse_xtype) (xtype, C2, cm) ;
A2 = C2 ;
}
if (A3 != NULL && A3->xtype != xtype)
{
C3 = CHOLMOD(copy_sparse) (A3, cm) ;
CHOLMOD(sparse_xtype) (xtype, C3, cm) ;
A3 = C3 ;
}
if (X->xtype != xtype)
{
X2 = CHOLMOD(copy_dense) (X, cm) ;
CHOLMOD(dense_xtype) (xtype, X2, cm) ;
X = X2 ;
}
if (B->xtype != xtype)
{
B2 = CHOLMOD(copy_dense) (B, cm) ;
CHOLMOD(dense_xtype) (xtype, B2, cm) ;
B = B2 ;
}
if (cm->status < CHOLMOD_OK)
{
ERROR (CHOLMOD_INVALID, "cannot compute resid3") ;
CHOLMOD(free_sparse) (&C1, cm) ;
CHOLMOD(free_sparse) (&C2, cm) ;
CHOLMOD(free_sparse) (&C3, cm) ;
CHOLMOD(free_dense) (&B2, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
return (1) ;
}
/* ---------------------------------------------------------------------- */
/* get B and its norm */
/* ---------------------------------------------------------------------- */
nrhs = B->ncol ;
d = B->d ;
bnorm = CHOLMOD(norm_dense) (B, 1, cm) ;
/* ---------------------------------------------------------------------- */
/* compute the residual */
/* ---------------------------------------------------------------------- */
if (A3 != NULL)
{
/* W1 = A3*X */
W1 = CHOLMOD(zeros) (n, nrhs, xtype, cm) ;
CHOLMOD(sdmult) (A3, FALSE, one, zero, X, W1, cm) ;
}
else
{
W1 = X ;
}
if (A2 != NULL)
{
/* W2 = A2*W1 */
W2 = CHOLMOD(eye) (n, nrhs, xtype, cm) ;
CHOLMOD(sdmult) (A2, FALSE, one, zero, W1, W2, cm) ;
}
else
{
W2 = W1 ;
}
/* R = B */
R = CHOLMOD(copy_dense) (B, cm) ;
/* R = A1*W2 - R */
CHOLMOD(sdmult) (A1, FALSE, one, minusone, W2, R, cm) ;
/* ---------------------------------------------------------------------- */
/* r = norm (R) / norm (B) */
/* ---------------------------------------------------------------------- */
r = CHOLMOD(norm_dense) (R, 1, cm) ;
CHOLMOD(free_dense) (&R, cm) ;
CHOLMOD(free_sparse) (&C1, cm) ;
CHOLMOD(free_sparse) (&C2, cm) ;
CHOLMOD(free_sparse) (&C3, cm) ;
CHOLMOD(free_dense) (&B2, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
if (A3 != NULL)
{
CHOLMOD(free_dense) (&W1, cm) ;
}
if (A2 != NULL)
{
CHOLMOD(free_dense) (&W2, cm) ;
}
if (bnorm > 0)
{
r /= bnorm ;
}
return (r) ;
}
/* ========================================================================== */
/* === pnorm ================================================================ */
/* ========================================================================== */
/* r = norm (x-Pb) or r = norm(x-P'b). This is lengthy because CHOLMOD does
* not provide any operations on dense matrices, and because it used to test
* the sparse-to-dense conversion routine. Multiple methods are used to compute
* the same thing.
*/
double pnorm (cholmod_dense *X, Int *P, cholmod_dense *B, Int inv)
{
cholmod_dense *R, *X2, *B2 ;
cholmod_factor *L ;
double *xx, *xz, *bx, *bz, *rx, *rz ;
Int *Pinv, *Perm ;
double rnorm, r ;
Int i, j, k, n, nrhs, xtype, ok, save, lxtype ;
if (X == NULL || P == NULL || B == NULL)
{
ERROR (CHOLMOD_INVALID, "cannot compute pnorm") ;
return (1) ;
}
save = cm->prefer_zomplex ;
n = X->nrow ;
nrhs = X->ncol ;
rnorm = 0 ;
Pinv = CHOLMOD(malloc) (n, sizeof (Int), cm) ;
if (Pinv != NULL)
{
for (k = 0 ; k < n ; k++)
{
Pinv [P [k]] = k ;
}
}
xtype = MAX (X->xtype, B->xtype) ;
R = CHOLMOD(zeros) (n, nrhs, CHOLMOD_ZOMPLEX, cm) ;
B2 = CHOLMOD(copy_dense) (B, cm) ;
ok = R != NULL && B2 != NULL ;
ok = ok && CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, B2, cm) ;
for (lxtype = CHOLMOD_REAL ; ok && lxtype <= CHOLMOD_ZOMPLEX ; lxtype++)
{
/* create a fake factor object */
L = CHOLMOD(allocate_factor) (n, cm) ;
CHOLMOD(change_factor) (lxtype, TRUE, FALSE, TRUE, TRUE, L, cm) ;
ok = ok && (L != NULL && L->Perm != NULL && Pinv != NULL) ;
if (ok)
{
L->ordering = CHOLMOD_GIVEN ;
Perm = L->Perm ;
for (k = 0 ; k < n ; k++)
{
Perm [k] = Pinv [k] ;
}
}
for (k = 0 ; k <= 1 ; k++)
{
/* solve the inverse permutation system, X2 = P*X or X2 = P'*X */
cm->prefer_zomplex = k ;
X2 = CHOLMOD(solve) (inv ? CHOLMOD_Pt : CHOLMOD_P, L, X, cm) ;
ok = ok && CHOLMOD(dense_xtype) (CHOLMOD_ZOMPLEX, X2, cm) ;
if (ok && X2 != NULL)
{
rx = R->x ;
rz = R->z ;
xx = X2->x ;
xz = X2->z ;
bx = B2->x ;
bz = B2->z ;
for (j = 0 ; j < nrhs ; j++)
{
for (i = 0 ; i < n ; i++)
{
rx [i+j*n] = xx [i+j*n] - bx [i+j*n] ;
rz [i+j*n] = xz [i+j*n] - bz [i+j*n] ;
}
}
}
r = CHOLMOD(norm_dense) (R, 0, cm) ;
rnorm = MAX (r, rnorm) ;
CHOLMOD(free_dense) (&X2, cm) ;
}
CHOLMOD(free_factor) (&L, cm) ;
}
CHOLMOD(free_dense) (&B2, cm) ;
CHOLMOD(free_dense) (&R, cm) ;
if (xtype == CHOLMOD_REAL)
{
/* X and B are both real */
cholmod_sparse *Bs, *Pb ;
Bs = CHOLMOD(dense_to_sparse) (B, TRUE, cm) ;
Pb = CHOLMOD(submatrix) (Bs, inv ? Pinv : P, n, NULL, -1, TRUE, TRUE,cm);
X2 = CHOLMOD(sparse_to_dense) (Pb, cm) ;
R = CHOLMOD(zeros) (n, nrhs, CHOLMOD_REAL, cm) ;
if (R != NULL && X != NULL && X2 != NULL)
{
rx = R->x ;
xx = X->x ;
bx = X2->x ;
for (j = 0 ; j < nrhs ; j++)
{
for (i = 0 ; i < n ; i++)
{
rx [i+j*n] = xx [i+j*n] - bx [i+j*n] ;
}
}
}
CHOLMOD(free_sparse) (&Bs, cm) ;
CHOLMOD(free_sparse) (&Pb, cm) ;
CHOLMOD(free_dense) (&X2, cm) ;
r = CHOLMOD(norm_dense) (R, 1, cm) ;
rnorm = MAX (rnorm, r) ;
CHOLMOD(free_dense) (&R, cm) ;
}
CHOLMOD(free) (n, sizeof (Int), Pinv, cm) ;
cm->prefer_zomplex = save ;
return (rnorm) ;
}
/* ========================================================================== */
/* === prune_row ============================================================ */
/* ========================================================================== */
/* Set row k and column k of a packed matrix A to zero. Set A(k,k) to 1
* if space is available. */
void prune_row (cholmod_sparse *A, Int k)
{
double *Ax ;
Int *Ap, *Ai ;
Int ncol, p, i, j, nz ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "nothing to prune") ;
return ;
}
Ap = A->p ;
Ai = A->i ;
Ax = A->x ;
nz = 0 ;
ncol = A->ncol ;
for (j = 0 ; j < ncol ; j++)
{
p = Ap [j] ;
Ap [j] = nz ;
if (j == k && nz < Ap [j+1])
{
Ai [nz] = k ;
Ax [nz] = 1 ;
nz++ ;
}
else
{
for ( ; p < Ap [j+1] ; p++)
{
i = Ai [p] ;
if (i != k)
{
Ai [nz] = i ;
Ax [nz] = Ax [p] ;
nz++ ;
}
}
}
}
Ap [ncol] = nz ;
}
/* ========================================================================== */
/* === do_matrix =========================================================== */
/* ========================================================================== */
double do_matrix (cholmod_sparse *A)
{
double err, maxerr = 0 ;
Int print, precise, maxprint, minprint, nmethods ;
if (A == NULL)
{
ERROR (CHOLMOD_INVALID, "no matrix") ;
return (1) ;
}
/* ---------------------------------------------------------------------- */
/* determine print level, based on matrix size */
/* ---------------------------------------------------------------------- */
if (A->nrow <= 4)
{
minprint = 5 ;
maxprint = 5 ;
}
else if (A->nrow <= 8)
{
minprint = 4 ;
maxprint = 4 ;
}
else
{
minprint = 1 ;
maxprint = 1 ;
}
/* ---------------------------------------------------------------------- */
/* for all print levels and precisions, do: */
/* ---------------------------------------------------------------------- */
for (print = minprint ; print <= maxprint ; print++)
{
for (precise = 0 ; precise <= (print >= 4) ? 1 : 0 ; precise++)
{
Int save1, save2 ;
maxerr = 0 ;
my_srand (42) ; /* RAND reset */
cm->print = print ;
cm->precise = precise ;
printf ("\n----------Print level %d precise: %d\n",
cm->print, cm->precise) ;
save1 = cm->final_asis ;
save2 = cm->final_super ;
cm->final_asis = FALSE ;
cm->final_super = TRUE ;
OK (CHOLMOD(print_common) ("cm", cm)) ;
cm->final_asis = save1 ;
cm->final_super = save2 ;
/* -------------------------------------------------------------- */
/* test various matrix operations */
/* -------------------------------------------------------------- */
err = test_ops (A) ; /* RAND */
MAXERR (maxerr, err, 1) ;
/* -------------------------------------------------------------- */
/* solve the augmented system */
/* -------------------------------------------------------------- */
err = aug (A) ; /* no random number use */
MAXERR (maxerr, err, 1) ;
/* -------------------------------------------------------------- */
/* solve using different methods */
/* -------------------------------------------------------------- */
printf ("test_solver (1)\n") ;
cm->nmethods = 9 ;
cm->final_asis = TRUE ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (2)\n") ;
cm->final_asis = TRUE ;
cm->method [0].ordering = CHOLMOD_NATURAL ;
for (nmethods = 0 ; nmethods < 7 ; nmethods++)
{
cm->nmethods = nmethods ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
}
printf ("test_solver (3)\n") ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_NESDIS ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (3b)\n") ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_NESDIS ;
cm->method [0].nd_camd = 2 ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (3c)\n") ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_NATURAL ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (4)\n") ;
cm->nmethods = 1 ;
cm->method[0].ordering = CHOLMOD_METIS ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (5)\n") ;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_AMD ;
CHOLMOD(free_work) (cm) ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
printf ("test_solver (6)\n") ;
cm->nmethods = 1 ;
cm->method[0].ordering = CHOLMOD_COLAMD ;
err = test_solver (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
/* -------------------------------------------------------------- */
/* restore default control parameters */
/* -------------------------------------------------------------- */
OK (CHOLMOD(print_common) ("cm", cm)) ;
CHOLMOD(defaults) (cm) ;
}
}
printf ("do_matrix max error %.1g\n", maxerr) ;
return (maxerr) ;
}
/* ========================================================================== */
/* === main ================================================================= */
/* ========================================================================== */
/* Usage:
* cm < matrix do not perform intensive memory-failure tests
* cm -m < matrix do perform memory tests
* cm -s < matrix matrix is singular, nan error expected
*
* (The memory tests are performed if any argument is given to cm).
*/
int main (int argc, char **argv)
{
cholmod_triplet *T ;
cholmod_sparse *A, *C, *AT ;
char *s ;
double err = 0, maxerr = 0 ;
Int n = 0, nmin = 0, nrow = 0, ncol = 0, save ;
int singular, do_memory, i, do_nantests ;
double v = CHOLMOD_VERSION ;
printf ("Testing CHOLMOD (%g)\n", v) ;
printf ("%s: argc: %d\n", argv [0], argc) ;
my_srand (42) ; /* RAND */
fflush (stdout) ;
singular = FALSE ;
do_memory = FALSE ;
do_nantests = FALSE ;
for (i = 1 ; i < argc ; i++)
{
s = argv [i] ;
if (s [0] == '-' && s [1] == 'm') do_memory = TRUE ;
if (s [0] == '-' && s [1] == 's') singular = TRUE ;
if (s [0] == '-' && s [1] == 'n') do_nantests = TRUE ;
}
printf ("do_memory: %d singular: %d\n", do_memory, singular) ;
/* ---------------------------------------------------------------------- */
/* initialize CHOLMOD */
/* ---------------------------------------------------------------------- */
cm = &Common ;
OK (CHOLMOD(start) (cm)) ;
/* ---------------------------------------------------------------------- */
/* test all methods with NULL common */
/* ---------------------------------------------------------------------- */
/* no user error handler, since lots of errors will be raised here */
cm->error_handler = NULL ;
null_test (NULL) ;
save = cm->itype ;
cm->itype = -999 ;
null_test (cm) ;
cm->itype = save ;
null_test2 ( ) ;
CHOLMOD(finish) (cm) ;
OK (cm->malloc_count == 0) ;
OK (CHOLMOD(start) (cm)) ;
/* ---------------------------------------------------------------------- */
/* create basic scalars */
/* ---------------------------------------------------------------------- */
Zero [0] = 0 ;
Zero [1] = 0 ;
zero [0] = 0 ;
zero [1] = 0 ;
one [0] = 1 ;
one [1] = 0 ;
minusone [0] = -1 ;
minusone [1] = 0 ;
M1 = CHOLMOD(ones) (1, 1, CHOLMOD_REAL, cm) ;
if (M1 != NULL)
{
((double *) (M1->x)) [0] = -1 ;
}
/* ---------------------------------------------------------------------- */
/* read in a triplet matrix and use it to test CHOLMOD */
/* ---------------------------------------------------------------------- */
for ( ; (T = read_triplet (stdin)) != NULL ; ) /* RAND */
{
if (T->nrow > 1000000)
{
/* do huge-problem tests only */
huge ( ) ;
CHOLMOD(free_triplet) (&T, cm) ;
continue ;
}
maxerr = 0 ;
CHOLMOD(defaults) (cm) ;
cm->error_handler = my_handler ;
cm->print = 4 ;
cm->precise = FALSE ;
cm->nmethods = 8 ;
/* ------------------------------------------------------------------ */
/* convert triplet to sparse matrix */
/* ------------------------------------------------------------------ */
A = CHOLMOD(triplet_to_sparse) (T, 0, cm) ;
AT = CHOLMOD(transpose) (A, 0, cm) ;
OK (CHOLMOD(print_sparse) (A, "Test matrix, A", cm)) ;
C = unpack (A) ; /* RAND */
OK (CHOLMOD(print_sparse) (C, "Unpacked/unsorted version of A", cm)) ;
cm->print = 1 ;
if (T != NULL)
{
nrow = T->nrow ;
ncol = T->ncol ;
n = MAX (nrow, ncol) ;
nmin = MIN (nrow, ncol) ;
}
/* ------------------------------------------------------------------ */
/* basic error tests */
/* ------------------------------------------------------------------ */
null2 (T, do_nantests) ; /* RAND */
printf ("Null2 OK : no error\n") ;
if (do_nantests)
{
maxerr = 0 ;
goto done ; /* yes, this is ugly */
}
/* ------------------------------------------------------------------ */
/* raw factorization tests */
/* ------------------------------------------------------------------ */
cm->error_handler = NULL ;
err = raw_factor (A, 2) ; /* RAND */
cm->error_handler = my_handler ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error %.1g\n", err) ;
err = raw_factor2 (A, 0., 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error2 %.1g\n", err) ;
err = raw_factor2 (A, 1e-16, 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error3 %.1g\n", err) ;
if (n < 1000 && A && T && A->stype == 1)
{
/* factorize a symmetric matrix (upper part stored), possibly
* with ignored entries in lower triangular part. */
cholmod_sparse *F ;
int save = T->stype ;
T->stype = 0 ;
F = CHOLMOD(triplet_to_sparse) (T, 0, cm) ;
T->stype = save ;
/*
ET = CHOLMOD(transpose) (E, 2, cm) ;
if (E) E->stype = 0 ;
if (ET) ET->stype = 0 ;
F = CHOLMOD(add) (E, ET, one, one, 1, 1, cm) ;
*/
if (F) F->stype = 1 ;
err = raw_factor2 (F, 0., 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error4 %.1g\n", err) ;
err = raw_factor2 (F, 1e-16, 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error5 %.1g\n", err) ;
cm->dbound = 1e-15 ;
err = raw_factor2 (F, 0., 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error6 %.1g\n", err) ;
err = raw_factor2 (F, 1e-16, 0) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error7 %.1g\n", err) ;
err = raw_factor2 (F, 1e-16, 1) ;
MAXERR (maxerr, err, 1) ;
printf ("raw factorization error8 %.1g\n", err) ;
cm->dbound = 0 ;
/*
CHOLMOD(free_sparse) (&E, cm) ;
CHOLMOD(free_sparse) (&ET, cm) ;
*/
CHOLMOD(free_sparse) (&F, cm) ;
}
/* ------------------------------------------------------------------ */
/* matrix ops */
/* ------------------------------------------------------------------ */
err = test_ops (A) ; /* RAND */
MAXERR (maxerr, err, 1) ;
printf ("initial testops error %.1g\n", err) ;
/* ------------------------------------------------------------------ */
/* analyze, factorize, solve */
/* ------------------------------------------------------------------ */
err = solve (A) ; /* RAND */
MAXERR (maxerr, err, 1) ;
printf ("initial solve error %.1g\n", err) ;
/* ------------------------------------------------------------------ */
/* CCOLAMD tests */
/* ------------------------------------------------------------------ */
cctest (A) ; /* RAND reset */
cctest (AT) ; /* RAND reset */
/* ------------------------------------------------------------------ */
/* COLAMD tests */
/* ------------------------------------------------------------------ */
ctest (A) ;
ctest (AT) ;
/* ------------------------------------------------------------------ */
/* AMD tests */
/* ------------------------------------------------------------------ */
if (n < NLARGE || A->stype)
{
/* for unsymmetric matrices, this forms A*A' and A'*A, which can
* fail if A has a dense row or column. So it is only done for
* modest-sized unsymmetric matrices. */
amdtest (A) ;
amdtest (AT) ;
camdtest (A) ; /* RAND */
camdtest (AT) ; /* RAND */
}
if (n < NSMALL)
{
/* -------------------------------------------------------------- */
/* do_matrix with an unpacked matrix */
/* -------------------------------------------------------------- */
/* try with an unpacked matrix, and a non-default dbound */
cm->dbound = 1e-15 ;
err = do_matrix (C) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
cm->dbound = 0 ;
/* -------------------------------------------------------------- */
/* do_matrix: analyze, factorize, and solve, with many options */
/* -------------------------------------------------------------- */
err = do_matrix (A) ; /* RAND reset */
MAXERR (maxerr, err, 1) ;
/* -------------------------------------------------------------- */
/* pretend to solve an LP */
/* -------------------------------------------------------------- */
if (nrow != ncol)
{
cm->print = 2 ;
err = lpdemo (T) ; /* RAND */
cm->print = 1 ;
MAXERR (maxerr, err, 1) ;
cm->print = 5; CHOLMOD(print_common) ("Common", cm);cm->print=1;
cm->nmethods = 1 ;
cm->method [0].ordering = CHOLMOD_COLAMD ;
err = lpdemo (T) ; /* RAND */
MAXERR (maxerr, err, 1) ;
printf ("initial lp error %.1g, dbound %g\n", err, cm->dbound) ;
cm->nmethods = 0 ;
cm->method [0].ordering = CHOLMOD_GIVEN ;
}
}
progress (1, '|') ;
if (n < NSMALL && do_memory)
{
/* -------------------------------------------------------------- */
/* Exhaustive memory-error handling */
/* -------------------------------------------------------------- */
memory_tests (T) ; /* RAND */
}
/* ------------------------------------------------------------------ */
/* free matrices and print results */
/* ------------------------------------------------------------------ */
done: /* an ugly "goto" target; added to minimize code changes
* when added "do_nantests", for version 1.0.2 */
CHOLMOD(free_sparse) (&C, cm) ;
CHOLMOD(free_sparse) (&AT, cm) ;
CHOLMOD(free_sparse) (&A, cm) ;
CHOLMOD(free_triplet) (&T, cm) ;
fprintf (stderr, "\n "
" Test OK") ;
if (nrow <= ncol && !singular)
{
/* maxerr should be NaN if nrow > ncol, so don't print it */
fprintf (stderr, ", maxerr %.1g", maxerr) ;
}
fprintf (stderr, "\n") ;
my_srand (42) ; /* RAND reset */
}
/* ---------------------------------------------------------------------- */
/* finalize CHOLMOD */
/* ---------------------------------------------------------------------- */
CHOLMOD(free_dense) (&M1, cm) ;
CHOLMOD(finish) (cm) ;
cm->print = 5 ; OK (CHOLMOD(print_common) ("cm", cm)) ;
printf ("malloc count "ID" inuse "ID"\n",
(Int) cm->malloc_count,
(Int) cm->memory_inuse) ;
OK (cm->malloc_count == 0) ;
OK (cm->memory_inuse == 0) ;
if (nrow > ncol)
{
/* maxerr should be NaN, so don't print it */
printf ("All tests successful\n") ;
}
else
{
printf ("All tests successful: max error %.1g\n", maxerr) ;
}
return (0) ;
}
syntax highlighted by Code2HTML, v. 0.9.1