00001 /* zget02.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 zget02_(char *trans, integer *m, integer *n, integer * 00022 nrhs, doublecomplex *a, integer *lda, doublecomplex *x, integer *ldx, 00023 doublecomplex *b, integer *ldb, doublereal *rwork, doublereal *resid) 00024 { 00025 /* System generated locals */ 00026 integer a_dim1, a_offset, 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, n1, n2; 00032 doublereal eps; 00033 extern logical lsame_(char *, char *); 00034 doublereal anorm, bnorm; 00035 extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 00036 integer *, doublecomplex *, doublecomplex *, integer *, 00037 doublecomplex *, integer *, doublecomplex *, doublecomplex *, 00038 integer *); 00039 doublereal xnorm; 00040 extern doublereal dlamch_(char *), zlange_(char *, integer *, 00041 integer *, doublecomplex *, integer *, doublereal *), 00042 dzasum_(integer *, doublecomplex *, integer *); 00043 00044 00045 /* -- LAPACK test routine (version 3.1) -- */ 00046 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00047 /* November 2006 */ 00048 00049 /* .. Scalar Arguments .. */ 00050 /* .. */ 00051 /* .. Array Arguments .. */ 00052 /* .. */ 00053 00054 /* Purpose */ 00055 /* ======= */ 00056 00057 /* ZGET02 computes the residual for a solution of a system of linear */ 00058 /* equations A*x = b or A'*x = b: */ 00059 /* RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), */ 00060 /* where EPS is the machine epsilon. */ 00061 00062 /* Arguments */ 00063 /* ========= */ 00064 00065 /* TRANS (input) CHARACTER*1 */ 00066 /* Specifies the form of the system of equations: */ 00067 /* = 'N': A *x = b */ 00068 /* = 'T': A^T*x = b, where A^T is the transpose of A */ 00069 /* = 'C': A^H*x = b, where A^H is the conjugate transpose of A */ 00070 00071 /* M (input) INTEGER */ 00072 /* The number of rows of the matrix A. M >= 0. */ 00073 00074 /* N (input) INTEGER */ 00075 /* The number of columns of the matrix A. N >= 0. */ 00076 00077 /* NRHS (input) INTEGER */ 00078 /* The number of columns of B, the matrix of right hand sides. */ 00079 /* NRHS >= 0. */ 00080 00081 /* A (input) COMPLEX*16 array, dimension (LDA,N) */ 00082 /* The original M x N matrix A. */ 00083 00084 /* LDA (input) INTEGER */ 00085 /* The leading dimension of the array A. LDA >= max(1,M). */ 00086 00087 /* X (input) COMPLEX*16 array, dimension (LDX,NRHS) */ 00088 /* The computed solution vectors for the system of linear */ 00089 /* equations. */ 00090 00091 /* LDX (input) INTEGER */ 00092 /* The leading dimension of the array X. If TRANS = 'N', */ 00093 /* LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). */ 00094 00095 /* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) */ 00096 /* On entry, the right hand side vectors for the system of */ 00097 /* linear equations. */ 00098 /* On exit, B is overwritten with the difference B - A*X. */ 00099 00100 /* LDB (input) INTEGER */ 00101 /* The leading dimension of the array B. IF TRANS = 'N', */ 00102 /* LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). */ 00103 00104 /* RWORK (workspace) DOUBLE PRECISION array, dimension (M) */ 00105 00106 /* RESID (output) DOUBLE PRECISION */ 00107 /* The maximum over the number of right hand sides of */ 00108 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ). */ 00109 00110 /* ===================================================================== */ 00111 00112 /* .. Parameters .. */ 00113 /* .. */ 00114 /* .. Local Scalars .. */ 00115 /* .. */ 00116 /* .. External Functions .. */ 00117 /* .. */ 00118 /* .. External Subroutines .. */ 00119 /* .. */ 00120 /* .. Intrinsic Functions .. */ 00121 /* .. */ 00122 /* .. Executable Statements .. */ 00123 00124 /* Quick exit if M = 0 or N = 0 or NRHS = 0 */ 00125 00126 /* Parameter adjustments */ 00127 a_dim1 = *lda; 00128 a_offset = 1 + a_dim1; 00129 a -= a_offset; 00130 x_dim1 = *ldx; 00131 x_offset = 1 + x_dim1; 00132 x -= x_offset; 00133 b_dim1 = *ldb; 00134 b_offset = 1 + b_dim1; 00135 b -= b_offset; 00136 --rwork; 00137 00138 /* Function Body */ 00139 if (*m <= 0 || *n <= 0 || *nrhs == 0) { 00140 *resid = 0.; 00141 return 0; 00142 } 00143 00144 if (lsame_(trans, "T") || lsame_(trans, "C")) { 00145 n1 = *n; 00146 n2 = *m; 00147 } else { 00148 n1 = *m; 00149 n2 = *n; 00150 } 00151 00152 /* Exit with RESID = 1/EPS if ANORM = 0. */ 00153 00154 eps = dlamch_("Epsilon"); 00155 anorm = zlange_("1", &n1, &n2, &a[a_offset], lda, &rwork[1]); 00156 if (anorm <= 0.) { 00157 *resid = 1. / eps; 00158 return 0; 00159 } 00160 00161 /* Compute B - A*X (or B - A'*X ) and store in B. */ 00162 00163 z__1.r = -1., z__1.i = -0.; 00164 zgemm_(trans, "No transpose", &n1, nrhs, &n2, &z__1, &a[a_offset], lda, & 00165 x[x_offset], ldx, &c_b1, &b[b_offset], ldb) 00166 ; 00167 00168 /* Compute the maximum over the number of right hand sides of */ 00169 /* norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . */ 00170 00171 *resid = 0.; 00172 i__1 = *nrhs; 00173 for (j = 1; j <= i__1; ++j) { 00174 bnorm = dzasum_(&n1, &b[j * b_dim1 + 1], &c__1); 00175 xnorm = dzasum_(&n2, &x[j * x_dim1 + 1], &c__1); 00176 if (xnorm <= 0.) { 00177 *resid = 1. / eps; 00178 } else { 00179 /* Computing MAX */ 00180 d__1 = *resid, d__2 = bnorm / anorm / xnorm / eps; 00181 *resid = max(d__1,d__2); 00182 } 00183 /* L10: */ 00184 } 00185 00186 return 0; 00187 00188 /* End of ZGET02 */ 00189 00190 } /* zget02_ */