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