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.

chpr.f 8.1 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279
  1. *> \brief \b CHPR
  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 CHPR(UPLO,N,ALPHA,X,INCX,AP)
  12. *
  13. * .. Scalar Arguments ..
  14. * REAL ALPHA
  15. * INTEGER INCX,N
  16. * CHARACTER UPLO
  17. * ..
  18. * .. Array Arguments ..
  19. * COMPLEX AP(*),X(*)
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> CHPR performs the hermitian rank 1 operation
  29. *>
  30. *> A := alpha*x*x**H + A,
  31. *>
  32. *> where alpha is a real scalar, x is an n element vector and A is an
  33. *> n by n hermitian matrix, supplied in packed form.
  34. *> \endverbatim
  35. *
  36. * Arguments:
  37. * ==========
  38. *
  39. *> \param[in] UPLO
  40. *> \verbatim
  41. *> UPLO is CHARACTER*1
  42. *> On entry, UPLO specifies whether the upper or lower
  43. *> triangular part of the matrix A is supplied in the packed
  44. *> array AP as follows:
  45. *>
  46. *> UPLO = 'U' or 'u' The upper triangular part of A is
  47. *> supplied in AP.
  48. *>
  49. *> UPLO = 'L' or 'l' The lower triangular part of A is
  50. *> supplied in AP.
  51. *> \endverbatim
  52. *>
  53. *> \param[in] N
  54. *> \verbatim
  55. *> N is INTEGER
  56. *> On entry, N specifies the order of the matrix A.
  57. *> N must be at least zero.
  58. *> \endverbatim
  59. *>
  60. *> \param[in] ALPHA
  61. *> \verbatim
  62. *> ALPHA is REAL
  63. *> On entry, ALPHA specifies the scalar alpha.
  64. *> \endverbatim
  65. *>
  66. *> \param[in] X
  67. *> \verbatim
  68. *> X is COMPLEX array, dimension at least
  69. *> ( 1 + ( n - 1 )*abs( INCX ) ).
  70. *> Before entry, the incremented array X must contain the n
  71. *> element vector x.
  72. *> \endverbatim
  73. *>
  74. *> \param[in] INCX
  75. *> \verbatim
  76. *> INCX is INTEGER
  77. *> On entry, INCX specifies the increment for the elements of
  78. *> X. INCX must not be zero.
  79. *> \endverbatim
  80. *>
  81. *> \param[in,out] AP
  82. *> \verbatim
  83. *> AP is COMPLEX array, dimension at least
  84. *> ( ( n*( n + 1 ) )/2 ).
  85. *> Before entry with UPLO = 'U' or 'u', the array AP must
  86. *> contain the upper triangular part of the hermitian matrix
  87. *> packed sequentially, column by column, so that AP( 1 )
  88. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  89. *> and a( 2, 2 ) respectively, and so on. On exit, the array
  90. *> AP is overwritten by the upper triangular part of the
  91. *> updated matrix.
  92. *> Before entry with UPLO = 'L' or 'l', the array AP must
  93. *> contain the lower triangular part of the hermitian matrix
  94. *> packed sequentially, column by column, so that AP( 1 )
  95. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  96. *> and a( 3, 1 ) respectively, and so on. On exit, the array
  97. *> AP is overwritten by the lower triangular part of the
  98. *> updated matrix.
  99. *> Note that the imaginary parts of the diagonal elements need
  100. *> not be set, they are assumed to be zero, and on exit they
  101. *> are set to zero.
  102. *> \endverbatim
  103. *
  104. * Authors:
  105. * ========
  106. *
  107. *> \author Univ. of Tennessee
  108. *> \author Univ. of California Berkeley
  109. *> \author Univ. of Colorado Denver
  110. *> \author NAG Ltd.
  111. *
  112. *> \date December 2016
  113. *
  114. *> \ingroup complex_blas_level2
  115. *
  116. *> \par Further Details:
  117. * =====================
  118. *>
  119. *> \verbatim
  120. *>
  121. *> Level 2 Blas routine.
  122. *>
  123. *> -- Written on 22-October-1986.
  124. *> Jack Dongarra, Argonne National Lab.
  125. *> Jeremy Du Croz, Nag Central Office.
  126. *> Sven Hammarling, Nag Central Office.
  127. *> Richard Hanson, Sandia National Labs.
  128. *> \endverbatim
  129. *>
  130. * =====================================================================
  131. SUBROUTINE CHPR(UPLO,N,ALPHA,X,INCX,AP)
  132. *
  133. * -- Reference BLAS level2 routine (version 3.7.0) --
  134. * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
  135. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  136. * December 2016
  137. *
  138. * .. Scalar Arguments ..
  139. REAL ALPHA
  140. INTEGER INCX,N
  141. CHARACTER UPLO
  142. * ..
  143. * .. Array Arguments ..
  144. COMPLEX AP(*),X(*)
  145. * ..
  146. *
  147. * =====================================================================
  148. *
  149. * .. Parameters ..
  150. COMPLEX ZERO
  151. PARAMETER (ZERO= (0.0E+0,0.0E+0))
  152. * ..
  153. * .. Local Scalars ..
  154. COMPLEX TEMP
  155. INTEGER I,INFO,IX,J,JX,K,KK,KX
  156. * ..
  157. * .. External Functions ..
  158. LOGICAL LSAME
  159. EXTERNAL LSAME
  160. * ..
  161. * .. External Subroutines ..
  162. EXTERNAL XERBLA
  163. * ..
  164. * .. Intrinsic Functions ..
  165. INTRINSIC CONJG,REAL
  166. * ..
  167. *
  168. * Test the input parameters.
  169. *
  170. INFO = 0
  171. IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  172. INFO = 1
  173. ELSE IF (N.LT.0) THEN
  174. INFO = 2
  175. ELSE IF (INCX.EQ.0) THEN
  176. INFO = 5
  177. END IF
  178. IF (INFO.NE.0) THEN
  179. CALL XERBLA('CHPR ',INFO)
  180. RETURN
  181. END IF
  182. *
  183. * Quick return if possible.
  184. *
  185. IF ((N.EQ.0) .OR. (ALPHA.EQ.REAL(ZERO))) RETURN
  186. *
  187. * Set the start point in X if the increment is not unity.
  188. *
  189. IF (INCX.LE.0) THEN
  190. KX = 1 - (N-1)*INCX
  191. ELSE IF (INCX.NE.1) THEN
  192. KX = 1
  193. END IF
  194. *
  195. * Start the operations. In this version the elements of the array AP
  196. * are accessed sequentially with one pass through AP.
  197. *
  198. KK = 1
  199. IF (LSAME(UPLO,'U')) THEN
  200. *
  201. * Form A when upper triangle is stored in AP.
  202. *
  203. IF (INCX.EQ.1) THEN
  204. DO 20 J = 1,N
  205. IF (X(J).NE.ZERO) THEN
  206. TEMP = ALPHA*CONJG(X(J))
  207. K = KK
  208. DO 10 I = 1,J - 1
  209. AP(K) = AP(K) + X(I)*TEMP
  210. K = K + 1
  211. 10 CONTINUE
  212. AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(J)*TEMP)
  213. ELSE
  214. AP(KK+J-1) = REAL(AP(KK+J-1))
  215. END IF
  216. KK = KK + J
  217. 20 CONTINUE
  218. ELSE
  219. JX = KX
  220. DO 40 J = 1,N
  221. IF (X(JX).NE.ZERO) THEN
  222. TEMP = ALPHA*CONJG(X(JX))
  223. IX = KX
  224. DO 30 K = KK,KK + J - 2
  225. AP(K) = AP(K) + X(IX)*TEMP
  226. IX = IX + INCX
  227. 30 CONTINUE
  228. AP(KK+J-1) = REAL(AP(KK+J-1)) + REAL(X(JX)*TEMP)
  229. ELSE
  230. AP(KK+J-1) = REAL(AP(KK+J-1))
  231. END IF
  232. JX = JX + INCX
  233. KK = KK + J
  234. 40 CONTINUE
  235. END IF
  236. ELSE
  237. *
  238. * Form A when lower triangle is stored in AP.
  239. *
  240. IF (INCX.EQ.1) THEN
  241. DO 60 J = 1,N
  242. IF (X(J).NE.ZERO) THEN
  243. TEMP = ALPHA*CONJG(X(J))
  244. AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(J))
  245. K = KK + 1
  246. DO 50 I = J + 1,N
  247. AP(K) = AP(K) + X(I)*TEMP
  248. K = K + 1
  249. 50 CONTINUE
  250. ELSE
  251. AP(KK) = REAL(AP(KK))
  252. END IF
  253. KK = KK + N - J + 1
  254. 60 CONTINUE
  255. ELSE
  256. JX = KX
  257. DO 80 J = 1,N
  258. IF (X(JX).NE.ZERO) THEN
  259. TEMP = ALPHA*CONJG(X(JX))
  260. AP(KK) = REAL(AP(KK)) + REAL(TEMP*X(JX))
  261. IX = JX
  262. DO 70 K = KK + 1,KK + N - J
  263. IX = IX + INCX
  264. AP(K) = AP(K) + X(IX)*TEMP
  265. 70 CONTINUE
  266. ELSE
  267. AP(KK) = REAL(AP(KK))
  268. END IF
  269. JX = JX + INCX
  270. KK = KK + N - J + 1
  271. 80 CONTINUE
  272. END IF
  273. END IF
  274. *
  275. RETURN
  276. *
  277. * End of CHPR .
  278. *
  279. END