1 SUBROUTINE ZGET02( 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 * November 2006
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 RWORK( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGET02 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) / ( norm(A) * norm(X) * 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^T*x = b, where A^T is the transpose of A
33 * = 'C': A^H*x = b, where A^H is the conjugate 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 COMPLEX*16 CONE
80 PARAMETER ( CONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER J, N1, N2
84 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
89 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANGE
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZGEMM
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Quick exit if M = 0 or N = 0 or NRHS = 0
100 *
101 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
102 RESID = ZERO
103 RETURN
104 END IF
105 *
106 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
107 N1 = N
108 N2 = M
109 ELSE
110 N1 = M
111 N2 = N
112 END IF
113 *
114 * Exit with RESID = 1/EPS if ANORM = 0.
115 *
116 EPS = DLAMCH( 'Epsilon' )
117 ANORM = ZLANGE( '1', N1, N2, A, LDA, RWORK )
118 IF( ANORM.LE.ZERO ) THEN
119 RESID = ONE / EPS
120 RETURN
121 END IF
122 *
123 * Compute B - A*X (or B - A'*X ) and store in B.
124 *
125 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
126 $ LDX, CONE, B, LDB )
127 *
128 * Compute the maximum over the number of right hand sides of
129 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
130 *
131 RESID = ZERO
132 DO 10 J = 1, NRHS
133 BNORM = DZASUM( N1, B( 1, J ), 1 )
134 XNORM = DZASUM( N2, X( 1, J ), 1 )
135 IF( XNORM.LE.ZERO ) THEN
136 RESID = ONE / EPS
137 ELSE
138 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
139 END IF
140 10 CONTINUE
141 *
142 RETURN
143 *
144 * End of ZGET02
145 *
146 END
2 $ RWORK, RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
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 RWORK( * )
15 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGET02 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) / ( norm(A) * norm(X) * 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^T*x = b, where A^T is the transpose of A
33 * = 'C': A^H*x = b, where A^H is the conjugate 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 COMPLEX*16 CONE
80 PARAMETER ( CONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER J, N1, N2
84 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
89 EXTERNAL LSAME, DLAMCH, DZASUM, ZLANGE
90 * ..
91 * .. External Subroutines ..
92 EXTERNAL ZGEMM
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 * Quick exit if M = 0 or N = 0 or NRHS = 0
100 *
101 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.EQ.0 ) THEN
102 RESID = ZERO
103 RETURN
104 END IF
105 *
106 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
107 N1 = N
108 N2 = M
109 ELSE
110 N1 = M
111 N2 = N
112 END IF
113 *
114 * Exit with RESID = 1/EPS if ANORM = 0.
115 *
116 EPS = DLAMCH( 'Epsilon' )
117 ANORM = ZLANGE( '1', N1, N2, A, LDA, RWORK )
118 IF( ANORM.LE.ZERO ) THEN
119 RESID = ONE / EPS
120 RETURN
121 END IF
122 *
123 * Compute B - A*X (or B - A'*X ) and store in B.
124 *
125 CALL ZGEMM( TRANS, 'No transpose', N1, NRHS, N2, -CONE, A, LDA, X,
126 $ LDX, CONE, B, LDB )
127 *
128 * Compute the maximum over the number of right hand sides of
129 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ) .
130 *
131 RESID = ZERO
132 DO 10 J = 1, NRHS
133 BNORM = DZASUM( N1, B( 1, J ), 1 )
134 XNORM = DZASUM( N2, X( 1, J ), 1 )
135 IF( XNORM.LE.ZERO ) THEN
136 RESID = ONE / EPS
137 ELSE
138 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
139 END IF
140 10 CONTINUE
141 *
142 RETURN
143 *
144 * End of ZGET02
145 *
146 END