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