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.

cpotrf.f 6.7 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243
  1. C> \brief \b CPOTRF VARIANT: right looking block version of the algorithm, calling Level 3 BLAS.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. * Definition:
  9. * ===========
  10. *
  11. * SUBROUTINE CPOTRF ( UPLO, N, A, LDA, INFO )
  12. *
  13. * .. Scalar Arguments ..
  14. * CHARACTER UPLO
  15. * INTEGER INFO, LDA, N
  16. * ..
  17. * .. Array Arguments ..
  18. * COMPLEX A( LDA, * )
  19. * ..
  20. *
  21. * Purpose
  22. * =======
  23. *
  24. C>\details \b Purpose:
  25. C>\verbatim
  26. C>
  27. C> CPOTRF computes the Cholesky factorization of a real Hermitian
  28. C> positive definite matrix A.
  29. C>
  30. C> The factorization has the form
  31. C> A = U**H * U, if UPLO = 'U', or
  32. C> A = L * L**H, if UPLO = 'L',
  33. C> where U is an upper triangular matrix and L is lower triangular.
  34. C>
  35. C> This is the right looking block version of the algorithm, calling Level 3 BLAS.
  36. C>
  37. C>\endverbatim
  38. *
  39. * Arguments:
  40. * ==========
  41. *
  42. C> \param[in] UPLO
  43. C> \verbatim
  44. C> UPLO is CHARACTER*1
  45. C> = 'U': Upper triangle of A is stored;
  46. C> = 'L': Lower triangle of A is stored.
  47. C> \endverbatim
  48. C>
  49. C> \param[in] N
  50. C> \verbatim
  51. C> N is INTEGER
  52. C> The order of the matrix A. N >= 0.
  53. C> \endverbatim
  54. C>
  55. C> \param[in,out] A
  56. C> \verbatim
  57. C> A is COMPLEX array, dimension (LDA,N)
  58. C> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
  59. C> N-by-N upper triangular part of A contains the upper
  60. C> triangular part of the matrix A, and the strictly lower
  61. C> triangular part of A is not referenced. If UPLO = 'L', the
  62. C> leading N-by-N lower triangular part of A contains the lower
  63. C> triangular part of the matrix A, and the strictly upper
  64. C> triangular part of A is not referenced.
  65. C> \endverbatim
  66. C> \verbatim
  67. C> On exit, if INFO = 0, the factor U or L from the Cholesky
  68. C> factorization A = U**H*U or A = L*L**H.
  69. C> \endverbatim
  70. C>
  71. C> \param[in] LDA
  72. C> \verbatim
  73. C> LDA is INTEGER
  74. C> The leading dimension of the array A. LDA >= max(1,N).
  75. C> \endverbatim
  76. C>
  77. C> \param[out] INFO
  78. C> \verbatim
  79. C> INFO is INTEGER
  80. C> = 0: successful exit
  81. C> < 0: if INFO = -i, the i-th argument had an illegal value
  82. C> > 0: if INFO = i, the leading minor of order i is not
  83. C> positive definite, and the factorization could not be
  84. C> completed.
  85. C> \endverbatim
  86. C>
  87. *
  88. * Authors:
  89. * ========
  90. *
  91. C> \author Univ. of Tennessee
  92. C> \author Univ. of California Berkeley
  93. C> \author Univ. of Colorado Denver
  94. C> \author NAG Ltd.
  95. *
  96. C> \date November 2011
  97. *
  98. C> \ingroup variantsPOcomputational
  99. *
  100. * =====================================================================
  101. SUBROUTINE CPOTRF ( UPLO, N, A, LDA, INFO )
  102. *
  103. * -- LAPACK computational routine (version 3.1) --
  104. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  105. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  106. * November 2011
  107. *
  108. * .. Scalar Arguments ..
  109. CHARACTER UPLO
  110. INTEGER INFO, LDA, N
  111. * ..
  112. * .. Array Arguments ..
  113. COMPLEX A( LDA, * )
  114. * ..
  115. *
  116. * =====================================================================
  117. *
  118. * .. Parameters ..
  119. REAL ONE
  120. COMPLEX CONE
  121. PARAMETER ( ONE = 1.0E+0, CONE = ( 1.0E+0, 0.0E+0 ) )
  122. * ..
  123. * .. Local Scalars ..
  124. LOGICAL UPPER
  125. INTEGER J, JB, NB
  126. * ..
  127. * .. External Functions ..
  128. LOGICAL LSAME
  129. INTEGER ILAENV
  130. EXTERNAL LSAME, ILAENV
  131. * ..
  132. * .. External Subroutines ..
  133. EXTERNAL CGEMM, CPOTF2, CHERK, CTRSM, XERBLA
  134. * ..
  135. * .. Intrinsic Functions ..
  136. INTRINSIC MAX, MIN
  137. * ..
  138. * .. Executable Statements ..
  139. *
  140. * Test the input parameters.
  141. *
  142. INFO = 0
  143. UPPER = LSAME( UPLO, 'U' )
  144. IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  145. INFO = -1
  146. ELSE IF( N.LT.0 ) THEN
  147. INFO = -2
  148. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  149. INFO = -4
  150. END IF
  151. IF( INFO.NE.0 ) THEN
  152. CALL XERBLA( 'CPOTRF', -INFO )
  153. RETURN
  154. END IF
  155. *
  156. * Quick return if possible
  157. *
  158. IF( N.EQ.0 )
  159. $ RETURN
  160. *
  161. * Determine the block size for this environment.
  162. *
  163. NB = ILAENV( 1, 'CPOTRF', UPLO, N, -1, -1, -1 )
  164. IF( NB.LE.1 .OR. NB.GE.N ) THEN
  165. *
  166. * Use unblocked code.
  167. *
  168. CALL CPOTF2( UPLO, N, A, LDA, INFO )
  169. ELSE
  170. *
  171. * Use blocked code.
  172. *
  173. IF( UPPER ) THEN
  174. *
  175. * Compute the Cholesky factorization A = U'*U.
  176. *
  177. DO 10 J = 1, N, NB
  178. *
  179. * Update and factorize the current diagonal block and test
  180. * for non-positive-definiteness.
  181. *
  182. JB = MIN( NB, N-J+1 )
  183. CALL CPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
  184. IF( INFO.NE.0 )
  185. $ GO TO 30
  186. IF( J+JB.LE.N ) THEN
  187. *
  188. * Updating the trailing submatrix.
  189. *
  190. CALL CTRSM( 'Left', 'Upper', 'Conjugate Transpose',
  191. $ 'Non-unit', JB, N-J-JB+1, CONE, A( J, J ),
  192. $ LDA, A( J, J+JB ), LDA )
  193. CALL CHERK( 'Upper', 'Conjugate transpose', N-J-JB+1,
  194. $ JB, -ONE, A( J, J+JB ), LDA,
  195. $ ONE, A( J+JB, J+JB ), LDA )
  196. END IF
  197. 10 CONTINUE
  198. *
  199. ELSE
  200. *
  201. * Compute the Cholesky factorization A = L*L'.
  202. *
  203. DO 20 J = 1, N, NB
  204. *
  205. * Update and factorize the current diagonal block and test
  206. * for non-positive-definiteness.
  207. *
  208. JB = MIN( NB, N-J+1 )
  209. CALL CPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
  210. IF( INFO.NE.0 )
  211. $ GO TO 30
  212. IF( J+JB.LE.N ) THEN
  213. *
  214. * Updating the trailing submatrix.
  215. *
  216. CALL CTRSM( 'Right', 'Lower', 'Conjugate Transpose',
  217. $ 'Non-unit', N-J-JB+1, JB, CONE, A( J, J ),
  218. $ LDA, A( J+JB, J ), LDA )
  219. CALL CHERK( 'Lower', 'No Transpose', N-J-JB+1, JB,
  220. $ -ONE, A( J+JB, J ), LDA,
  221. $ ONE, A( J+JB, J+JB ), LDA )
  222. END IF
  223. 20 CONTINUE
  224. END IF
  225. END IF
  226. GO TO 40
  227. *
  228. 30 CONTINUE
  229. INFO = INFO + J - 1
  230. *
  231. 40 CONTINUE
  232. RETURN
  233. *
  234. * End of CPOTRF
  235. *
  236. END