1: SUBROUTINE DPTTRF( 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( * ), E( * )
13: * ..
14: *
15: * Purpose
16: * =======
17: *
18: * DPTTRF computes the L*D*L' factorization of a real symmetric
19: * positive definite tridiagonal matrix A. The factorization may also
20: * be regarded as having the form A = U'*D*U.
21: *
22: * Arguments
23: * =========
24: *
25: * N (input) INTEGER
26: * The order of the matrix A. N >= 0.
27: *
28: * D (input/output) DOUBLE PRECISION array, dimension (N)
29: * On entry, the n diagonal elements of the tridiagonal matrix
30: * A. On exit, the n diagonal elements of the diagonal matrix
31: * D from the L*D*L' factorization of A.
32: *
33: * E (input/output) DOUBLE PRECISION array, dimension (N-1)
34: * On entry, the (n-1) subdiagonal elements of the tridiagonal
35: * matrix A. On exit, the (n-1) subdiagonal elements of the
36: * unit bidiagonal factor L from the L*D*L' factorization of A.
37: * E can also be regarded as the superdiagonal of the unit
38: * bidiagonal factor U from the U'*D*U factorization of A.
39: *
40: * INFO (output) INTEGER
41: * = 0: successful exit
42: * < 0: if INFO = -k, the k-th argument had an illegal value
43: * > 0: if INFO = k, the leading minor of order k is not
44: * positive definite; if k < N, the factorization could not
45: * be completed, while if k = N, the factorization was
46: * completed, but D(N) <= 0.
47: *
48: * =====================================================================
49: *
50: * .. Parameters ..
51: DOUBLE PRECISION ZERO
52: PARAMETER ( ZERO = 0.0D+0 )
53: * ..
54: * .. Local Scalars ..
55: INTEGER I, I4
56: DOUBLE PRECISION EI
57: * ..
58: * .. External Subroutines ..
59: EXTERNAL XERBLA
60: * ..
61: * .. Intrinsic Functions ..
62: INTRINSIC MOD
63: * ..
64: * .. Executable Statements ..
65: *
66: * Test the input parameters.
67: *
68: INFO = 0
69: IF( N.LT.0 ) THEN
70: INFO = -1
71: CALL XERBLA( 'DPTTRF', -INFO )
72: RETURN
73: END IF
74: *
75: * Quick return if possible
76: *
77: IF( N.EQ.0 )
78: $ RETURN
79: *
80: * Compute the L*D*L' (or U'*D*U) factorization of A.
81: *
82: I4 = MOD( N-1, 4 )
83: DO 10 I = 1, I4
84: IF( D( I ).LE.ZERO ) THEN
85: INFO = I
86: GO TO 30
87: END IF
88: EI = E( I )
89: E( I ) = EI / D( I )
90: D( I+1 ) = D( I+1 ) - E( I )*EI
91: 10 CONTINUE
92: *
93: DO 20 I = I4 + 1, N - 4, 4
94: *
95: * Drop out of the loop if d(i) <= 0: the matrix is not positive
96: * definite.
97: *
98: IF( D( I ).LE.ZERO ) THEN
99: INFO = I
100: GO TO 30
101: END IF
102: *
103: * Solve for e(i) and d(i+1).
104: *
105: EI = E( I )
106: E( I ) = EI / D( I )
107: D( I+1 ) = D( I+1 ) - E( I )*EI
108: *
109: IF( D( I+1 ).LE.ZERO ) THEN
110: INFO = I + 1
111: GO TO 30
112: END IF
113: *
114: * Solve for e(i+1) and d(i+2).
115: *
116: EI = E( I+1 )
117: E( I+1 ) = EI / D( I+1 )
118: D( I+2 ) = D( I+2 ) - E( I+1 )*EI
119: *
120: IF( D( I+2 ).LE.ZERO ) THEN
121: INFO = I + 2
122: GO TO 30
123: END IF
124: *
125: * Solve for e(i+2) and d(i+3).
126: *
127: EI = E( I+2 )
128: E( I+2 ) = EI / D( I+2 )
129: D( I+3 ) = D( I+3 ) - E( I+2 )*EI
130: *
131: IF( D( I+3 ).LE.ZERO ) THEN
132: INFO = I + 3
133: GO TO 30
134: END IF
135: *
136: * Solve for e(i+3) and d(i+4).
137: *
138: EI = E( I+3 )
139: E( I+3 ) = EI / D( I+3 )
140: D( I+4 ) = D( I+4 ) - E( I+3 )*EI
141: 20 CONTINUE
142: *
143: * Check d(n) for positive definiteness.
144: *
145: IF( D( N ).LE.ZERO )
146: $ INFO = N
147: *
148: 30 CONTINUE
149: RETURN
150: *
151: * End of DPTTRF
152: *
153: END
CVSweb interface <joel.bertrand@systella.fr>