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