1 SUBROUTINE ZGBT02( TRANS, M, N, KL, KU, NRHS, A, LDA, X, LDX, B,
2 $ LDB, 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 KL, KU, LDA, LDB, LDX, M, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGBT02 computes the residual for a solution of a banded system of
21 * equations A*x = b or A'*x = b:
22 * RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS).
23 * where EPS is the machine precision.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANS (input) CHARACTER*1
29 * Specifies the form of the system of equations:
30 * = 'N': A *x = b
31 * = 'T': A'*x = b, where A' is the transpose of A
32 * = 'C': A'*x = b, where A' is the transpose of A
33 *
34 * M (input) INTEGER
35 * The number of rows of the matrix A. M >= 0.
36 *
37 * N (input) INTEGER
38 * The number of columns of the matrix A. N >= 0.
39 *
40 * KL (input) INTEGER
41 * The number of subdiagonals within the band of A. KL >= 0.
42 *
43 * KU (input) INTEGER
44 * The number of superdiagonals within the band of A. KU >= 0.
45 *
46 * NRHS (input) INTEGER
47 * The number of columns of B. NRHS >= 0.
48 *
49 * A (input) COMPLEX*16 array, dimension (LDA,N)
50 * The original matrix A in band storage, stored in rows 1 to
51 * KL+KU+1.
52 *
53 * LDA (input) INTEGER
54 * The leading dimension of the array A. LDA >= max(1,KL+KU+1).
55 *
56 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
57 * The computed solution vectors for the system of linear
58 * equations.
59 *
60 * LDX (input) INTEGER
61 * The leading dimension of the array X. If TRANS = 'N',
62 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
63 *
64 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
65 * On entry, the right hand side vectors for the system of
66 * linear equations.
67 * On exit, B is overwritten with the difference B - A*X.
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of the array B. IF TRANS = 'N',
71 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
72 *
73 * RESID (output) DOUBLE PRECISION
74 * The maximum over the number of right hand sides of
75 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE
81 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
82 COMPLEX*16 CONE
83 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
84 * ..
85 * .. Local Scalars ..
86 INTEGER I1, I2, J, KD, N1
87 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
88 * ..
89 * .. External Functions ..
90 LOGICAL LSAME
91 DOUBLE PRECISION DLAMCH, DZASUM
92 EXTERNAL LSAME, DLAMCH, DZASUM
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL ZGBMV
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX, MIN
99 * ..
100 * .. Executable Statements ..
101 *
102 * Quick return if N = 0 pr NRHS = 0
103 *
104 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.LE.0 ) THEN
105 RESID = ZERO
106 RETURN
107 END IF
108 *
109 * Exit with RESID = 1/EPS if ANORM = 0.
110 *
111 EPS = DLAMCH( 'Epsilon' )
112 KD = KU + 1
113 ANORM = ZERO
114 DO 10 J = 1, N
115 I1 = MAX( KD+1-J, 1 )
116 I2 = MIN( KD+M-J, KL+KD )
117 ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
118 10 CONTINUE
119 IF( ANORM.LE.ZERO ) THEN
120 RESID = ONE / EPS
121 RETURN
122 END IF
123 *
124 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
125 N1 = N
126 ELSE
127 N1 = M
128 END IF
129 *
130 * Compute B - A*X (or B - A'*X )
131 *
132 DO 20 J = 1, NRHS
133 CALL ZGBMV( TRANS, M, N, KL, KU, -CONE, A, LDA, X( 1, J ), 1,
134 $ CONE, B( 1, J ), 1 )
135 20 CONTINUE
136 *
137 * Compute the maximum over the number of right hand sides of
138 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
139 *
140 RESID = ZERO
141 DO 30 J = 1, NRHS
142 BNORM = DZASUM( N1, B( 1, J ), 1 )
143 XNORM = DZASUM( N1, X( 1, J ), 1 )
144 IF( XNORM.LE.ZERO ) THEN
145 RESID = ONE / EPS
146 ELSE
147 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
148 END IF
149 30 CONTINUE
150 *
151 RETURN
152 *
153 * End of ZGBT02
154 *
155 END
2 $ LDB, 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 KL, KU, LDA, LDB, LDX, M, N, NRHS
11 DOUBLE PRECISION RESID
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), B( LDB, * ), X( LDX, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZGBT02 computes the residual for a solution of a banded system of
21 * equations A*x = b or A'*x = b:
22 * RESID = norm( B - A*X ) / ( norm(A) * norm(X) * EPS).
23 * where EPS is the machine precision.
24 *
25 * Arguments
26 * =========
27 *
28 * TRANS (input) CHARACTER*1
29 * Specifies the form of the system of equations:
30 * = 'N': A *x = b
31 * = 'T': A'*x = b, where A' is the transpose of A
32 * = 'C': A'*x = b, where A' is the transpose of A
33 *
34 * M (input) INTEGER
35 * The number of rows of the matrix A. M >= 0.
36 *
37 * N (input) INTEGER
38 * The number of columns of the matrix A. N >= 0.
39 *
40 * KL (input) INTEGER
41 * The number of subdiagonals within the band of A. KL >= 0.
42 *
43 * KU (input) INTEGER
44 * The number of superdiagonals within the band of A. KU >= 0.
45 *
46 * NRHS (input) INTEGER
47 * The number of columns of B. NRHS >= 0.
48 *
49 * A (input) COMPLEX*16 array, dimension (LDA,N)
50 * The original matrix A in band storage, stored in rows 1 to
51 * KL+KU+1.
52 *
53 * LDA (input) INTEGER
54 * The leading dimension of the array A. LDA >= max(1,KL+KU+1).
55 *
56 * X (input) COMPLEX*16 array, dimension (LDX,NRHS)
57 * The computed solution vectors for the system of linear
58 * equations.
59 *
60 * LDX (input) INTEGER
61 * The leading dimension of the array X. If TRANS = 'N',
62 * LDX >= max(1,N); if TRANS = 'T' or 'C', LDX >= max(1,M).
63 *
64 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
65 * On entry, the right hand side vectors for the system of
66 * linear equations.
67 * On exit, B is overwritten with the difference B - A*X.
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of the array B. IF TRANS = 'N',
71 * LDB >= max(1,M); if TRANS = 'T' or 'C', LDB >= max(1,N).
72 *
73 * RESID (output) DOUBLE PRECISION
74 * The maximum over the number of right hand sides of
75 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
76 *
77 * =====================================================================
78 *
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO, ONE
81 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
82 COMPLEX*16 CONE
83 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
84 * ..
85 * .. Local Scalars ..
86 INTEGER I1, I2, J, KD, N1
87 DOUBLE PRECISION ANORM, BNORM, EPS, XNORM
88 * ..
89 * .. External Functions ..
90 LOGICAL LSAME
91 DOUBLE PRECISION DLAMCH, DZASUM
92 EXTERNAL LSAME, DLAMCH, DZASUM
93 * ..
94 * .. External Subroutines ..
95 EXTERNAL ZGBMV
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX, MIN
99 * ..
100 * .. Executable Statements ..
101 *
102 * Quick return if N = 0 pr NRHS = 0
103 *
104 IF( M.LE.0 .OR. N.LE.0 .OR. NRHS.LE.0 ) THEN
105 RESID = ZERO
106 RETURN
107 END IF
108 *
109 * Exit with RESID = 1/EPS if ANORM = 0.
110 *
111 EPS = DLAMCH( 'Epsilon' )
112 KD = KU + 1
113 ANORM = ZERO
114 DO 10 J = 1, N
115 I1 = MAX( KD+1-J, 1 )
116 I2 = MIN( KD+M-J, KL+KD )
117 ANORM = MAX( ANORM, DZASUM( I2-I1+1, A( I1, J ), 1 ) )
118 10 CONTINUE
119 IF( ANORM.LE.ZERO ) THEN
120 RESID = ONE / EPS
121 RETURN
122 END IF
123 *
124 IF( LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' ) ) THEN
125 N1 = N
126 ELSE
127 N1 = M
128 END IF
129 *
130 * Compute B - A*X (or B - A'*X )
131 *
132 DO 20 J = 1, NRHS
133 CALL ZGBMV( TRANS, M, N, KL, KU, -CONE, A, LDA, X( 1, J ), 1,
134 $ CONE, B( 1, J ), 1 )
135 20 CONTINUE
136 *
137 * Compute the maximum over the number of right hand sides of
138 * norm(B - A*X) / ( norm(A) * norm(X) * EPS ).
139 *
140 RESID = ZERO
141 DO 30 J = 1, NRHS
142 BNORM = DZASUM( N1, B( 1, J ), 1 )
143 XNORM = DZASUM( N1, X( 1, J ), 1 )
144 IF( XNORM.LE.ZERO ) THEN
145 RESID = ONE / EPS
146 ELSE
147 RESID = MAX( RESID, ( ( BNORM / ANORM ) / XNORM ) / EPS )
148 END IF
149 30 CONTINUE
150 *
151 RETURN
152 *
153 * End of ZGBT02
154 *
155 END