LAPACK 3.3.0

dpptri.f

Go to the documentation of this file.
00001       SUBROUTINE DPPTRI( UPLO, N, AP, 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       CHARACTER          UPLO
00010       INTEGER            INFO, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DPPTRI computes the inverse of a real symmetric positive definite
00020 *  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
00021 *  computed by DPPTRF.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  UPLO    (input) CHARACTER*1
00027 *          = 'U':  Upper triangular factor is stored in AP;
00028 *          = 'L':  Lower triangular factor is stored in AP.
00029 *
00030 *  N       (input) INTEGER
00031 *          The order of the matrix A.  N >= 0.
00032 *
00033 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00034 *          On entry, the triangular factor U or L from the Cholesky
00035 *          factorization A = U**T*U or A = L*L**T, packed columnwise as
00036 *          a linear array.  The j-th column of U or L is stored in the
00037 *          array AP as follows:
00038 *          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
00039 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
00040 *
00041 *          On exit, the upper or lower triangle of the (symmetric)
00042 *          inverse of A, overwriting the input factor U or L.
00043 *
00044 *  INFO    (output) INTEGER
00045 *          = 0:  successful exit
00046 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00047 *          > 0:  if INFO = i, the (i,i) element of the factor U or L is
00048 *                zero, and the inverse could not be computed.
00049 *
00050 *  =====================================================================
00051 *
00052 *     .. Parameters ..
00053       DOUBLE PRECISION   ONE
00054       PARAMETER          ( ONE = 1.0D+0 )
00055 *     ..
00056 *     .. Local Scalars ..
00057       LOGICAL            UPPER
00058       INTEGER            J, JC, JJ, JJN
00059       DOUBLE PRECISION   AJJ
00060 *     ..
00061 *     .. External Functions ..
00062       LOGICAL            LSAME
00063       DOUBLE PRECISION   DDOT
00064       EXTERNAL           LSAME, DDOT
00065 *     ..
00066 *     .. External Subroutines ..
00067       EXTERNAL           DSCAL, DSPR, DTPMV, DTPTRI, XERBLA
00068 *     ..
00069 *     .. Executable Statements ..
00070 *
00071 *     Test the input parameters.
00072 *
00073       INFO = 0
00074       UPPER = LSAME( UPLO, 'U' )
00075       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00076          INFO = -1
00077       ELSE IF( N.LT.0 ) THEN
00078          INFO = -2
00079       END IF
00080       IF( INFO.NE.0 ) THEN
00081          CALL XERBLA( 'DPPTRI', -INFO )
00082          RETURN
00083       END IF
00084 *
00085 *     Quick return if possible
00086 *
00087       IF( N.EQ.0 )
00088      $   RETURN
00089 *
00090 *     Invert the triangular Cholesky factor U or L.
00091 *
00092       CALL DTPTRI( UPLO, 'Non-unit', N, AP, INFO )
00093       IF( INFO.GT.0 )
00094      $   RETURN
00095 *
00096       IF( UPPER ) THEN
00097 *
00098 *        Compute the product inv(U) * inv(U)'.
00099 *
00100          JJ = 0
00101          DO 10 J = 1, N
00102             JC = JJ + 1
00103             JJ = JJ + J
00104             IF( J.GT.1 )
00105      $         CALL DSPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
00106             AJJ = AP( JJ )
00107             CALL DSCAL( J, AJJ, AP( JC ), 1 )
00108    10    CONTINUE
00109 *
00110       ELSE
00111 *
00112 *        Compute the product inv(L)' * inv(L).
00113 *
00114          JJ = 1
00115          DO 20 J = 1, N
00116             JJN = JJ + N - J + 1
00117             AP( JJ ) = DDOT( N-J+1, AP( JJ ), 1, AP( JJ ), 1 )
00118             IF( J.LT.N )
00119      $         CALL DTPMV( 'Lower', 'Transpose', 'Non-unit', N-J,
00120      $                     AP( JJN ), AP( JJ+1 ), 1 )
00121             JJ = JJN
00122    20    CONTINUE
00123       END IF
00124 *
00125       RETURN
00126 *
00127 *     End of DPPTRI
00128 *
00129       END
 All Files Functions