00001 /* cgtt02.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_b6 = -1.f; 00019 static real c_b7 = 1.f; 00020 static integer c__1 = 1; 00021 00022 /* Subroutine */ int cgtt02_(char *trans, integer *n, integer *nrhs, complex * 00023 dl, complex *d__, complex *du, complex *x, integer *ldx, complex *b, 00024 integer *ldb, real *rwork, real *resid) 00025 { 00026 /* System generated locals */ 00027 integer 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; 00033 extern logical lsame_(char *, char *); 00034 real anorm, bnorm, xnorm; 00035 extern doublereal slamch_(char *), clangt_(char *, integer *, 00036 complex *, complex *, complex *); 00037 extern /* Subroutine */ int clagtm_(char *, integer *, integer *, real *, 00038 complex *, complex *, complex *, complex *, integer *, real *, 00039 complex *, integer *); 00040 extern doublereal scasum_(integer *, complex *, integer *); 00041 00042 00043 /* -- LAPACK test routine (version 3.1) -- */ 00044 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00045 /* November 2006 */ 00046 00047 /* .. Scalar Arguments .. */ 00048 /* .. */ 00049 /* .. Array Arguments .. */ 00050 /* .. */ 00051 00052 /* Purpose */ 00053 /* ======= */ 00054 00055 /* CGTT02 computes the residual for the solution to a tridiagonal */ 00056 /* system of equations: */ 00057 /* RESID = norm(B - op(A)*X) / (norm(A) * norm(X) * EPS), */ 00058 /* where EPS is the machine epsilon. */ 00059 00060 /* Arguments */ 00061 /* ========= */ 00062 00063 /* TRANS (input) CHARACTER */ 00064 /* Specifies the form of the residual. */ 00065 /* = 'N': B - A * X (No transpose) */ 00066 /* = 'T': B - A**T * X (Transpose) */ 00067 /* = 'C': B - A**H * X (Conjugate transpose) */ 00068 00069 /* N (input) INTEGTER */ 00070 /* The order of the matrix A. N >= 0. */ 00071 00072 /* NRHS (input) INTEGER */ 00073 /* The number of right hand sides, i.e., the number of columns */ 00074 /* of the matrices B and X. NRHS >= 0. */ 00075 00076 /* DL (input) COMPLEX array, dimension (N-1) */ 00077 /* The (n-1) sub-diagonal elements of A. */ 00078 00079 /* D (input) COMPLEX array, dimension (N) */ 00080 /* The diagonal elements of A. */ 00081 00082 /* DU (input) COMPLEX array, dimension (N-1) */ 00083 /* The (n-1) super-diagonal elements of A. */ 00084 00085 /* X (input) COMPLEX array, dimension (LDX,NRHS) */ 00086 /* The computed solution vectors X. */ 00087 00088 /* LDX (input) INTEGER */ 00089 /* The leading dimension of the array X. LDX >= max(1,N). */ 00090 00091 /* B (input/output) COMPLEX array, dimension (LDB,NRHS) */ 00092 /* On entry, the right hand side vectors for the system of */ 00093 /* linear equations. */ 00094 /* On exit, B is overwritten with the difference B - op(A)*X. */ 00095 00096 /* LDB (input) INTEGER */ 00097 /* The leading dimension of the array B. LDB >= max(1,N). */ 00098 00099 /* RWORK (workspace) REAL array, dimension (N) */ 00100 00101 /* RESID (output) REAL */ 00102 /* norm(B - op(A)*X) / (norm(A) * norm(X) * EPS) */ 00103 00104 /* ===================================================================== */ 00105 00106 /* .. Parameters .. */ 00107 /* .. */ 00108 /* .. Local Scalars .. */ 00109 /* .. */ 00110 /* .. External Functions .. */ 00111 /* .. */ 00112 /* .. External Subroutines .. */ 00113 /* .. */ 00114 /* .. Intrinsic Functions .. */ 00115 /* .. */ 00116 /* .. Executable Statements .. */ 00117 00118 /* Quick exit if N = 0 or NRHS = 0 */ 00119 00120 /* Parameter adjustments */ 00121 --dl; 00122 --d__; 00123 --du; 00124 x_dim1 = *ldx; 00125 x_offset = 1 + x_dim1; 00126 x -= x_offset; 00127 b_dim1 = *ldb; 00128 b_offset = 1 + b_dim1; 00129 b -= b_offset; 00130 --rwork; 00131 00132 /* Function Body */ 00133 *resid = 0.f; 00134 if (*n <= 0 || *nrhs == 0) { 00135 return 0; 00136 } 00137 00138 /* Compute the maximum over the number of right hand sides of */ 00139 /* norm(B - op(A)*X) / ( norm(A) * norm(X) * EPS ). */ 00140 00141 if (lsame_(trans, "N")) { 00142 anorm = clangt_("1", n, &dl[1], &d__[1], &du[1]); 00143 } else { 00144 anorm = clangt_("I", n, &dl[1], &d__[1], &du[1]); 00145 } 00146 00147 /* Exit with RESID = 1/EPS if ANORM = 0. */ 00148 00149 eps = slamch_("Epsilon"); 00150 if (anorm <= 0.f) { 00151 *resid = 1.f / eps; 00152 return 0; 00153 } 00154 00155 /* Compute B - op(A)*X. */ 00156 00157 clagtm_(trans, n, nrhs, &c_b6, &dl[1], &d__[1], &du[1], &x[x_offset], ldx, 00158 &c_b7, &b[b_offset], ldb); 00159 00160 i__1 = *nrhs; 00161 for (j = 1; j <= i__1; ++j) { 00162 bnorm = scasum_(n, &b[j * b_dim1 + 1], &c__1); 00163 xnorm = scasum_(n, &x[j * x_dim1 + 1], &c__1); 00164 if (xnorm <= 0.f) { 00165 *resid = 1.f / eps; 00166 } else { 00167 /* Computing MAX */ 00168 r__1 = *resid, r__2 = bnorm / anorm / xnorm / eps; 00169 *resid = dmax(r__1,r__2); 00170 } 00171 /* L10: */ 00172 } 00173 00174 return 0; 00175 00176 /* End of CGTT02 */ 00177 00178 } /* cgtt02_ */