/usr/share/doc/libsuperlu-dev/tests/sp_cget02.c is in libsuperlu3-dev 3.0+20070106-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 | /*
* -- SuperLU routine (version 3.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* October 15, 2003
*
*/
#include "slu_cdefs.h"
int sp_cget02(trans_t trans, int m, int n, int nrhs, SuperMatrix *A,
complex *x, int ldx, complex *b, int ldb, float *resid)
{
/*
Purpose
=======
SP_CGET02 computes the residual for a solution of a system of linear
equations A*x = b or A'*x = b:
RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
where EPS is the machine epsilon.
Arguments
=========
TRANS (input) trans_t
Specifies the form of the system of equations:
= NOTRANS: A *x = b
= TRANS : A'*x = b, where A' is the transpose of A
= CONJ : A'*x = b, where A' is the transpose of A
M (input) INTEGER
The number of rows of the matrix A. M >= 0.
N (input) INTEGER
The number of columns of the matrix A. N >= 0.
NRHS (input) INTEGER
The number of columns of B, the matrix of right hand sides.
NRHS >= 0.
A (input) SuperMatrix*, dimension (LDA,N)
The original M x N sparse matrix A.
X (input) COMPLEX PRECISION array, dimension (LDX,NRHS)
The computed solution vectors for the system of linear
equations.
LDX (input) INTEGER
The leading dimension of the array X. If TRANS = NOTRANS,
LDX >= max(1,N); if TRANS = TRANS or CONJ, LDX >= max(1,M).
B (input/output) COMPLEX PRECISION array, dimension (LDB,NRHS)
On entry, the right hand side vectors for the system of
linear equations.
On exit, B is overwritten with the difference B - A*X.
LDB (input) INTEGER
The leading dimension of the array B. IF TRANS = NOTRANS,
LDB >= max(1,M); if TRANS = TRANS or CONJ, LDB >= max(1,N).
RESID (output) FLOAT PRECISION
The maximum over the number of right hand sides of
norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
=====================================================================
*/
/* Table of constant values */
complex alpha = {-1., 0.0};
complex beta = {1., 0.0};
int c__1 = 1;
/* System generated locals */
float d__1, d__2;
/* Local variables */
int j;
int n1, n2;
float anorm, bnorm;
float xnorm;
float eps;
char transc[1];
/* Function prototypes */
extern int lsame_(char *, char *);
extern float clangs(char *, SuperMatrix *);
extern float scasum_(int *, complex *, int *);
extern double slamch_(char *);
/* Function Body */
if ( m <= 0 || n <= 0 || nrhs == 0) {
*resid = 0.;
return 0;
}
if ( (trans == TRANS) || (trans == CONJ) ) {
n1 = n;
n2 = m;
*transc = 'T';
} else {
n1 = m;
n2 = n;
*transc = 'N';
}
/* Exit with RESID = 1/EPS if ANORM = 0. */
eps = slamch_("Epsilon");
anorm = clangs("1", A);
if (anorm <= 0.) {
*resid = 1. / eps;
return 0;
}
/* Compute B - A*X (or B - A'*X ) and store in B. */
sp_cgemm(transc, "N", n1, nrhs, n2, alpha, A, x, ldx, beta, b, ldb);
/* Compute the maximum over the number of right hand sides of
norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . */
*resid = 0.;
for (j = 0; j < nrhs; ++j) {
bnorm = scasum_(&n1, &b[j*ldb], &c__1);
xnorm = scasum_(&n2, &x[j*ldx], &c__1);
if (xnorm <= 0.) {
*resid = 1. / eps;
} else {
/* Computing MAX */
d__1 = *resid, d__2 = bnorm / anorm / xnorm / eps;
*resid = SUPERLU_MAX(d__1, d__2);
}
}
return 0;
} /* sp_cget02 */
|