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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztprfs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack