LAPACK: SRC/zpttrf.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZPTTRF( N, D, E, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       DOUBLE PRECISION   D( * )
00013       COMPLEX*16         E( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZPTTRF computes the L*D*L' factorization of a complex Hermitian
00020 *  positive definite tridiagonal matrix A.  The factorization may also
00021 *  be regarded as having the form A = U'*D*U.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N       (input) INTEGER
00027 *          The order of the matrix A.  N >= 0.
00028 *
00029 *  D       (input/output) DOUBLE PRECISION array, dimension (N)
00030 *          On entry, the n diagonal elements of the tridiagonal matrix
00031 *          A.  On exit, the n diagonal elements of the diagonal matrix
00032 *          D from the L*D*L' factorization of A.
00033 *
00034 *  E       (input/output) COMPLEX*16 array, dimension (N-1)
00035 *          On entry, the (n-1) subdiagonal elements of the tridiagonal
00036 *          matrix A.  On exit, the (n-1) subdiagonal elements of the
00037 *          unit bidiagonal factor L from the L*D*L' factorization of A.
00038 *          E can also be regarded as the superdiagonal of the unit
00039 *          bidiagonal factor U from the U'*D*U factorization of A.
00040 *
00041 *  INFO    (output) INTEGER
00042 *          = 0: successful exit
00043 *          < 0: if INFO = -k, the k-th argument had an illegal value
00044 *          > 0: if INFO = k, the leading minor of order k is not
00045 *               positive definite; if k < N, the factorization could not
00046 *               be completed, while if k = N, the factorization was
00047 *               completed, but D(N) <= 0.
00048 *
00049 *  =====================================================================
00050 *
00051 *     .. Parameters ..
00052       DOUBLE PRECISION   ZERO
00053       PARAMETER          ( ZERO = 0.0D+0 )
00054 *     ..
00055 *     .. Local Scalars ..
00056       INTEGER            I, I4
00057       DOUBLE PRECISION   EII, EIR, F, G
00058 *     ..
00059 *     .. External Subroutines ..
00060       EXTERNAL           XERBLA
00061 *     ..
00062 *     .. Intrinsic Functions ..
00063       INTRINSIC          DBLE, DCMPLX, DIMAG, MOD
00064 *     ..
00065 *     .. Executable Statements ..
00066 *
00067 *     Test the input parameters.
00068 *
00069       INFO = 0
00070       IF( N.LT.0 ) THEN
00071          INFO = -1
00072          CALL XERBLA( 'ZPTTRF', -INFO )
00073          RETURN
00074       END IF
00075 *
00076 *     Quick return if possible
00077 *
00078       IF( N.EQ.0 )
00079      $   RETURN
00080 *
00081 *     Compute the L*D*L' (or U'*D*U) factorization of A.
00082 *
00083       I4 = MOD( N-1, 4 )
00084       DO 10 I = 1, I4
00085          IF( D( I ).LE.ZERO ) THEN
00086             INFO = I
00087             GO TO 30
00088          END IF
00089          EIR = DBLE( E( I ) )
00090          EII = DIMAG( E( I ) )
00091          F = EIR / D( I )
00092          G = EII / D( I )
00093          E( I ) = DCMPLX( F, G )
00094          D( I+1 ) = D( I+1 ) - F*EIR - G*EII
00095    10 CONTINUE
00096 *
00097       DO 20 I = I4 + 1, N - 4, 4
00098 *
00099 *        Drop out of the loop if d(i) <= 0: the matrix is not positive
00100 *        definite.
00101 *
00102          IF( D( I ).LE.ZERO ) THEN
00103             INFO = I
00104             GO TO 30
00105          END IF
00106 *
00107 *        Solve for e(i) and d(i+1).
00108 *
00109          EIR = DBLE( E( I ) )
00110          EII = DIMAG( E( I ) )
00111          F = EIR / D( I )
00112          G = EII / D( I )
00113          E( I ) = DCMPLX( F, G )
00114          D( I+1 ) = D( I+1 ) - F*EIR - G*EII
00115 *
00116          IF( D( I+1 ).LE.ZERO ) THEN
00117             INFO = I + 1
00118             GO TO 30
00119          END IF
00120 *
00121 *        Solve for e(i+1) and d(i+2).
00122 *
00123          EIR = DBLE( E( I+1 ) )
00124          EII = DIMAG( E( I+1 ) )
00125          F = EIR / D( I+1 )
00126          G = EII / D( I+1 )
00127          E( I+1 ) = DCMPLX( F, G )
00128          D( I+2 ) = D( I+2 ) - F*EIR - G*EII
00129 *
00130          IF( D( I+2 ).LE.ZERO ) THEN
00131             INFO = I + 2
00132             GO TO 30
00133          END IF
00134 *
00135 *        Solve for e(i+2) and d(i+3).
00136 *
00137          EIR = DBLE( E( I+2 ) )
00138          EII = DIMAG( E( I+2 ) )
00139          F = EIR / D( I+2 )
00140          G = EII / D( I+2 )
00141          E( I+2 ) = DCMPLX( F, G )
00142          D( I+3 ) = D( I+3 ) - F*EIR - G*EII
00143 *
00144          IF( D( I+3 ).LE.ZERO ) THEN
00145             INFO = I + 3
00146             GO TO 30
00147          END IF
00148 *
00149 *        Solve for e(i+3) and d(i+4).
00150 *
00151          EIR = DBLE( E( I+3 ) )
00152          EII = DIMAG( E( I+3 ) )
00153          F = EIR / D( I+3 )
00154          G = EII / D( I+3 )
00155          E( I+3 ) = DCMPLX( F, G )
00156          D( I+4 ) = D( I+4 ) - F*EIR - G*EII
00157    20 CONTINUE
00158 *
00159 *     Check d(n) for positive definiteness.
00160 *
00161       IF( D( N ).LE.ZERO )
00162      $   INFO = N
00163 *
00164    30 CONTINUE
00165       RETURN
00166 *
00167 *     End of ZPTTRF
00168 *
00169       END