1 SUBROUTINE DGTRFS( TRANS, N, NRHS, DL, D, DU, DLF, DF, DUF, DU2,
2 $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK,
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 DLACN2 in place of DLACON, 5 Feb 03, SJH.
11 *
12 * .. Scalar Arguments ..
13 CHARACTER TRANS
14 INTEGER INFO, LDB, LDX, N, NRHS
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * ), IWORK( * )
18 DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ),
19 $ DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
20 $ FERR( * ), WORK( * ), X( LDX, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DGTRFS improves the computed solution to a system of linear
27 * equations when the coefficient matrix is tridiagonal, 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 = Transpose)
38 *
39 * N (input) INTEGER
40 * The order of the matrix A. N >= 0.
41 *
42 * NRHS (input) INTEGER
43 * The number of right hand sides, i.e., the number of columns
44 * of the matrix B. NRHS >= 0.
45 *
46 * DL (input) DOUBLE PRECISION array, dimension (N-1)
47 * The (n-1) subdiagonal elements of A.
48 *
49 * D (input) DOUBLE PRECISION array, dimension (N)
50 * The diagonal elements of A.
51 *
52 * DU (input) DOUBLE PRECISION array, dimension (N-1)
53 * The (n-1) superdiagonal elements of A.
54 *
55 * DLF (input) DOUBLE PRECISION array, dimension (N-1)
56 * The (n-1) multipliers that define the matrix L from the
57 * LU factorization of A as computed by DGTTRF.
58 *
59 * DF (input) DOUBLE PRECISION array, dimension (N)
60 * The n diagonal elements of the upper triangular matrix U from
61 * the LU factorization of A.
62 *
63 * DUF (input) DOUBLE PRECISION array, dimension (N-1)
64 * The (n-1) elements of the first superdiagonal of U.
65 *
66 * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
67 * The (n-2) elements of the second superdiagonal of U.
68 *
69 * IPIV (input) INTEGER array, dimension (N)
70 * The pivot indices; for 1 <= i <= n, row i of the matrix was
71 * interchanged with row IPIV(i). IPIV(i) will always be either
72 * i or i+1; IPIV(i) = i indicates a row interchange was not
73 * required.
74 *
75 * B (input) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDX,NRHS)
82 * On entry, the solution matrix X, as computed by DGTTRS.
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) DOUBLE PRECISION array, dimension (3*N)
104 *
105 * IWORK (workspace) INTEGER 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, ONE
122 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
123 DOUBLE PRECISION TWO
124 PARAMETER ( TWO = 2.0D+0 )
125 DOUBLE PRECISION THREE
126 PARAMETER ( THREE = 3.0D+0 )
127 * ..
128 * .. Local Scalars ..
129 LOGICAL NOTRAN
130 CHARACTER TRANSN, TRANST
131 INTEGER COUNT, I, J, KASE, NZ
132 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
133 * ..
134 * .. Local Arrays ..
135 INTEGER ISAVE( 3 )
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DAXPY, DCOPY, DGTTRS, DLACN2, DLAGTM, XERBLA
139 * ..
140 * .. Intrinsic Functions ..
141 INTRINSIC ABS, MAX
142 * ..
143 * .. External Functions ..
144 LOGICAL LSAME
145 DOUBLE PRECISION DLAMCH
146 EXTERNAL LSAME, DLAMCH
147 * ..
148 * .. Executable Statements ..
149 *
150 * Test the input parameters.
151 *
152 INFO = 0
153 NOTRAN = LSAME( TRANS, 'N' )
154 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
155 $ LSAME( TRANS, 'C' ) ) THEN
156 INFO = -1
157 ELSE IF( N.LT.0 ) THEN
158 INFO = -2
159 ELSE IF( NRHS.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
162 INFO = -13
163 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
164 INFO = -15
165 END IF
166 IF( INFO.NE.0 ) THEN
167 CALL XERBLA( 'DGTRFS', -INFO )
168 RETURN
169 END IF
170 *
171 * Quick return if possible
172 *
173 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
174 DO 10 J = 1, NRHS
175 FERR( J ) = ZERO
176 BERR( J ) = ZERO
177 10 CONTINUE
178 RETURN
179 END IF
180 *
181 IF( NOTRAN ) THEN
182 TRANSN = 'N'
183 TRANST = 'T'
184 ELSE
185 TRANSN = 'T'
186 TRANST = 'N'
187 END IF
188 *
189 * NZ = maximum number of nonzero elements in each row of A, plus 1
190 *
191 NZ = 4
192 EPS = DLAMCH( 'Epsilon' )
193 SAFMIN = DLAMCH( 'Safe minimum' )
194 SAFE1 = NZ*SAFMIN
195 SAFE2 = SAFE1 / EPS
196 *
197 * Do for each right hand side
198 *
199 DO 110 J = 1, NRHS
200 *
201 COUNT = 1
202 LSTRES = THREE
203 20 CONTINUE
204 *
205 * Loop until stopping criterion is satisfied.
206 *
207 * Compute residual R = B - op(A) * X,
208 * where op(A) = A, A**T, or A**H, depending on TRANS.
209 *
210 CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
211 CALL DLAGTM( TRANS, N, 1, -ONE, DL, D, DU, X( 1, J ), LDX, ONE,
212 $ WORK( N+1 ), N )
213 *
214 * Compute abs(op(A))*abs(x) + abs(b) for use in the backward
215 * error bound.
216 *
217 IF( NOTRAN ) THEN
218 IF( N.EQ.1 ) THEN
219 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) )
220 ELSE
221 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) ) +
222 $ ABS( DU( 1 )*X( 2, J ) )
223 DO 30 I = 2, N - 1
224 WORK( I ) = ABS( B( I, J ) ) +
225 $ ABS( DL( I-1 )*X( I-1, J ) ) +
226 $ ABS( D( I )*X( I, J ) ) +
227 $ ABS( DU( I )*X( I+1, J ) )
228 30 CONTINUE
229 WORK( N ) = ABS( B( N, J ) ) +
230 $ ABS( DL( N-1 )*X( N-1, J ) ) +
231 $ ABS( D( N )*X( N, J ) )
232 END IF
233 ELSE
234 IF( N.EQ.1 ) THEN
235 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) )
236 ELSE
237 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) ) +
238 $ ABS( DL( 1 )*X( 2, J ) )
239 DO 40 I = 2, N - 1
240 WORK( I ) = ABS( B( I, J ) ) +
241 $ ABS( DU( I-1 )*X( I-1, J ) ) +
242 $ ABS( D( I )*X( I, J ) ) +
243 $ ABS( DL( I )*X( I+1, J ) )
244 40 CONTINUE
245 WORK( N ) = ABS( B( N, J ) ) +
246 $ ABS( DU( N-1 )*X( N-1, J ) ) +
247 $ ABS( D( N )*X( N, J ) )
248 END IF
249 END IF
250 *
251 * Compute componentwise relative backward error from formula
252 *
253 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
254 *
255 * where abs(Z) is the componentwise absolute value of the matrix
256 * or vector Z. If the i-th component of the denominator is less
257 * than SAFE2, then SAFE1 is added to the i-th components of the
258 * numerator and denominator before dividing.
259 *
260 S = ZERO
261 DO 50 I = 1, N
262 IF( WORK( I ).GT.SAFE2 ) THEN
263 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
264 ELSE
265 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
266 $ ( WORK( I )+SAFE1 ) )
267 END IF
268 50 CONTINUE
269 BERR( J ) = S
270 *
271 * Test stopping criterion. Continue iterating if
272 * 1) The residual BERR(J) is larger than machine epsilon, and
273 * 2) BERR(J) decreased by at least a factor of 2 during the
274 * last iteration, and
275 * 3) At most ITMAX iterations tried.
276 *
277 IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
278 $ COUNT.LE.ITMAX ) THEN
279 *
280 * Update solution and try again.
281 *
282 CALL DGTTRS( TRANS, N, 1, DLF, DF, DUF, DU2, IPIV,
283 $ WORK( N+1 ), N, INFO )
284 CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
285 LSTRES = BERR( J )
286 COUNT = COUNT + 1
287 GO TO 20
288 END IF
289 *
290 * Bound error from formula
291 *
292 * norm(X - XTRUE) / norm(X) .le. FERR =
293 * norm( abs(inv(op(A)))*
294 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
295 *
296 * where
297 * norm(Z) is the magnitude of the largest component of Z
298 * inv(op(A)) is the inverse of op(A)
299 * abs(Z) is the componentwise absolute value of the matrix or
300 * vector Z
301 * NZ is the maximum number of nonzeros in any row of A, plus 1
302 * EPS is machine epsilon
303 *
304 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
305 * is incremented by SAFE1 if the i-th component of
306 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
307 *
308 * Use DLACN2 to estimate the infinity-norm of the matrix
309 * inv(op(A)) * diag(W),
310 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
311 *
312 DO 60 I = 1, N
313 IF( WORK( I ).GT.SAFE2 ) THEN
314 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
315 ELSE
316 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
317 END IF
318 60 CONTINUE
319 *
320 KASE = 0
321 70 CONTINUE
322 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
323 $ KASE, ISAVE )
324 IF( KASE.NE.0 ) THEN
325 IF( KASE.EQ.1 ) THEN
326 *
327 * Multiply by diag(W)*inv(op(A)**T).
328 *
329 CALL DGTTRS( TRANST, N, 1, DLF, DF, DUF, DU2, IPIV,
330 $ WORK( N+1 ), N, INFO )
331 DO 80 I = 1, N
332 WORK( N+I ) = WORK( I )*WORK( N+I )
333 80 CONTINUE
334 ELSE
335 *
336 * Multiply by inv(op(A))*diag(W).
337 *
338 DO 90 I = 1, N
339 WORK( N+I ) = WORK( I )*WORK( N+I )
340 90 CONTINUE
341 CALL DGTTRS( TRANSN, N, 1, DLF, DF, DUF, DU2, IPIV,
342 $ WORK( N+1 ), N, INFO )
343 END IF
344 GO TO 70
345 END IF
346 *
347 * Normalize error.
348 *
349 LSTRES = ZERO
350 DO 100 I = 1, N
351 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
352 100 CONTINUE
353 IF( LSTRES.NE.ZERO )
354 $ FERR( J ) = FERR( J ) / LSTRES
355 *
356 110 CONTINUE
357 *
358 RETURN
359 *
360 * End of DGTRFS
361 *
362 END
2 $ IPIV, B, LDB, X, LDX, FERR, BERR, WORK, IWORK,
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 DLACN2 in place of DLACON, 5 Feb 03, SJH.
11 *
12 * .. Scalar Arguments ..
13 CHARACTER TRANS
14 INTEGER INFO, LDB, LDX, N, NRHS
15 * ..
16 * .. Array Arguments ..
17 INTEGER IPIV( * ), IWORK( * )
18 DOUBLE PRECISION B( LDB, * ), BERR( * ), D( * ), DF( * ),
19 $ DL( * ), DLF( * ), DU( * ), DU2( * ), DUF( * ),
20 $ FERR( * ), WORK( * ), X( LDX, * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * DGTRFS improves the computed solution to a system of linear
27 * equations when the coefficient matrix is tridiagonal, 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 = Transpose)
38 *
39 * N (input) INTEGER
40 * The order of the matrix A. N >= 0.
41 *
42 * NRHS (input) INTEGER
43 * The number of right hand sides, i.e., the number of columns
44 * of the matrix B. NRHS >= 0.
45 *
46 * DL (input) DOUBLE PRECISION array, dimension (N-1)
47 * The (n-1) subdiagonal elements of A.
48 *
49 * D (input) DOUBLE PRECISION array, dimension (N)
50 * The diagonal elements of A.
51 *
52 * DU (input) DOUBLE PRECISION array, dimension (N-1)
53 * The (n-1) superdiagonal elements of A.
54 *
55 * DLF (input) DOUBLE PRECISION array, dimension (N-1)
56 * The (n-1) multipliers that define the matrix L from the
57 * LU factorization of A as computed by DGTTRF.
58 *
59 * DF (input) DOUBLE PRECISION array, dimension (N)
60 * The n diagonal elements of the upper triangular matrix U from
61 * the LU factorization of A.
62 *
63 * DUF (input) DOUBLE PRECISION array, dimension (N-1)
64 * The (n-1) elements of the first superdiagonal of U.
65 *
66 * DU2 (input) DOUBLE PRECISION array, dimension (N-2)
67 * The (n-2) elements of the second superdiagonal of U.
68 *
69 * IPIV (input) INTEGER array, dimension (N)
70 * The pivot indices; for 1 <= i <= n, row i of the matrix was
71 * interchanged with row IPIV(i). IPIV(i) will always be either
72 * i or i+1; IPIV(i) = i indicates a row interchange was not
73 * required.
74 *
75 * B (input) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDX,NRHS)
82 * On entry, the solution matrix X, as computed by DGTTRS.
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) DOUBLE PRECISION array, dimension (3*N)
104 *
105 * IWORK (workspace) INTEGER 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, ONE
122 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
123 DOUBLE PRECISION TWO
124 PARAMETER ( TWO = 2.0D+0 )
125 DOUBLE PRECISION THREE
126 PARAMETER ( THREE = 3.0D+0 )
127 * ..
128 * .. Local Scalars ..
129 LOGICAL NOTRAN
130 CHARACTER TRANSN, TRANST
131 INTEGER COUNT, I, J, KASE, NZ
132 DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
133 * ..
134 * .. Local Arrays ..
135 INTEGER ISAVE( 3 )
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DAXPY, DCOPY, DGTTRS, DLACN2, DLAGTM, XERBLA
139 * ..
140 * .. Intrinsic Functions ..
141 INTRINSIC ABS, MAX
142 * ..
143 * .. External Functions ..
144 LOGICAL LSAME
145 DOUBLE PRECISION DLAMCH
146 EXTERNAL LSAME, DLAMCH
147 * ..
148 * .. Executable Statements ..
149 *
150 * Test the input parameters.
151 *
152 INFO = 0
153 NOTRAN = LSAME( TRANS, 'N' )
154 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
155 $ LSAME( TRANS, 'C' ) ) THEN
156 INFO = -1
157 ELSE IF( N.LT.0 ) THEN
158 INFO = -2
159 ELSE IF( NRHS.LT.0 ) THEN
160 INFO = -3
161 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
162 INFO = -13
163 ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
164 INFO = -15
165 END IF
166 IF( INFO.NE.0 ) THEN
167 CALL XERBLA( 'DGTRFS', -INFO )
168 RETURN
169 END IF
170 *
171 * Quick return if possible
172 *
173 IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
174 DO 10 J = 1, NRHS
175 FERR( J ) = ZERO
176 BERR( J ) = ZERO
177 10 CONTINUE
178 RETURN
179 END IF
180 *
181 IF( NOTRAN ) THEN
182 TRANSN = 'N'
183 TRANST = 'T'
184 ELSE
185 TRANSN = 'T'
186 TRANST = 'N'
187 END IF
188 *
189 * NZ = maximum number of nonzero elements in each row of A, plus 1
190 *
191 NZ = 4
192 EPS = DLAMCH( 'Epsilon' )
193 SAFMIN = DLAMCH( 'Safe minimum' )
194 SAFE1 = NZ*SAFMIN
195 SAFE2 = SAFE1 / EPS
196 *
197 * Do for each right hand side
198 *
199 DO 110 J = 1, NRHS
200 *
201 COUNT = 1
202 LSTRES = THREE
203 20 CONTINUE
204 *
205 * Loop until stopping criterion is satisfied.
206 *
207 * Compute residual R = B - op(A) * X,
208 * where op(A) = A, A**T, or A**H, depending on TRANS.
209 *
210 CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
211 CALL DLAGTM( TRANS, N, 1, -ONE, DL, D, DU, X( 1, J ), LDX, ONE,
212 $ WORK( N+1 ), N )
213 *
214 * Compute abs(op(A))*abs(x) + abs(b) for use in the backward
215 * error bound.
216 *
217 IF( NOTRAN ) THEN
218 IF( N.EQ.1 ) THEN
219 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) )
220 ELSE
221 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) ) +
222 $ ABS( DU( 1 )*X( 2, J ) )
223 DO 30 I = 2, N - 1
224 WORK( I ) = ABS( B( I, J ) ) +
225 $ ABS( DL( I-1 )*X( I-1, J ) ) +
226 $ ABS( D( I )*X( I, J ) ) +
227 $ ABS( DU( I )*X( I+1, J ) )
228 30 CONTINUE
229 WORK( N ) = ABS( B( N, J ) ) +
230 $ ABS( DL( N-1 )*X( N-1, J ) ) +
231 $ ABS( D( N )*X( N, J ) )
232 END IF
233 ELSE
234 IF( N.EQ.1 ) THEN
235 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) )
236 ELSE
237 WORK( 1 ) = ABS( B( 1, J ) ) + ABS( D( 1 )*X( 1, J ) ) +
238 $ ABS( DL( 1 )*X( 2, J ) )
239 DO 40 I = 2, N - 1
240 WORK( I ) = ABS( B( I, J ) ) +
241 $ ABS( DU( I-1 )*X( I-1, J ) ) +
242 $ ABS( D( I )*X( I, J ) ) +
243 $ ABS( DL( I )*X( I+1, J ) )
244 40 CONTINUE
245 WORK( N ) = ABS( B( N, J ) ) +
246 $ ABS( DU( N-1 )*X( N-1, J ) ) +
247 $ ABS( D( N )*X( N, J ) )
248 END IF
249 END IF
250 *
251 * Compute componentwise relative backward error from formula
252 *
253 * max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
254 *
255 * where abs(Z) is the componentwise absolute value of the matrix
256 * or vector Z. If the i-th component of the denominator is less
257 * than SAFE2, then SAFE1 is added to the i-th components of the
258 * numerator and denominator before dividing.
259 *
260 S = ZERO
261 DO 50 I = 1, N
262 IF( WORK( I ).GT.SAFE2 ) THEN
263 S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
264 ELSE
265 S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
266 $ ( WORK( I )+SAFE1 ) )
267 END IF
268 50 CONTINUE
269 BERR( J ) = S
270 *
271 * Test stopping criterion. Continue iterating if
272 * 1) The residual BERR(J) is larger than machine epsilon, and
273 * 2) BERR(J) decreased by at least a factor of 2 during the
274 * last iteration, and
275 * 3) At most ITMAX iterations tried.
276 *
277 IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
278 $ COUNT.LE.ITMAX ) THEN
279 *
280 * Update solution and try again.
281 *
282 CALL DGTTRS( TRANS, N, 1, DLF, DF, DUF, DU2, IPIV,
283 $ WORK( N+1 ), N, INFO )
284 CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
285 LSTRES = BERR( J )
286 COUNT = COUNT + 1
287 GO TO 20
288 END IF
289 *
290 * Bound error from formula
291 *
292 * norm(X - XTRUE) / norm(X) .le. FERR =
293 * norm( abs(inv(op(A)))*
294 * ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
295 *
296 * where
297 * norm(Z) is the magnitude of the largest component of Z
298 * inv(op(A)) is the inverse of op(A)
299 * abs(Z) is the componentwise absolute value of the matrix or
300 * vector Z
301 * NZ is the maximum number of nonzeros in any row of A, plus 1
302 * EPS is machine epsilon
303 *
304 * The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
305 * is incremented by SAFE1 if the i-th component of
306 * abs(op(A))*abs(X) + abs(B) is less than SAFE2.
307 *
308 * Use DLACN2 to estimate the infinity-norm of the matrix
309 * inv(op(A)) * diag(W),
310 * where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
311 *
312 DO 60 I = 1, N
313 IF( WORK( I ).GT.SAFE2 ) THEN
314 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
315 ELSE
316 WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
317 END IF
318 60 CONTINUE
319 *
320 KASE = 0
321 70 CONTINUE
322 CALL DLACN2( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
323 $ KASE, ISAVE )
324 IF( KASE.NE.0 ) THEN
325 IF( KASE.EQ.1 ) THEN
326 *
327 * Multiply by diag(W)*inv(op(A)**T).
328 *
329 CALL DGTTRS( TRANST, N, 1, DLF, DF, DUF, DU2, IPIV,
330 $ WORK( N+1 ), N, INFO )
331 DO 80 I = 1, N
332 WORK( N+I ) = WORK( I )*WORK( N+I )
333 80 CONTINUE
334 ELSE
335 *
336 * Multiply by inv(op(A))*diag(W).
337 *
338 DO 90 I = 1, N
339 WORK( N+I ) = WORK( I )*WORK( N+I )
340 90 CONTINUE
341 CALL DGTTRS( TRANSN, N, 1, DLF, DF, DUF, DU2, IPIV,
342 $ WORK( N+1 ), N, INFO )
343 END IF
344 GO TO 70
345 END IF
346 *
347 * Normalize error.
348 *
349 LSTRES = ZERO
350 DO 100 I = 1, N
351 LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
352 100 CONTINUE
353 IF( LSTRES.NE.ZERO )
354 $ FERR( J ) = FERR( J ) / LSTRES
355 *
356 110 CONTINUE
357 *
358 RETURN
359 *
360 * End of DGTRFS
361 *
362 END