1 SUBROUTINE DTRRFS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB, X,
2 $ LDX, FERR, BERR, WORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER DIAG, TRANS, UPLO
13 INTEGER INFO, LDA, LDB, LDX, N, NRHS
14 * ..
15 * .. Array Arguments ..
16 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
18 $ WORK( * ), X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DTRRFS provides error bounds and backward error estimates for the
25 * solution to a system of linear equations with a triangular
26 * coefficient matrix.
27 *
28 * The solution matrix X must be computed by DTRTRS or some other
29 * means before entering this routine. DTRRFS does not do iterative
30 * refinement because doing so cannot improve the backward error.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * = 'U': A is upper triangular;
37 * = 'L': A is lower triangular.
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the form of the system of equations:
41 * = 'N': A * X = B (No transpose)
42 * = 'T': A**T * X = B (Transpose)
43 * = 'C': A**H * X = B (Conjugate transpose = Transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * = 'N': A is non-unit triangular;
47 * = 'U': A is unit triangular.
48 *
49 * N (input) INTEGER
50 * The order of the matrix A. N >= 0.
51 *
52 * NRHS (input) INTEGER
53 * The number of right hand sides, i.e., the number of columns
54 * of the matrices B and X. NRHS >= 0.
55 *
56 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
57 * The triangular matrix A. If UPLO = 'U', the leading N-by-N
58 * upper triangular part of the array A contains the upper
59 * triangular matrix, and the strictly lower triangular part of
60 * A is not referenced. If UPLO = 'L', the leading N-by-N lower
61 * triangular part of the array A contains the lower triangular
62 * matrix, and the strictly upper triangular part of A is not
63 * referenced. If DIAG = 'U', the diagonal elements of A are
64 * also not referenced and are assumed to be 1.
65 *
66 * LDA (input) INTEGER
67 * The leading dimension of the array A. LDA >= max(1,N).
68 *
69 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
70 * The right hand side matrix B.
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the array B. LDB >= max(1,N).
74 *
75 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
76 * The solution matrix X.
77 *
78 * LDX (input) INTEGER
79 * The leading dimension of the array X. LDX >= max(1,N).
80 *
81 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
82 * The estimated forward error bound for each solution vector
83 * X(j) (the j-th column of the solution matrix X).
84 * If XTRUE is the true solution corresponding to X(j), FERR(j)
85 * is an estimated upper bound for the magnitude of the largest
86 * element in (X(j) - XTRUE) divided by the magnitude of the
87 * largest element in X(j). The estimate is as reliable as
88 * the estimate for RCOND, and is almost always a slight
89 * overestimate of the true error.
90 *
91 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
92 * The componentwise relative backward error of each solution
93 * vector X(j) (i.e., the smallest relative change in
94 * any element of A or B that makes X(j) an exact solution).
95 *
96 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
97 *
98 * IWORK (workspace) INTEGER array, dimension (N)
99 *
100 * INFO (output) INTEGER
101 * = 0: successful exit
102 * < 0: if INFO = -i, the i-th argument had an illegal value
103 *
104 * =====================================================================
105 *
106 * .. Parameters ..
107 DOUBLE PRECISION ZERO
108 PARAMETER ( ZERO = 0.0D+0 )
109 DOUBLE PRECISION ONE
110 PARAMETER ( ONE = 1.0D+0 )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOTRAN, NOUNIT, UPPER
114 CHARACTER TRANST
115 INTEGER I, J, K, KASE, NZ
116 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
117 * ..
118 * .. Local Arrays ..
119 INTEGER ISAVE( 3 )
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL DAXPY, DCOPY, DLACN2, DTRMV, DTRSV, XERBLA
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, MAX
126 * ..
127 * .. External Functions ..
128 LOGICAL LSAME
129 DOUBLE PRECISION DLAMCH
130 EXTERNAL LSAME, DLAMCH
131 * ..
132 * .. Executable Statements ..
133 *
134 * Test the input parameters.
135 *
136 INFO = 0
137 UPPER = LSAME( UPLO, 'U' )
138 NOTRAN = LSAME( TRANS, 'N' )
139 NOUNIT = LSAME( DIAG, 'N' )
140 *
141 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
142 INFO = -1
143 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
144 $ LSAME( TRANS, 'C' ) ) THEN
145 INFO = -2
146 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
147 INFO = -3
148 ELSE IF( N.LT.0 ) THEN
149 INFO = -4
150 ELSE IF( NRHS.LT.0 ) THEN
151 INFO = -5
152 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
153 INFO = -7
154 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
155 INFO = -9
156 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
157 INFO = -11
158 END IF
159 IF( INFO.NE.0 ) THEN
160 CALL XERBLA( 'DTRRFS', -INFO )
161 RETURN
162 END IF
163 *
164 * Quick return if possible
165 *
166 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
167 DO 10 J = 1, NRHS
168 FERR( J ) = ZERO
169 BERR( J ) = ZERO
170 10 CONTINUE
171 RETURN
172 END IF
173 *
174 IF( NOTRAN ) THEN
175 TRANST = 'T'
176 ELSE
177 TRANST = 'N'
178 END IF
179 *
180 * NZ = maximum number of nonzero elements in each row of A, plus 1
181 *
182 NZ = N + 1
183 EPS = DLAMCH( 'Epsilon' )
184 SAFMIN = DLAMCH( 'Safe minimum' )
185 SAFE1 = NZ*SAFMIN
186 SAFE2 = SAFE1 / EPS
187 *
188 * Do for each right hand side
189 *
190 DO 250 J = 1, NRHS
191 *
192 * Compute residual R = B - op(A) * X,
193 * where op(A) = A or A**T, depending on TRANS.
194 *
195 CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
196 CALL DTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ), 1 )
197 CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
198 *
199 * Compute componentwise relative backward error from formula
200 *
201 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
202 *
203 * where abs(Z) is the componentwise absolute value of the matrix
204 * or vector Z. If the i-th component of the denominator is less
205 * than SAFE2, then SAFE1 is added to the i-th components of the
206 * numerator and denominator before dividing.
207 *
208 DO 20 I = 1, N
209 WORK( I ) = ABS( B( I, J ) )
210 20 CONTINUE
211 *
212 IF( NOTRAN ) THEN
213 *
214 * Compute abs(A)*abs(X) + abs(B).
215 *
216 IF( UPPER ) THEN
217 IF( NOUNIT ) THEN
218 DO 40 K = 1, N
219 XK = ABS( X( K, J ) )
220 DO 30 I = 1, K
221 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
222 30 CONTINUE
223 40 CONTINUE
224 ELSE
225 DO 60 K = 1, N
226 XK = ABS( X( K, J ) )
227 DO 50 I = 1, K - 1
228 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
229 50 CONTINUE
230 WORK( K ) = WORK( K ) + XK
231 60 CONTINUE
232 END IF
233 ELSE
234 IF( NOUNIT ) THEN
235 DO 80 K = 1, N
236 XK = ABS( X( K, J ) )
237 DO 70 I = K, N
238 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
239 70 CONTINUE
240 80 CONTINUE
241 ELSE
242 DO 100 K = 1, N
243 XK = ABS( X( K, J ) )
244 DO 90 I = K + 1, N
245 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
246 90 CONTINUE
247 WORK( K ) = WORK( K ) + XK
248 100 CONTINUE
249 END IF
250 END IF
251 ELSE
252 *
253 * Compute abs(A**T)*abs(X) + abs(B).
254 *
255 IF( UPPER ) THEN
256 IF( NOUNIT ) THEN
257 DO 120 K = 1, N
258 S = ZERO
259 DO 110 I = 1, K
260 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
261 110 CONTINUE
262 WORK( K ) = WORK( K ) + S
263 120 CONTINUE
264 ELSE
265 DO 140 K = 1, N
266 S = ABS( X( K, J ) )
267 DO 130 I = 1, K - 1
268 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
269 130 CONTINUE
270 WORK( K ) = WORK( K ) + S
271 140 CONTINUE
272 END IF
273 ELSE
274 IF( NOUNIT ) THEN
275 DO 160 K = 1, N
276 S = ZERO
277 DO 150 I = K, N
278 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
279 150 CONTINUE
280 WORK( K ) = WORK( K ) + S
281 160 CONTINUE
282 ELSE
283 DO 180 K = 1, N
284 S = ABS( X( K, J ) )
285 DO 170 I = K + 1, N
286 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
287 170 CONTINUE
288 WORK( K ) = WORK( K ) + S
289 180 CONTINUE
290 END IF
291 END IF
292 END IF
293 S = ZERO
294 DO 190 I = 1, N
295 IF( WORK( I ).GT.SAFE2 ) THEN
296 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
297 ELSE
298 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
299 $ ( WORK( I )+SAFE1 ) )
300 END IF
301 190 CONTINUE
302 BERR( J ) = S
303 *
304 * Bound error from formula
305 *
306 * norm(X - XTRUE) / norm(X) .le. FERR =
307 * norm( abs(inv(op(A)))*
308 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
309 *
310 * where
311 * norm(Z) is the magnitude of the largest component of Z
312 * inv(op(A)) is the inverse of op(A)
313 * abs(Z) is the componentwise absolute value of the matrix or
314 * vector Z
315 * NZ is the maximum number of nonzeros in any row of A, plus 1
316 * EPS is machine epsilon
317 *
318 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
319 * is incremented by SAFE1 if the i-th component of
320 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
321 *
322 * Use DLACN2 to estimate the infinity-norm of the matrix
323 * inv(op(A)) * diag(W),
324 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
325 *
326 DO 200 I = 1, N
327 IF( WORK( I ).GT.SAFE2 ) THEN
328 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
329 ELSE
330 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
331 END IF
332 200 CONTINUE
333 *
334 KASE = 0
335 210 CONTINUE
336 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
337 $ KASE, ISAVE )
338 IF( KASE.NE.0 ) THEN
339 IF( KASE.EQ.1 ) THEN
340 *
341 * Multiply by diag(W)*inv(op(A)**T).
342 *
343 CALL DTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK( N+1 ),
344 $ 1 )
345 DO 220 I = 1, N
346 WORK( N+I ) = WORK( I )*WORK( N+I )
347 220 CONTINUE
348 ELSE
349 *
350 * Multiply by inv(op(A))*diag(W).
351 *
352 DO 230 I = 1, N
353 WORK( N+I ) = WORK( I )*WORK( N+I )
354 230 CONTINUE
355 CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ),
356 $ 1 )
357 END IF
358 GO TO 210
359 END IF
360 *
361 * Normalize error.
362 *
363 LSTRES = ZERO
364 DO 240 I = 1, N
365 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
366 240 CONTINUE
367 IF( LSTRES.NE.ZERO )
368 $ FERR( J ) = FERR( J ) / LSTRES
369 *
370 250 CONTINUE
371 *
372 RETURN
373 *
374 * End of DTRRFS
375 *
376 END
2 $ LDX, FERR, BERR, WORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
10 *
11 * .. Scalar Arguments ..
12 CHARACTER DIAG, TRANS, UPLO
13 INTEGER INFO, LDA, LDB, LDX, N, NRHS
14 * ..
15 * .. Array Arguments ..
16 INTEGER IWORK( * )
17 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ),
18 $ WORK( * ), X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DTRRFS provides error bounds and backward error estimates for the
25 * solution to a system of linear equations with a triangular
26 * coefficient matrix.
27 *
28 * The solution matrix X must be computed by DTRTRS or some other
29 * means before entering this routine. DTRRFS does not do iterative
30 * refinement because doing so cannot improve the backward error.
31 *
32 * Arguments
33 * =========
34 *
35 * UPLO (input) CHARACTER*1
36 * = 'U': A is upper triangular;
37 * = 'L': A is lower triangular.
38 *
39 * TRANS (input) CHARACTER*1
40 * Specifies the form of the system of equations:
41 * = 'N': A * X = B (No transpose)
42 * = 'T': A**T * X = B (Transpose)
43 * = 'C': A**H * X = B (Conjugate transpose = Transpose)
44 *
45 * DIAG (input) CHARACTER*1
46 * = 'N': A is non-unit triangular;
47 * = 'U': A is unit triangular.
48 *
49 * N (input) INTEGER
50 * The order of the matrix A. N >= 0.
51 *
52 * NRHS (input) INTEGER
53 * The number of right hand sides, i.e., the number of columns
54 * of the matrices B and X. NRHS >= 0.
55 *
56 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
57 * The triangular matrix A. If UPLO = 'U', the leading N-by-N
58 * upper triangular part of the array A contains the upper
59 * triangular matrix, and the strictly lower triangular part of
60 * A is not referenced. If UPLO = 'L', the leading N-by-N lower
61 * triangular part of the array A contains the lower triangular
62 * matrix, and the strictly upper triangular part of A is not
63 * referenced. If DIAG = 'U', the diagonal elements of A are
64 * also not referenced and are assumed to be 1.
65 *
66 * LDA (input) INTEGER
67 * The leading dimension of the array A. LDA >= max(1,N).
68 *
69 * B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
70 * The right hand side matrix B.
71 *
72 * LDB (input) INTEGER
73 * The leading dimension of the array B. LDB >= max(1,N).
74 *
75 * X (input) DOUBLE PRECISION array, dimension (LDX,NRHS)
76 * The solution matrix X.
77 *
78 * LDX (input) INTEGER
79 * The leading dimension of the array X. LDX >= max(1,N).
80 *
81 * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
82 * The estimated forward error bound for each solution vector
83 * X(j) (the j-th column of the solution matrix X).
84 * If XTRUE is the true solution corresponding to X(j), FERR(j)
85 * is an estimated upper bound for the magnitude of the largest
86 * element in (X(j) - XTRUE) divided by the magnitude of the
87 * largest element in X(j). The estimate is as reliable as
88 * the estimate for RCOND, and is almost always a slight
89 * overestimate of the true error.
90 *
91 * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
92 * The componentwise relative backward error of each solution
93 * vector X(j) (i.e., the smallest relative change in
94 * any element of A or B that makes X(j) an exact solution).
95 *
96 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
97 *
98 * IWORK (workspace) INTEGER array, dimension (N)
99 *
100 * INFO (output) INTEGER
101 * = 0: successful exit
102 * < 0: if INFO = -i, the i-th argument had an illegal value
103 *
104 * =====================================================================
105 *
106 * .. Parameters ..
107 DOUBLE PRECISION ZERO
108 PARAMETER ( ZERO = 0.0D+0 )
109 DOUBLE PRECISION ONE
110 PARAMETER ( ONE = 1.0D+0 )
111 * ..
112 * .. Local Scalars ..
113 LOGICAL NOTRAN, NOUNIT, UPPER
114 CHARACTER TRANST
115 INTEGER I, J, K, KASE, NZ
116 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
117 * ..
118 * .. Local Arrays ..
119 INTEGER ISAVE( 3 )
120 * ..
121 * .. External Subroutines ..
122 EXTERNAL DAXPY, DCOPY, DLACN2, DTRMV, DTRSV, XERBLA
123 * ..
124 * .. Intrinsic Functions ..
125 INTRINSIC ABS, MAX
126 * ..
127 * .. External Functions ..
128 LOGICAL LSAME
129 DOUBLE PRECISION DLAMCH
130 EXTERNAL LSAME, DLAMCH
131 * ..
132 * .. Executable Statements ..
133 *
134 * Test the input parameters.
135 *
136 INFO = 0
137 UPPER = LSAME( UPLO, 'U' )
138 NOTRAN = LSAME( TRANS, 'N' )
139 NOUNIT = LSAME( DIAG, 'N' )
140 *
141 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
142 INFO = -1
143 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
144 $ LSAME( TRANS, 'C' ) ) THEN
145 INFO = -2
146 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
147 INFO = -3
148 ELSE IF( N.LT.0 ) THEN
149 INFO = -4
150 ELSE IF( NRHS.LT.0 ) THEN
151 INFO = -5
152 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
153 INFO = -7
154 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
155 INFO = -9
156 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
157 INFO = -11
158 END IF
159 IF( INFO.NE.0 ) THEN
160 CALL XERBLA( 'DTRRFS', -INFO )
161 RETURN
162 END IF
163 *
164 * Quick return if possible
165 *
166 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
167 DO 10 J = 1, NRHS
168 FERR( J ) = ZERO
169 BERR( J ) = ZERO
170 10 CONTINUE
171 RETURN
172 END IF
173 *
174 IF( NOTRAN ) THEN
175 TRANST = 'T'
176 ELSE
177 TRANST = 'N'
178 END IF
179 *
180 * NZ = maximum number of nonzero elements in each row of A, plus 1
181 *
182 NZ = N + 1
183 EPS = DLAMCH( 'Epsilon' )
184 SAFMIN = DLAMCH( 'Safe minimum' )
185 SAFE1 = NZ*SAFMIN
186 SAFE2 = SAFE1 / EPS
187 *
188 * Do for each right hand side
189 *
190 DO 250 J = 1, NRHS
191 *
192 * Compute residual R = B - op(A) * X,
193 * where op(A) = A or A**T, depending on TRANS.
194 *
195 CALL DCOPY( N, X( 1, J ), 1, WORK( N+1 ), 1 )
196 CALL DTRMV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ), 1 )
197 CALL DAXPY( N, -ONE, B( 1, J ), 1, WORK( N+1 ), 1 )
198 *
199 * Compute componentwise relative backward error from formula
200 *
201 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
202 *
203 * where abs(Z) is the componentwise absolute value of the matrix
204 * or vector Z. If the i-th component of the denominator is less
205 * than SAFE2, then SAFE1 is added to the i-th components of the
206 * numerator and denominator before dividing.
207 *
208 DO 20 I = 1, N
209 WORK( I ) = ABS( B( I, J ) )
210 20 CONTINUE
211 *
212 IF( NOTRAN ) THEN
213 *
214 * Compute abs(A)*abs(X) + abs(B).
215 *
216 IF( UPPER ) THEN
217 IF( NOUNIT ) THEN
218 DO 40 K = 1, N
219 XK = ABS( X( K, J ) )
220 DO 30 I = 1, K
221 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
222 30 CONTINUE
223 40 CONTINUE
224 ELSE
225 DO 60 K = 1, N
226 XK = ABS( X( K, J ) )
227 DO 50 I = 1, K - 1
228 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
229 50 CONTINUE
230 WORK( K ) = WORK( K ) + XK
231 60 CONTINUE
232 END IF
233 ELSE
234 IF( NOUNIT ) THEN
235 DO 80 K = 1, N
236 XK = ABS( X( K, J ) )
237 DO 70 I = K, N
238 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
239 70 CONTINUE
240 80 CONTINUE
241 ELSE
242 DO 100 K = 1, N
243 XK = ABS( X( K, J ) )
244 DO 90 I = K + 1, N
245 WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
246 90 CONTINUE
247 WORK( K ) = WORK( K ) + XK
248 100 CONTINUE
249 END IF
250 END IF
251 ELSE
252 *
253 * Compute abs(A**T)*abs(X) + abs(B).
254 *
255 IF( UPPER ) THEN
256 IF( NOUNIT ) THEN
257 DO 120 K = 1, N
258 S = ZERO
259 DO 110 I = 1, K
260 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
261 110 CONTINUE
262 WORK( K ) = WORK( K ) + S
263 120 CONTINUE
264 ELSE
265 DO 140 K = 1, N
266 S = ABS( X( K, J ) )
267 DO 130 I = 1, K - 1
268 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
269 130 CONTINUE
270 WORK( K ) = WORK( K ) + S
271 140 CONTINUE
272 END IF
273 ELSE
274 IF( NOUNIT ) THEN
275 DO 160 K = 1, N
276 S = ZERO
277 DO 150 I = K, N
278 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
279 150 CONTINUE
280 WORK( K ) = WORK( K ) + S
281 160 CONTINUE
282 ELSE
283 DO 180 K = 1, N
284 S = ABS( X( K, J ) )
285 DO 170 I = K + 1, N
286 S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
287 170 CONTINUE
288 WORK( K ) = WORK( K ) + S
289 180 CONTINUE
290 END IF
291 END IF
292 END IF
293 S = ZERO
294 DO 190 I = 1, N
295 IF( WORK( I ).GT.SAFE2 ) THEN
296 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
297 ELSE
298 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
299 $ ( WORK( I )+SAFE1 ) )
300 END IF
301 190 CONTINUE
302 BERR( J ) = S
303 *
304 * Bound error from formula
305 *
306 * norm(X - XTRUE) / norm(X) .le. FERR =
307 * norm( abs(inv(op(A)))*
308 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
309 *
310 * where
311 * norm(Z) is the magnitude of the largest component of Z
312 * inv(op(A)) is the inverse of op(A)
313 * abs(Z) is the componentwise absolute value of the matrix or
314 * vector Z
315 * NZ is the maximum number of nonzeros in any row of A, plus 1
316 * EPS is machine epsilon
317 *
318 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
319 * is incremented by SAFE1 if the i-th component of
320 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
321 *
322 * Use DLACN2 to estimate the infinity-norm of the matrix
323 * inv(op(A)) * diag(W),
324 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
325 *
326 DO 200 I = 1, N
327 IF( WORK( I ).GT.SAFE2 ) THEN
328 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
329 ELSE
330 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
331 END IF
332 200 CONTINUE
333 *
334 KASE = 0
335 210 CONTINUE
336 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
337 $ KASE, ISAVE )
338 IF( KASE.NE.0 ) THEN
339 IF( KASE.EQ.1 ) THEN
340 *
341 * Multiply by diag(W)*inv(op(A)**T).
342 *
343 CALL DTRSV( UPLO, TRANST, DIAG, N, A, LDA, WORK( N+1 ),
344 $ 1 )
345 DO 220 I = 1, N
346 WORK( N+I ) = WORK( I )*WORK( N+I )
347 220 CONTINUE
348 ELSE
349 *
350 * Multiply by inv(op(A))*diag(W).
351 *
352 DO 230 I = 1, N
353 WORK( N+I ) = WORK( I )*WORK( N+I )
354 230 CONTINUE
355 CALL DTRSV( UPLO, TRANS, DIAG, N, A, LDA, WORK( N+1 ),
356 $ 1 )
357 END IF
358 GO TO 210
359 END IF
360 *
361 * Normalize error.
362 *
363 LSTRES = ZERO
364 DO 240 I = 1, N
365 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
366 240 CONTINUE
367 IF( LSTRES.NE.ZERO )
368 $ FERR( J ) = FERR( J ) / LSTRES
369 *
370 250 CONTINUE
371 *
372 RETURN
373 *
374 * End of DTRRFS
375 *
376 END