Mercurial > octave-nkf

comparison libcruft/lapack/zpttrf.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 ZPTTRF( 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( * )
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