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.

chpr2.f 9.6 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318
  1. *> \brief \b CHPR2
  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 CHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  12. *
  13. * .. Scalar Arguments ..
  14. * COMPLEX ALPHA
  15. * INTEGER INCX,INCY,N
  16. * CHARACTER UPLO
  17. * ..
  18. * .. Array Arguments ..
  19. * COMPLEX AP(*),X(*),Y(*)
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> CHPR2 performs the hermitian rank 2 operation
  29. *>
  30. *> A := alpha*x*y**H + conjg( alpha )*y*x**H + A,
  31. *>
  32. *> where alpha is a scalar, x and y are n element vectors 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 COMPLEX
  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] Y
  82. *> \verbatim
  83. *> Y is COMPLEX array, dimension at least
  84. *> ( 1 + ( n - 1 )*abs( INCY ) ).
  85. *> Before entry, the incremented array Y must contain the n
  86. *> element vector y.
  87. *> \endverbatim
  88. *>
  89. *> \param[in] INCY
  90. *> \verbatim
  91. *> INCY is INTEGER
  92. *> On entry, INCY specifies the increment for the elements of
  93. *> Y. INCY must not be zero.
  94. *> \endverbatim
  95. *>
  96. *> \param[in,out] AP
  97. *> \verbatim
  98. *> AP is COMPLEX array, dimension at least
  99. *> ( ( n*( n + 1 ) )/2 ).
  100. *> Before entry with UPLO = 'U' or 'u', the array AP must
  101. *> contain the upper triangular part of the hermitian matrix
  102. *> packed sequentially, column by column, so that AP( 1 )
  103. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  104. *> and a( 2, 2 ) respectively, and so on. On exit, the array
  105. *> AP is overwritten by the upper triangular part of the
  106. *> updated matrix.
  107. *> Before entry with UPLO = 'L' or 'l', the array AP must
  108. *> contain the lower triangular part of the hermitian matrix
  109. *> packed sequentially, column by column, so that AP( 1 )
  110. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  111. *> and a( 3, 1 ) respectively, and so on. On exit, the array
  112. *> AP is overwritten by the lower triangular part of the
  113. *> updated matrix.
  114. *> Note that the imaginary parts of the diagonal elements need
  115. *> not be set, they are assumed to be zero, and on exit they
  116. *> are set to zero.
  117. *> \endverbatim
  118. *
  119. * Authors:
  120. * ========
  121. *
  122. *> \author Univ. of Tennessee
  123. *> \author Univ. of California Berkeley
  124. *> \author Univ. of Colorado Denver
  125. *> \author NAG Ltd.
  126. *
  127. *> \date December 2016
  128. *
  129. *> \ingroup complex_blas_level2
  130. *
  131. *> \par Further Details:
  132. * =====================
  133. *>
  134. *> \verbatim
  135. *>
  136. *> Level 2 Blas routine.
  137. *>
  138. *> -- Written on 22-October-1986.
  139. *> Jack Dongarra, Argonne National Lab.
  140. *> Jeremy Du Croz, Nag Central Office.
  141. *> Sven Hammarling, Nag Central Office.
  142. *> Richard Hanson, Sandia National Labs.
  143. *> \endverbatim
  144. *>
  145. * =====================================================================
  146. SUBROUTINE CHPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  147. *
  148. * -- Reference BLAS level2 routine (version 3.7.0) --
  149. * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
  150. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  151. * December 2016
  152. *
  153. * .. Scalar Arguments ..
  154. COMPLEX ALPHA
  155. INTEGER INCX,INCY,N
  156. CHARACTER UPLO
  157. * ..
  158. * .. Array Arguments ..
  159. COMPLEX AP(*),X(*),Y(*)
  160. * ..
  161. *
  162. * =====================================================================
  163. *
  164. * .. Parameters ..
  165. COMPLEX ZERO
  166. PARAMETER (ZERO= (0.0E+0,0.0E+0))
  167. * ..
  168. * .. Local Scalars ..
  169. COMPLEX TEMP1,TEMP2
  170. INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
  171. * ..
  172. * .. External Functions ..
  173. LOGICAL LSAME
  174. EXTERNAL LSAME
  175. * ..
  176. * .. External Subroutines ..
  177. EXTERNAL XERBLA
  178. * ..
  179. * .. Intrinsic Functions ..
  180. INTRINSIC CONJG,REAL
  181. * ..
  182. *
  183. * Test the input parameters.
  184. *
  185. INFO = 0
  186. IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  187. INFO = 1
  188. ELSE IF (N.LT.0) THEN
  189. INFO = 2
  190. ELSE IF (INCX.EQ.0) THEN
  191. INFO = 5
  192. ELSE IF (INCY.EQ.0) THEN
  193. INFO = 7
  194. END IF
  195. IF (INFO.NE.0) THEN
  196. CALL XERBLA('CHPR2 ',INFO)
  197. RETURN
  198. END IF
  199. *
  200. * Quick return if possible.
  201. *
  202. IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  203. *
  204. * Set up the start points in X and Y if the increments are not both
  205. * unity.
  206. *
  207. IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
  208. IF (INCX.GT.0) THEN
  209. KX = 1
  210. ELSE
  211. KX = 1 - (N-1)*INCX
  212. END IF
  213. IF (INCY.GT.0) THEN
  214. KY = 1
  215. ELSE
  216. KY = 1 - (N-1)*INCY
  217. END IF
  218. JX = KX
  219. JY = KY
  220. END IF
  221. *
  222. * Start the operations. In this version the elements of the array AP
  223. * are accessed sequentially with one pass through AP.
  224. *
  225. KK = 1
  226. IF (LSAME(UPLO,'U')) THEN
  227. *
  228. * Form A when upper triangle is stored in AP.
  229. *
  230. IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  231. DO 20 J = 1,N
  232. IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  233. TEMP1 = ALPHA*CONJG(Y(J))
  234. TEMP2 = CONJG(ALPHA*X(J))
  235. K = KK
  236. DO 10 I = 1,J - 1
  237. AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  238. K = K + 1
  239. 10 CONTINUE
  240. AP(KK+J-1) = REAL(AP(KK+J-1)) +
  241. + REAL(X(J)*TEMP1+Y(J)*TEMP2)
  242. ELSE
  243. AP(KK+J-1) = REAL(AP(KK+J-1))
  244. END IF
  245. KK = KK + J
  246. 20 CONTINUE
  247. ELSE
  248. DO 40 J = 1,N
  249. IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  250. TEMP1 = ALPHA*CONJG(Y(JY))
  251. TEMP2 = CONJG(ALPHA*X(JX))
  252. IX = KX
  253. IY = KY
  254. DO 30 K = KK,KK + J - 2
  255. AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  256. IX = IX + INCX
  257. IY = IY + INCY
  258. 30 CONTINUE
  259. AP(KK+J-1) = REAL(AP(KK+J-1)) +
  260. + REAL(X(JX)*TEMP1+Y(JY)*TEMP2)
  261. ELSE
  262. AP(KK+J-1) = REAL(AP(KK+J-1))
  263. END IF
  264. JX = JX + INCX
  265. JY = JY + INCY
  266. KK = KK + J
  267. 40 CONTINUE
  268. END IF
  269. ELSE
  270. *
  271. * Form A when lower triangle is stored in AP.
  272. *
  273. IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  274. DO 60 J = 1,N
  275. IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  276. TEMP1 = ALPHA*CONJG(Y(J))
  277. TEMP2 = CONJG(ALPHA*X(J))
  278. AP(KK) = REAL(AP(KK)) +
  279. + REAL(X(J)*TEMP1+Y(J)*TEMP2)
  280. K = KK + 1
  281. DO 50 I = J + 1,N
  282. AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  283. K = K + 1
  284. 50 CONTINUE
  285. ELSE
  286. AP(KK) = REAL(AP(KK))
  287. END IF
  288. KK = KK + N - J + 1
  289. 60 CONTINUE
  290. ELSE
  291. DO 80 J = 1,N
  292. IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  293. TEMP1 = ALPHA*CONJG(Y(JY))
  294. TEMP2 = CONJG(ALPHA*X(JX))
  295. AP(KK) = REAL(AP(KK)) +
  296. + REAL(X(JX)*TEMP1+Y(JY)*TEMP2)
  297. IX = JX
  298. IY = JY
  299. DO 70 K = KK + 1,KK + N - J
  300. IX = IX + INCX
  301. IY = IY + INCY
  302. AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  303. 70 CONTINUE
  304. ELSE
  305. AP(KK) = REAL(AP(KK))
  306. END IF
  307. JX = JX + INCX
  308. JY = JY + INCY
  309. KK = KK + N - J + 1
  310. 80 CONTINUE
  311. END IF
  312. END IF
  313. *
  314. RETURN
  315. *
  316. * End of CHPR2 .
  317. *
  318. END