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.

cher2k.f 15 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442
  1. *> \brief \b CHER2K
  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 CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  12. *
  13. * .. Scalar Arguments ..
  14. * COMPLEX ALPHA
  15. * REAL BETA
  16. * INTEGER K,LDA,LDB,LDC,N
  17. * CHARACTER TRANS,UPLO
  18. * ..
  19. * .. Array Arguments ..
  20. * COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
  21. * ..
  22. *
  23. *
  24. *> \par Purpose:
  25. * =============
  26. *>
  27. *> \verbatim
  28. *>
  29. *> CHER2K performs one of the hermitian rank 2k operations
  30. *>
  31. *> C := alpha*A*B**H + conjg( alpha )*B*A**H + beta*C,
  32. *>
  33. *> or
  34. *>
  35. *> C := alpha*A**H*B + conjg( alpha )*B**H*A + beta*C,
  36. *>
  37. *> where alpha and beta are scalars with beta real, C is an n by n
  38. *> hermitian matrix and A and B are n by k matrices in the first case
  39. *> and k by n matrices in the second case.
  40. *> \endverbatim
  41. *
  42. * Arguments:
  43. * ==========
  44. *
  45. *> \param[in] UPLO
  46. *> \verbatim
  47. *> UPLO is CHARACTER*1
  48. *> On entry, UPLO specifies whether the upper or lower
  49. *> triangular part of the array C is to be referenced as
  50. *> follows:
  51. *>
  52. *> UPLO = 'U' or 'u' Only the upper triangular part of C
  53. *> is to be referenced.
  54. *>
  55. *> UPLO = 'L' or 'l' Only the lower triangular part of C
  56. *> is to be referenced.
  57. *> \endverbatim
  58. *>
  59. *> \param[in] TRANS
  60. *> \verbatim
  61. *> TRANS is CHARACTER*1
  62. *> On entry, TRANS specifies the operation to be performed as
  63. *> follows:
  64. *>
  65. *> TRANS = 'N' or 'n' C := alpha*A*B**H +
  66. *> conjg( alpha )*B*A**H +
  67. *> beta*C.
  68. *>
  69. *> TRANS = 'C' or 'c' C := alpha*A**H*B +
  70. *> conjg( alpha )*B**H*A +
  71. *> beta*C.
  72. *> \endverbatim
  73. *>
  74. *> \param[in] N
  75. *> \verbatim
  76. *> N is INTEGER
  77. *> On entry, N specifies the order of the matrix C. N must be
  78. *> at least zero.
  79. *> \endverbatim
  80. *>
  81. *> \param[in] K
  82. *> \verbatim
  83. *> K is INTEGER
  84. *> On entry with TRANS = 'N' or 'n', K specifies the number
  85. *> of columns of the matrices A and B, and on entry with
  86. *> TRANS = 'C' or 'c', K specifies the number of rows of the
  87. *> matrices A and B. K must be at least zero.
  88. *> \endverbatim
  89. *>
  90. *> \param[in] ALPHA
  91. *> \verbatim
  92. *> ALPHA is COMPLEX
  93. *> On entry, ALPHA specifies the scalar alpha.
  94. *> \endverbatim
  95. *>
  96. *> \param[in] A
  97. *> \verbatim
  98. *> A is COMPLEX array of DIMENSION ( LDA, ka ), where ka is
  99. *> k when TRANS = 'N' or 'n', and is n otherwise.
  100. *> Before entry with TRANS = 'N' or 'n', the leading n by k
  101. *> part of the array A must contain the matrix A, otherwise
  102. *> the leading k by n part of the array A must contain the
  103. *> matrix A.
  104. *> \endverbatim
  105. *>
  106. *> \param[in] LDA
  107. *> \verbatim
  108. *> LDA is INTEGER
  109. *> On entry, LDA specifies the first dimension of A as declared
  110. *> in the calling (sub) program. When TRANS = 'N' or 'n'
  111. *> then LDA must be at least max( 1, n ), otherwise LDA must
  112. *> be at least max( 1, k ).
  113. *> \endverbatim
  114. *>
  115. *> \param[in] B
  116. *> \verbatim
  117. *> B is COMPLEX array of DIMENSION ( LDB, kb ), where kb is
  118. *> k when TRANS = 'N' or 'n', and is n otherwise.
  119. *> Before entry with TRANS = 'N' or 'n', the leading n by k
  120. *> part of the array B must contain the matrix B, otherwise
  121. *> the leading k by n part of the array B must contain the
  122. *> matrix B.
  123. *> \endverbatim
  124. *>
  125. *> \param[in] LDB
  126. *> \verbatim
  127. *> LDB is INTEGER
  128. *> On entry, LDB specifies the first dimension of B as declared
  129. *> in the calling (sub) program. When TRANS = 'N' or 'n'
  130. *> then LDB must be at least max( 1, n ), otherwise LDB must
  131. *> be at least max( 1, k ).
  132. *> \endverbatim
  133. *>
  134. *> \param[in] BETA
  135. *> \verbatim
  136. *> BETA is REAL
  137. *> On entry, BETA specifies the scalar beta.
  138. *> \endverbatim
  139. *>
  140. *> \param[in,out] C
  141. *> \verbatim
  142. *> C is COMPLEX array of DIMENSION ( LDC, n ).
  143. *> Before entry with UPLO = 'U' or 'u', the leading n by n
  144. *> upper triangular part of the array C must contain the upper
  145. *> triangular part of the hermitian matrix and the strictly
  146. *> lower triangular part of C is not referenced. On exit, the
  147. *> upper triangular part of the array C is overwritten by the
  148. *> upper triangular part of the updated matrix.
  149. *> Before entry with UPLO = 'L' or 'l', the leading n by n
  150. *> lower triangular part of the array C must contain the lower
  151. *> triangular part of the hermitian matrix and the strictly
  152. *> upper triangular part of C is not referenced. On exit, the
  153. *> lower triangular part of the array C is overwritten by the
  154. *> lower triangular part of the updated matrix.
  155. *> Note that the imaginary parts of the diagonal elements need
  156. *> not be set, they are assumed to be zero, and on exit they
  157. *> are set to zero.
  158. *> \endverbatim
  159. *>
  160. *> \param[in] LDC
  161. *> \verbatim
  162. *> LDC is INTEGER
  163. *> On entry, LDC specifies the first dimension of C as declared
  164. *> in the calling (sub) program. LDC must be at least
  165. *> max( 1, n ).
  166. *> \endverbatim
  167. *
  168. * Authors:
  169. * ========
  170. *
  171. *> \author Univ. of Tennessee
  172. *> \author Univ. of California Berkeley
  173. *> \author Univ. of Colorado Denver
  174. *> \author NAG Ltd.
  175. *
  176. *> \date December 2016
  177. *
  178. *> \ingroup complex_blas_level3
  179. *
  180. *> \par Further Details:
  181. * =====================
  182. *>
  183. *> \verbatim
  184. *>
  185. *> Level 3 Blas routine.
  186. *>
  187. *> -- Written on 8-February-1989.
  188. *> Jack Dongarra, Argonne National Laboratory.
  189. *> Iain Duff, AERE Harwell.
  190. *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
  191. *> Sven Hammarling, Numerical Algorithms Group Ltd.
  192. *>
  193. *> -- Modified 8-Nov-93 to set C(J,J) to REAL( C(J,J) ) when BETA = 1.
  194. *> Ed Anderson, Cray Research Inc.
  195. *> \endverbatim
  196. *>
  197. * =====================================================================
  198. SUBROUTINE CHER2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  199. *
  200. * -- Reference BLAS level3 routine (version 3.7.0) --
  201. * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
  202. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  203. * December 2016
  204. *
  205. * .. Scalar Arguments ..
  206. COMPLEX ALPHA
  207. REAL BETA
  208. INTEGER K,LDA,LDB,LDC,N
  209. CHARACTER TRANS,UPLO
  210. * ..
  211. * .. Array Arguments ..
  212. COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
  213. * ..
  214. *
  215. * =====================================================================
  216. *
  217. * .. External Functions ..
  218. LOGICAL LSAME
  219. EXTERNAL LSAME
  220. * ..
  221. * .. External Subroutines ..
  222. EXTERNAL XERBLA
  223. * ..
  224. * .. Intrinsic Functions ..
  225. INTRINSIC CONJG,MAX,REAL
  226. * ..
  227. * .. Local Scalars ..
  228. COMPLEX TEMP1,TEMP2
  229. INTEGER I,INFO,J,L,NROWA
  230. LOGICAL UPPER
  231. * ..
  232. * .. Parameters ..
  233. REAL ONE
  234. PARAMETER (ONE=1.0E+0)
  235. COMPLEX ZERO
  236. PARAMETER (ZERO= (0.0E+0,0.0E+0))
  237. * ..
  238. *
  239. * Test the input parameters.
  240. *
  241. IF (LSAME(TRANS,'N')) THEN
  242. NROWA = N
  243. ELSE
  244. NROWA = K
  245. END IF
  246. UPPER = LSAME(UPLO,'U')
  247. *
  248. INFO = 0
  249. IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
  250. INFO = 1
  251. ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
  252. + (.NOT.LSAME(TRANS,'C'))) THEN
  253. INFO = 2
  254. ELSE IF (N.LT.0) THEN
  255. INFO = 3
  256. ELSE IF (K.LT.0) THEN
  257. INFO = 4
  258. ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
  259. INFO = 7
  260. ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
  261. INFO = 9
  262. ELSE IF (LDC.LT.MAX(1,N)) THEN
  263. INFO = 12
  264. END IF
  265. IF (INFO.NE.0) THEN
  266. CALL XERBLA('CHER2K',INFO)
  267. RETURN
  268. END IF
  269. *
  270. * Quick return if possible.
  271. *
  272. IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
  273. + (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
  274. *
  275. * And when alpha.eq.zero.
  276. *
  277. IF (ALPHA.EQ.ZERO) THEN
  278. IF (UPPER) THEN
  279. IF (BETA.EQ.REAL(ZERO)) THEN
  280. DO 20 J = 1,N
  281. DO 10 I = 1,J
  282. C(I,J) = ZERO
  283. 10 CONTINUE
  284. 20 CONTINUE
  285. ELSE
  286. DO 40 J = 1,N
  287. DO 30 I = 1,J - 1
  288. C(I,J) = BETA*C(I,J)
  289. 30 CONTINUE
  290. C(J,J) = BETA*REAL(C(J,J))
  291. 40 CONTINUE
  292. END IF
  293. ELSE
  294. IF (BETA.EQ.REAL(ZERO)) THEN
  295. DO 60 J = 1,N
  296. DO 50 I = J,N
  297. C(I,J) = ZERO
  298. 50 CONTINUE
  299. 60 CONTINUE
  300. ELSE
  301. DO 80 J = 1,N
  302. C(J,J) = BETA*REAL(C(J,J))
  303. DO 70 I = J + 1,N
  304. C(I,J) = BETA*C(I,J)
  305. 70 CONTINUE
  306. 80 CONTINUE
  307. END IF
  308. END IF
  309. RETURN
  310. END IF
  311. *
  312. * Start the operations.
  313. *
  314. IF (LSAME(TRANS,'N')) THEN
  315. *
  316. * Form C := alpha*A*B**H + conjg( alpha )*B*A**H +
  317. * C.
  318. *
  319. IF (UPPER) THEN
  320. DO 130 J = 1,N
  321. IF (BETA.EQ.REAL(ZERO)) THEN
  322. DO 90 I = 1,J
  323. C(I,J) = ZERO
  324. 90 CONTINUE
  325. ELSE IF (BETA.NE.ONE) THEN
  326. DO 100 I = 1,J - 1
  327. C(I,J) = BETA*C(I,J)
  328. 100 CONTINUE
  329. C(J,J) = BETA*REAL(C(J,J))
  330. ELSE
  331. C(J,J) = REAL(C(J,J))
  332. END IF
  333. DO 120 L = 1,K
  334. IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
  335. TEMP1 = ALPHA*CONJG(B(J,L))
  336. TEMP2 = CONJG(ALPHA*A(J,L))
  337. DO 110 I = 1,J - 1
  338. C(I,J) = C(I,J) + A(I,L)*TEMP1 +
  339. + B(I,L)*TEMP2
  340. 110 CONTINUE
  341. C(J,J) = REAL(C(J,J)) +
  342. + REAL(A(J,L)*TEMP1+B(J,L)*TEMP2)
  343. END IF
  344. 120 CONTINUE
  345. 130 CONTINUE
  346. ELSE
  347. DO 180 J = 1,N
  348. IF (BETA.EQ.REAL(ZERO)) THEN
  349. DO 140 I = J,N
  350. C(I,J) = ZERO
  351. 140 CONTINUE
  352. ELSE IF (BETA.NE.ONE) THEN
  353. DO 150 I = J + 1,N
  354. C(I,J) = BETA*C(I,J)
  355. 150 CONTINUE
  356. C(J,J) = BETA*REAL(C(J,J))
  357. ELSE
  358. C(J,J) = REAL(C(J,J))
  359. END IF
  360. DO 170 L = 1,K
  361. IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
  362. TEMP1 = ALPHA*CONJG(B(J,L))
  363. TEMP2 = CONJG(ALPHA*A(J,L))
  364. DO 160 I = J + 1,N
  365. C(I,J) = C(I,J) + A(I,L)*TEMP1 +
  366. + B(I,L)*TEMP2
  367. 160 CONTINUE
  368. C(J,J) = REAL(C(J,J)) +
  369. + REAL(A(J,L)*TEMP1+B(J,L)*TEMP2)
  370. END IF
  371. 170 CONTINUE
  372. 180 CONTINUE
  373. END IF
  374. ELSE
  375. *
  376. * Form C := alpha*A**H*B + conjg( alpha )*B**H*A +
  377. * C.
  378. *
  379. IF (UPPER) THEN
  380. DO 210 J = 1,N
  381. DO 200 I = 1,J
  382. TEMP1 = ZERO
  383. TEMP2 = ZERO
  384. DO 190 L = 1,K
  385. TEMP1 = TEMP1 + CONJG(A(L,I))*B(L,J)
  386. TEMP2 = TEMP2 + CONJG(B(L,I))*A(L,J)
  387. 190 CONTINUE
  388. IF (I.EQ.J) THEN
  389. IF (BETA.EQ.REAL(ZERO)) THEN
  390. C(J,J) = REAL(ALPHA*TEMP1+
  391. + CONJG(ALPHA)*TEMP2)
  392. ELSE
  393. C(J,J) = BETA*REAL(C(J,J)) +
  394. + REAL(ALPHA*TEMP1+
  395. + CONJG(ALPHA)*TEMP2)
  396. END IF
  397. ELSE
  398. IF (BETA.EQ.REAL(ZERO)) THEN
  399. C(I,J) = ALPHA*TEMP1 + CONJG(ALPHA)*TEMP2
  400. ELSE
  401. C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
  402. + CONJG(ALPHA)*TEMP2
  403. END IF
  404. END IF
  405. 200 CONTINUE
  406. 210 CONTINUE
  407. ELSE
  408. DO 240 J = 1,N
  409. DO 230 I = J,N
  410. TEMP1 = ZERO
  411. TEMP2 = ZERO
  412. DO 220 L = 1,K
  413. TEMP1 = TEMP1 + CONJG(A(L,I))*B(L,J)
  414. TEMP2 = TEMP2 + CONJG(B(L,I))*A(L,J)
  415. 220 CONTINUE
  416. IF (I.EQ.J) THEN
  417. IF (BETA.EQ.REAL(ZERO)) THEN
  418. C(J,J) = REAL(ALPHA*TEMP1+
  419. + CONJG(ALPHA)*TEMP2)
  420. ELSE
  421. C(J,J) = BETA*REAL(C(J,J)) +
  422. + REAL(ALPHA*TEMP1+
  423. + CONJG(ALPHA)*TEMP2)
  424. END IF
  425. ELSE
  426. IF (BETA.EQ.REAL(ZERO)) THEN
  427. C(I,J) = ALPHA*TEMP1 + CONJG(ALPHA)*TEMP2
  428. ELSE
  429. C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
  430. + CONJG(ALPHA)*TEMP2
  431. END IF
  432. END IF
  433. 230 CONTINUE
  434. 240 CONTINUE
  435. END IF
  436. END IF
  437. *
  438. RETURN
  439. *
  440. * End of CHER2K.
  441. *
  442. END