00001 /* spbt02.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 real c_b5 = -1.f; 00019 static integer c__1 = 1; 00020 static real c_b7 = 1.f; 00021 00022 /* Subroutine */ int spbt02_(char *uplo, integer *n, integer *kd, integer * 00023 nrhs, real *a, integer *lda, real *x, integer *ldx, real *b, integer * 00024 ldb, real *rwork, real *resid) 00025 { 00026 /* System generated locals */ 00027 integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, i__1; 00028 real r__1, r__2; 00029 00030 /* Local variables */ 00031 integer j; 00032 real eps, anorm, bnorm; 00033 extern doublereal sasum_(integer *, real *, integer *); 00034 extern /* Subroutine */ int ssbmv_(char *, integer *, integer *, real *, 00035 real *, integer *, real *, integer *, real *, real *, integer *); 00036 real xnorm; 00037 extern doublereal slamch_(char *), slansb_(char *, char *, 00038 integer *, integer *, real *, integer *, real *); 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 /* SPBT02 computes the residual for a solution of a symmetric banded */ 00054 /* system of equations A*x = b: */ 00055 /* RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS) */ 00056 /* where EPS is the machine precision. */ 00057 00058 /* Arguments */ 00059 /* ========= */ 00060 00061 /* UPLO (input) CHARACTER*1 */ 00062 /* Specifies whether the upper or lower triangular part of the */ 00063 /* symmetric matrix A is stored: */ 00064 /* = 'U': Upper triangular */ 00065 /* = 'L': Lower triangular */ 00066 00067 /* N (input) INTEGER */ 00068 /* The number of rows and columns of the matrix A. N >= 0. */ 00069 00070 /* KD (input) INTEGER */ 00071 /* The number of super-diagonals of the matrix A if UPLO = 'U', */ 00072 /* or the number of sub-diagonals if UPLO = 'L'. KD >= 0. */ 00073 00074 /* A (input) REAL array, dimension (LDA,N) */ 00075 /* The original symmetric band matrix A. If UPLO = 'U', the */ 00076 /* upper triangular part of A is stored as a band matrix; if */ 00077 /* UPLO = 'L', the lower triangular part of A is stored. The */ 00078 /* columns of the appropriate triangle are stored in the columns */ 00079 /* of A and the diagonals of the triangle are stored in the rows */ 00080 /* of A. See SPBTRF for further details. */ 00081 00082 /* LDA (input) INTEGER. */ 00083 /* The leading dimension of the array A. LDA >= max(1,KD+1). */ 00084 00085 /* X (input) REAL array, dimension (LDX,NRHS) */ 00086 /* The computed solution vectors for the system of linear */ 00087 /* equations. */ 00088 00089 /* LDX (input) INTEGER */ 00090 /* The leading dimension of the array X. LDX >= max(1,N). */ 00091 00092 /* B (input/output) REAL array, dimension (LDB,NRHS) */ 00093 /* On entry, the right hand side vectors for the system of */ 00094 /* linear equations. */ 00095 /* On exit, B is overwritten with the difference B - A*X. */ 00096 00097 /* LDB (input) INTEGER */ 00098 /* The leading dimension of the array B. LDB >= max(1,N). */ 00099 00100 /* RWORK (workspace) REAL array, dimension (N) */ 00101 00102 /* RESID (output) REAL */ 00103 /* The maximum over the number of right hand sides of */ 00104 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). */ 00105 00106 /* ===================================================================== */ 00107 00108 /* .. Parameters .. */ 00109 /* .. */ 00110 /* .. Local Scalars .. */ 00111 /* .. */ 00112 /* .. External Functions .. */ 00113 /* .. */ 00114 /* .. External Subroutines .. */ 00115 /* .. */ 00116 /* .. Intrinsic Functions .. */ 00117 /* .. */ 00118 /* .. Executable Statements .. */ 00119 00120 /* Quick exit if N = 0 or NRHS = 0. */ 00121 00122 /* Parameter adjustments */ 00123 a_dim1 = *lda; 00124 a_offset = 1 + a_dim1; 00125 a -= a_offset; 00126 x_dim1 = *ldx; 00127 x_offset = 1 + x_dim1; 00128 x -= x_offset; 00129 b_dim1 = *ldb; 00130 b_offset = 1 + b_dim1; 00131 b -= b_offset; 00132 --rwork; 00133 00134 /* Function Body */ 00135 if (*n <= 0 || *nrhs <= 0) { 00136 *resid = 0.f; 00137 return 0; 00138 } 00139 00140 /* Exit with RESID = 1/EPS if ANORM = 0. */ 00141 00142 eps = slamch_("Epsilon"); 00143 anorm = slansb_("1", uplo, n, kd, &a[a_offset], lda, &rwork[1]); 00144 if (anorm <= 0.f) { 00145 *resid = 1.f / eps; 00146 return 0; 00147 } 00148 00149 /* Compute B - A*X */ 00150 00151 i__1 = *nrhs; 00152 for (j = 1; j <= i__1; ++j) { 00153 ssbmv_(uplo, n, kd, &c_b5, &a[a_offset], lda, &x[j * x_dim1 + 1], & 00154 c__1, &c_b7, &b[j * b_dim1 + 1], &c__1); 00155 /* L10: */ 00156 } 00157 00158 /* Compute the maximum over the number of right hand sides of */ 00159 /* norm( B - A*X ) / ( norm(A) * norm(X) * EPS ) */ 00160 00161 *resid = 0.f; 00162 i__1 = *nrhs; 00163 for (j = 1; j <= i__1; ++j) { 00164 bnorm = sasum_(n, &b[j * b_dim1 + 1], &c__1); 00165 xnorm = sasum_(n, &x[j * x_dim1 + 1], &c__1); 00166 if (xnorm <= 0.f) { 00167 *resid = 1.f / eps; 00168 } else { 00169 /* Computing MAX */ 00170 r__1 = *resid, r__2 = bnorm / anorm / xnorm / eps; 00171 *resid = dmax(r__1,r__2); 00172 } 00173 /* L20: */ 00174 } 00175 00176 return 0; 00177 00178 /* End of SPBT02 */ 00179 00180 } /* spbt02_ */