1 SUBROUTINE SPBT05( 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 AB( LDAB, * ), B( LDB, * ), BERR( * ),
14 $ FERR( * ), RESLTS( * ), X( LDX, * ),
15 $ XACT( LDXACT, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * SPBT05 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 * symmetric 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 * symmetric 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) REAL array, dimension (LDAB,N)
57 * The upper or lower triangle of the symmetric 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) REAL 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) REAL 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) REAL 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 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 INTEGER ISAMAX
117 REAL SLAMCH
118 EXTERNAL LSAME, ISAMAX, SLAMCH
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, MAX, MIN
122 * ..
123 * .. Executable Statements ..
124 *
125 * Quick exit if N = 0 or NRHS = 0.
126 *
127 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
128 RESLTS( 1 ) = ZERO
129 RESLTS( 2 ) = ZERO
130 RETURN
131 END IF
132 *
133 EPS = SLAMCH( 'Epsilon' )
134 UNFL = SLAMCH( 'Safe minimum' )
135 OVFL = ONE / UNFL
136 UPPER = LSAME( UPLO, 'U' )
137 NZ = 2*MAX( KD, N-1 ) + 1
138 *
139 * Test 1: Compute the maximum of
140 * norm(X - XACT) / ( norm(X) * FERR )
141 * over all the vectors X and XACT using the infinity-norm.
142 *
143 ERRBND = ZERO
144 DO 30 J = 1, NRHS
145 IMAX = ISAMAX( N, X( 1, J ), 1 )
146 XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
147 DIFF = ZERO
148 DO 10 I = 1, N
149 DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
150 10 CONTINUE
151 *
152 IF( XNORM.GT.ONE ) THEN
153 GO TO 20
154 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
155 GO TO 20
156 ELSE
157 ERRBND = ONE / EPS
158 GO TO 30
159 END IF
160 *
161 20 CONTINUE
162 IF( DIFF / XNORM.LE.FERR( J ) ) THEN
163 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
164 ELSE
165 ERRBND = ONE / EPS
166 END IF
167 30 CONTINUE
168 RESLTS( 1 ) = ERRBND
169 *
170 * Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
171 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
172 *
173 DO 90 K = 1, NRHS
174 DO 80 I = 1, N
175 TMP = ABS( B( I, K ) )
176 IF( UPPER ) THEN
177 DO 40 J = MAX( I-KD, 1 ), I
178 TMP = TMP + ABS( AB( KD+1-I+J, I ) )*ABS( X( J, K ) )
179 40 CONTINUE
180 DO 50 J = I + 1, MIN( I+KD, N )
181 TMP = TMP + ABS( AB( KD+1+I-J, J ) )*ABS( X( J, K ) )
182 50 CONTINUE
183 ELSE
184 DO 60 J = MAX( I-KD, 1 ), I - 1
185 TMP = TMP + ABS( AB( 1+I-J, J ) )*ABS( X( J, K ) )
186 60 CONTINUE
187 DO 70 J = I, MIN( I+KD, N )
188 TMP = TMP + ABS( AB( 1+J-I, I ) )*ABS( X( J, K ) )
189 70 CONTINUE
190 END IF
191 IF( I.EQ.1 ) THEN
192 AXBI = TMP
193 ELSE
194 AXBI = MIN( AXBI, TMP )
195 END IF
196 80 CONTINUE
197 TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
198 IF( K.EQ.1 ) THEN
199 RESLTS( 2 ) = TMP
200 ELSE
201 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
202 END IF
203 90 CONTINUE
204 *
205 RETURN
206 *
207 * End of SPBT05
208 *
209 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 AB( LDAB, * ), B( LDB, * ), BERR( * ),
14 $ FERR( * ), RESLTS( * ), X( LDX, * ),
15 $ XACT( LDXACT, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * SPBT05 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 * symmetric 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 * symmetric 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) REAL array, dimension (LDAB,N)
57 * The upper or lower triangle of the symmetric 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) REAL 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) REAL 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) REAL 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 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 INTEGER ISAMAX
117 REAL SLAMCH
118 EXTERNAL LSAME, ISAMAX, SLAMCH
119 * ..
120 * .. Intrinsic Functions ..
121 INTRINSIC ABS, MAX, MIN
122 * ..
123 * .. Executable Statements ..
124 *
125 * Quick exit if N = 0 or NRHS = 0.
126 *
127 IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
128 RESLTS( 1 ) = ZERO
129 RESLTS( 2 ) = ZERO
130 RETURN
131 END IF
132 *
133 EPS = SLAMCH( 'Epsilon' )
134 UNFL = SLAMCH( 'Safe minimum' )
135 OVFL = ONE / UNFL
136 UPPER = LSAME( UPLO, 'U' )
137 NZ = 2*MAX( KD, N-1 ) + 1
138 *
139 * Test 1: Compute the maximum of
140 * norm(X - XACT) / ( norm(X) * FERR )
141 * over all the vectors X and XACT using the infinity-norm.
142 *
143 ERRBND = ZERO
144 DO 30 J = 1, NRHS
145 IMAX = ISAMAX( N, X( 1, J ), 1 )
146 XNORM = MAX( ABS( X( IMAX, J ) ), UNFL )
147 DIFF = ZERO
148 DO 10 I = 1, N
149 DIFF = MAX( DIFF, ABS( X( I, J )-XACT( I, J ) ) )
150 10 CONTINUE
151 *
152 IF( XNORM.GT.ONE ) THEN
153 GO TO 20
154 ELSE IF( DIFF.LE.OVFL*XNORM ) THEN
155 GO TO 20
156 ELSE
157 ERRBND = ONE / EPS
158 GO TO 30
159 END IF
160 *
161 20 CONTINUE
162 IF( DIFF / XNORM.LE.FERR( J ) ) THEN
163 ERRBND = MAX( ERRBND, ( DIFF / XNORM ) / FERR( J ) )
164 ELSE
165 ERRBND = ONE / EPS
166 END IF
167 30 CONTINUE
168 RESLTS( 1 ) = ERRBND
169 *
170 * Test 2: Compute the maximum of BERR / ( NZ*EPS + (*) ), where
171 * (*) = NZ*UNFL / (min_i (abs(A)*abs(X) +abs(b))_i )
172 *
173 DO 90 K = 1, NRHS
174 DO 80 I = 1, N
175 TMP = ABS( B( I, K ) )
176 IF( UPPER ) THEN
177 DO 40 J = MAX( I-KD, 1 ), I
178 TMP = TMP + ABS( AB( KD+1-I+J, I ) )*ABS( X( J, K ) )
179 40 CONTINUE
180 DO 50 J = I + 1, MIN( I+KD, N )
181 TMP = TMP + ABS( AB( KD+1+I-J, J ) )*ABS( X( J, K ) )
182 50 CONTINUE
183 ELSE
184 DO 60 J = MAX( I-KD, 1 ), I - 1
185 TMP = TMP + ABS( AB( 1+I-J, J ) )*ABS( X( J, K ) )
186 60 CONTINUE
187 DO 70 J = I, MIN( I+KD, N )
188 TMP = TMP + ABS( AB( 1+J-I, I ) )*ABS( X( J, K ) )
189 70 CONTINUE
190 END IF
191 IF( I.EQ.1 ) THEN
192 AXBI = TMP
193 ELSE
194 AXBI = MIN( AXBI, TMP )
195 END IF
196 80 CONTINUE
197 TMP = BERR( K ) / ( NZ*EPS+NZ*UNFL / MAX( AXBI, NZ*UNFL ) )
198 IF( K.EQ.1 ) THEN
199 RESLTS( 2 ) = TMP
200 ELSE
201 RESLTS( 2 ) = MAX( RESLTS( 2 ), TMP )
202 END IF
203 90 CONTINUE
204 *
205 RETURN
206 *
207 * End of SPBT05
208 *
209 END