zher2k.c 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671
  1. /* zher2k.f -- translated by f2c (version 20061008).
  2. You must link the resulting object file with libf2c:
  3. on Microsoft Windows system, link with libf2c.lib;
  4. on Linux or Unix systems, link with .../path/to/libf2c.a -lm
  5. or, if you install libf2c.a in a standard place, with -lf2c -lm
  6. -- in that order, at the end of the command line, as in
  7. cc *.o -lf2c -lm
  8. Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
  9. http://www.netlib.org/f2c/libf2c.zip
  10. */
  11. #include "f2c.h"
  12. #include "blaswrap.h"
  13. /* Subroutine */ int zher2k_(char *uplo, char *trans, integer *n, integer *k,
  14. doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *
  15. b, integer *ldb, doublereal *beta, doublecomplex *c__, integer *ldc)
  16. {
  17. /* System generated locals */
  18. integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2,
  19. i__3, i__4, i__5, i__6, i__7;
  20. doublereal d__1;
  21. doublecomplex z__1, z__2, z__3, z__4, z__5, z__6;
  22. /* Builtin functions */
  23. void d_cnjg(doublecomplex *, doublecomplex *);
  24. /* Local variables */
  25. integer i__, j, l, info;
  26. doublecomplex temp1, temp2;
  27. extern logical lsame_(char *, char *);
  28. integer nrowa;
  29. logical upper;
  30. extern /* Subroutine */ int xerbla_(char *, integer *);
  31. /* .. Scalar Arguments .. */
  32. /* .. */
  33. /* .. Array Arguments .. */
  34. /* .. */
  35. /* Purpose */
  36. /* ======= */
  37. /* ZHER2K performs one of the hermitian rank 2k operations */
  38. /* C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C, */
  39. /* or */
  40. /* C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C, */
  41. /* where alpha and beta are scalars with beta real, C is an n by n */
  42. /* hermitian matrix and A and B are n by k matrices in the first case */
  43. /* and k by n matrices in the second case. */
  44. /* Arguments */
  45. /* ========== */
  46. /* UPLO - CHARACTER*1. */
  47. /* On entry, UPLO specifies whether the upper or lower */
  48. /* triangular part of the array C is to be referenced as */
  49. /* follows: */
  50. /* UPLO = 'U' or 'u' Only the upper triangular part of C */
  51. /* is to be referenced. */
  52. /* UPLO = 'L' or 'l' Only the lower triangular part of C */
  53. /* is to be referenced. */
  54. /* Unchanged on exit. */
  55. /* TRANS - CHARACTER*1. */
  56. /* On entry, TRANS specifies the operation to be performed as */
  57. /* follows: */
  58. /* TRANS = 'N' or 'n' C := alpha*A*conjg( B' ) + */
  59. /* conjg( alpha )*B*conjg( A' ) + */
  60. /* beta*C. */
  61. /* TRANS = 'C' or 'c' C := alpha*conjg( A' )*B + */
  62. /* conjg( alpha )*conjg( B' )*A + */
  63. /* beta*C. */
  64. /* Unchanged on exit. */
  65. /* N - INTEGER. */
  66. /* On entry, N specifies the order of the matrix C. N must be */
  67. /* at least zero. */
  68. /* Unchanged on exit. */
  69. /* K - INTEGER. */
  70. /* On entry with TRANS = 'N' or 'n', K specifies the number */
  71. /* of columns of the matrices A and B, and on entry with */
  72. /* TRANS = 'C' or 'c', K specifies the number of rows of the */
  73. /* matrices A and B. K must be at least zero. */
  74. /* Unchanged on exit. */
  75. /* ALPHA - COMPLEX*16 . */
  76. /* On entry, ALPHA specifies the scalar alpha. */
  77. /* Unchanged on exit. */
  78. /* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is */
  79. /* k when TRANS = 'N' or 'n', and is n otherwise. */
  80. /* Before entry with TRANS = 'N' or 'n', the leading n by k */
  81. /* part of the array A must contain the matrix A, otherwise */
  82. /* the leading k by n part of the array A must contain the */
  83. /* matrix A. */
  84. /* Unchanged on exit. */
  85. /* LDA - INTEGER. */
  86. /* On entry, LDA specifies the first dimension of A as declared */
  87. /* in the calling (sub) program. When TRANS = 'N' or 'n' */
  88. /* then LDA must be at least max( 1, n ), otherwise LDA must */
  89. /* be at least max( 1, k ). */
  90. /* Unchanged on exit. */
  91. /* B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is */
  92. /* k when TRANS = 'N' or 'n', and is n otherwise. */
  93. /* Before entry with TRANS = 'N' or 'n', the leading n by k */
  94. /* part of the array B must contain the matrix B, otherwise */
  95. /* the leading k by n part of the array B must contain the */
  96. /* matrix B. */
  97. /* Unchanged on exit. */
  98. /* LDB - INTEGER. */
  99. /* On entry, LDB specifies the first dimension of B as declared */
  100. /* in the calling (sub) program. When TRANS = 'N' or 'n' */
  101. /* then LDB must be at least max( 1, n ), otherwise LDB must */
  102. /* be at least max( 1, k ). */
  103. /* Unchanged on exit. */
  104. /* BETA - DOUBLE PRECISION . */
  105. /* On entry, BETA specifies the scalar beta. */
  106. /* Unchanged on exit. */
  107. /* C - COMPLEX*16 array of DIMENSION ( LDC, n ). */
  108. /* Before entry with UPLO = 'U' or 'u', the leading n by n */
  109. /* upper triangular part of the array C must contain the upper */
  110. /* triangular part of the hermitian matrix and the strictly */
  111. /* lower triangular part of C is not referenced. On exit, the */
  112. /* upper triangular part of the array C is overwritten by the */
  113. /* upper triangular part of the updated matrix. */
  114. /* Before entry with UPLO = 'L' or 'l', the leading n by n */
  115. /* lower triangular part of the array C must contain the lower */
  116. /* triangular part of the hermitian matrix and the strictly */
  117. /* upper triangular part of C is not referenced. On exit, the */
  118. /* lower triangular part of the array C is overwritten by the */
  119. /* lower triangular part of the updated matrix. */
  120. /* Note that the imaginary parts of the diagonal elements need */
  121. /* not be set, they are assumed to be zero, and on exit they */
  122. /* are set to zero. */
  123. /* LDC - INTEGER. */
  124. /* On entry, LDC specifies the first dimension of C as declared */
  125. /* in the calling (sub) program. LDC must be at least */
  126. /* max( 1, n ). */
  127. /* Unchanged on exit. */
  128. /* Level 3 Blas routine. */
  129. /* -- Written on 8-February-1989. */
  130. /* Jack Dongarra, Argonne National Laboratory. */
  131. /* Iain Duff, AERE Harwell. */
  132. /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
  133. /* Sven Hammarling, Numerical Algorithms Group Ltd. */
  134. /* -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. */
  135. /* Ed Anderson, Cray Research Inc. */
  136. /* .. External Functions .. */
  137. /* .. */
  138. /* .. External Subroutines .. */
  139. /* .. */
  140. /* .. Intrinsic Functions .. */
  141. /* .. */
  142. /* .. Local Scalars .. */
  143. /* .. */
  144. /* .. Parameters .. */
  145. /* .. */
  146. /* Test the input parameters. */
  147. /* Parameter adjustments */
  148. a_dim1 = *lda;
  149. a_offset = 1 + a_dim1;
  150. a -= a_offset;
  151. b_dim1 = *ldb;
  152. b_offset = 1 + b_dim1;
  153. b -= b_offset;
  154. c_dim1 = *ldc;
  155. c_offset = 1 + c_dim1;
  156. c__ -= c_offset;
  157. /* Function Body */
  158. if (lsame_(trans, "N")) {
  159. nrowa = *n;
  160. } else {
  161. nrowa = *k;
  162. }
  163. upper = lsame_(uplo, "U");
  164. info = 0;
  165. if (! upper && ! lsame_(uplo, "L")) {
  166. info = 1;
  167. } else if (! lsame_(trans, "N") && ! lsame_(trans,
  168. "C")) {
  169. info = 2;
  170. } else if (*n < 0) {
  171. info = 3;
  172. } else if (*k < 0) {
  173. info = 4;
  174. } else if (*lda < max(1,nrowa)) {
  175. info = 7;
  176. } else if (*ldb < max(1,nrowa)) {
  177. info = 9;
  178. } else if (*ldc < max(1,*n)) {
  179. info = 12;
  180. }
  181. if (info != 0) {
  182. xerbla_("ZHER2K", &info);
  183. return 0;
  184. }
  185. /* Quick return if possible. */
  186. if (*n == 0 || (alpha->r == 0. && alpha->i == 0. || *k == 0) && *beta ==
  187. 1.) {
  188. return 0;
  189. }
  190. /* And when alpha.eq.zero. */
  191. if (alpha->r == 0. && alpha->i == 0.) {
  192. if (upper) {
  193. if (*beta == 0.) {
  194. i__1 = *n;
  195. for (j = 1; j <= i__1; ++j) {
  196. i__2 = j;
  197. for (i__ = 1; i__ <= i__2; ++i__) {
  198. i__3 = i__ + j * c_dim1;
  199. c__[i__3].r = 0., c__[i__3].i = 0.;
  200. /* L10: */
  201. }
  202. /* L20: */
  203. }
  204. } else {
  205. i__1 = *n;
  206. for (j = 1; j <= i__1; ++j) {
  207. i__2 = j - 1;
  208. for (i__ = 1; i__ <= i__2; ++i__) {
  209. i__3 = i__ + j * c_dim1;
  210. i__4 = i__ + j * c_dim1;
  211. z__1.r = *beta * c__[i__4].r, z__1.i = *beta * c__[
  212. i__4].i;
  213. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  214. /* L30: */
  215. }
  216. i__2 = j + j * c_dim1;
  217. i__3 = j + j * c_dim1;
  218. d__1 = *beta * c__[i__3].r;
  219. c__[i__2].r = d__1, c__[i__2].i = 0.;
  220. /* L40: */
  221. }
  222. }
  223. } else {
  224. if (*beta == 0.) {
  225. i__1 = *n;
  226. for (j = 1; j <= i__1; ++j) {
  227. i__2 = *n;
  228. for (i__ = j; i__ <= i__2; ++i__) {
  229. i__3 = i__ + j * c_dim1;
  230. c__[i__3].r = 0., c__[i__3].i = 0.;
  231. /* L50: */
  232. }
  233. /* L60: */
  234. }
  235. } else {
  236. i__1 = *n;
  237. for (j = 1; j <= i__1; ++j) {
  238. i__2 = j + j * c_dim1;
  239. i__3 = j + j * c_dim1;
  240. d__1 = *beta * c__[i__3].r;
  241. c__[i__2].r = d__1, c__[i__2].i = 0.;
  242. i__2 = *n;
  243. for (i__ = j + 1; i__ <= i__2; ++i__) {
  244. i__3 = i__ + j * c_dim1;
  245. i__4 = i__ + j * c_dim1;
  246. z__1.r = *beta * c__[i__4].r, z__1.i = *beta * c__[
  247. i__4].i;
  248. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  249. /* L70: */
  250. }
  251. /* L80: */
  252. }
  253. }
  254. }
  255. return 0;
  256. }
  257. /* Start the operations. */
  258. if (lsame_(trans, "N")) {
  259. /* Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + */
  260. /* C. */
  261. if (upper) {
  262. i__1 = *n;
  263. for (j = 1; j <= i__1; ++j) {
  264. if (*beta == 0.) {
  265. i__2 = j;
  266. for (i__ = 1; i__ <= i__2; ++i__) {
  267. i__3 = i__ + j * c_dim1;
  268. c__[i__3].r = 0., c__[i__3].i = 0.;
  269. /* L90: */
  270. }
  271. } else if (*beta != 1.) {
  272. i__2 = j - 1;
  273. for (i__ = 1; i__ <= i__2; ++i__) {
  274. i__3 = i__ + j * c_dim1;
  275. i__4 = i__ + j * c_dim1;
  276. z__1.r = *beta * c__[i__4].r, z__1.i = *beta * c__[
  277. i__4].i;
  278. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  279. /* L100: */
  280. }
  281. i__2 = j + j * c_dim1;
  282. i__3 = j + j * c_dim1;
  283. d__1 = *beta * c__[i__3].r;
  284. c__[i__2].r = d__1, c__[i__2].i = 0.;
  285. } else {
  286. i__2 = j + j * c_dim1;
  287. i__3 = j + j * c_dim1;
  288. d__1 = c__[i__3].r;
  289. c__[i__2].r = d__1, c__[i__2].i = 0.;
  290. }
  291. i__2 = *k;
  292. for (l = 1; l <= i__2; ++l) {
  293. i__3 = j + l * a_dim1;
  294. i__4 = j + l * b_dim1;
  295. if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r !=
  296. 0. || b[i__4].i != 0.)) {
  297. d_cnjg(&z__2, &b[j + l * b_dim1]);
  298. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i,
  299. z__1.i = alpha->r * z__2.i + alpha->i *
  300. z__2.r;
  301. temp1.r = z__1.r, temp1.i = z__1.i;
  302. i__3 = j + l * a_dim1;
  303. z__2.r = alpha->r * a[i__3].r - alpha->i * a[i__3].i,
  304. z__2.i = alpha->r * a[i__3].i + alpha->i * a[
  305. i__3].r;
  306. d_cnjg(&z__1, &z__2);
  307. temp2.r = z__1.r, temp2.i = z__1.i;
  308. i__3 = j - 1;
  309. for (i__ = 1; i__ <= i__3; ++i__) {
  310. i__4 = i__ + j * c_dim1;
  311. i__5 = i__ + j * c_dim1;
  312. i__6 = i__ + l * a_dim1;
  313. z__3.r = a[i__6].r * temp1.r - a[i__6].i *
  314. temp1.i, z__3.i = a[i__6].r * temp1.i + a[
  315. i__6].i * temp1.r;
  316. z__2.r = c__[i__5].r + z__3.r, z__2.i = c__[i__5]
  317. .i + z__3.i;
  318. i__7 = i__ + l * b_dim1;
  319. z__4.r = b[i__7].r * temp2.r - b[i__7].i *
  320. temp2.i, z__4.i = b[i__7].r * temp2.i + b[
  321. i__7].i * temp2.r;
  322. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i +
  323. z__4.i;
  324. c__[i__4].r = z__1.r, c__[i__4].i = z__1.i;
  325. /* L110: */
  326. }
  327. i__3 = j + j * c_dim1;
  328. i__4 = j + j * c_dim1;
  329. i__5 = j + l * a_dim1;
  330. z__2.r = a[i__5].r * temp1.r - a[i__5].i * temp1.i,
  331. z__2.i = a[i__5].r * temp1.i + a[i__5].i *
  332. temp1.r;
  333. i__6 = j + l * b_dim1;
  334. z__3.r = b[i__6].r * temp2.r - b[i__6].i * temp2.i,
  335. z__3.i = b[i__6].r * temp2.i + b[i__6].i *
  336. temp2.r;
  337. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  338. d__1 = c__[i__4].r + z__1.r;
  339. c__[i__3].r = d__1, c__[i__3].i = 0.;
  340. }
  341. /* L120: */
  342. }
  343. /* L130: */
  344. }
  345. } else {
  346. i__1 = *n;
  347. for (j = 1; j <= i__1; ++j) {
  348. if (*beta == 0.) {
  349. i__2 = *n;
  350. for (i__ = j; i__ <= i__2; ++i__) {
  351. i__3 = i__ + j * c_dim1;
  352. c__[i__3].r = 0., c__[i__3].i = 0.;
  353. /* L140: */
  354. }
  355. } else if (*beta != 1.) {
  356. i__2 = *n;
  357. for (i__ = j + 1; i__ <= i__2; ++i__) {
  358. i__3 = i__ + j * c_dim1;
  359. i__4 = i__ + j * c_dim1;
  360. z__1.r = *beta * c__[i__4].r, z__1.i = *beta * c__[
  361. i__4].i;
  362. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  363. /* L150: */
  364. }
  365. i__2 = j + j * c_dim1;
  366. i__3 = j + j * c_dim1;
  367. d__1 = *beta * c__[i__3].r;
  368. c__[i__2].r = d__1, c__[i__2].i = 0.;
  369. } else {
  370. i__2 = j + j * c_dim1;
  371. i__3 = j + j * c_dim1;
  372. d__1 = c__[i__3].r;
  373. c__[i__2].r = d__1, c__[i__2].i = 0.;
  374. }
  375. i__2 = *k;
  376. for (l = 1; l <= i__2; ++l) {
  377. i__3 = j + l * a_dim1;
  378. i__4 = j + l * b_dim1;
  379. if (a[i__3].r != 0. || a[i__3].i != 0. || (b[i__4].r !=
  380. 0. || b[i__4].i != 0.)) {
  381. d_cnjg(&z__2, &b[j + l * b_dim1]);
  382. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i,
  383. z__1.i = alpha->r * z__2.i + alpha->i *
  384. z__2.r;
  385. temp1.r = z__1.r, temp1.i = z__1.i;
  386. i__3 = j + l * a_dim1;
  387. z__2.r = alpha->r * a[i__3].r - alpha->i * a[i__3].i,
  388. z__2.i = alpha->r * a[i__3].i + alpha->i * a[
  389. i__3].r;
  390. d_cnjg(&z__1, &z__2);
  391. temp2.r = z__1.r, temp2.i = z__1.i;
  392. i__3 = *n;
  393. for (i__ = j + 1; i__ <= i__3; ++i__) {
  394. i__4 = i__ + j * c_dim1;
  395. i__5 = i__ + j * c_dim1;
  396. i__6 = i__ + l * a_dim1;
  397. z__3.r = a[i__6].r * temp1.r - a[i__6].i *
  398. temp1.i, z__3.i = a[i__6].r * temp1.i + a[
  399. i__6].i * temp1.r;
  400. z__2.r = c__[i__5].r + z__3.r, z__2.i = c__[i__5]
  401. .i + z__3.i;
  402. i__7 = i__ + l * b_dim1;
  403. z__4.r = b[i__7].r * temp2.r - b[i__7].i *
  404. temp2.i, z__4.i = b[i__7].r * temp2.i + b[
  405. i__7].i * temp2.r;
  406. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i +
  407. z__4.i;
  408. c__[i__4].r = z__1.r, c__[i__4].i = z__1.i;
  409. /* L160: */
  410. }
  411. i__3 = j + j * c_dim1;
  412. i__4 = j + j * c_dim1;
  413. i__5 = j + l * a_dim1;
  414. z__2.r = a[i__5].r * temp1.r - a[i__5].i * temp1.i,
  415. z__2.i = a[i__5].r * temp1.i + a[i__5].i *
  416. temp1.r;
  417. i__6 = j + l * b_dim1;
  418. z__3.r = b[i__6].r * temp2.r - b[i__6].i * temp2.i,
  419. z__3.i = b[i__6].r * temp2.i + b[i__6].i *
  420. temp2.r;
  421. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  422. d__1 = c__[i__4].r + z__1.r;
  423. c__[i__3].r = d__1, c__[i__3].i = 0.;
  424. }
  425. /* L170: */
  426. }
  427. /* L180: */
  428. }
  429. }
  430. } else {
  431. /* Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + */
  432. /* C. */
  433. if (upper) {
  434. i__1 = *n;
  435. for (j = 1; j <= i__1; ++j) {
  436. i__2 = j;
  437. for (i__ = 1; i__ <= i__2; ++i__) {
  438. temp1.r = 0., temp1.i = 0.;
  439. temp2.r = 0., temp2.i = 0.;
  440. i__3 = *k;
  441. for (l = 1; l <= i__3; ++l) {
  442. d_cnjg(&z__3, &a[l + i__ * a_dim1]);
  443. i__4 = l + j * b_dim1;
  444. z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
  445. z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
  446. .r;
  447. z__1.r = temp1.r + z__2.r, z__1.i = temp1.i + z__2.i;
  448. temp1.r = z__1.r, temp1.i = z__1.i;
  449. d_cnjg(&z__3, &b[l + i__ * b_dim1]);
  450. i__4 = l + j * a_dim1;
  451. z__2.r = z__3.r * a[i__4].r - z__3.i * a[i__4].i,
  452. z__2.i = z__3.r * a[i__4].i + z__3.i * a[i__4]
  453. .r;
  454. z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
  455. temp2.r = z__1.r, temp2.i = z__1.i;
  456. /* L190: */
  457. }
  458. if (i__ == j) {
  459. if (*beta == 0.) {
  460. i__3 = j + j * c_dim1;
  461. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  462. z__2.i = alpha->r * temp1.i + alpha->i *
  463. temp1.r;
  464. d_cnjg(&z__4, alpha);
  465. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  466. z__3.i = z__4.r * temp2.i + z__4.i *
  467. temp2.r;
  468. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  469. z__3.i;
  470. d__1 = z__1.r;
  471. c__[i__3].r = d__1, c__[i__3].i = 0.;
  472. } else {
  473. i__3 = j + j * c_dim1;
  474. i__4 = j + j * c_dim1;
  475. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  476. z__2.i = alpha->r * temp1.i + alpha->i *
  477. temp1.r;
  478. d_cnjg(&z__4, alpha);
  479. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  480. z__3.i = z__4.r * temp2.i + z__4.i *
  481. temp2.r;
  482. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  483. z__3.i;
  484. d__1 = *beta * c__[i__4].r + z__1.r;
  485. c__[i__3].r = d__1, c__[i__3].i = 0.;
  486. }
  487. } else {
  488. if (*beta == 0.) {
  489. i__3 = i__ + j * c_dim1;
  490. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  491. z__2.i = alpha->r * temp1.i + alpha->i *
  492. temp1.r;
  493. d_cnjg(&z__4, alpha);
  494. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  495. z__3.i = z__4.r * temp2.i + z__4.i *
  496. temp2.r;
  497. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  498. z__3.i;
  499. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  500. } else {
  501. i__3 = i__ + j * c_dim1;
  502. i__4 = i__ + j * c_dim1;
  503. z__3.r = *beta * c__[i__4].r, z__3.i = *beta *
  504. c__[i__4].i;
  505. z__4.r = alpha->r * temp1.r - alpha->i * temp1.i,
  506. z__4.i = alpha->r * temp1.i + alpha->i *
  507. temp1.r;
  508. z__2.r = z__3.r + z__4.r, z__2.i = z__3.i +
  509. z__4.i;
  510. d_cnjg(&z__6, alpha);
  511. z__5.r = z__6.r * temp2.r - z__6.i * temp2.i,
  512. z__5.i = z__6.r * temp2.i + z__6.i *
  513. temp2.r;
  514. z__1.r = z__2.r + z__5.r, z__1.i = z__2.i +
  515. z__5.i;
  516. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  517. }
  518. }
  519. /* L200: */
  520. }
  521. /* L210: */
  522. }
  523. } else {
  524. i__1 = *n;
  525. for (j = 1; j <= i__1; ++j) {
  526. i__2 = *n;
  527. for (i__ = j; i__ <= i__2; ++i__) {
  528. temp1.r = 0., temp1.i = 0.;
  529. temp2.r = 0., temp2.i = 0.;
  530. i__3 = *k;
  531. for (l = 1; l <= i__3; ++l) {
  532. d_cnjg(&z__3, &a[l + i__ * a_dim1]);
  533. i__4 = l + j * b_dim1;
  534. z__2.r = z__3.r * b[i__4].r - z__3.i * b[i__4].i,
  535. z__2.i = z__3.r * b[i__4].i + z__3.i * b[i__4]
  536. .r;
  537. z__1.r = temp1.r + z__2.r, z__1.i = temp1.i + z__2.i;
  538. temp1.r = z__1.r, temp1.i = z__1.i;
  539. d_cnjg(&z__3, &b[l + i__ * b_dim1]);
  540. i__4 = l + j * a_dim1;
  541. z__2.r = z__3.r * a[i__4].r - z__3.i * a[i__4].i,
  542. z__2.i = z__3.r * a[i__4].i + z__3.i * a[i__4]
  543. .r;
  544. z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
  545. temp2.r = z__1.r, temp2.i = z__1.i;
  546. /* L220: */
  547. }
  548. if (i__ == j) {
  549. if (*beta == 0.) {
  550. i__3 = j + j * c_dim1;
  551. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  552. z__2.i = alpha->r * temp1.i + alpha->i *
  553. temp1.r;
  554. d_cnjg(&z__4, alpha);
  555. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  556. z__3.i = z__4.r * temp2.i + z__4.i *
  557. temp2.r;
  558. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  559. z__3.i;
  560. d__1 = z__1.r;
  561. c__[i__3].r = d__1, c__[i__3].i = 0.;
  562. } else {
  563. i__3 = j + j * c_dim1;
  564. i__4 = j + j * c_dim1;
  565. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  566. z__2.i = alpha->r * temp1.i + alpha->i *
  567. temp1.r;
  568. d_cnjg(&z__4, alpha);
  569. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  570. z__3.i = z__4.r * temp2.i + z__4.i *
  571. temp2.r;
  572. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  573. z__3.i;
  574. d__1 = *beta * c__[i__4].r + z__1.r;
  575. c__[i__3].r = d__1, c__[i__3].i = 0.;
  576. }
  577. } else {
  578. if (*beta == 0.) {
  579. i__3 = i__ + j * c_dim1;
  580. z__2.r = alpha->r * temp1.r - alpha->i * temp1.i,
  581. z__2.i = alpha->r * temp1.i + alpha->i *
  582. temp1.r;
  583. d_cnjg(&z__4, alpha);
  584. z__3.r = z__4.r * temp2.r - z__4.i * temp2.i,
  585. z__3.i = z__4.r * temp2.i + z__4.i *
  586. temp2.r;
  587. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i +
  588. z__3.i;
  589. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  590. } else {
  591. i__3 = i__ + j * c_dim1;
  592. i__4 = i__ + j * c_dim1;
  593. z__3.r = *beta * c__[i__4].r, z__3.i = *beta *
  594. c__[i__4].i;
  595. z__4.r = alpha->r * temp1.r - alpha->i * temp1.i,
  596. z__4.i = alpha->r * temp1.i + alpha->i *
  597. temp1.r;
  598. z__2.r = z__3.r + z__4.r, z__2.i = z__3.i +
  599. z__4.i;
  600. d_cnjg(&z__6, alpha);
  601. z__5.r = z__6.r * temp2.r - z__6.i * temp2.i,
  602. z__5.i = z__6.r * temp2.i + z__6.i *
  603. temp2.r;
  604. z__1.r = z__2.r + z__5.r, z__1.i = z__2.i +
  605. z__5.i;
  606. c__[i__3].r = z__1.r, c__[i__3].i = z__1.i;
  607. }
  608. }
  609. /* L230: */
  610. }
  611. /* L240: */
  612. }
  613. }
  614. }
  615. return 0;
  616. /* End of ZHER2K. */
  617. } /* zher2k_ */