1: *> \brief \b DPTTRF
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DPTTRF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpttrf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpttrf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpttrf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DPTTRF( N, D, E, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION D( * ), E( * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DPTTRF computes the L*D*L**T factorization of a real symmetric
37: *> positive definite tridiagonal matrix A. The factorization may also
38: *> be regarded as having the form A = U**T*D*U.
39: *> \endverbatim
40: *
41: * Arguments:
42: * ==========
43: *
44: *> \param[in] N
45: *> \verbatim
46: *> N is INTEGER
47: *> The order of the matrix A. N >= 0.
48: *> \endverbatim
49: *>
50: *> \param[in,out] D
51: *> \verbatim
52: *> D is DOUBLE PRECISION array, dimension (N)
53: *> On entry, the n diagonal elements of the tridiagonal matrix
54: *> A. On exit, the n diagonal elements of the diagonal matrix
55: *> D from the L*D*L**T factorization of A.
56: *> \endverbatim
57: *>
58: *> \param[in,out] E
59: *> \verbatim
60: *> E is DOUBLE PRECISION array, dimension (N-1)
61: *> On entry, the (n-1) subdiagonal elements of the tridiagonal
62: *> matrix A. On exit, the (n-1) subdiagonal elements of the
63: *> unit bidiagonal factor L from the L*D*L**T factorization of A.
64: *> E can also be regarded as the superdiagonal of the unit
65: *> bidiagonal factor U from the U**T*D*U factorization of A.
66: *> \endverbatim
67: *>
68: *> \param[out] INFO
69: *> \verbatim
70: *> INFO is INTEGER
71: *> = 0: successful exit
72: *> < 0: if INFO = -k, the k-th argument had an illegal value
73: *> > 0: if INFO = k, the leading minor of order k is not
74: *> positive definite; if k < N, the factorization could not
75: *> be completed, while if k = N, the factorization was
76: *> completed, but D(N) <= 0.
77: *> \endverbatim
78: *
79: * Authors:
80: * ========
81: *
82: *> \author Univ. of Tennessee
83: *> \author Univ. of California Berkeley
84: *> \author Univ. of Colorado Denver
85: *> \author NAG Ltd.
86: *
87: *> \date November 2011
88: *
89: *> \ingroup auxOTHERcomputational
90: *
91: * =====================================================================
92: SUBROUTINE DPTTRF( N, D, E, INFO )
93: *
94: * -- LAPACK computational routine (version 3.4.0) --
95: * -- LAPACK is a software package provided by Univ. of Tennessee, --
96: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
97: * November 2011
98: *
99: * .. Scalar Arguments ..
100: INTEGER INFO, N
101: * ..
102: * .. Array Arguments ..
103: DOUBLE PRECISION D( * ), E( * )
104: * ..
105: *
106: * =====================================================================
107: *
108: * .. Parameters ..
109: DOUBLE PRECISION ZERO
110: PARAMETER ( ZERO = 0.0D+0 )
111: * ..
112: * .. Local Scalars ..
113: INTEGER I, I4
114: DOUBLE PRECISION EI
115: * ..
116: * .. External Subroutines ..
117: EXTERNAL XERBLA
118: * ..
119: * .. Intrinsic Functions ..
120: INTRINSIC MOD
121: * ..
122: * .. Executable Statements ..
123: *
124: * Test the input parameters.
125: *
126: INFO = 0
127: IF( N.LT.0 ) THEN
128: INFO = -1
129: CALL XERBLA( 'DPTTRF', -INFO )
130: RETURN
131: END IF
132: *
133: * Quick return if possible
134: *
135: IF( N.EQ.0 )
136: $ RETURN
137: *
138: * Compute the L*D*L**T (or U**T*D*U) factorization of A.
139: *
140: I4 = MOD( N-1, 4 )
141: DO 10 I = 1, I4
142: IF( D( I ).LE.ZERO ) THEN
143: INFO = I
144: GO TO 30
145: END IF
146: EI = E( I )
147: E( I ) = EI / D( I )
148: D( I+1 ) = D( I+1 ) - E( I )*EI
149: 10 CONTINUE
150: *
151: DO 20 I = I4 + 1, N - 4, 4
152: *
153: * Drop out of the loop if d(i) <= 0: the matrix is not positive
154: * definite.
155: *
156: IF( D( I ).LE.ZERO ) THEN
157: INFO = I
158: GO TO 30
159: END IF
160: *
161: * Solve for e(i) and d(i+1).
162: *
163: EI = E( I )
164: E( I ) = EI / D( I )
165: D( I+1 ) = D( I+1 ) - E( I )*EI
166: *
167: IF( D( I+1 ).LE.ZERO ) THEN
168: INFO = I + 1
169: GO TO 30
170: END IF
171: *
172: * Solve for e(i+1) and d(i+2).
173: *
174: EI = E( I+1 )
175: E( I+1 ) = EI / D( I+1 )
176: D( I+2 ) = D( I+2 ) - E( I+1 )*EI
177: *
178: IF( D( I+2 ).LE.ZERO ) THEN
179: INFO = I + 2
180: GO TO 30
181: END IF
182: *
183: * Solve for e(i+2) and d(i+3).
184: *
185: EI = E( I+2 )
186: E( I+2 ) = EI / D( I+2 )
187: D( I+3 ) = D( I+3 ) - E( I+2 )*EI
188: *
189: IF( D( I+3 ).LE.ZERO ) THEN
190: INFO = I + 3
191: GO TO 30
192: END IF
193: *
194: * Solve for e(i+3) and d(i+4).
195: *
196: EI = E( I+3 )
197: E( I+3 ) = EI / D( I+3 )
198: D( I+4 ) = D( I+4 ) - E( I+3 )*EI
199: 20 CONTINUE
200: *
201: * Check d(n) for positive definiteness.
202: *
203: IF( D( N ).LE.ZERO )
204: $ INFO = N
205: *
206: 30 CONTINUE
207: RETURN
208: *
209: * End of DPTTRF
210: *
211: END
CVSweb interface <joel.bertrand@systella.fr>