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