LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGET02( 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 A( LDA, * ), B( LDB, * ), RWORK( * ), 00015 $ X( LDX, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DGET02 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'*x = b, where A' is the transpose of A 00033 * = 'C': A'*x = b, where A' is the 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 * .. 00080 * .. Local Scalars .. 00081 INTEGER J, N1, N2 00082 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM 00083 * .. 00084 * .. External Functions .. 00085 LOGICAL LSAME 00086 DOUBLE PRECISION DASUM, DLAMCH, DLANGE 00087 EXTERNAL LSAME, DASUM, DLAMCH, DLANGE 00088 * .. 00089 * .. External Subroutines .. 00090 EXTERNAL DGEMM 00091 * .. 00092 * .. Intrinsic Functions .. 00093 INTRINSIC MAX 00094 * .. 00095 * .. Executable Statements .. 00096 * 00097 * Quick exit if M = 0 or N = 0 or NRHS = 0 00098 * 00099 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN 00100 RESID = ZERO 00101 RETURN 00102 END IF 00103 * 00104 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN 00105 N1 = N 00106 N2 = M 00107 ELSE 00108 N1 = M 00109 N2 = N 00110 END IF 00111 * 00112 * Exit with RESID = 1/EPS if ANORM = 0. 00113 * 00114 EPS = DLAMCH( 'Epsilon' ) 00115 ANORM = DLANGE( '1', N1, N2, A, LDA, RWORK ) 00116 IF( ANORM.LE.ZERO ) THEN 00117 RESID = ONE / EPS 00118 RETURN 00119 END IF 00120 * 00121 * Compute B - A*X (or B - A'*X ) and store in B. 00122 * 00123 CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X, 00124 $ LDX, ONE, B, LDB ) 00125 * 00126 * Compute the maximum over the number of right hand sides of 00127 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . 00128 * 00129 RESID = ZERO 00130 DO 10 J = 1, NRHS 00131 BNORM = DASUM( N1, B( 1, J ), 1 ) 00132 XNORM = DASUM( N2, X( 1, J ), 1 ) 00133 IF( XNORM.LE.ZERO ) THEN 00134 RESID = ONE / EPS 00135 ELSE 00136 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS ) 00137 END IF 00138 10 CONTINUE 00139 * 00140 RETURN 00141 * 00142 * End of DGET02 00143 * 00144 END