diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/zpttrf.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/zpttrf.f')
-rw-r--r-- | SRC/zpttrf.f | 168 |
1 files changed, 168 insertions, 0 deletions
diff --git a/SRC/zpttrf.f b/SRC/zpttrf.f new file mode 100644 index 00000000..6deda45f --- /dev/null +++ b/SRC/zpttrf.f @@ -0,0 +1,168 @@ + SUBROUTINE ZPTTRF( N, D, E, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION D( * ) + COMPLEX*16 E( * ) +* .. +* +* Purpose +* ======= +* +* ZPTTRF computes the L*D*L' factorization of a complex Hermitian +* positive definite tridiagonal matrix A. The factorization may also +* be regarded as having the form A = U'*D*U. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. N >= 0. +* +* D (input/output) DOUBLE PRECISION array, dimension (N) +* On entry, the n diagonal elements of the tridiagonal matrix +* A. On exit, the n diagonal elements of the diagonal matrix +* D from the L*D*L' factorization of A. +* +* E (input/output) COMPLEX*16 array, dimension (N-1) +* On entry, the (n-1) subdiagonal elements of the tridiagonal +* matrix A. On exit, the (n-1) subdiagonal elements of the +* unit bidiagonal factor L from the L*D*L' factorization of A. +* E can also be regarded as the superdiagonal of the unit +* bidiagonal factor U from the U'*D*U factorization of A. +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -k, the k-th argument had an illegal value +* > 0: if INFO = k, the leading minor of order k is not +* positive definite; if k < N, the factorization could not +* be completed, while if k = N, the factorization was +* completed, but D(N) <= 0. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, I4 + DOUBLE PRECISION EII, EIR, F, G +* .. +* .. External Subroutines .. + EXTERNAL XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC DBLE, DCMPLX, DIMAG, MOD +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + INFO = 0 + IF( N.LT.0 ) THEN + INFO = -1 + CALL XERBLA( 'ZPTTRF', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Compute the L*D*L' (or U'*D*U) factorization of A. +* + I4 = MOD( N-1, 4 ) + DO 10 I = 1, I4 + IF( D( I ).LE.ZERO ) THEN + INFO = I + GO TO 30 + END IF + EIR = DBLE( E( I ) ) + EII = DIMAG( E( I ) ) + F = EIR / D( I ) + G = EII / D( I ) + E( I ) = DCMPLX( F, G ) + D( I+1 ) = D( I+1 ) - F*EIR - G*EII + 10 CONTINUE +* + DO 20 I = I4 + 1, N - 4, 4 +* +* Drop out of the loop if d(i) <= 0: the matrix is not positive +* definite. +* + IF( D( I ).LE.ZERO ) THEN + INFO = I + GO TO 30 + END IF +* +* Solve for e(i) and d(i+1). +* + EIR = DBLE( E( I ) ) + EII = DIMAG( E( I ) ) + F = EIR / D( I ) + G = EII / D( I ) + E( I ) = DCMPLX( F, G ) + D( I+1 ) = D( I+1 ) - F*EIR - G*EII +* + IF( D( I+1 ).LE.ZERO ) THEN + INFO = I + 1 + GO TO 30 + END IF +* +* Solve for e(i+1) and d(i+2). +* + EIR = DBLE( E( I+1 ) ) + EII = DIMAG( E( I+1 ) ) + F = EIR / D( I+1 ) + G = EII / D( I+1 ) + E( I+1 ) = DCMPLX( F, G ) + D( I+2 ) = D( I+2 ) - F*EIR - G*EII +* + IF( D( I+2 ).LE.ZERO ) THEN + INFO = I + 2 + GO TO 30 + END IF +* +* Solve for e(i+2) and d(i+3). +* + EIR = DBLE( E( I+2 ) ) + EII = DIMAG( E( I+2 ) ) + F = EIR / D( I+2 ) + G = EII / D( I+2 ) + E( I+2 ) = DCMPLX( F, G ) + D( I+3 ) = D( I+3 ) - F*EIR - G*EII +* + IF( D( I+3 ).LE.ZERO ) THEN + INFO = I + 3 + GO TO 30 + END IF +* +* Solve for e(i+3) and d(i+4). +* + EIR = DBLE( E( I+3 ) ) + EII = DIMAG( E( I+3 ) ) + F = EIR / D( I+3 ) + G = EII / D( I+3 ) + E( I+3 ) = DCMPLX( F, G ) + D( I+4 ) = D( I+4 ) - F*EIR - G*EII + 20 CONTINUE +* +* Check d(n) for positive definiteness. +* + IF( D( N ).LE.ZERO ) + $ INFO = N +* + 30 CONTINUE + RETURN +* +* End of ZPTTRF +* + END |