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
1.7.2