Return to dpptri.f CVS log

Up to [local] / rpl / lapack / lapack

Annotation of rpl/lapack/lapack/dpptri.f, revision 1.8

1.1       bertrand    1:       SUBROUTINE DPPTRI( UPLO, N, AP, INFO )
                      2: *
1.8     ! bertrand    3: *  -- LAPACK routine (version 3.3.1) --
1.1       bertrand    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                      5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8     ! bertrand    6: *  -- April 2011                                                      --
1.1       bertrand    7: *
                      8: *     .. Scalar Arguments ..
                      9:       CHARACTER          UPLO
                     10:       INTEGER            INFO, N
                     11: *     ..
                     12: *     .. Array Arguments ..
                     13:       DOUBLE PRECISION   AP( * )
                     14: *     ..
                     15: *
                     16: *  Purpose
                     17: *  =======
                     18: *
                     19: *  DPPTRI computes the inverse of a real symmetric positive definite
                     20: *  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
                     21: *  computed by DPPTRF.
                     22: *
                     23: *  Arguments
                     24: *  =========
                     25: *
                     26: *  UPLO    (input) CHARACTER*1
                     27: *          = 'U':  Upper triangular factor is stored in AP;
                     28: *          = 'L':  Lower triangular factor is stored in AP.
                     29: *
                     30: *  N       (input) INTEGER
                     31: *          The order of the matrix A.  N >= 0.
                     32: *
                     33: *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
                     34: *          On entry, the triangular factor U or L from the Cholesky
                     35: *          factorization A = U**T*U or A = L*L**T, packed columnwise as
                     36: *          a linear array.  The j-th column of U or L is stored in the
                     37: *          array AP as follows:
                     38: *          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
                     39: *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
                     40: *
                     41: *          On exit, the upper or lower triangle of the (symmetric)
                     42: *          inverse of A, overwriting the input factor U or L.
                     43: *
                     44: *  INFO    (output) INTEGER
                     45: *          = 0:  successful exit
                     46: *          < 0:  if INFO = -i, the i-th argument had an illegal value
                     47: *          > 0:  if INFO = i, the (i,i) element of the factor U or L is
                     48: *                zero, and the inverse could not be computed.
                     49: *
                     50: *  =====================================================================
                     51: *
                     52: *     .. Parameters ..
                     53:       DOUBLE PRECISION   ONE
                     54:       PARAMETER          ( ONE = 1.0D+0 )
                     55: *     ..
                     56: *     .. Local Scalars ..
                     57:       LOGICAL            UPPER
                     58:       INTEGER            J, JC, JJ, JJN
                     59:       DOUBLE PRECISION   AJJ
                     60: *     ..
                     61: *     .. External Functions ..
                     62:       LOGICAL            LSAME
                     63:       DOUBLE PRECISION   DDOT
                     64:       EXTERNAL           LSAME, DDOT
                     65: *     ..
                     66: *     .. External Subroutines ..
                     67:       EXTERNAL           DSCAL, DSPR, DTPMV, DTPTRI, XERBLA
                     68: *     ..
                     69: *     .. Executable Statements ..
                     70: *
                     71: *     Test the input parameters.
                     72: *
                     73:       INFO = 0
                     74:       UPPER = LSAME( UPLO, 'U' )
                     75:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
                     76:          INFO = -1
                     77:       ELSE IF( N.LT.0 ) THEN
                     78:          INFO = -2
                     79:       END IF
                     80:       IF( INFO.NE.0 ) THEN
                     81:          CALL XERBLA( 'DPPTRI', -INFO )
                     82:          RETURN
                     83:       END IF
                     84: *
                     85: *     Quick return if possible
                     86: *
                     87:       IF( N.EQ.0 )
                     88:      $   RETURN
                     89: *
                     90: *     Invert the triangular Cholesky factor U or L.
                     91: *
                     92:       CALL DTPTRI( UPLO, 'Non-unit', N, AP, INFO )
                     93:       IF( INFO.GT.0 )
                     94:      $   RETURN
                     95: *
                     96:       IF( UPPER ) THEN
                     97: *
1.8     ! bertrand   98: *        Compute the product inv(U) * inv(U)**T.
1.1       bertrand   99: *
                    100:          JJ = 0
                    101:          DO 10 J = 1, N
                    102:             JC = JJ + 1
                    103:             JJ = JJ + J
                    104:             IF( J.GT.1 )
                    105:      $         CALL DSPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
                    106:             AJJ = AP( JJ )
                    107:             CALL DSCAL( J, AJJ, AP( JC ), 1 )
                    108:    10    CONTINUE
                    109: *
                    110:       ELSE
                    111: *
1.8     ! bertrand  112: *        Compute the product inv(L)**T * inv(L).
1.1       bertrand  113: *
                    114:          JJ = 1
                    115:          DO 20 J = 1, N
                    116:             JJN = JJ + N - J + 1
                    117:             AP( JJ ) = DDOT( N-J+1, AP( JJ ), 1, AP( JJ ), 1 )
                    118:             IF( J.LT.N )
                    119:      $         CALL DTPMV( 'Lower', 'Transpose', 'Non-unit', N-J,
                    120:      $                     AP( JJN ), AP( JJ+1 ), 1 )
                    121:             JJ = JJN
                    122:    20    CONTINUE
                    123:       END IF
                    124: *
                    125:       RETURN
                    126: *
                    127: *     End of DPPTRI
                    128: *
                    129:       END