1 SUBROUTINE CPBT05( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, X, LDX,
2 $ XACT, LDXACT, FERR, BERR, RESLTS )
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 UPLO
10 INTEGER KD, LDAB, LDB, LDX, LDXACT, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 REAL BERR( * ), FERR( * ), RESLTS( * )
14 COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
15 $ XACT( LDXACT, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CPBT05 tests the error bounds from iterative refinement for the
22 * computed solution to a system of equations A*X = B, where A is a
23 * Hermitian band matrix.
24 *
25 * RESLTS(1) = test of the error bound
26 * = norm(X - XACT) / ( norm(X) * FERR )
27 *
28 * A large value is returned if this ratio is not less than one.
29 *
30 * RESLTS(2) = residual from the iterative refinement routine
31 * = the maximum of BERR / ( NZ*EPS + (*) ), where
32 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
33 * and NZ = max. number of nonzeros in any row of A, plus 1
34 *
35 * Arguments
36 * =========
37 *
38 * UPLO (input) CHARACTER*1
39 * Specifies whether the upper or lower triangular part of the
40 * Hermitian matrix A is stored.
41 * = 'U': Upper triangular
42 * = 'L': Lower triangular
43 *
44 * N (input) INTEGER
45 * The number of rows of the matrices X, B, and XACT, and the
46 * order of the matrix A. N >= 0.
47 *
48 * KD (input) INTEGER
49 * The number of super-diagonals of the matrix A if UPLO = 'U',
50 * or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
51 *
52 * NRHS (input) INTEGER
53 * The number of columns of the matrices X, B, and XACT.
54 * NRHS >= 0.
55 *
56 * AB (input) COMPLEX array, dimension (LDAB,N)
57 * The upper or lower triangle of the Hermitian band matrix A,
58 * stored in the first KD+1 rows of the array. The j-th column
59 * of A is stored in the j-th column of the array AB as follows:
60 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
61 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
62 *
63 * LDAB (input) INTEGER
64 * The leading dimension of the array AB. LDAB >= KD+1.
65 *
66 * B (input) COMPLEX array, dimension (LDB,NRHS)
67 * The right hand side vectors for the system of linear
68 * equations.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(1,N).
72 *
73 * X (input) COMPLEX array, dimension (LDX,NRHS)
74 * The computed solution vectors. Each vector is stored as a
75 * column of the matrix X.
76 *
77 * LDX (input) INTEGER
78 * The leading dimension of the array X. LDX >= max(1,N).
79 *
80 * XACT (input) COMPLEX array, dimension (LDX,NRHS)
81 * The exact solution vectors. Each vector is stored as a
82 * column of the matrix XACT.
83 *
84 * LDXACT (input) INTEGER
85 * The leading dimension of the array XACT. LDXACT >= max(1,N).
86 *
87 * FERR (input) REAL array, dimension (NRHS)
88 * The estimated forward error bounds for each solution vector
89 * X. If XTRUE is the true solution, FERR bounds the magnitude
90 * of the largest entry in (X - XTRUE) divided by the magnitude
91 * of the largest entry in X.
92 *
93 * BERR (input) REAL array, dimension (NRHS)
94 * The componentwise relative backward error of each solution
95 * vector (i.e., the smallest relative change in any entry of A
96 * or B that makes X an exact solution).
97 *
98 * RESLTS (output) REAL array, dimension (2)
99 * The maximum over the NRHS solution vectors of the ratios:
100 * RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
101 * RESLTS(2) = BERR / ( NZ*EPS + (*) )
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106 REAL ZERO, ONE
107 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
108 * ..
109 * .. Local Scalars ..
110 LOGICAL UPPER
111 INTEGER I, IMAX, J, K, NZ
112 REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
113 COMPLEX ZDUM
114 * ..
115 * .. External Functions ..
116 LOGICAL LSAME
117 INTEGER ICAMAX
118 REAL SLAMCH
119 EXTERNAL LSAME, ICAMAX, SLAMCH
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC ABS, AIMAG, MAX, MIN, REAL
123 * ..
124 * .. Statement Functions ..
125 REAL CABS1
126 * ..
127 * .. Statement Function definitions ..
128 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
129 * ..
130 * .. Executable Statements ..
131 *
132 * Quick exit if N = 0 or NRHS = 0.
133 *
134 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
135 RESLTS( 1 ) = ZERO
136 RESLTS( 2 ) = ZERO
137 RETURN
138 END IF
139 *
140 EPS = SLAMCH( 'Epsilon' )
141 UNFL = SLAMCH( 'Safe minimum' )
142 OVFL = ONE / UNFL
143 UPPER = LSAME( UPLO, 'U' )
144 NZ = 2*MAX( KD, N-1 ) + 1
145 *
146 * Test 1: Compute the maximum of
147 * norm(X - XACT) / ( norm(X) * FERR )
148 * over all the vectors X and XACT using the infinity-norm.
149 *
150 ERRBND = ZERO
151 DO 30 J = 1, NRHS
152 IMAX = ICAMAX( N, X( 1, J ), 1 )
153 XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
154 DIFF = ZERO
155 DO 10 I = 1, N
156 DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
157 10 CONTINUE
158 *
159 IF( XNORM.GT.ONE ) THEN
160 GO TO 20
161 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
162 GO TO 20
163 ELSE
164 ERRBND = ONE / EPS
165 GO TO 30
166 END IF
167 *
168 20 CONTINUE
169 IF( DIFF / XNORM.LE.FERR( J ) ) THEN
170 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
171 ELSE
172 ERRBND = ONE / EPS
173 END IF
174 30 CONTINUE
175 RESLTS( 1 ) = ERRBND
176 *
177 * Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
178 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
179 *
180 DO 90 K = 1, NRHS
181 DO 80 I = 1, N
182 TMP = CABS1( B( I, K ) )
183 IF( UPPER ) THEN
184 DO 40 J = MAX( I-KD, 1 ), I - 1
185 TMP = TMP + CABS1( AB( KD+1-I+J, I ) )*
186 $ CABS1( X( J, K ) )
187 40 CONTINUE
188 TMP = TMP + ABS( REAL( AB( KD+1, I ) ) )*
189 $ CABS1( X( I, K ) )
190 DO 50 J = I + 1, MIN( I+KD, N )
191 TMP = TMP + CABS1( AB( KD+1+I-J, J ) )*
192 $ CABS1( X( J, K ) )
193 50 CONTINUE
194 ELSE
195 DO 60 J = MAX( I-KD, 1 ), I - 1
196 TMP = TMP + CABS1( AB( 1+I-J, J ) )*CABS1( X( J, K ) )
197 60 CONTINUE
198 TMP = TMP + ABS( REAL( AB( 1, I ) ) )*CABS1( X( I, K ) )
199 DO 70 J = I + 1, MIN( I+KD, N )
200 TMP = TMP + CABS1( AB( 1+J-I, I ) )*CABS1( X( J, K ) )
201 70 CONTINUE
202 END IF
203 IF( I.EQ.1 ) THEN
204 AXBI = TMP
205 ELSE
206 AXBI = MIN( AXBI, TMP )
207 END IF
208 80 CONTINUE
209 TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
210 IF( K.EQ.1 ) THEN
211 RESLTS( 2 ) = TMP
212 ELSE
213 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
214 END IF
215 90 CONTINUE
216 *
217 RETURN
218 *
219 * End of CPBT05
220 *
221 END
2 $ XACT, LDXACT, FERR, BERR, RESLTS )
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 UPLO
10 INTEGER KD, LDAB, LDB, LDX, LDXACT, N, NRHS
11 * ..
12 * .. Array Arguments ..
13 REAL BERR( * ), FERR( * ), RESLTS( * )
14 COMPLEX AB( LDAB, * ), B( LDB, * ), X( LDX, * ),
15 $ XACT( LDXACT, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * CPBT05 tests the error bounds from iterative refinement for the
22 * computed solution to a system of equations A*X = B, where A is a
23 * Hermitian band matrix.
24 *
25 * RESLTS(1) = test of the error bound
26 * = norm(X - XACT) / ( norm(X) * FERR )
27 *
28 * A large value is returned if this ratio is not less than one.
29 *
30 * RESLTS(2) = residual from the iterative refinement routine
31 * = the maximum of BERR / ( NZ*EPS + (*) ), where
32 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
33 * and NZ = max. number of nonzeros in any row of A, plus 1
34 *
35 * Arguments
36 * =========
37 *
38 * UPLO (input) CHARACTER*1
39 * Specifies whether the upper or lower triangular part of the
40 * Hermitian matrix A is stored.
41 * = 'U': Upper triangular
42 * = 'L': Lower triangular
43 *
44 * N (input) INTEGER
45 * The number of rows of the matrices X, B, and XACT, and the
46 * order of the matrix A. N >= 0.
47 *
48 * KD (input) INTEGER
49 * The number of super-diagonals of the matrix A if UPLO = 'U',
50 * or the number of sub-diagonals if UPLO = 'L'. KD >= 0.
51 *
52 * NRHS (input) INTEGER
53 * The number of columns of the matrices X, B, and XACT.
54 * NRHS >= 0.
55 *
56 * AB (input) COMPLEX array, dimension (LDAB,N)
57 * The upper or lower triangle of the Hermitian band matrix A,
58 * stored in the first KD+1 rows of the array. The j-th column
59 * of A is stored in the j-th column of the array AB as follows:
60 * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
61 * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
62 *
63 * LDAB (input) INTEGER
64 * The leading dimension of the array AB. LDAB >= KD+1.
65 *
66 * B (input) COMPLEX array, dimension (LDB,NRHS)
67 * The right hand side vectors for the system of linear
68 * equations.
69 *
70 * LDB (input) INTEGER
71 * The leading dimension of the array B. LDB >= max(1,N).
72 *
73 * X (input) COMPLEX array, dimension (LDX,NRHS)
74 * The computed solution vectors. Each vector is stored as a
75 * column of the matrix X.
76 *
77 * LDX (input) INTEGER
78 * The leading dimension of the array X. LDX >= max(1,N).
79 *
80 * XACT (input) COMPLEX array, dimension (LDX,NRHS)
81 * The exact solution vectors. Each vector is stored as a
82 * column of the matrix XACT.
83 *
84 * LDXACT (input) INTEGER
85 * The leading dimension of the array XACT. LDXACT >= max(1,N).
86 *
87 * FERR (input) REAL array, dimension (NRHS)
88 * The estimated forward error bounds for each solution vector
89 * X. If XTRUE is the true solution, FERR bounds the magnitude
90 * of the largest entry in (X - XTRUE) divided by the magnitude
91 * of the largest entry in X.
92 *
93 * BERR (input) REAL array, dimension (NRHS)
94 * The componentwise relative backward error of each solution
95 * vector (i.e., the smallest relative change in any entry of A
96 * or B that makes X an exact solution).
97 *
98 * RESLTS (output) REAL array, dimension (2)
99 * The maximum over the NRHS solution vectors of the ratios:
100 * RESLTS(1) = norm(X - XACT) / ( norm(X) * FERR )
101 * RESLTS(2) = BERR / ( NZ*EPS + (*) )
102 *
103 * =====================================================================
104 *
105 * .. Parameters ..
106 REAL ZERO, ONE
107 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
108 * ..
109 * .. Local Scalars ..
110 LOGICAL UPPER
111 INTEGER I, IMAX, J, K, NZ
112 REAL AXBI, DIFF, EPS, ERRBND, OVFL, TMP, UNFL, XNORM
113 COMPLEX ZDUM
114 * ..
115 * .. External Functions ..
116 LOGICAL LSAME
117 INTEGER ICAMAX
118 REAL SLAMCH
119 EXTERNAL LSAME, ICAMAX, SLAMCH
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC ABS, AIMAG, MAX, MIN, REAL
123 * ..
124 * .. Statement Functions ..
125 REAL CABS1
126 * ..
127 * .. Statement Function definitions ..
128 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
129 * ..
130 * .. Executable Statements ..
131 *
132 * Quick exit if N = 0 or NRHS = 0.
133 *
134 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
135 RESLTS( 1 ) = ZERO
136 RESLTS( 2 ) = ZERO
137 RETURN
138 END IF
139 *
140 EPS = SLAMCH( 'Epsilon' )
141 UNFL = SLAMCH( 'Safe minimum' )
142 OVFL = ONE / UNFL
143 UPPER = LSAME( UPLO, 'U' )
144 NZ = 2*MAX( KD, N-1 ) + 1
145 *
146 * Test 1: Compute the maximum of
147 * norm(X - XACT) / ( norm(X) * FERR )
148 * over all the vectors X and XACT using the infinity-norm.
149 *
150 ERRBND = ZERO
151 DO 30 J = 1, NRHS
152 IMAX = ICAMAX( N, X( 1, J ), 1 )
153 XNORM = MAX( CABS1( X( IMAX, J ) ), UNFL )
154 DIFF = ZERO
155 DO 10 I = 1, N
156 DIFF = MAX( DIFF, CABS1( X( I, J )-XACT( I, J ) ) )
157 10 CONTINUE
158 *
159 IF( XNORM.GT.ONE ) THEN
160 GO TO 20
161 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
162 GO TO 20
163 ELSE
164 ERRBND = ONE / EPS
165 GO TO 30
166 END IF
167 *
168 20 CONTINUE
169 IF( DIFF / XNORM.LE.FERR( J ) ) THEN
170 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
171 ELSE
172 ERRBND = ONE / EPS
173 END IF
174 30 CONTINUE
175 RESLTS( 1 ) = ERRBND
176 *
177 * Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
178 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
179 *
180 DO 90 K = 1, NRHS
181 DO 80 I = 1, N
182 TMP = CABS1( B( I, K ) )
183 IF( UPPER ) THEN
184 DO 40 J = MAX( I-KD, 1 ), I - 1
185 TMP = TMP + CABS1( AB( KD+1-I+J, I ) )*
186 $ CABS1( X( J, K ) )
187 40 CONTINUE
188 TMP = TMP + ABS( REAL( AB( KD+1, I ) ) )*
189 $ CABS1( X( I, K ) )
190 DO 50 J = I + 1, MIN( I+KD, N )
191 TMP = TMP + CABS1( AB( KD+1+I-J, J ) )*
192 $ CABS1( X( J, K ) )
193 50 CONTINUE
194 ELSE
195 DO 60 J = MAX( I-KD, 1 ), I - 1
196 TMP = TMP + CABS1( AB( 1+I-J, J ) )*CABS1( X( J, K ) )
197 60 CONTINUE
198 TMP = TMP + ABS( REAL( AB( 1, I ) ) )*CABS1( X( I, K ) )
199 DO 70 J = I + 1, MIN( I+KD, N )
200 TMP = TMP + CABS1( AB( 1+J-I, I ) )*CABS1( X( J, K ) )
201 70 CONTINUE
202 END IF
203 IF( I.EQ.1 ) THEN
204 AXBI = TMP
205 ELSE
206 AXBI = MIN( AXBI, TMP )
207 END IF
208 80 CONTINUE
209 TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
210 IF( K.EQ.1 ) THEN
211 RESLTS( 2 ) = TMP
212 ELSE
213 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
214 END IF
215 90 CONTINUE
216 *
217 RETURN
218 *
219 * End of CPBT05
220 *
221 END