1 SUBROUTINE ZTBRFS( UPLO, TRANS, DIAG, N, KD, NRHS, AB, LDAB, B,
2 $ LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * Modified to call ZLACN2 in place of ZLACON, 10 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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
17 COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
18 $ X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZTBRFS 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 ZTBTRS or some other
29 * means before entering this routine. ZTBRFS 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)
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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 array, dimension (2*N)
100 *
101 * RWORK (workspace) DOUBLE PRECISION 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 COMPLEX*16 ONE
113 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
114 * ..
115 * .. Local Scalars ..
116 LOGICAL NOTRAN, NOUNIT, UPPER
117 CHARACTER TRANSN, TRANST
118 INTEGER I, J, K, KASE, NZ
119 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
120 COMPLEX*16 ZDUM
121 * ..
122 * .. Local Arrays ..
123 INTEGER ISAVE( 3 )
124 * ..
125 * .. External Subroutines ..
126 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZLACN2, ZTBMV, ZTBSV
127 * ..
128 * .. Intrinsic Functions ..
129 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
130 * ..
131 * .. External Functions ..
132 LOGICAL LSAME
133 DOUBLE PRECISION DLAMCH
134 EXTERNAL LSAME, DLAMCH
135 * ..
136 * .. Statement Functions ..
137 DOUBLE PRECISION CABS1
138 * ..
139 * .. Statement Function definitions ..
140 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
141 * ..
142 * .. Executable Statements ..
143 *
144 * Test the input parameters.
145 *
146 INFO = 0
147 UPPER = LSAME( UPLO, 'U' )
148 NOTRAN = LSAME( TRANS, 'N' )
149 NOUNIT = LSAME( DIAG, 'N' )
150 *
151 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
152 INFO = -1
153 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
154 $ LSAME( TRANS, 'C' ) ) THEN
155 INFO = -2
156 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
157 INFO = -3
158 ELSE IF( N.LT.0 ) THEN
159 INFO = -4
160 ELSE IF( KD.LT.0 ) THEN
161 INFO = -5
162 ELSE IF( NRHS.LT.0 ) THEN
163 INFO = -6
164 ELSE IF( LDAB.LT.KD+1 ) THEN
165 INFO = -8
166 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
167 INFO = -10
168 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
169 INFO = -12
170 END IF
171 IF( INFO.NE.0 ) THEN
172 CALL XERBLA( 'ZTBRFS', -INFO )
173 RETURN
174 END IF
175 *
176 * Quick return if possible
177 *
178 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
179 DO 10 J = 1, NRHS
180 FERR( J ) = ZERO
181 BERR( J ) = ZERO
182 10 CONTINUE
183 RETURN
184 END IF
185 *
186 IF( NOTRAN ) THEN
187 TRANSN = 'N'
188 TRANST = 'C'
189 ELSE
190 TRANSN = 'C'
191 TRANST = 'N'
192 END IF
193 *
194 * NZ = maximum number of nonzero elements in each row of A, plus 1
195 *
196 NZ = KD + 2
197 EPS = DLAMCH( 'Epsilon' )
198 SAFMIN = DLAMCH( 'Safe minimum' )
199 SAFE1 = NZ*SAFMIN
200 SAFE2 = SAFE1 / EPS
201 *
202 * Do for each right hand side
203 *
204 DO 250 J = 1, NRHS
205 *
206 * Compute residual R = B - op(A) * X,
207 * where op(A) = A, A**T, or A**H, depending on TRANS.
208 *
209 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
210 CALL ZTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
211 CALL ZAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
212 *
213 * Compute componentwise relative backward error from formula
214 *
215 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
216 *
217 * where abs(Z) is the componentwise absolute value of the matrix
218 * or vector Z. If the i-th component of the denominator is less
219 * than SAFE2, then SAFE1 is added to the i-th components of the
220 * numerator and denominator before dividing.
221 *
222 DO 20 I = 1, N
223 RWORK( I ) = CABS1( B( I, J ) )
224 20 CONTINUE
225 *
226 IF( NOTRAN ) THEN
227 *
228 * Compute abs(A)*abs(X) + abs(B).
229 *
230 IF( UPPER ) THEN
231 IF( NOUNIT ) THEN
232 DO 40 K = 1, N
233 XK = CABS1( X( K, J ) )
234 DO 30 I = MAX( 1, K-KD ), K
235 RWORK( I ) = RWORK( I ) +
236 $ CABS1( AB( KD+1+I-K, K ) )*XK
237 30 CONTINUE
238 40 CONTINUE
239 ELSE
240 DO 60 K = 1, N
241 XK = CABS1( X( K, J ) )
242 DO 50 I = MAX( 1, K-KD ), K - 1
243 RWORK( I ) = RWORK( I ) +
244 $ CABS1( AB( KD+1+I-K, K ) )*XK
245 50 CONTINUE
246 RWORK( K ) = RWORK( K ) + XK
247 60 CONTINUE
248 END IF
249 ELSE
250 IF( NOUNIT ) THEN
251 DO 80 K = 1, N
252 XK = CABS1( X( K, J ) )
253 DO 70 I = K, MIN( N, K+KD )
254 RWORK( I ) = RWORK( I ) +
255 $ CABS1( AB( 1+I-K, K ) )*XK
256 70 CONTINUE
257 80 CONTINUE
258 ELSE
259 DO 100 K = 1, N
260 XK = CABS1( X( K, J ) )
261 DO 90 I = K + 1, MIN( N, K+KD )
262 RWORK( I ) = RWORK( I ) +
263 $ CABS1( AB( 1+I-K, K ) )*XK
264 90 CONTINUE
265 RWORK( K ) = RWORK( K ) + XK
266 100 CONTINUE
267 END IF
268 END IF
269 ELSE
270 *
271 * Compute abs(A**H)*abs(X) + abs(B).
272 *
273 IF( UPPER ) THEN
274 IF( NOUNIT ) THEN
275 DO 120 K = 1, N
276 S = ZERO
277 DO 110 I = MAX( 1, K-KD ), K
278 S = S + CABS1( AB( KD+1+I-K, K ) )*
279 $ CABS1( X( I, J ) )
280 110 CONTINUE
281 RWORK( K ) = RWORK( K ) + S
282 120 CONTINUE
283 ELSE
284 DO 140 K = 1, N
285 S = CABS1( X( K, J ) )
286 DO 130 I = MAX( 1, K-KD ), K - 1
287 S = S + CABS1( AB( KD+1+I-K, K ) )*
288 $ CABS1( X( I, J ) )
289 130 CONTINUE
290 RWORK( K ) = RWORK( K ) + S
291 140 CONTINUE
292 END IF
293 ELSE
294 IF( NOUNIT ) THEN
295 DO 160 K = 1, N
296 S = ZERO
297 DO 150 I = K, MIN( N, K+KD )
298 S = S + CABS1( AB( 1+I-K, K ) )*
299 $ CABS1( X( I, J ) )
300 150 CONTINUE
301 RWORK( K ) = RWORK( K ) + S
302 160 CONTINUE
303 ELSE
304 DO 180 K = 1, N
305 S = CABS1( X( K, J ) )
306 DO 170 I = K + 1, MIN( N, K+KD )
307 S = S + CABS1( AB( 1+I-K, K ) )*
308 $ CABS1( X( I, J ) )
309 170 CONTINUE
310 RWORK( K ) = RWORK( K ) + S
311 180 CONTINUE
312 END IF
313 END IF
314 END IF
315 S = ZERO
316 DO 190 I = 1, N
317 IF( RWORK( I ).GT.SAFE2 ) THEN
318 S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
319 ELSE
320 S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
321 $ ( RWORK( I )+SAFE1 ) )
322 END IF
323 190 CONTINUE
324 BERR( J ) = S
325 *
326 * Bound error from formula
327 *
328 * norm(X - XTRUE) / norm(X) .le. FERR =
329 * norm( abs(inv(op(A)))*
330 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
331 *
332 * where
333 * norm(Z) is the magnitude of the largest component of Z
334 * inv(op(A)) is the inverse of op(A)
335 * abs(Z) is the componentwise absolute value of the matrix or
336 * vector Z
337 * NZ is the maximum number of nonzeros in any row of A, plus 1
338 * EPS is machine epsilon
339 *
340 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
341 * is incremented by SAFE1 if the i-th component of
342 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
343 *
344 * Use ZLACN2 to estimate the infinity-norm of the matrix
345 * inv(op(A)) * diag(W),
346 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
347 *
348 DO 200 I = 1, N
349 IF( RWORK( I ).GT.SAFE2 ) THEN
350 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
351 ELSE
352 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
353 $ SAFE1
354 END IF
355 200 CONTINUE
356 *
357 KASE = 0
358 210 CONTINUE
359 CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
360 IF( KASE.NE.0 ) THEN
361 IF( KASE.EQ.1 ) THEN
362 *
363 * Multiply by diag(W)*inv(op(A)**H).
364 *
365 CALL ZTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, WORK,
366 $ 1 )
367 DO 220 I = 1, N
368 WORK( I ) = RWORK( I )*WORK( I )
369 220 CONTINUE
370 ELSE
371 *
372 * Multiply by inv(op(A))*diag(W).
373 *
374 DO 230 I = 1, N
375 WORK( I ) = RWORK( I )*WORK( I )
376 230 CONTINUE
377 CALL ZTBSV( UPLO, TRANSN, DIAG, N, KD, AB, LDAB, WORK,
378 $ 1 )
379 END IF
380 GO TO 210
381 END IF
382 *
383 * Normalize error.
384 *
385 LSTRES = ZERO
386 DO 240 I = 1, N
387 LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
388 240 CONTINUE
389 IF( LSTRES.NE.ZERO )
390 $ FERR( J ) = FERR( J ) / LSTRES
391 *
392 250 CONTINUE
393 *
394 RETURN
395 *
396 * End of ZTBRFS
397 *
398 END
2 $ LDB, X, LDX, FERR, BERR, WORK, RWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * Modified to call ZLACN2 in place of ZLACON, 10 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 DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * )
17 COMPLEX*16 AB( LDAB, * ), B( LDB, * ), WORK( * ),
18 $ X( LDX, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * ZTBRFS 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 ZTBTRS or some other
29 * means before entering this routine. ZTBRFS 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)
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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 array, dimension (2*N)
100 *
101 * RWORK (workspace) DOUBLE PRECISION 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 COMPLEX*16 ONE
113 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
114 * ..
115 * .. Local Scalars ..
116 LOGICAL NOTRAN, NOUNIT, UPPER
117 CHARACTER TRANSN, TRANST
118 INTEGER I, J, K, KASE, NZ
119 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
120 COMPLEX*16 ZDUM
121 * ..
122 * .. Local Arrays ..
123 INTEGER ISAVE( 3 )
124 * ..
125 * .. External Subroutines ..
126 EXTERNAL XERBLA, ZAXPY, ZCOPY, ZLACN2, ZTBMV, ZTBSV
127 * ..
128 * .. Intrinsic Functions ..
129 INTRINSIC ABS, DBLE, DIMAG, MAX, MIN
130 * ..
131 * .. External Functions ..
132 LOGICAL LSAME
133 DOUBLE PRECISION DLAMCH
134 EXTERNAL LSAME, DLAMCH
135 * ..
136 * .. Statement Functions ..
137 DOUBLE PRECISION CABS1
138 * ..
139 * .. Statement Function definitions ..
140 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
141 * ..
142 * .. Executable Statements ..
143 *
144 * Test the input parameters.
145 *
146 INFO = 0
147 UPPER = LSAME( UPLO, 'U' )
148 NOTRAN = LSAME( TRANS, 'N' )
149 NOUNIT = LSAME( DIAG, 'N' )
150 *
151 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
152 INFO = -1
153 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
154 $ LSAME( TRANS, 'C' ) ) THEN
155 INFO = -2
156 ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
157 INFO = -3
158 ELSE IF( N.LT.0 ) THEN
159 INFO = -4
160 ELSE IF( KD.LT.0 ) THEN
161 INFO = -5
162 ELSE IF( NRHS.LT.0 ) THEN
163 INFO = -6
164 ELSE IF( LDAB.LT.KD+1 ) THEN
165 INFO = -8
166 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
167 INFO = -10
168 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
169 INFO = -12
170 END IF
171 IF( INFO.NE.0 ) THEN
172 CALL XERBLA( 'ZTBRFS', -INFO )
173 RETURN
174 END IF
175 *
176 * Quick return if possible
177 *
178 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
179 DO 10 J = 1, NRHS
180 FERR( J ) = ZERO
181 BERR( J ) = ZERO
182 10 CONTINUE
183 RETURN
184 END IF
185 *
186 IF( NOTRAN ) THEN
187 TRANSN = 'N'
188 TRANST = 'C'
189 ELSE
190 TRANSN = 'C'
191 TRANST = 'N'
192 END IF
193 *
194 * NZ = maximum number of nonzero elements in each row of A, plus 1
195 *
196 NZ = KD + 2
197 EPS = DLAMCH( 'Epsilon' )
198 SAFMIN = DLAMCH( 'Safe minimum' )
199 SAFE1 = NZ*SAFMIN
200 SAFE2 = SAFE1 / EPS
201 *
202 * Do for each right hand side
203 *
204 DO 250 J = 1, NRHS
205 *
206 * Compute residual R = B - op(A) * X,
207 * where op(A) = A, A**T, or A**H, depending on TRANS.
208 *
209 CALL ZCOPY( N, X( 1, J ), 1, WORK, 1 )
210 CALL ZTBMV( UPLO, TRANS, DIAG, N, KD, AB, LDAB, WORK, 1 )
211 CALL ZAXPY( N, -ONE, B( 1, J ), 1, WORK, 1 )
212 *
213 * Compute componentwise relative backward error from formula
214 *
215 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
216 *
217 * where abs(Z) is the componentwise absolute value of the matrix
218 * or vector Z. If the i-th component of the denominator is less
219 * than SAFE2, then SAFE1 is added to the i-th components of the
220 * numerator and denominator before dividing.
221 *
222 DO 20 I = 1, N
223 RWORK( I ) = CABS1( B( I, J ) )
224 20 CONTINUE
225 *
226 IF( NOTRAN ) THEN
227 *
228 * Compute abs(A)*abs(X) + abs(B).
229 *
230 IF( UPPER ) THEN
231 IF( NOUNIT ) THEN
232 DO 40 K = 1, N
233 XK = CABS1( X( K, J ) )
234 DO 30 I = MAX( 1, K-KD ), K
235 RWORK( I ) = RWORK( I ) +
236 $ CABS1( AB( KD+1+I-K, K ) )*XK
237 30 CONTINUE
238 40 CONTINUE
239 ELSE
240 DO 60 K = 1, N
241 XK = CABS1( X( K, J ) )
242 DO 50 I = MAX( 1, K-KD ), K - 1
243 RWORK( I ) = RWORK( I ) +
244 $ CABS1( AB( KD+1+I-K, K ) )*XK
245 50 CONTINUE
246 RWORK( K ) = RWORK( K ) + XK
247 60 CONTINUE
248 END IF
249 ELSE
250 IF( NOUNIT ) THEN
251 DO 80 K = 1, N
252 XK = CABS1( X( K, J ) )
253 DO 70 I = K, MIN( N, K+KD )
254 RWORK( I ) = RWORK( I ) +
255 $ CABS1( AB( 1+I-K, K ) )*XK
256 70 CONTINUE
257 80 CONTINUE
258 ELSE
259 DO 100 K = 1, N
260 XK = CABS1( X( K, J ) )
261 DO 90 I = K + 1, MIN( N, K+KD )
262 RWORK( I ) = RWORK( I ) +
263 $ CABS1( AB( 1+I-K, K ) )*XK
264 90 CONTINUE
265 RWORK( K ) = RWORK( K ) + XK
266 100 CONTINUE
267 END IF
268 END IF
269 ELSE
270 *
271 * Compute abs(A**H)*abs(X) + abs(B).
272 *
273 IF( UPPER ) THEN
274 IF( NOUNIT ) THEN
275 DO 120 K = 1, N
276 S = ZERO
277 DO 110 I = MAX( 1, K-KD ), K
278 S = S + CABS1( AB( KD+1+I-K, K ) )*
279 $ CABS1( X( I, J ) )
280 110 CONTINUE
281 RWORK( K ) = RWORK( K ) + S
282 120 CONTINUE
283 ELSE
284 DO 140 K = 1, N
285 S = CABS1( X( K, J ) )
286 DO 130 I = MAX( 1, K-KD ), K - 1
287 S = S + CABS1( AB( KD+1+I-K, K ) )*
288 $ CABS1( X( I, J ) )
289 130 CONTINUE
290 RWORK( K ) = RWORK( K ) + S
291 140 CONTINUE
292 END IF
293 ELSE
294 IF( NOUNIT ) THEN
295 DO 160 K = 1, N
296 S = ZERO
297 DO 150 I = K, MIN( N, K+KD )
298 S = S + CABS1( AB( 1+I-K, K ) )*
299 $ CABS1( X( I, J ) )
300 150 CONTINUE
301 RWORK( K ) = RWORK( K ) + S
302 160 CONTINUE
303 ELSE
304 DO 180 K = 1, N
305 S = CABS1( X( K, J ) )
306 DO 170 I = K + 1, MIN( N, K+KD )
307 S = S + CABS1( AB( 1+I-K, K ) )*
308 $ CABS1( X( I, J ) )
309 170 CONTINUE
310 RWORK( K ) = RWORK( K ) + S
311 180 CONTINUE
312 END IF
313 END IF
314 END IF
315 S = ZERO
316 DO 190 I = 1, N
317 IF( RWORK( I ).GT.SAFE2 ) THEN
318 S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
319 ELSE
320 S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
321 $ ( RWORK( I )+SAFE1 ) )
322 END IF
323 190 CONTINUE
324 BERR( J ) = S
325 *
326 * Bound error from formula
327 *
328 * norm(X - XTRUE) / norm(X) .le. FERR =
329 * norm( abs(inv(op(A)))*
330 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
331 *
332 * where
333 * norm(Z) is the magnitude of the largest component of Z
334 * inv(op(A)) is the inverse of op(A)
335 * abs(Z) is the componentwise absolute value of the matrix or
336 * vector Z
337 * NZ is the maximum number of nonzeros in any row of A, plus 1
338 * EPS is machine epsilon
339 *
340 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
341 * is incremented by SAFE1 if the i-th component of
342 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
343 *
344 * Use ZLACN2 to estimate the infinity-norm of the matrix
345 * inv(op(A)) * diag(W),
346 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
347 *
348 DO 200 I = 1, N
349 IF( RWORK( I ).GT.SAFE2 ) THEN
350 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
351 ELSE
352 RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
353 $ SAFE1
354 END IF
355 200 CONTINUE
356 *
357 KASE = 0
358 210 CONTINUE
359 CALL ZLACN2( N, WORK( N+1 ), WORK, FERR( J ), KASE, ISAVE )
360 IF( KASE.NE.0 ) THEN
361 IF( KASE.EQ.1 ) THEN
362 *
363 * Multiply by diag(W)*inv(op(A)**H).
364 *
365 CALL ZTBSV( UPLO, TRANST, DIAG, N, KD, AB, LDAB, WORK,
366 $ 1 )
367 DO 220 I = 1, N
368 WORK( I ) = RWORK( I )*WORK( I )
369 220 CONTINUE
370 ELSE
371 *
372 * Multiply by inv(op(A))*diag(W).
373 *
374 DO 230 I = 1, N
375 WORK( I ) = RWORK( I )*WORK( I )
376 230 CONTINUE
377 CALL ZTBSV( UPLO, TRANSN, DIAG, N, KD, AB, LDAB, WORK,
378 $ 1 )
379 END IF
380 GO TO 210
381 END IF
382 *
383 * Normalize error.
384 *
385 LSTRES = ZERO
386 DO 240 I = 1, N
387 LSTRES = MAX( LSTRES, CABS1( X( I, J ) ) )
388 240 CONTINUE
389 IF( LSTRES.NE.ZERO )
390 $ FERR( J ) = FERR( J ) / LSTRES
391 *
392 250 CONTINUE
393 *
394 RETURN
395 *
396 * End of ZTBRFS
397 *
398 END