1 SUBROUTINE DGET08( TRANS, M, N, NRHS, A, LDA, X, LDX, B, LDB,
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 CHARACTER TRANS
10 INTEGER LDA, LDB, LDX, M, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), RWORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DGET08 computes the residual for a solution of a system of linear
22 * equations A*x = b or A'*x = b:
23 * RESID = norm(B - A*X,inf) / ( norm(A,inf) * norm(X,inf) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANS (input) CHARACTER*1
30 * Specifies the form of the system of equations:
31 * = 'N': A *x = b
32 * = 'T': A'*x = b, where A' is the transpose of A
33 * = 'C': A'*x = b, where A' is the transpose of A
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of columns of the matrix A. N >= 0.
40 *
41 * NRHS (input) INTEGER
42 * The number of columns of B, the matrix of right hand sides.
43 * NRHS >= 0.
44 *
45 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
46 * The original M x N matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
52 * The computed solution vectors for the system of linear
53 * equations.
54 *
55 * LDX (input) INTEGER
56 * The leading dimension of the array X. If TRANS = 'N',
57 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
58 *
59 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
60 * On entry, the right hand side vectors for the system of
61 * linear equations.
62 * On exit, B is overwritten with the difference B - A*X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. IF TRANS = 'N',
66 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * The maximum over the number of right hand sides of
72 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 * ..
80 * .. Local Scalars ..
81 INTEGER J, N1, N2
82 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 INTEGER IDAMAX
87 DOUBLE PRECISION DLAMCH, DLANGE
88 EXTERNAL LSAME, IDAMAX, DLAMCH, DLANGE
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL DGEMM
92 * ..
93 * .. Intrinsic Functions ..
94 INTRINSIC MAX, ABS
95 * ..
96 * .. Executable Statements ..
97 *
98 * Quick exit if M = 0 or N = 0 or NRHS = 0
99 *
100 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
101 RESID = ZERO
102 RETURN
103 END IF
104 *
105 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
106 N1 = N
107 N2 = M
108 ELSE
109 N1 = M
110 N2 = N
111 END IF
112 *
113 * Exit with RESID = 1/EPS if ANORM = 0.
114 *
115 EPS = DLAMCH( 'Epsilon' )
116 ANORM = DLANGE( 'I', N1, N2, A, LDA, RWORK )
117 IF( ANORM.LE.ZERO ) THEN
118 RESID = ONE / EPS
119 RETURN
120 END IF
121 *
122 * Compute B - A*X (or B - A'*X ) and store in B.
123 *
124 CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X,
125 $ LDX, ONE, B, LDB )
126 *
127 * Compute the maximum over the number of right hand sides of
128 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
129 *
130 RESID = ZERO
131 DO 10 J = 1, NRHS
132 BNORM = ABS(B(IDAMAX( N1, B( 1, J ), 1 ),J))
133 XNORM = ABS(X(IDAMAX( N2, X( 1, J ), 1 ),J))
134 IF( XNORM.LE.ZERO ) THEN
135 RESID = ONE / EPS
136 ELSE
137 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
138 END IF
139 10 CONTINUE
140 *
141 RETURN
142 *
143 * End of DGET02
144 *
145 END
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * June 2010
7 *
8 * .. Scalar Arguments ..
9 CHARACTER TRANS
10 INTEGER LDA, LDB, LDX, M, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), RWORK( * ),
15 $ X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DGET08 computes the residual for a solution of a system of linear
22 * equations A*x = b or A'*x = b:
23 * RESID = norm(B - A*X,inf) / ( norm(A,inf) * norm(X,inf) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * =========
28 *
29 * TRANS (input) CHARACTER*1
30 * Specifies the form of the system of equations:
31 * = 'N': A *x = b
32 * = 'T': A'*x = b, where A' is the transpose of A
33 * = 'C': A'*x = b, where A' is the transpose of A
34 *
35 * M (input) INTEGER
36 * The number of rows of the matrix A. M >= 0.
37 *
38 * N (input) INTEGER
39 * The number of columns of the matrix A. N >= 0.
40 *
41 * NRHS (input) INTEGER
42 * The number of columns of B, the matrix of right hand sides.
43 * NRHS >= 0.
44 *
45 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
46 * The original M x N matrix A.
47 *
48 * LDA (input) INTEGER
49 * The leading dimension of the array A. LDA >= max(1,M).
50 *
51 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
52 * The computed solution vectors for the system of linear
53 * equations.
54 *
55 * LDX (input) INTEGER
56 * The leading dimension of the array X. If TRANS = 'N',
57 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
58 *
59 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
60 * On entry, the right hand side vectors for the system of
61 * linear equations.
62 * On exit, B is overwritten with the difference B - A*X.
63 *
64 * LDB (input) INTEGER
65 * The leading dimension of the array B. IF TRANS = 'N',
66 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
67 *
68 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
69 *
70 * RESID (output) DOUBLE PRECISION
71 * The maximum over the number of right hand sides of
72 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
73 *
74 * =====================================================================
75 *
76 * .. Parameters ..
77 DOUBLE PRECISION ZERO, ONE
78 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
79 * ..
80 * .. Local Scalars ..
81 INTEGER J, N1, N2
82 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 INTEGER IDAMAX
87 DOUBLE PRECISION DLAMCH, DLANGE
88 EXTERNAL LSAME, IDAMAX, DLAMCH, DLANGE
89 * ..
90 * .. External Subroutines ..
91 EXTERNAL DGEMM
92 * ..
93 * .. Intrinsic Functions ..
94 INTRINSIC MAX, ABS
95 * ..
96 * .. Executable Statements ..
97 *
98 * Quick exit if M = 0 or N = 0 or NRHS = 0
99 *
100 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
101 RESID = ZERO
102 RETURN
103 END IF
104 *
105 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
106 N1 = N
107 N2 = M
108 ELSE
109 N1 = M
110 N2 = N
111 END IF
112 *
113 * Exit with RESID = 1/EPS if ANORM = 0.
114 *
115 EPS = DLAMCH( 'Epsilon' )
116 ANORM = DLANGE( 'I', N1, N2, A, LDA, RWORK )
117 IF( ANORM.LE.ZERO ) THEN
118 RESID = ONE / EPS
119 RETURN
120 END IF
121 *
122 * Compute B - A*X (or B - A'*X ) and store in B.
123 *
124 CALL DGEMM( TRANS, 'No transpose', N1, NRHS, N2, -ONE, A, LDA, X,
125 $ LDX, ONE, B, LDB )
126 *
127 * Compute the maximum over the number of right hand sides of
128 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
129 *
130 RESID = ZERO
131 DO 10 J = 1, NRHS
132 BNORM = ABS(B(IDAMAX( N1, B( 1, J ), 1 ),J))
133 XNORM = ABS(X(IDAMAX( N2, X( 1, J ), 1 ),J))
134 IF( XNORM.LE.ZERO ) THEN
135 RESID = ONE / EPS
136 ELSE
137 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
138 END IF
139 10 CONTINUE
140 *
141 RETURN
142 *
143 * End of DGET02
144 *
145 END