zher2.c 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447
  1. /* zher2.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 zher2_(char *uplo, integer *n, doublecomplex *alpha,
  14. doublecomplex *x, integer *incx, doublecomplex *y, integer *incy,
  15. doublecomplex *a, integer *lda)
  16. {
  17. /* System generated locals */
  18. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  19. doublereal d__1;
  20. doublecomplex z__1, z__2, z__3, z__4;
  21. /* Builtin functions */
  22. void d_cnjg(doublecomplex *, doublecomplex *);
  23. /* Local variables */
  24. integer i__, j, ix, iy, jx, jy, kx, ky, info;
  25. doublecomplex temp1, temp2;
  26. extern logical lsame_(char *, char *);
  27. extern /* Subroutine */ int xerbla_(char *, integer *);
  28. /* .. Scalar Arguments .. */
  29. /* .. */
  30. /* .. Array Arguments .. */
  31. /* .. */
  32. /* Purpose */
  33. /* ======= */
  34. /* ZHER2 performs the hermitian rank 2 operation */
  35. /* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
  36. /* where alpha is a scalar, x and y are n element vectors and A is an n */
  37. /* by n hermitian matrix. */
  38. /* Arguments */
  39. /* ========== */
  40. /* UPLO - CHARACTER*1. */
  41. /* On entry, UPLO specifies whether the upper or lower */
  42. /* triangular part of the array A is to be referenced as */
  43. /* follows: */
  44. /* UPLO = 'U' or 'u' Only the upper triangular part of A */
  45. /* is to be referenced. */
  46. /* UPLO = 'L' or 'l' Only the lower triangular part of A */
  47. /* is to be referenced. */
  48. /* Unchanged on exit. */
  49. /* N - INTEGER. */
  50. /* On entry, N specifies the order of the matrix A. */
  51. /* N must be at least zero. */
  52. /* Unchanged on exit. */
  53. /* ALPHA - COMPLEX*16 . */
  54. /* On entry, ALPHA specifies the scalar alpha. */
  55. /* Unchanged on exit. */
  56. /* X - COMPLEX*16 array of dimension at least */
  57. /* ( 1 + ( n - 1 )*abs( INCX ) ). */
  58. /* Before entry, the incremented array X must contain the n */
  59. /* element vector x. */
  60. /* Unchanged on exit. */
  61. /* INCX - INTEGER. */
  62. /* On entry, INCX specifies the increment for the elements of */
  63. /* X. INCX must not be zero. */
  64. /* Unchanged on exit. */
  65. /* Y - COMPLEX*16 array of dimension at least */
  66. /* ( 1 + ( n - 1 )*abs( INCY ) ). */
  67. /* Before entry, the incremented array Y must contain the n */
  68. /* element vector y. */
  69. /* Unchanged on exit. */
  70. /* INCY - INTEGER. */
  71. /* On entry, INCY specifies the increment for the elements of */
  72. /* Y. INCY must not be zero. */
  73. /* Unchanged on exit. */
  74. /* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
  75. /* Before entry with UPLO = 'U' or 'u', the leading n by n */
  76. /* upper triangular part of the array A must contain the upper */
  77. /* triangular part of the hermitian matrix and the strictly */
  78. /* lower triangular part of A is not referenced. On exit, the */
  79. /* upper triangular part of the array A is overwritten by the */
  80. /* upper triangular part of the updated matrix. */
  81. /* Before entry with UPLO = 'L' or 'l', the leading n by n */
  82. /* lower triangular part of the array A must contain the lower */
  83. /* triangular part of the hermitian matrix and the strictly */
  84. /* upper triangular part of A is not referenced. On exit, the */
  85. /* lower triangular part of the array A is overwritten by the */
  86. /* lower triangular part of the updated matrix. */
  87. /* Note that the imaginary parts of the diagonal elements need */
  88. /* not be set, they are assumed to be zero, and on exit they */
  89. /* are set to zero. */
  90. /* LDA - INTEGER. */
  91. /* On entry, LDA specifies the first dimension of A as declared */
  92. /* in the calling (sub) program. LDA must be at least */
  93. /* max( 1, n ). */
  94. /* Unchanged on exit. */
  95. /* Level 2 Blas routine. */
  96. /* -- Written on 22-October-1986. */
  97. /* Jack Dongarra, Argonne National Lab. */
  98. /* Jeremy Du Croz, Nag Central Office. */
  99. /* Sven Hammarling, Nag Central Office. */
  100. /* Richard Hanson, Sandia National Labs. */
  101. /* .. Parameters .. */
  102. /* .. */
  103. /* .. Local Scalars .. */
  104. /* .. */
  105. /* .. External Functions .. */
  106. /* .. */
  107. /* .. External Subroutines .. */
  108. /* .. */
  109. /* .. Intrinsic Functions .. */
  110. /* .. */
  111. /* Test the input parameters. */
  112. /* Parameter adjustments */
  113. --x;
  114. --y;
  115. a_dim1 = *lda;
  116. a_offset = 1 + a_dim1;
  117. a -= a_offset;
  118. /* Function Body */
  119. info = 0;
  120. if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
  121. info = 1;
  122. } else if (*n < 0) {
  123. info = 2;
  124. } else if (*incx == 0) {
  125. info = 5;
  126. } else if (*incy == 0) {
  127. info = 7;
  128. } else if (*lda < max(1,*n)) {
  129. info = 9;
  130. }
  131. if (info != 0) {
  132. xerbla_("ZHER2 ", &info);
  133. return 0;
  134. }
  135. /* Quick return if possible. */
  136. if (*n == 0 || alpha->r == 0. && alpha->i == 0.) {
  137. return 0;
  138. }
  139. /* Set up the start points in X and Y if the increments are not both */
  140. /* unity. */
  141. if (*incx != 1 || *incy != 1) {
  142. if (*incx > 0) {
  143. kx = 1;
  144. } else {
  145. kx = 1 - (*n - 1) * *incx;
  146. }
  147. if (*incy > 0) {
  148. ky = 1;
  149. } else {
  150. ky = 1 - (*n - 1) * *incy;
  151. }
  152. jx = kx;
  153. jy = ky;
  154. }
  155. /* Start the operations. In this version the elements of A are */
  156. /* accessed sequentially with one pass through the triangular part */
  157. /* of A. */
  158. if (lsame_(uplo, "U")) {
  159. /* Form A when A is stored in the upper triangle. */
  160. if (*incx == 1 && *incy == 1) {
  161. i__1 = *n;
  162. for (j = 1; j <= i__1; ++j) {
  163. i__2 = j;
  164. i__3 = j;
  165. if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. ||
  166. y[i__3].i != 0.)) {
  167. d_cnjg(&z__2, &y[j]);
  168. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
  169. alpha->r * z__2.i + alpha->i * z__2.r;
  170. temp1.r = z__1.r, temp1.i = z__1.i;
  171. i__2 = j;
  172. z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
  173. z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
  174. .r;
  175. d_cnjg(&z__1, &z__2);
  176. temp2.r = z__1.r, temp2.i = z__1.i;
  177. i__2 = j - 1;
  178. for (i__ = 1; i__ <= i__2; ++i__) {
  179. i__3 = i__ + j * a_dim1;
  180. i__4 = i__ + j * a_dim1;
  181. i__5 = i__;
  182. z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i,
  183. z__3.i = x[i__5].r * temp1.i + x[i__5].i *
  184. temp1.r;
  185. z__2.r = a[i__4].r + z__3.r, z__2.i = a[i__4].i +
  186. z__3.i;
  187. i__6 = i__;
  188. z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i,
  189. z__4.i = y[i__6].r * temp2.i + y[i__6].i *
  190. temp2.r;
  191. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
  192. a[i__3].r = z__1.r, a[i__3].i = z__1.i;
  193. /* L10: */
  194. }
  195. i__2 = j + j * a_dim1;
  196. i__3 = j + j * a_dim1;
  197. i__4 = j;
  198. z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i,
  199. z__2.i = x[i__4].r * temp1.i + x[i__4].i *
  200. temp1.r;
  201. i__5 = j;
  202. z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i,
  203. z__3.i = y[i__5].r * temp2.i + y[i__5].i *
  204. temp2.r;
  205. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  206. d__1 = a[i__3].r + z__1.r;
  207. a[i__2].r = d__1, a[i__2].i = 0.;
  208. } else {
  209. i__2 = j + j * a_dim1;
  210. i__3 = j + j * a_dim1;
  211. d__1 = a[i__3].r;
  212. a[i__2].r = d__1, a[i__2].i = 0.;
  213. }
  214. /* L20: */
  215. }
  216. } else {
  217. i__1 = *n;
  218. for (j = 1; j <= i__1; ++j) {
  219. i__2 = jx;
  220. i__3 = jy;
  221. if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. ||
  222. y[i__3].i != 0.)) {
  223. d_cnjg(&z__2, &y[jy]);
  224. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
  225. alpha->r * z__2.i + alpha->i * z__2.r;
  226. temp1.r = z__1.r, temp1.i = z__1.i;
  227. i__2 = jx;
  228. z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
  229. z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
  230. .r;
  231. d_cnjg(&z__1, &z__2);
  232. temp2.r = z__1.r, temp2.i = z__1.i;
  233. ix = kx;
  234. iy = ky;
  235. i__2 = j - 1;
  236. for (i__ = 1; i__ <= i__2; ++i__) {
  237. i__3 = i__ + j * a_dim1;
  238. i__4 = i__ + j * a_dim1;
  239. i__5 = ix;
  240. z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i,
  241. z__3.i = x[i__5].r * temp1.i + x[i__5].i *
  242. temp1.r;
  243. z__2.r = a[i__4].r + z__3.r, z__2.i = a[i__4].i +
  244. z__3.i;
  245. i__6 = iy;
  246. z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i,
  247. z__4.i = y[i__6].r * temp2.i + y[i__6].i *
  248. temp2.r;
  249. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
  250. a[i__3].r = z__1.r, a[i__3].i = z__1.i;
  251. ix += *incx;
  252. iy += *incy;
  253. /* L30: */
  254. }
  255. i__2 = j + j * a_dim1;
  256. i__3 = j + j * a_dim1;
  257. i__4 = jx;
  258. z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i,
  259. z__2.i = x[i__4].r * temp1.i + x[i__4].i *
  260. temp1.r;
  261. i__5 = jy;
  262. z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i,
  263. z__3.i = y[i__5].r * temp2.i + y[i__5].i *
  264. temp2.r;
  265. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  266. d__1 = a[i__3].r + z__1.r;
  267. a[i__2].r = d__1, a[i__2].i = 0.;
  268. } else {
  269. i__2 = j + j * a_dim1;
  270. i__3 = j + j * a_dim1;
  271. d__1 = a[i__3].r;
  272. a[i__2].r = d__1, a[i__2].i = 0.;
  273. }
  274. jx += *incx;
  275. jy += *incy;
  276. /* L40: */
  277. }
  278. }
  279. } else {
  280. /* Form A when A is stored in the lower triangle. */
  281. if (*incx == 1 && *incy == 1) {
  282. i__1 = *n;
  283. for (j = 1; j <= i__1; ++j) {
  284. i__2 = j;
  285. i__3 = j;
  286. if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. ||
  287. y[i__3].i != 0.)) {
  288. d_cnjg(&z__2, &y[j]);
  289. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
  290. alpha->r * z__2.i + alpha->i * z__2.r;
  291. temp1.r = z__1.r, temp1.i = z__1.i;
  292. i__2 = j;
  293. z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
  294. z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
  295. .r;
  296. d_cnjg(&z__1, &z__2);
  297. temp2.r = z__1.r, temp2.i = z__1.i;
  298. i__2 = j + j * a_dim1;
  299. i__3 = j + j * a_dim1;
  300. i__4 = j;
  301. z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i,
  302. z__2.i = x[i__4].r * temp1.i + x[i__4].i *
  303. temp1.r;
  304. i__5 = j;
  305. z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i,
  306. z__3.i = y[i__5].r * temp2.i + y[i__5].i *
  307. temp2.r;
  308. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  309. d__1 = a[i__3].r + z__1.r;
  310. a[i__2].r = d__1, a[i__2].i = 0.;
  311. i__2 = *n;
  312. for (i__ = j + 1; i__ <= i__2; ++i__) {
  313. i__3 = i__ + j * a_dim1;
  314. i__4 = i__ + j * a_dim1;
  315. i__5 = i__;
  316. z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i,
  317. z__3.i = x[i__5].r * temp1.i + x[i__5].i *
  318. temp1.r;
  319. z__2.r = a[i__4].r + z__3.r, z__2.i = a[i__4].i +
  320. z__3.i;
  321. i__6 = i__;
  322. z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i,
  323. z__4.i = y[i__6].r * temp2.i + y[i__6].i *
  324. temp2.r;
  325. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
  326. a[i__3].r = z__1.r, a[i__3].i = z__1.i;
  327. /* L50: */
  328. }
  329. } else {
  330. i__2 = j + j * a_dim1;
  331. i__3 = j + j * a_dim1;
  332. d__1 = a[i__3].r;
  333. a[i__2].r = d__1, a[i__2].i = 0.;
  334. }
  335. /* L60: */
  336. }
  337. } else {
  338. i__1 = *n;
  339. for (j = 1; j <= i__1; ++j) {
  340. i__2 = jx;
  341. i__3 = jy;
  342. if (x[i__2].r != 0. || x[i__2].i != 0. || (y[i__3].r != 0. ||
  343. y[i__3].i != 0.)) {
  344. d_cnjg(&z__2, &y[jy]);
  345. z__1.r = alpha->r * z__2.r - alpha->i * z__2.i, z__1.i =
  346. alpha->r * z__2.i + alpha->i * z__2.r;
  347. temp1.r = z__1.r, temp1.i = z__1.i;
  348. i__2 = jx;
  349. z__2.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i,
  350. z__2.i = alpha->r * x[i__2].i + alpha->i * x[i__2]
  351. .r;
  352. d_cnjg(&z__1, &z__2);
  353. temp2.r = z__1.r, temp2.i = z__1.i;
  354. i__2 = j + j * a_dim1;
  355. i__3 = j + j * a_dim1;
  356. i__4 = jx;
  357. z__2.r = x[i__4].r * temp1.r - x[i__4].i * temp1.i,
  358. z__2.i = x[i__4].r * temp1.i + x[i__4].i *
  359. temp1.r;
  360. i__5 = jy;
  361. z__3.r = y[i__5].r * temp2.r - y[i__5].i * temp2.i,
  362. z__3.i = y[i__5].r * temp2.i + y[i__5].i *
  363. temp2.r;
  364. z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
  365. d__1 = a[i__3].r + z__1.r;
  366. a[i__2].r = d__1, a[i__2].i = 0.;
  367. ix = jx;
  368. iy = jy;
  369. i__2 = *n;
  370. for (i__ = j + 1; i__ <= i__2; ++i__) {
  371. ix += *incx;
  372. iy += *incy;
  373. i__3 = i__ + j * a_dim1;
  374. i__4 = i__ + j * a_dim1;
  375. i__5 = ix;
  376. z__3.r = x[i__5].r * temp1.r - x[i__5].i * temp1.i,
  377. z__3.i = x[i__5].r * temp1.i + x[i__5].i *
  378. temp1.r;
  379. z__2.r = a[i__4].r + z__3.r, z__2.i = a[i__4].i +
  380. z__3.i;
  381. i__6 = iy;
  382. z__4.r = y[i__6].r * temp2.r - y[i__6].i * temp2.i,
  383. z__4.i = y[i__6].r * temp2.i + y[i__6].i *
  384. temp2.r;
  385. z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
  386. a[i__3].r = z__1.r, a[i__3].i = z__1.i;
  387. /* L70: */
  388. }
  389. } else {
  390. i__2 = j + j * a_dim1;
  391. i__3 = j + j * a_dim1;
  392. d__1 = a[i__3].r;
  393. a[i__2].r = d__1, a[i__2].i = 0.;
  394. }
  395. jx += *incx;
  396. jy += *incy;
  397. /* L80: */
  398. }
  399. }
  400. }
  401. return 0;
  402. /* End of ZHER2 . */
  403. } /* zher2_ */