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