00001 /* cget02.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 complex c_b1 = {1.f,0.f}; 00019 static integer c__1 = 1; 00020 00021 /* Subroutine */ int cget02_(char *trans, integer *m, integer *n, integer * 00022 nrhs, complex *a, integer *lda, complex *x, integer *ldx, complex *b, 00023 integer *ldb, real *rwork, real *resid) 00024 { 00025 /* System generated locals */ 00026 integer a_dim1, a_offset, b_dim1, b_offset, x_dim1, x_offset, i__1; 00027 real r__1, r__2; 00028 complex q__1; 00029 00030 /* Local variables */ 00031 integer j, n1, n2; 00032 real eps; 00033 extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, 00034 integer *, complex *, complex *, integer *, complex *, integer *, 00035 complex *, complex *, integer *); 00036 extern logical lsame_(char *, char *); 00037 real anorm, bnorm, xnorm; 00038 extern doublereal clange_(char *, integer *, integer *, complex *, 00039 integer *, real *), slamch_(char *), scasum_( 00040 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 /* CGET02 computes the residual for a solution of a system of linear */ 00056 /* equations A*x = b or A'*x = b: */ 00057 /* RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), */ 00058 /* where EPS is the machine epsilon. */ 00059 00060 /* Arguments */ 00061 /* ========= */ 00062 00063 /* TRANS (input) CHARACTER*1 */ 00064 /* Specifies the form of the system of equations: */ 00065 /* = 'N': A *x = b */ 00066 /* = 'T': A^T*x = b, where A^T is the transpose of A */ 00067 /* = 'C': A^H*x = b, where A^H is the conjugate transpose of A */ 00068 00069 /* M (input) INTEGER */ 00070 /* The number of rows of the matrix A. M >= 0. */ 00071 00072 /* N (input) INTEGER */ 00073 /* The number of columns of the matrix A. N >= 0. */ 00074 00075 /* NRHS (input) INTEGER */ 00076 /* The number of columns of B, the matrix of right hand sides. */ 00077 /* NRHS >= 0. */ 00078 00079 /* A (input) COMPLEX array, dimension (LDA,N) */ 00080 /* The original M x N matrix A. */ 00081 00082 /* LDA (input) INTEGER */ 00083 /* The leading dimension of the array A. LDA >= max(1,M). */ 00084 00085 /* X (input) COMPLEX 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. If TRANS = 'N', */ 00091 /* LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). */ 00092 00093 /* B (input/output) COMPLEX array, dimension (LDB,NRHS) */ 00094 /* On entry, the right hand side vectors for the system of */ 00095 /* linear equations. */ 00096 /* On exit, B is overwritten with the difference B - A*X. */ 00097 00098 /* LDB (input) INTEGER */ 00099 /* The leading dimension of the array B. IF TRANS = 'N', */ 00100 /* LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). */ 00101 00102 /* RWORK (workspace) REAL array, dimension (M) */ 00103 00104 /* RESID (output) REAL */ 00105 /* The maximum over the number of right hand sides of */ 00106 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). */ 00107 00108 /* ===================================================================== */ 00109 00110 /* .. Parameters .. */ 00111 /* .. */ 00112 /* .. Local Scalars .. */ 00113 /* .. */ 00114 /* .. External Functions .. */ 00115 /* .. */ 00116 /* .. External Subroutines .. */ 00117 /* .. */ 00118 /* .. Intrinsic Functions .. */ 00119 /* .. */ 00120 /* .. Executable Statements .. */ 00121 00122 /* Quick exit if M = 0 or N = 0 or NRHS = 0 */ 00123 00124 /* Parameter adjustments */ 00125 a_dim1 = *lda; 00126 a_offset = 1 + a_dim1; 00127 a -= a_offset; 00128 x_dim1 = *ldx; 00129 x_offset = 1 + x_dim1; 00130 x -= x_offset; 00131 b_dim1 = *ldb; 00132 b_offset = 1 + b_dim1; 00133 b -= b_offset; 00134 --rwork; 00135 00136 /* Function Body */ 00137 if (*m <= 0 || *n <= 0 || *nrhs == 0) { 00138 *resid = 0.f; 00139 return 0; 00140 } 00141 00142 if (lsame_(trans, "T") || lsame_(trans, "C")) { 00143 n1 = *n; 00144 n2 = *m; 00145 } else { 00146 n1 = *m; 00147 n2 = *n; 00148 } 00149 00150 /* Exit with RESID = 1/EPS if ANORM = 0. */ 00151 00152 eps = slamch_("Epsilon"); 00153 anorm = clange_("1", &n1, &n2, &a[a_offset], lda, &rwork[1]); 00154 if (anorm <= 0.f) { 00155 *resid = 1.f / eps; 00156 return 0; 00157 } 00158 00159 /* Compute B - A*X (or B - A'*X ) and store in B. */ 00160 00161 q__1.r = -1.f, q__1.i = -0.f; 00162 cgemm_(trans, "No transpose", &n1, nrhs, &n2, &q__1, &a[a_offset], lda, & 00163 x[x_offset], ldx, &c_b1, &b[b_offset], ldb) 00164 ; 00165 00166 /* Compute the maximum over the number of right hand sides of */ 00167 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . */ 00168 00169 *resid = 0.f; 00170 i__1 = *nrhs; 00171 for (j = 1; j <= i__1; ++j) { 00172 bnorm = scasum_(&n1, &b[j * b_dim1 + 1], &c__1); 00173 xnorm = scasum_(&n2, &x[j * x_dim1 + 1], &c__1); 00174 if (xnorm <= 0.f) { 00175 *resid = 1.f / eps; 00176 } else { 00177 /* Computing MAX */ 00178 r__1 = *resid, r__2 = bnorm / anorm / xnorm / eps; 00179 *resid = dmax(r__1,r__2); 00180 } 00181 /* L10: */ 00182 } 00183 00184 return 0; 00185 00186 /* End of CGET02 */ 00187 00188 } /* cget02_ */