Mercurial > octave-nkf

comparison libcruft/lapack/dpttrf.f @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children 68db500cb558
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 SUBROUTINE DPTTRF( N, D, E, INFO )
2 *
3 * -- LAPACK routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * June 30, 1999
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