1: SUBROUTINE ZPTRFS( UPLO, N, NRHS, D, E, DF, EF, B, LDB, X, LDX,
2: $ FERR, BERR, WORK, RWORK, 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: * .. Scalar Arguments ..
10: CHARACTER UPLO
11: INTEGER INFO, LDB, LDX, N, NRHS
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION BERR( * ), D( * ), DF( * ), FERR( * ),
15: $ RWORK( * )
16: COMPLEX*16 B( LDB, * ), E( * ), EF( * ), WORK( * ),
17: $ X( LDX, * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * ZPTRFS improves the computed solution to a system of linear
24: * equations when the coefficient matrix is Hermitian positive definite
25: * and tridiagonal, and provides error bounds and backward error
26: * estimates for the solution.
27: *
28: * Arguments
29: * =========
30: *
31: * UPLO (input) CHARACTER*1
32: * Specifies whether the superdiagonal or the subdiagonal of the
33: * tridiagonal matrix A is stored and the form of the
34: * factorization:
35: * = 'U': E is the superdiagonal of A, and A = U**H*D*U;
36: * = 'L': E is the subdiagonal of A, and A = L*D*L**H.
37: * (The two forms are equivalent if A is real.)
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: * D (input) DOUBLE PRECISION array, dimension (N)
47: * The n real diagonal elements of the tridiagonal matrix A.
48: *
49: * E (input) COMPLEX*16 array, dimension (N-1)
50: * The (n-1) off-diagonal elements of the tridiagonal matrix A
51: * (see UPLO).
52: *
53: * DF (input) DOUBLE PRECISION array, dimension (N)
54: * The n diagonal elements of the diagonal matrix D from
55: * the factorization computed by ZPTTRF.
56: *
57: * EF (input) COMPLEX*16 array, dimension (N-1)
58: * The (n-1) off-diagonal elements of the unit bidiagonal
59: * factor U or L from the factorization computed by ZPTTRF
60: * (see UPLO).
61: *
62: * B (input) COMPLEX*16 array, dimension (LDB,NRHS)
63: * The right hand side matrix B.
64: *
65: * LDB (input) INTEGER
66: * The leading dimension of the array B. LDB >= max(1,N).
67: *
68: * X (input/output) COMPLEX*16 array, dimension (LDX,NRHS)
69: * On entry, the solution matrix X, as computed by ZPTTRS.
70: * On exit, the improved solution matrix X.
71: *
72: * LDX (input) INTEGER
73: * The leading dimension of the array X. LDX >= max(1,N).
74: *
75: * FERR (output) DOUBLE PRECISION array, dimension (NRHS)
76: * The forward error bound for each solution vector
77: * X(j) (the j-th column of the solution matrix X).
78: * If XTRUE is the true solution corresponding to X(j), FERR(j)
79: * is an estimated upper bound for the magnitude of the largest
80: * element in (X(j) - XTRUE) divided by the magnitude of the
81: * largest element in X(j).
82: *
83: * BERR (output) DOUBLE PRECISION array, dimension (NRHS)
84: * The componentwise relative backward error of each solution
85: * vector X(j) (i.e., the smallest relative change in
86: * any element of A or B that makes X(j) an exact solution).
87: *
88: * WORK (workspace) COMPLEX*16 array, dimension (N)
89: *
90: * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
91: *
92: * INFO (output) INTEGER
93: * = 0: successful exit
94: * < 0: if INFO = -i, the i-th argument had an illegal value
95: *
96: * Internal Parameters
97: * ===================
98: *
99: * ITMAX is the maximum number of steps of iterative refinement.
100: *
101: * =====================================================================
102: *
103: * .. Parameters ..
104: INTEGER ITMAX
105: PARAMETER ( ITMAX = 5 )
106: DOUBLE PRECISION ZERO
107: PARAMETER ( ZERO = 0.0D+0 )
108: DOUBLE PRECISION ONE
109: PARAMETER ( ONE = 1.0D+0 )
110: DOUBLE PRECISION TWO
111: PARAMETER ( TWO = 2.0D+0 )
112: DOUBLE PRECISION THREE
113: PARAMETER ( THREE = 3.0D+0 )
114: * ..
115: * .. Local Scalars ..
116: LOGICAL UPPER
117: INTEGER COUNT, I, IX, J, NZ
118: DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN
119: COMPLEX*16 BI, CX, DX, EX, ZDUM
120: * ..
121: * .. External Functions ..
122: LOGICAL LSAME
123: INTEGER IDAMAX
124: DOUBLE PRECISION DLAMCH
125: EXTERNAL LSAME, IDAMAX, DLAMCH
126: * ..
127: * .. External Subroutines ..
128: EXTERNAL XERBLA, ZAXPY, ZPTTRS
129: * ..
130: * .. Intrinsic Functions ..
131: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
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: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
146: INFO = -1
147: ELSE IF( N.LT.0 ) THEN
148: INFO = -2
149: ELSE IF( NRHS.LT.0 ) THEN
150: INFO = -3
151: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
152: INFO = -9
153: ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
154: INFO = -11
155: END IF
156: IF( INFO.NE.0 ) THEN
157: CALL XERBLA( 'ZPTRFS', -INFO )
158: RETURN
159: END IF
160: *
161: * Quick return if possible
162: *
163: IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
164: DO 10 J = 1, NRHS
165: FERR( J ) = ZERO
166: BERR( J ) = ZERO
167: 10 CONTINUE
168: RETURN
169: END IF
170: *
171: * NZ = maximum number of nonzero elements in each row of A, plus 1
172: *
173: NZ = 4
174: EPS = DLAMCH( 'Epsilon' )
175: SAFMIN = DLAMCH( 'Safe minimum' )
176: SAFE1 = NZ*SAFMIN
177: SAFE2 = SAFE1 / EPS
178: *
179: * Do for each right hand side
180: *
181: DO 100 J = 1, NRHS
182: *
183: COUNT = 1
184: LSTRES = THREE
185: 20 CONTINUE
186: *
187: * Loop until stopping criterion is satisfied.
188: *
189: * Compute residual R = B - A * X. Also compute
190: * abs(A)*abs(x) + abs(b) for use in the backward error bound.
191: *
192: IF( UPPER ) THEN
193: IF( N.EQ.1 ) THEN
194: BI = B( 1, J )
195: DX = D( 1 )*X( 1, J )
196: WORK( 1 ) = BI - DX
197: RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
198: ELSE
199: BI = B( 1, J )
200: DX = D( 1 )*X( 1, J )
201: EX = E( 1 )*X( 2, J )
202: WORK( 1 ) = BI - DX - EX
203: RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
204: $ CABS1( E( 1 ) )*CABS1( X( 2, J ) )
205: DO 30 I = 2, N - 1
206: BI = B( I, J )
207: CX = DCONJG( E( I-1 ) )*X( I-1, J )
208: DX = D( I )*X( I, J )
209: EX = E( I )*X( I+1, J )
210: WORK( I ) = BI - CX - DX - EX
211: RWORK( I ) = CABS1( BI ) +
212: $ CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
213: $ CABS1( DX ) + CABS1( E( I ) )*
214: $ CABS1( X( I+1, J ) )
215: 30 CONTINUE
216: BI = B( N, J )
217: CX = DCONJG( E( N-1 ) )*X( N-1, J )
218: DX = D( N )*X( N, J )
219: WORK( N ) = BI - CX - DX
220: RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
221: $ CABS1( X( N-1, J ) ) + CABS1( DX )
222: END IF
223: ELSE
224: IF( N.EQ.1 ) THEN
225: BI = B( 1, J )
226: DX = D( 1 )*X( 1, J )
227: WORK( 1 ) = BI - DX
228: RWORK( 1 ) = CABS1( BI ) + CABS1( DX )
229: ELSE
230: BI = B( 1, J )
231: DX = D( 1 )*X( 1, J )
232: EX = DCONJG( E( 1 ) )*X( 2, J )
233: WORK( 1 ) = BI - DX - EX
234: RWORK( 1 ) = CABS1( BI ) + CABS1( DX ) +
235: $ CABS1( E( 1 ) )*CABS1( X( 2, J ) )
236: DO 40 I = 2, N - 1
237: BI = B( I, J )
238: CX = E( I-1 )*X( I-1, J )
239: DX = D( I )*X( I, J )
240: EX = DCONJG( E( I ) )*X( I+1, J )
241: WORK( I ) = BI - CX - DX - EX
242: RWORK( I ) = CABS1( BI ) +
243: $ CABS1( E( I-1 ) )*CABS1( X( I-1, J ) ) +
244: $ CABS1( DX ) + CABS1( E( I ) )*
245: $ CABS1( X( I+1, J ) )
246: 40 CONTINUE
247: BI = B( N, J )
248: CX = E( N-1 )*X( N-1, J )
249: DX = D( N )*X( N, J )
250: WORK( N ) = BI - CX - DX
251: RWORK( N ) = CABS1( BI ) + CABS1( E( N-1 ) )*
252: $ CABS1( X( N-1, J ) ) + CABS1( DX )
253: END IF
254: END IF
255: *
256: * Compute componentwise relative backward error from formula
257: *
258: * max(i) ( abs(R(i)) / ( abs(A)*abs(X) + abs(B) )(i) )
259: *
260: * where abs(Z) is the componentwise absolute value of the matrix
261: * or vector Z. If the i-th component of the denominator is less
262: * than SAFE2, then SAFE1 is added to the i-th components of the
263: * numerator and denominator before dividing.
264: *
265: S = ZERO
266: DO 50 I = 1, N
267: IF( RWORK( I ).GT.SAFE2 ) THEN
268: S = MAX( S, CABS1( WORK( I ) ) / RWORK( I ) )
269: ELSE
270: S = MAX( S, ( CABS1( WORK( I ) )+SAFE1 ) /
271: $ ( RWORK( I )+SAFE1 ) )
272: END IF
273: 50 CONTINUE
274: BERR( J ) = S
275: *
276: * Test stopping criterion. Continue iterating if
277: * 1) The residual BERR(J) is larger than machine epsilon, and
278: * 2) BERR(J) decreased by at least a factor of 2 during the
279: * last iteration, and
280: * 3) At most ITMAX iterations tried.
281: *
282: IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
283: $ COUNT.LE.ITMAX ) THEN
284: *
285: * Update solution and try again.
286: *
287: CALL ZPTTRS( UPLO, N, 1, DF, EF, WORK, N, INFO )
288: CALL ZAXPY( N, DCMPLX( ONE ), WORK, 1, X( 1, J ), 1 )
289: LSTRES = BERR( J )
290: COUNT = COUNT + 1
291: GO TO 20
292: END IF
293: *
294: * Bound error from formula
295: *
296: * norm(X - XTRUE) / norm(X) .le. FERR =
297: * norm( abs(inv(A))*
298: * ( abs(R) + NZ*EPS*( abs(A)*abs(X)+abs(B) ))) / norm(X)
299: *
300: * where
301: * norm(Z) is the magnitude of the largest component of Z
302: * inv(A) is the inverse of A
303: * abs(Z) is the componentwise absolute value of the matrix or
304: * vector Z
305: * NZ is the maximum number of nonzeros in any row of A, plus 1
306: * EPS is machine epsilon
307: *
308: * The i-th component of abs(R)+NZ*EPS*(abs(A)*abs(X)+abs(B))
309: * is incremented by SAFE1 if the i-th component of
310: * abs(A)*abs(X) + abs(B) is less than SAFE2.
311: *
312: DO 60 I = 1, N
313: IF( RWORK( I ).GT.SAFE2 ) THEN
314: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I )
315: ELSE
316: RWORK( I ) = CABS1( WORK( I ) ) + NZ*EPS*RWORK( I ) +
317: $ SAFE1
318: END IF
319: 60 CONTINUE
320: IX = IDAMAX( N, RWORK, 1 )
321: FERR( J ) = RWORK( IX )
322: *
323: * Estimate the norm of inv(A).
324: *
325: * Solve M(A) * x = e, where M(A) = (m(i,j)) is given by
326: *
327: * m(i,j) = abs(A(i,j)), i = j,
328: * m(i,j) = -abs(A(i,j)), i .ne. j,
329: *
330: * and e = [ 1, 1, ..., 1 ]**T. Note M(A) = M(L)*D*M(L)**H.
331: *
332: * Solve M(L) * x = e.
333: *
334: RWORK( 1 ) = ONE
335: DO 70 I = 2, N
336: RWORK( I ) = ONE + RWORK( I-1 )*ABS( EF( I-1 ) )
337: 70 CONTINUE
338: *
339: * Solve D * M(L)**H * x = b.
340: *
341: RWORK( N ) = RWORK( N ) / DF( N )
342: DO 80 I = N - 1, 1, -1
343: RWORK( I ) = RWORK( I ) / DF( I ) +
344: $ RWORK( I+1 )*ABS( EF( I ) )
345: 80 CONTINUE
346: *
347: * Compute norm(inv(A)) = max(x(i)), 1<=i<=n.
348: *
349: IX = IDAMAX( N, RWORK, 1 )
350: FERR( J ) = FERR( J )*ABS( RWORK( IX ) )
351: *
352: * Normalize error.
353: *
354: LSTRES = ZERO
355: DO 90 I = 1, N
356: LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
357: 90 CONTINUE
358: IF( LSTRES.NE.ZERO )
359: $ FERR( J ) = FERR( J ) / LSTRES
360: *
361: 100 CONTINUE
362: *
363: RETURN
364: *
365: * End of ZPTRFS
366: *
367: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zptrfs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack