00001 SUBROUTINE ZGET02( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB, 00002 $ RWORK, RESID ) 00003 * 00004 * -- LAPACK test routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER TRANS 00010 INTEGER LDA, LDB, LDX, M, N, NRHS 00011 DOUBLE PRECISION RESID 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION RWORK( * ) 00015 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * ZGET02 computes the residual for a solution of a system of linear 00022 * equations A*x = b or A'*x = b: 00023 * RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ), 00024 * where EPS is the machine epsilon. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * TRANS (input) CHARACTER*1 00030 * Specifies the form of the system of equations: 00031 * = 'N': A *x = b 00032 * = 'T': A^T*x = b, where A^T is the transpose of A 00033 * = 'C': A^H*x = b, where A^H is the conjugate transpose of A 00034 * 00035 * M (input) INTEGER 00036 * The number of rows of the matrix A. M >= 0. 00037 * 00038 * N (input) INTEGER 00039 * The number of columns of the matrix A. N >= 0. 00040 * 00041 * NRHS (input) INTEGER 00042 * The number of columns of B, the matrix of right hand sides. 00043 * NRHS >= 0. 00044 * 00045 * A (input) COMPLEX*16 array, dimension (LDA,N) 00046 * The original M x N matrix A. 00047 * 00048 * LDA (input) INTEGER 00049 * The leading dimension of the array A. LDA >= max(1,M). 00050 * 00051 * X (input) COMPLEX*16 array, dimension (LDX,NRHS) 00052 * The computed solution vectors for the system of linear 00053 * equations. 00054 * 00055 * LDX (input) INTEGER 00056 * The leading dimension of the array X. If TRANS = 'N', 00057 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M). 00058 * 00059 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) 00060 * On entry, the right hand side vectors for the system of 00061 * linear equations. 00062 * On exit, B is overwritten with the difference B - A*X. 00063 * 00064 * LDB (input) INTEGER 00065 * The leading dimension of the array B. IF TRANS = 'N', 00066 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N). 00067 * 00068 * RWORK (workspace) DOUBLE PRECISION array, dimension (M) 00069 * 00070 * RESID (output) DOUBLE PRECISION 00071 * The maximum over the number of right hand sides of 00072 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ). 00073 * 00074 * ===================================================================== 00075 * 00076 * .. Parameters .. 00077 DOUBLE PRECISION ZERO, ONE 00078 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00079 COMPLEX*16 CONE 00080 PARAMETER ( CONE = 1.0D+0 ) 00081 * .. 00082 * .. Local Scalars .. 00083 INTEGER J, N1, N2 00084 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 00085 * .. 00086 * .. External Functions .. 00087 LOGICAL LSAME 00088 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE 00089 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANGE 00090 * .. 00091 * .. External Subroutines .. 00092 EXTERNAL ZGEMM 00093 * .. 00094 * .. Intrinsic Functions .. 00095 INTRINSIC MAX 00096 * .. 00097 * .. Executable Statements .. 00098 * 00099 * Quick exit if M = 0 or N = 0 or NRHS = 0 00100 * 00101 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN 00102 RESID = ZERO 00103 RETURN 00104 END IF 00105 * 00106 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN 00107 N1 = N 00108 N2 = M 00109 ELSE 00110 N1 = M 00111 N2 = N 00112 END IF 00113 * 00114 * Exit with RESID = 1/EPS if ANORM = 0. 00115 * 00116 EPS = DLAMCH( 'Epsilon' ) 00117 ANORM = ZLANGE( '1', N1, N2, A, LDA, RWORK ) 00118 IF( ANORM.LE.ZERO ) THEN 00119 RESID = ONE / EPS 00120 RETURN 00121 END IF 00122 * 00123 * Compute B - A*X (or B - A'*X ) and store in B. 00124 * 00125 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, 00126 $ LDX, CONE, B, LDB ) 00127 * 00128 * Compute the maximum over the number of right hand sides of 00129 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . 00130 * 00131 RESID = ZERO 00132 DO 10 J = 1, NRHS 00133 BNORM = DZASUM( N1, B( 1, J ), 1 ) 00134 XNORM = DZASUM( N2, X( 1, J ), 1 ) 00135 IF( XNORM.LE.ZERO ) THEN 00136 RESID = ONE / EPS 00137 ELSE 00138 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00139 END IF 00140 10 CONTINUE 00141 * 00142 RETURN 00143 * 00144 * End of ZGET02 00145 * 00146 END
1.7.2