/usr/share/doc/libsuperlu-dev/tests/sp_dget04.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 | /*
* -- SuperLU routine (version 2.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November 15, 1997
*
*/
#include <math.h>
#include "slu_ddefs.h"
int sp_dget04(int n, int nrhs, double *x, int ldx, double *xact,
int ldxact, double rcond, double *resid)
{
/*
Purpose
=======
SP_DGET04 computes the difference between a computed solution and the
true solution to a system of linear equations.
RESID = ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
where RCOND is the reciprocal of the condition number and EPS is the
machine epsilon.
Arguments
=========
N (input) INT
The number of rows of the matrices X and XACT. N >= 0.
NRHS (input) INT
The number of columns of the matrices X and XACT. NRHS >= 0.
X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
The computed solution vectors. Each vector is stored as a
column of the matrix X.
LDX (input) INT
The leading dimension of the array X. LDX >= max(1,N).
XACT (input) DOUBLE PRECISION array, dimension( LDX, NRHS )
The exact solution vectors. Each vector is stored as a
column of the matrix XACT.
LDXACT (input) INT
The leading dimension of the array XACT. LDXACT >= max(1,N).
RCOND (input) DOUBLE PRECISION
The reciprocal of the condition number of the coefficient
matrix in the system of equations.
RESID (output) DOUBLE PRECISION
The maximum over the NRHS solution vectors of
( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS )
=====================================================================
*/
/* Table of constant values */
int c__1 = 1;
/* System generated locals */
double d__1, d__2, d__3, d__4;
/* Local variables */
int i, j, n__1;
int ix;
double xnorm;
double eps;
double diffnm;
/* Function prototypes */
extern int idamax_(int *, double *, int *);
extern double dlamch_(char *);
/* Quick exit if N = 0 or NRHS = 0. */
if ( n <= 0 || nrhs <= 0 ) {
*resid = 0.;
return 0;
}
/* Exit with RESID = 1/EPS if RCOND is invalid. */
eps = dlamch_("Epsilon");
if ( rcond < 0. ) {
*resid = 1. / eps;
return 0;
}
/* Compute the maximum of norm(X - XACT) / ( norm(XACT) * EPS )
over all the vectors X and XACT . */
*resid = 0.;
for (j = 0; j < nrhs; ++j) {
n__1 = n;
ix = idamax_(&n__1, &xact[j*ldxact], &c__1);
xnorm = (d__1 = xact[ix-1 + j*ldxact], fabs(d__1));
diffnm = 0.;
for (i = 0; i < n; ++i) {
/* Computing MAX */
d__3 = diffnm;
d__4 = (d__1 = x[i+j*ldx]-xact[i+j*ldxact], fabs(d__1));
diffnm = SUPERLU_MAX(d__3,d__4);
}
if (xnorm <= 0.) {
if (diffnm > 0.) {
*resid = 1. / eps;
}
} else {
/* Computing MAX */
d__1 = *resid, d__2 = diffnm / xnorm * rcond;
*resid = SUPERLU_MAX(d__1,d__2);
}
}
if (*resid * eps < 1.) {
*resid /= eps;
}
return 0;
} /* sp_dget04_ */
|