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