SRC/cpptri.f

   1 *> \brief \b CPPTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CPPTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpptri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpptri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpptri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CPPTRI( UPLO, N, AP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX            AP( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> CPPTRI computes the inverse of a complex Hermitian positive definite
  38 *> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
  39 *> computed by CPPTRF.
  40 *> \endverbatim
  41 *
  42 *  Arguments:
  43 *  ==========
  44 *
  45 *> \param[in] UPLO
  46 *> \verbatim
  47 *>          UPLO is CHARACTER*1
  48 *>          = 'U':  Upper triangular factor is stored in AP;
  49 *>          = 'L':  Lower triangular factor is stored in AP.
  50 *> \endverbatim
  51 *>
  52 *> \param[in] N
  53 *> \verbatim
  54 *>          N is INTEGER
  55 *>          The order of the matrix A.  N >= 0.
  56 *> \endverbatim
  57 *>
  58 *> \param[in,out] AP
  59 *> \verbatim
  60 *>          AP is COMPLEX array, dimension (N*(N+1)/2)
  61 *>          On entry, the triangular factor U or L from the Cholesky
  62 *>          factorization A = U**H*U or A = L*L**H, packed columnwise as
  63 *>          a linear array.  The j-th column of U or L is stored in the
  64 *>          array AP as follows:
  65 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
  66 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
  67 *>
  68 *>          On exit, the upper or lower triangle of the (Hermitian)
  69 *>          inverse of A, overwriting the input factor U or L.
  70 *> \endverbatim
  71 *>
  72 *> \param[out] INFO
  73 *> \verbatim
  74 *>          INFO is INTEGER
  75 *>          = 0:  successful exit
  76 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  77 *>          > 0:  if INFO = i, the (i,i) element of the factor U or L is
  78 *>                zero, and the inverse could not be computed.
  79 *> \endverbatim
  80 *
  81 *  Authors:
  82 *  ========
  83 *
  84 *> \author Univ. of Tennessee
  85 *> \author Univ. of California Berkeley
  86 *> \author Univ. of Colorado Denver
  87 *> \author NAG Ltd.
  88 *
  89 *> \date November 2011
  90 *
  91 *> \ingroup complexOTHERcomputational
  92 *
  93 *  =====================================================================
  94       SUBROUTINE CPPTRI( UPLO, N, AP, INFO )
  95 *
  96 *  -- LAPACK computational routine (version 3.4.0) --
  97 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  98 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  99 *     November 2011
 100 *
 101 *     .. Scalar Arguments ..
 102       CHARACTER          UPLO
 103       INTEGER            INFO, N
 104 *     ..
 105 *     .. Array Arguments ..
 106       COMPLEX            AP( * )
 107 *     ..
 108 *
 109 *  =====================================================================
 110 *
 111 *     .. Parameters ..
 112       REAL               ONE
 113       PARAMETER          ( ONE = 1.0E+0 )
 114 *     ..
 115 *     .. Local Scalars ..
 116       LOGICAL            UPPER
 117       INTEGER            J, JC, JJ, JJN
 118       REAL               AJJ
 119 *     ..
 120 *     .. External Functions ..
 121       LOGICAL            LSAME
 122       COMPLEX            CDOTC
 123       EXTERNAL           LSAME, CDOTC
 124 *     ..
 125 *     .. External Subroutines ..
 126       EXTERNAL           CHPR, CSSCAL, CTPMV, CTPTRI, XERBLA
 127 *     ..
 128 *     .. Intrinsic Functions ..
 129       INTRINSIC          REAL
 130 *     ..
 131 *     .. Executable Statements ..
 132 *
 133 *     Test the input parameters.
 134 *
 135       INFO = 0
 136       UPPER = LSAME( UPLO, 'U' )
 137       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 138          INFO = -1
 139       ELSE IF( N.LT.0 ) THEN
 140          INFO = -2
 141       END IF
 142       IF( INFO.NE.0 ) THEN
 143          CALL XERBLA( 'CPPTRI', -INFO )
 144          RETURN
 145       END IF
 146 *
 147 *     Quick return if possible
 148 *
 149       IF( N.EQ.0 )
 150      $   RETURN
 151 *
 152 *     Invert the triangular Cholesky factor U or L.
 153 *
 154       CALL CTPTRI( UPLO, 'Non-unit', N, AP, INFO )
 155       IF( INFO.GT.0 )
 156      $   RETURN
 157       IF( UPPER ) THEN
 158 *
 159 *        Compute the product inv(U) * inv(U)**H.
 160 *
 161          JJ = 0
 162          DO 10 J = 1, N
 163             JC = JJ + 1
 164             JJ = JJ + J
 165             IF( J.GT.1 )
 166      $         CALL CHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
 167             AJJ = AP( JJ )
 168             CALL CSSCAL( J, AJJ, AP( JC ), 1 )
 169    10    CONTINUE
 170 *
 171       ELSE
 172 *
 173 *        Compute the product inv(L)**H * inv(L).
 174 *
 175          JJ = 1
 176          DO 20 J = 1, N
 177             JJN = JJ + N - J + 1
 178             AP( JJ ) = REAL( CDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) )
 179             IF( J.LT.N )
 180      $         CALL CTPMV( 'Lower', 'Conjugate transpose', 'Non-unit',
 181      $                     N-J, AP( JJN ), AP( JJ+1 ), 1 )
 182             JJ = JJN
 183    20    CONTINUE
 184       END IF
 185 *
 186       RETURN
 187 *
 188 *     End of CPPTRI
 189 *
 190       END