LAPACK 3.3.0

dpttrf.f

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