|
5164
|
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 |
