00001 SUBROUTINE ZGET08( 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 * June 2010 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 * ZGET08 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, 0.0D+0 ) ) 00081 * .. 00082 * .. Local Scalars .. 00083 INTEGER J, N1, N2 00084 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 00085 COMPLEX*16 ZDUM 00086 * .. 00087 * .. External Functions .. 00088 LOGICAL LSAME 00089 INTEGER IZAMAX 00090 DOUBLE PRECISION DLAMCH, ZLANGE 00091 EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANGE 00092 * .. 00093 * .. External Subroutines .. 00094 EXTERNAL ZGEMM 00095 * .. 00096 * .. Intrinsic Functions .. 00097 INTRINSIC ABS, DBLE, DIMAG, MAX 00098 * .. 00099 * .. Statement Functions .. 00100 DOUBLE PRECISION CABS1 00101 * .. 00102 * .. Statement Function definitions .. 00103 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00104 * .. 00105 * .. Executable Statements .. 00106 * 00107 * Quick exit if M = 0 or N = 0 or NRHS = 0 00108 * 00109 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN 00110 RESID = ZERO 00111 RETURN 00112 END IF 00113 * 00114 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN 00115 N1 = N 00116 N2 = M 00117 ELSE 00118 N1 = M 00119 N2 = N 00120 END IF 00121 * 00122 * Exit with RESID = 1/EPS if ANORM = 0. 00123 * 00124 EPS = DLAMCH( 'Epsilon' ) 00125 ANORM = ZLANGE( 'I', N1, N2, A, LDA, RWORK ) 00126 IF( ANORM.LE.ZERO ) THEN 00127 RESID = ONE / EPS 00128 RETURN 00129 END IF 00130 * 00131 * Compute B - A*X (or B - A'*X ) and store in B. 00132 * 00133 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X, 00134 $ LDX, CONE, B, LDB ) 00135 * 00136 * Compute the maximum over the number of right hand sides of 00137 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . 00138 * 00139 RESID = ZERO 00140 DO 10 J = 1, NRHS 00141 BNORM = CABS1( B( IZAMAX( N1, B( 1, J ), 1 ), J ) ) 00142 XNORM = CABS1( X( IZAMAX( N2, X( 1, J ), 1 ), J ) ) 00143 IF( XNORM.LE.ZERO ) THEN 00144 RESID = ONE / EPS 00145 ELSE 00146 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00147 END IF 00148 10 CONTINUE 00149 * 00150 RETURN 00151 * 00152 * End of ZGET02 00153 * 00154 END
1.7.2