Mercurial > octave

comparison libcruft/lapack/spttrf.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SPTTRF( N, D, E, INFO )
2 *
3 * -- LAPACK routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 INTEGER INFO, N
9 * ..
10 * .. Array Arguments ..
11 REAL D( * ), E( * )
12 * ..
13 *
14 * Purpose
15 * =======
16 *
17 * SPTTRF computes the L*D*L' factorization of a real symmetric
18 * positive definite tridiagonal matrix A. The factorization may also
19 * be regarded as having the form A = U'*D*U.
20 *
21 * Arguments
22 * =========
23 *
24 * N (input) INTEGER
25 * The order of the matrix A. N >= 0.
26 *
27 * D (input/output) REAL array, dimension (N)
28 * On entry, the n diagonal elements of the tridiagonal matrix
29 * A. On exit, the n diagonal elements of the diagonal matrix
30 * D from the L*D*L' factorization of A.
31 *
32 * E (input/output) REAL array, dimension (N-1)
33 * On entry, the (n-1) subdiagonal elements of the tridiagonal
34 * matrix A. On exit, the (n-1) subdiagonal elements of the
35 * unit bidiagonal factor L from the L*D*L' factorization of A.
36 * E can also be regarded as the superdiagonal of the unit
37 * bidiagonal factor U from the U'*D*U factorization of A.
38 *
39 * INFO (output) INTEGER
40 * = 0: successful exit
41 * < 0: if INFO = -k, the k-th argument had an illegal value
42 * > 0: if INFO = k, the leading minor of order k is not
43 * positive definite; if k < N, the factorization could not
44 * be completed, while if k = N, the factorization was
45 * completed, but D(N) <= 0.
46 *
47 * =====================================================================
48 *
49 * .. Parameters ..
50 REAL ZERO
51 PARAMETER ( ZERO = 0.0E+0 )
52 * ..
53 * .. Local Scalars ..
54 INTEGER I, I4
55 REAL EI
56 * ..
57 * .. External Subroutines ..
58 EXTERNAL XERBLA
59 * ..
60 * .. Intrinsic Functions ..
61 INTRINSIC MOD
62 * ..
63 * .. Executable Statements ..
64 *
65 * Test the input parameters.
66 *
67 INFO = 0
68 IF( N.LT.0 ) THEN
69 INFO = -1
70 CALL XERBLA( 'SPTTRF', -INFO )
71 RETURN
72 END IF
73 *
74 * Quick return if possible
75 *
76 IF( N.EQ.0 )
77 $ RETURN
78 *
79 * Compute the L*D*L' (or U'*D*U) factorization of A.
80 *
81 I4 = MOD( N-1, 4 )
82 DO 10 I = 1, I4
83 IF( D( I ).LE.ZERO ) THEN
84 INFO = I
85 GO TO 30
86 END IF
87 EI = E( I )
88 E( I ) = EI / D( I )
89 D( I+1 ) = D( I+1 ) - E( I )*EI
90 10 CONTINUE
91 *
92 DO 20 I = I4 + 1, N - 4, 4
93 *
94 * Drop out of the loop if d(i) <= 0: the matrix is not positive
95 * definite.
96 *
97 IF( D( I ).LE.ZERO ) THEN
98 INFO = I
99 GO TO 30
100 END IF
101 *
102 * Solve for e(i) and d(i+1).
103 *
104 EI = E( I )
105 E( I ) = EI / D( I )
106 D( I+1 ) = D( I+1 ) - E( I )*EI
107 *
108 IF( D( I+1 ).LE.ZERO ) THEN
109 INFO = I + 1
110 GO TO 30
111 END IF
112 *
113 * Solve for e(i+1) and d(i+2).
114 *
115 EI = E( I+1 )
116 E( I+1 ) = EI / D( I+1 )
117 D( I+2 ) = D( I+2 ) - E( I+1 )*EI
118 *
119 IF( D( I+2 ).LE.ZERO ) THEN
120 INFO = I + 2
121 GO TO 30
122 END IF
123 *
124 * Solve for e(i+2) and d(i+3).
125 *
126 EI = E( I+2 )
127 E( I+2 ) = EI / D( I+2 )
128 D( I+3 ) = D( I+3 ) - E( I+2 )*EI
129 *
130 IF( D( I+3 ).LE.ZERO ) THEN
131 INFO = I + 3
132 GO TO 30
133 END IF
134 *
135 * Solve for e(i+3) and d(i+4).
136 *
137 EI = E( I+3 )
138 E( I+3 ) = EI / D( I+3 )
139 D( I+4 ) = D( I+4 ) - E( I+3 )*EI
140 20 CONTINUE
141 *
142 * Check d(n) for positive definiteness.
143 *
144 IF( D( N ).LE.ZERO )
145 $ INFO = N
146 *
147 30 CONTINUE
148 RETURN
149 *
150 * End of SPTTRF
151 *
152 END