![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpttrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to zpttrf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE ZPTTRF( N, D, E, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER INFO, N
10: * ..
11: * .. Array Arguments ..
12: DOUBLE PRECISION D( * )
13: COMPLEX*16 E( * )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * ZPTTRF computes the L*D*L' factorization of a complex Hermitian
20: * positive definite tridiagonal matrix A. The factorization may also
21: * be regarded as having the form A = U'*D*U.
22: *
23: * Arguments
24: * =========
25: *
26: * N (input) INTEGER
27: * The order of the matrix A. N >= 0.
28: *
29: * D (input/output) DOUBLE PRECISION array, dimension (N)
30: * On entry, the n diagonal elements of the tridiagonal matrix
31: * A. On exit, the n diagonal elements of the diagonal matrix
32: * D from the L*D*L' factorization of A.
33: *
34: * E (input/output) COMPLEX*16 array, dimension (N-1)
35: * On entry, the (n-1) subdiagonal elements of the tridiagonal
36: * matrix A. On exit, the (n-1) subdiagonal elements of the
37: * unit bidiagonal factor L from the L*D*L' factorization of A.
38: * E can also be regarded as the superdiagonal of the unit
39: * bidiagonal factor U from the U'*D*U factorization of A.
40: *
41: * INFO (output) INTEGER
42: * = 0: successful exit
43: * < 0: if INFO = -k, the k-th argument had an illegal value
44: * > 0: if INFO = k, the leading minor of order k is not
45: * positive definite; if k < N, the factorization could not
46: * be completed, while if k = N, the factorization was
47: * completed, but D(N) <= 0.
48: *
49: * =====================================================================
50: *
51: * .. Parameters ..
52: DOUBLE PRECISION ZERO
53: PARAMETER ( ZERO = 0.0D+0 )
54: * ..
55: * .. Local Scalars ..
56: INTEGER I, I4
57: DOUBLE PRECISION EII, EIR, F, G
58: * ..
59: * .. External Subroutines ..
60: EXTERNAL XERBLA
61: * ..
62: * .. Intrinsic Functions ..
63: INTRINSIC DBLE, DCMPLX, DIMAG, MOD
64: * ..
65: * .. Executable Statements ..
66: *
67: * Test the input parameters.
68: *
69: INFO = 0
70: IF( N.LT.0 ) THEN
71: INFO = -1
72: CALL XERBLA( 'ZPTTRF', -INFO )
73: RETURN
74: END IF
75: *
76: * Quick return if possible
77: *
78: IF( N.EQ.0 )
79: $ RETURN
80: *
81: * Compute the L*D*L' (or U'*D*U) factorization of A.
82: *
83: I4 = MOD( N-1, 4 )
84: DO 10 I = 1, I4
85: IF( D( I ).LE.ZERO ) THEN
86: INFO = I
87: GO TO 30
88: END IF
89: EIR = DBLE( E( I ) )
90: EII = DIMAG( E( I ) )
91: F = EIR / D( I )
92: G = EII / D( I )
93: E( I ) = DCMPLX( F, G )
94: D( I+1 ) = D( I+1 ) - F*EIR - G*EII
95: 10 CONTINUE
96: *
97: DO 20 I = I4 + 1, N - 4, 4
98: *
99: * Drop out of the loop if d(i) <= 0: the matrix is not positive
100: * definite.
101: *
102: IF( D( I ).LE.ZERO ) THEN
103: INFO = I
104: GO TO 30
105: END IF
106: *
107: * Solve for e(i) and d(i+1).
108: *
109: EIR = DBLE( E( I ) )
110: EII = DIMAG( E( I ) )
111: F = EIR / D( I )
112: G = EII / D( I )
113: E( I ) = DCMPLX( F, G )
114: D( I+1 ) = D( I+1 ) - F*EIR - G*EII
115: *
116: IF( D( I+1 ).LE.ZERO ) THEN
117: INFO = I + 1
118: GO TO 30
119: END IF
120: *
121: * Solve for e(i+1) and d(i+2).
122: *
123: EIR = DBLE( E( I+1 ) )
124: EII = DIMAG( E( I+1 ) )
125: F = EIR / D( I+1 )
126: G = EII / D( I+1 )
127: E( I+1 ) = DCMPLX( F, G )
128: D( I+2 ) = D( I+2 ) - F*EIR - G*EII
129: *
130: IF( D( I+2 ).LE.ZERO ) THEN
131: INFO = I + 2
132: GO TO 30
133: END IF
134: *
135: * Solve for e(i+2) and d(i+3).
136: *
137: EIR = DBLE( E( I+2 ) )
138: EII = DIMAG( E( I+2 ) )
139: F = EIR / D( I+2 )
140: G = EII / D( I+2 )
141: E( I+2 ) = DCMPLX( F, G )
142: D( I+3 ) = D( I+3 ) - F*EIR - G*EII
143: *
144: IF( D( I+3 ).LE.ZERO ) THEN
145: INFO = I + 3
146: GO TO 30
147: END IF
148: *
149: * Solve for e(i+3) and d(i+4).
150: *
151: EIR = DBLE( E( I+3 ) )
152: EII = DIMAG( E( I+3 ) )
153: F = EIR / D( I+3 )
154: G = EII / D( I+3 )
155: E( I+3 ) = DCMPLX( F, G )
156: D( I+4 ) = D( I+4 ) - F*EIR - G*EII
157: 20 CONTINUE
158: *
159: * Check d(n) for positive definiteness.
160: *
161: IF( D( N ).LE.ZERO )
162: $ INFO = N
163: *
164: 30 CONTINUE
165: RETURN
166: *
167: * End of ZPTTRF
168: *
169: END