001:       SUBROUTINE ZPTTRF( N, D, E, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, N
010: *     ..
011: *     .. Array Arguments ..
012:       DOUBLE PRECISION   D( * )
013:       COMPLEX*16         E( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZPTTRF computes the L*D*L' factorization of a complex Hermitian
020: *  positive definite tridiagonal matrix A.  The factorization may also
021: *  be regarded as having the form A = U'*D*U.
022: *
023: *  Arguments
024: *  =========
025: *
026: *  N       (input) INTEGER
027: *          The order of the matrix A.  N >= 0.
028: *
029: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
030: *          On entry, the n diagonal elements of the tridiagonal matrix
031: *          A.  On exit, the n diagonal elements of the diagonal matrix
032: *          D from the L*D*L' factorization of A.
033: *
034: *  E       (input/output) COMPLEX*16 array, dimension (N-1)
035: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
036: *          matrix A.  On exit, the (n-1) subdiagonal elements of the
037: *          unit bidiagonal factor L from the L*D*L' factorization of A.
038: *          E can also be regarded as the superdiagonal of the unit
039: *          bidiagonal factor U from the U'*D*U factorization of A.
040: *
041: *  INFO    (output) INTEGER
042: *          = 0: successful exit
043: *          < 0: if INFO = -k, the k-th argument had an illegal value
044: *          > 0: if INFO = k, the leading minor of order k is not
045: *               positive definite; if k < N, the factorization could not
046: *               be completed, while if k = N, the factorization was
047: *               completed, but D(N) <= 0.
048: *
049: *  =====================================================================
050: *
051: *     .. Parameters ..
052:       DOUBLE PRECISION   ZERO
053:       PARAMETER          ( ZERO = 0.0D+0 )
054: *     ..
055: *     .. Local Scalars ..
056:       INTEGER            I, I4
057:       DOUBLE PRECISION   EII, EIR, F, G
058: *     ..
059: *     .. External Subroutines ..
060:       EXTERNAL           XERBLA
061: *     ..
062: *     .. Intrinsic Functions ..
063:       INTRINSIC          DBLE, DCMPLX, DIMAG, MOD
064: *     ..
065: *     .. Executable Statements ..
066: *
067: *     Test the input parameters.
068: *
069:       INFO = 0
070:       IF( N.LT.0 ) THEN
071:          INFO = -1
072:          CALL XERBLA( 'ZPTTRF', -INFO )
073:          RETURN
074:       END IF
075: *
076: *     Quick return if possible
077: *
078:       IF( N.EQ.0 )
079:      $   RETURN
080: *
081: *     Compute the L*D*L' (or U'*D*U) factorization of A.
082: *
083:       I4 = MOD( N-1, 4 )
084:       DO 10 I = 1, I4
085:          IF( D( I ).LE.ZERO ) THEN
086:             INFO = I
087:             GO TO 30
088:          END IF
089:          EIR = DBLE( E( I ) )
090:          EII = DIMAG( E( I ) )
091:          F = EIR / D( I )
092:          G = EII / D( I )
093:          E( I ) = DCMPLX( F, G )
094:          D( I+1 ) = D( I+1 ) - F*EIR - G*EII
095:    10 CONTINUE
096: *
097:       DO 20 I = I4 + 1, N - 4, 4
098: *
099: *        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
170: