SRC/spptri.f

   1 *> \brief \b SPPTRI
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SPPTRI + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spptri.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spptri.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spptri.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SPPTRI( UPLO, N, AP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               AP( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SPPTRI computes the inverse of a real symmetric positive definite
  38 *> matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
  39 *> computed by SPPTRF.
  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 REAL array, dimension (N*(N+1)/2)
  61 *>          On entry, the triangular factor U or L from the Cholesky
  62 *>          factorization A = U**T*U or A = L*L**T, 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 (symmetric)
  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 realOTHERcomputational
  92 *
  93 *  =====================================================================
  94       SUBROUTINE SPPTRI( 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       REAL               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       REAL               SDOT
 123       EXTERNAL           LSAME, SDOT
 124 *     ..
 125 *     .. External Subroutines ..
 126       EXTERNAL           SSCAL, SSPR, STPMV, STPTRI, XERBLA
 127 *     ..
 128 *     .. Executable Statements ..
 129 *
 130 *     Test the input parameters.
 131 *
 132       INFO = 0
 133       UPPER = LSAME( UPLO, 'U' )
 134       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 135          INFO = -1
 136       ELSE IF( N.LT.0 ) THEN
 137          INFO = -2
 138       END IF
 139       IF( INFO.NE.0 ) THEN
 140          CALL XERBLA( 'SPPTRI', -INFO )
 141          RETURN
 142       END IF
 143 *
 144 *     Quick return if possible
 145 *
 146       IF( N.EQ.0 )
 147      $   RETURN
 148 *
 149 *     Invert the triangular Cholesky factor U or L.
 150 *
 151       CALL STPTRI( UPLO, 'Non-unit', N, AP, INFO )
 152       IF( INFO.GT.0 )
 153      $   RETURN
 154 *
 155       IF( UPPER ) THEN
 156 *
 157 *        Compute the product inv(U) * inv(U)**T.
 158 *
 159          JJ = 0
 160          DO 10 J = 1, N
 161             JC = JJ + 1
 162             JJ = JJ + J
 163             IF( J.GT.1 )
 164      $         CALL SSPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
 165             AJJ = AP( JJ )
 166             CALL SSCAL( J, AJJ, AP( JC ), 1 )
 167    10    CONTINUE
 168 *
 169       ELSE
 170 *
 171 *        Compute the product inv(L)**T * inv(L).
 172 *
 173          JJ = 1
 174          DO 20 J = 1, N
 175             JJN = JJ + N - J + 1
 176             AP( JJ ) = SDOT( N-J+1, AP( JJ ), 1, AP( JJ ), 1 )
 177             IF( J.LT.N )
 178      $         CALL STPMV( 'Lower', 'Transpose', 'Non-unit', N-J,
 179      $                     AP( JJN ), AP( JJ+1 ), 1 )
 180             JJ = JJN
 181    20    CONTINUE
 182       END IF
 183 *
 184       RETURN
 185 *
 186 *     End of SPPTRI
 187 *
 188       END