LAPACK 3.3.1
Linear Algebra PACKage

zpptri.f

Go to the documentation of this file.
00001       SUBROUTINE ZPPTRI( 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       COMPLEX*16         AP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZPPTRI computes the inverse of a complex Hermitian positive definite
00020 *  matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
00021 *  computed by ZPPTRF.
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) COMPLEX*16 array, dimension (N*(N+1)/2)
00034 *          On entry, the triangular factor U or L from the Cholesky
00035 *          factorization A = U**H*U or A = L*L**H, 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 (Hermitian)
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       COMPLEX*16         ZDOTC
00064       EXTERNAL           LSAME, ZDOTC
00065 *     ..
00066 *     .. External Subroutines ..
00067       EXTERNAL           XERBLA, ZDSCAL, ZHPR, ZTPMV, ZTPTRI
00068 *     ..
00069 *     .. Intrinsic Functions ..
00070       INTRINSIC          DBLE
00071 *     ..
00072 *     .. Executable Statements ..
00073 *
00074 *     Test the input parameters.
00075 *
00076       INFO = 0
00077       UPPER = LSAME( UPLO, 'U' )
00078       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00079          INFO = -1
00080       ELSE IF( N.LT.0 ) THEN
00081          INFO = -2
00082       END IF
00083       IF( INFO.NE.0 ) THEN
00084          CALL XERBLA( 'ZPPTRI', -INFO )
00085          RETURN
00086       END IF
00087 *
00088 *     Quick return if possible
00089 *
00090       IF( N.EQ.0 )
00091      $   RETURN
00092 *
00093 *     Invert the triangular Cholesky factor U or L.
00094 *
00095       CALL ZTPTRI( UPLO, 'Non-unit', N, AP, INFO )
00096       IF( INFO.GT.0 )
00097      $   RETURN
00098       IF( UPPER ) THEN
00099 *
00100 *        Compute the product inv(U) * inv(U)**H.
00101 *
00102          JJ = 0
00103          DO 10 J = 1, N
00104             JC = JJ + 1
00105             JJ = JJ + J
00106             IF( J.GT.1 )
00107      $         CALL ZHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
00108             AJJ = AP( JJ )
00109             CALL ZDSCAL( J, AJJ, AP( JC ), 1 )
00110    10    CONTINUE
00111 *
00112       ELSE
00113 *
00114 *        Compute the product inv(L)**H * inv(L).
00115 *
00116          JJ = 1
00117          DO 20 J = 1, N
00118             JJN = JJ + N - J + 1
00119             AP( JJ ) = DBLE( ZDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) )
00120             IF( J.LT.N )
00121      $         CALL ZTPMV( 'Lower', 'Conjugate transpose', 'Non-unit',
00122      $                     N-J, AP( JJN ), AP( JJ+1 ), 1 )
00123             JJ = JJN
00124    20    CONTINUE
00125       END IF
00126 *
00127       RETURN
00128 *
00129 *     End of ZPPTRI
00130 *
00131       END
 All Files Functions