LAPACK 3.3.1 Linear Algebra PACKage

# dpptrf.f

Go to the documentation of this file.
```00001       SUBROUTINE DPPTRF( UPLO, N, AP, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DPPTRF computes the Cholesky factorization of a real symmetric
00020 *  positive definite matrix A stored in packed format.
00021 *
00022 *  The factorization has the form
00023 *     A = U**T * U,  if UPLO = 'U', or
00024 *     A = L  * L**T,  if UPLO = 'L',
00025 *  where U is an upper triangular matrix and L is lower triangular.
00026 *
00027 *  Arguments
00028 *  =========
00029 *
00030 *  UPLO    (input) CHARACTER*1
00031 *          = 'U':  Upper triangle of A is stored;
00032 *          = 'L':  Lower triangle of A is stored.
00033 *
00034 *  N       (input) INTEGER
00035 *          The order of the matrix A.  N >= 0.
00036 *
00037 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00038 *          On entry, the upper or lower triangle of the symmetric matrix
00039 *          A, packed columnwise in a linear array.  The j-th column of A
00040 *          is stored in the array AP as follows:
00041 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00042 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00043 *          See below for further details.
00044 *
00045 *          On exit, if INFO = 0, the triangular factor U or L from the
00046 *          Cholesky factorization A = U**T*U or A = L*L**T, in the same
00047 *          storage format as A.
00048 *
00049 *  INFO    (output) INTEGER
00050 *          = 0:  successful exit
00051 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00052 *          > 0:  if INFO = i, the leading minor of order i is not
00053 *                positive definite, and the factorization could not be
00054 *                completed.
00055 *
00056 *  Further Details
00057 *  ======= =======
00058 *
00059 *  The packed storage scheme is illustrated by the following example
00060 *  when N = 4, UPLO = 'U':
00061 *
00062 *  Two-dimensional storage of the symmetric matrix A:
00063 *
00064 *     a11 a12 a13 a14
00065 *         a22 a23 a24
00066 *             a33 a34     (aij = aji)
00067 *                 a44
00068 *
00069 *  Packed storage of the upper triangle of A:
00070 *
00071 *  AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
00072 *
00073 *  =====================================================================
00074 *
00075 *     .. Parameters ..
00076       DOUBLE PRECISION   ONE, ZERO
00077       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00078 *     ..
00079 *     .. Local Scalars ..
00080       LOGICAL            UPPER
00081       INTEGER            J, JC, JJ
00082       DOUBLE PRECISION   AJJ
00083 *     ..
00084 *     .. External Functions ..
00085       LOGICAL            LSAME
00086       DOUBLE PRECISION   DDOT
00087       EXTERNAL           LSAME, DDOT
00088 *     ..
00089 *     .. External Subroutines ..
00090       EXTERNAL           DSCAL, DSPR, DTPSV, XERBLA
00091 *     ..
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          SQRT
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097 *     Test the input parameters.
00098 *
00099       INFO = 0
00100       UPPER = LSAME( UPLO, 'U' )
00101       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00102          INFO = -1
00103       ELSE IF( N.LT.0 ) THEN
00104          INFO = -2
00105       END IF
00106       IF( INFO.NE.0 ) THEN
00107          CALL XERBLA( 'DPPTRF', -INFO )
00108          RETURN
00109       END IF
00110 *
00111 *     Quick return if possible
00112 *
00113       IF( N.EQ.0 )
00114      \$   RETURN
00115 *
00116       IF( UPPER ) THEN
00117 *
00118 *        Compute the Cholesky factorization A = U**T*U.
00119 *
00120          JJ = 0
00121          DO 10 J = 1, N
00122             JC = JJ + 1
00123             JJ = JJ + J
00124 *
00125 *           Compute elements 1:J-1 of column J.
00126 *
00127             IF( J.GT.1 )
00128      \$         CALL DTPSV( 'Upper', 'Transpose', 'Non-unit', J-1, AP,
00129      \$                     AP( JC ), 1 )
00130 *
00131 *           Compute U(J,J) and test for non-positive-definiteness.
00132 *
00133             AJJ = AP( JJ ) - DDOT( J-1, AP( JC ), 1, AP( JC ), 1 )
00134             IF( AJJ.LE.ZERO ) THEN
00135                AP( JJ ) = AJJ
00136                GO TO 30
00137             END IF
00138             AP( JJ ) = SQRT( AJJ )
00139    10    CONTINUE
00140       ELSE
00141 *
00142 *        Compute the Cholesky factorization A = L*L**T.
00143 *
00144          JJ = 1
00145          DO 20 J = 1, N
00146 *
00147 *           Compute L(J,J) and test for non-positive-definiteness.
00148 *
00149             AJJ = AP( JJ )
00150             IF( AJJ.LE.ZERO ) THEN
00151                AP( JJ ) = AJJ
00152                GO TO 30
00153             END IF
00154             AJJ = SQRT( AJJ )
00155             AP( JJ ) = AJJ
00156 *
00157 *           Compute elements J+1:N of column J and update the trailing
00158 *           submatrix.
00159 *
00160             IF( J.LT.N ) THEN
00161                CALL DSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
00162                CALL DSPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
00163      \$                    AP( JJ+N-J+1 ) )
00164                JJ = JJ + N - J + 1
00165             END IF
00166    20    CONTINUE
00167       END IF
00168       GO TO 40
00169 *
00170    30 CONTINUE
00171       INFO = J
00172 *
00173    40 CONTINUE
00174       RETURN
00175 *
00176 *     End of DPPTRF
00177 *
00178       END
```