00001 /* zspt02.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static doublecomplex c_b1 = {1.,0.}; 00019 static integer c__1 = 1; 00020 00021 /* Subroutine */ int zspt02_(char *uplo, integer *n, integer *nrhs, 00022 doublecomplex *a, doublecomplex *x, integer *ldx, doublecomplex *b, 00023 integer *ldb, doublereal *rwork, doublereal *resid) 00024 { 00025 /* System generated locals */ 00026 integer b_dim1, b_offset, x_dim1, x_offset, i__1; 00027 doublereal d__1, d__2; 00028 doublecomplex z__1; 00029 00030 /* Local variables */ 00031 integer j; 00032 doublereal eps, anorm, bnorm, xnorm; 00033 extern /* Subroutine */ int zspmv_(char *, integer *, doublecomplex *, 00034 doublecomplex *, doublecomplex *, integer *, doublecomplex *, 00035 doublecomplex *, integer *); 00036 extern doublereal dlamch_(char *), dzasum_(integer *, 00037 doublecomplex *, integer *), zlansp_(char *, char *, integer *, 00038 doublecomplex *, doublereal *); 00039 00040 00041 /* -- LAPACK test routine (version 3.1) -- */ 00042 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00043 /* November 2006 */ 00044 00045 /* .. Scalar Arguments .. */ 00046 /* .. */ 00047 /* .. Array Arguments .. */ 00048 /* .. */ 00049 00050 /* Purpose */ 00051 /* ======= */ 00052 00053 /* ZSPT02 computes the residual in the solution of a complex symmetric */ 00054 /* system of linear equations A*x = b when packed storage is used for */ 00055 /* the coefficient matrix. The ratio computed is */ 00056 00057 /* RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS). */ 00058 00059 /* where EPS is the machine precision. */ 00060 00061 /* Arguments */ 00062 /* ========= */ 00063 00064 /* UPLO (input) CHARACTER*1 */ 00065 /* Specifies whether the upper or lower triangular part of the */ 00066 /* complex symmetric matrix A is stored: */ 00067 /* = 'U': Upper triangular */ 00068 /* = 'L': Lower triangular */ 00069 00070 /* N (input) INTEGER */ 00071 /* The number of rows and columns of the matrix A. N >= 0. */ 00072 00073 /* NRHS (input) INTEGER */ 00074 /* The number of columns of B, the matrix of right hand sides. */ 00075 /* NRHS >= 0. */ 00076 00077 /* A (input) COMPLEX*16 array, dimension (N*(N+1)/2) */ 00078 /* The original complex symmetric matrix A, stored as a packed */ 00079 /* triangular matrix. */ 00080 00081 /* X (input) COMPLEX*16 array, dimension (LDX,NRHS) */ 00082 /* The computed solution vectors for the system of linear */ 00083 /* equations. */ 00084 00085 /* LDX (input) INTEGER */ 00086 /* The leading dimension of the array X. LDX >= max(1,N). */ 00087 00088 /* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) */ 00089 /* On entry, the right hand side vectors for the system of */ 00090 /* linear equations. */ 00091 /* On exit, B is overwritten with the difference B - A*X. */ 00092 00093 /* LDB (input) INTEGER */ 00094 /* The leading dimension of the array B. LDB >= max(1,N). */ 00095 00096 /* RWORK (workspace) DOUBLE PRECISION array, dimension (N) */ 00097 00098 /* RESID (output) DOUBLE PRECISION */ 00099 /* The maximum over the number of right hand sides of */ 00100 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). */ 00101 00102 /* ===================================================================== */ 00103 00104 /* .. Parameters .. */ 00105 /* .. */ 00106 /* .. Local Scalars .. */ 00107 /* .. */ 00108 /* .. External Functions .. */ 00109 /* .. */ 00110 /* .. External Subroutines .. */ 00111 /* .. */ 00112 /* .. Intrinsic Functions .. */ 00113 /* .. */ 00114 /* .. Executable Statements .. */ 00115 00116 /* Quick exit if N = 0 or NRHS = 0 */ 00117 00118 /* Parameter adjustments */ 00119 --a; 00120 x_dim1 = *ldx; 00121 x_offset = 1 + x_dim1; 00122 x -= x_offset; 00123 b_dim1 = *ldb; 00124 b_offset = 1 + b_dim1; 00125 b -= b_offset; 00126 --rwork; 00127 00128 /* Function Body */ 00129 if (*n <= 0 || *nrhs <= 0) { 00130 *resid = 0.; 00131 return 0; 00132 } 00133 00134 /* Exit with RESID = 1/EPS if ANORM = 0. */ 00135 00136 eps = dlamch_("Epsilon"); 00137 anorm = zlansp_("1", uplo, n, &a[1], &rwork[1]); 00138 if (anorm <= 0.) { 00139 *resid = 1. / eps; 00140 return 0; 00141 } 00142 00143 /* Compute B - A*X for the matrix of right hand sides B. */ 00144 00145 i__1 = *nrhs; 00146 for (j = 1; j <= i__1; ++j) { 00147 z__1.r = -1., z__1.i = -0.; 00148 zspmv_(uplo, n, &z__1, &a[1], &x[j * x_dim1 + 1], &c__1, &c_b1, &b[j * 00149 b_dim1 + 1], &c__1); 00150 /* L10: */ 00151 } 00152 00153 /* Compute the maximum over the number of right hand sides of */ 00154 /* norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) . */ 00155 00156 *resid = 0.; 00157 i__1 = *nrhs; 00158 for (j = 1; j <= i__1; ++j) { 00159 bnorm = dzasum_(n, &b[j * b_dim1 + 1], &c__1); 00160 xnorm = dzasum_(n, &x[j * x_dim1 + 1], &c__1); 00161 if (xnorm <= 0.) { 00162 *resid = 1. / eps; 00163 } else { 00164 /* Computing MAX */ 00165 d__1 = *resid, d__2 = bnorm / anorm / xnorm / eps; 00166 *resid = max(d__1,d__2); 00167 } 00168 /* L20: */ 00169 } 00170 00171 return 0; 00172 00173 /* End of ZSPT02 */ 00174 00175 } /* zspt02_ */