You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

zpptri.f 5.2 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190
  1. *> \brief \b ZPPTRI
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download ZPPTRI + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpptri.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpptri.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpptri.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZPPTRI( UPLO, N, AP, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER UPLO
  25. * INTEGER INFO, N
  26. * ..
  27. * .. Array Arguments ..
  28. * COMPLEX*16 AP( * )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> ZPPTRI 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 ZPPTRF.
  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*16 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 complex16OTHERcomputational
  92. *
  93. * =====================================================================
  94. SUBROUTINE ZPPTRI( 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*16 AP( * )
  107. * ..
  108. *
  109. * =====================================================================
  110. *
  111. * .. Parameters ..
  112. DOUBLE PRECISION ONE
  113. PARAMETER ( ONE = 1.0D+0 )
  114. * ..
  115. * .. Local Scalars ..
  116. LOGICAL UPPER
  117. INTEGER J, JC, JJ, JJN
  118. DOUBLE PRECISION AJJ
  119. * ..
  120. * .. External Functions ..
  121. LOGICAL LSAME
  122. COMPLEX*16 ZDOTC
  123. EXTERNAL LSAME, ZDOTC
  124. * ..
  125. * .. External Subroutines ..
  126. EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPMV, ZTPTRI
  127. * ..
  128. * .. Intrinsic Functions ..
  129. INTRINSIC DBLE
  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( 'ZPPTRI', -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 ZTPTRI( 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 ZHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
  167. AJJ = AP( JJ )
  168. CALL ZDSCAL( 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 ) = DBLE( ZDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) )
  179. IF( J.LT.N )
  180. $ CALL ZTPMV( '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 ZPPTRI
  189. *
  190. END