ctrmv.c 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554
  1. /* ctrmv.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 ctrmv_(char *uplo, char *trans, char *diag, integer *n,
  14. complex *a, integer *lda, complex *x, integer *incx)
  15. {
  16. /* System generated locals */
  17. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
  18. complex q__1, q__2, q__3;
  19. /* Builtin functions */
  20. void r_cnjg(complex *, complex *);
  21. /* Local variables */
  22. integer i__, j, ix, jx, kx, info;
  23. complex temp;
  24. extern logical lsame_(char *, char *);
  25. extern /* Subroutine */ int xerbla_(char *, integer *);
  26. logical noconj, nounit;
  27. /* .. Scalar Arguments .. */
  28. /* .. */
  29. /* .. Array Arguments .. */
  30. /* .. */
  31. /* Purpose */
  32. /* ======= */
  33. /* CTRMV performs one of the matrix-vector operations */
  34. /* x := A*x, or x := A'*x, or x := conjg( A' )*x, */
  35. /* where x is an n element vector and A is an n by n unit, or non-unit, */
  36. /* upper or lower triangular matrix. */
  37. /* Arguments */
  38. /* ========== */
  39. /* UPLO - CHARACTER*1. */
  40. /* On entry, UPLO specifies whether the matrix is an upper or */
  41. /* lower triangular matrix as follows: */
  42. /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
  43. /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
  44. /* Unchanged on exit. */
  45. /* TRANS - CHARACTER*1. */
  46. /* On entry, TRANS specifies the operation to be performed as */
  47. /* follows: */
  48. /* TRANS = 'N' or 'n' x := A*x. */
  49. /* TRANS = 'T' or 't' x := A'*x. */
  50. /* TRANS = 'C' or 'c' x := conjg( A' )*x. */
  51. /* Unchanged on exit. */
  52. /* DIAG - CHARACTER*1. */
  53. /* On entry, DIAG specifies whether or not A is unit */
  54. /* triangular as follows: */
  55. /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
  56. /* DIAG = 'N' or 'n' A is not assumed to be unit */
  57. /* triangular. */
  58. /* Unchanged on exit. */
  59. /* N - INTEGER. */
  60. /* On entry, N specifies the order of the matrix A. */
  61. /* N must be at least zero. */
  62. /* Unchanged on exit. */
  63. /* A - COMPLEX array of DIMENSION ( LDA, n ). */
  64. /* Before entry with UPLO = 'U' or 'u', the leading n by n */
  65. /* upper triangular part of the array A must contain the upper */
  66. /* triangular matrix and the strictly lower triangular part of */
  67. /* A is not referenced. */
  68. /* Before entry with UPLO = 'L' or 'l', the leading n by n */
  69. /* lower triangular part of the array A must contain the lower */
  70. /* triangular matrix and the strictly upper triangular part of */
  71. /* A is not referenced. */
  72. /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
  73. /* A are not referenced either, but are assumed to be unity. */
  74. /* Unchanged on exit. */
  75. /* LDA - INTEGER. */
  76. /* On entry, LDA specifies the first dimension of A as declared */
  77. /* in the calling (sub) program. LDA must be at least */
  78. /* max( 1, n ). */
  79. /* Unchanged on exit. */
  80. /* X - COMPLEX array of dimension at least */
  81. /* ( 1 + ( n - 1 )*abs( INCX ) ). */
  82. /* Before entry, the incremented array X must contain the n */
  83. /* element vector x. On exit, X is overwritten with the */
  84. /* tranformed vector x. */
  85. /* INCX - INTEGER. */
  86. /* On entry, INCX specifies the increment for the elements of */
  87. /* X. INCX must not be zero. */
  88. /* Unchanged on exit. */
  89. /* Level 2 Blas routine. */
  90. /* -- Written on 22-October-1986. */
  91. /* Jack Dongarra, Argonne National Lab. */
  92. /* Jeremy Du Croz, Nag Central Office. */
  93. /* Sven Hammarling, Nag Central Office. */
  94. /* Richard Hanson, Sandia National Labs. */
  95. /* .. Parameters .. */
  96. /* .. */
  97. /* .. Local Scalars .. */
  98. /* .. */
  99. /* .. External Functions .. */
  100. /* .. */
  101. /* .. External Subroutines .. */
  102. /* .. */
  103. /* .. Intrinsic Functions .. */
  104. /* .. */
  105. /* Test the input parameters. */
  106. /* Parameter adjustments */
  107. a_dim1 = *lda;
  108. a_offset = 1 + a_dim1;
  109. a -= a_offset;
  110. --x;
  111. /* Function Body */
  112. info = 0;
  113. if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
  114. info = 1;
  115. } else if (! lsame_(trans, "N") && ! lsame_(trans,
  116. "T") && ! lsame_(trans, "C")) {
  117. info = 2;
  118. } else if (! lsame_(diag, "U") && ! lsame_(diag,
  119. "N")) {
  120. info = 3;
  121. } else if (*n < 0) {
  122. info = 4;
  123. } else if (*lda < max(1,*n)) {
  124. info = 6;
  125. } else if (*incx == 0) {
  126. info = 8;
  127. }
  128. if (info != 0) {
  129. xerbla_("CTRMV ", &info);
  130. return 0;
  131. }
  132. /* Quick return if possible. */
  133. if (*n == 0) {
  134. return 0;
  135. }
  136. noconj = lsame_(trans, "T");
  137. nounit = lsame_(diag, "N");
  138. /* Set up the start point in X if the increment is not unity. This */
  139. /* will be ( N - 1 )*INCX too small for descending loops. */
  140. if (*incx <= 0) {
  141. kx = 1 - (*n - 1) * *incx;
  142. } else if (*incx != 1) {
  143. kx = 1;
  144. }
  145. /* Start the operations. In this version the elements of A are */
  146. /* accessed sequentially with one pass through A. */
  147. if (lsame_(trans, "N")) {
  148. /* Form x := A*x. */
  149. if (lsame_(uplo, "U")) {
  150. if (*incx == 1) {
  151. i__1 = *n;
  152. for (j = 1; j <= i__1; ++j) {
  153. i__2 = j;
  154. if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
  155. i__2 = j;
  156. temp.r = x[i__2].r, temp.i = x[i__2].i;
  157. i__2 = j - 1;
  158. for (i__ = 1; i__ <= i__2; ++i__) {
  159. i__3 = i__;
  160. i__4 = i__;
  161. i__5 = i__ + j * a_dim1;
  162. q__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
  163. q__2.i = temp.r * a[i__5].i + temp.i * a[
  164. i__5].r;
  165. q__1.r = x[i__4].r + q__2.r, q__1.i = x[i__4].i +
  166. q__2.i;
  167. x[i__3].r = q__1.r, x[i__3].i = q__1.i;
  168. /* L10: */
  169. }
  170. if (nounit) {
  171. i__2 = j;
  172. i__3 = j;
  173. i__4 = j + j * a_dim1;
  174. q__1.r = x[i__3].r * a[i__4].r - x[i__3].i * a[
  175. i__4].i, q__1.i = x[i__3].r * a[i__4].i +
  176. x[i__3].i * a[i__4].r;
  177. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  178. }
  179. }
  180. /* L20: */
  181. }
  182. } else {
  183. jx = kx;
  184. i__1 = *n;
  185. for (j = 1; j <= i__1; ++j) {
  186. i__2 = jx;
  187. if (x[i__2].r != 0.f || x[i__2].i != 0.f) {
  188. i__2 = jx;
  189. temp.r = x[i__2].r, temp.i = x[i__2].i;
  190. ix = kx;
  191. i__2 = j - 1;
  192. for (i__ = 1; i__ <= i__2; ++i__) {
  193. i__3 = ix;
  194. i__4 = ix;
  195. i__5 = i__ + j * a_dim1;
  196. q__2.r = temp.r * a[i__5].r - temp.i * a[i__5].i,
  197. q__2.i = temp.r * a[i__5].i + temp.i * a[
  198. i__5].r;
  199. q__1.r = x[i__4].r + q__2.r, q__1.i = x[i__4].i +
  200. q__2.i;
  201. x[i__3].r = q__1.r, x[i__3].i = q__1.i;
  202. ix += *incx;
  203. /* L30: */
  204. }
  205. if (nounit) {
  206. i__2 = jx;
  207. i__3 = jx;
  208. i__4 = j + j * a_dim1;
  209. q__1.r = x[i__3].r * a[i__4].r - x[i__3].i * a[
  210. i__4].i, q__1.i = x[i__3].r * a[i__4].i +
  211. x[i__3].i * a[i__4].r;
  212. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  213. }
  214. }
  215. jx += *incx;
  216. /* L40: */
  217. }
  218. }
  219. } else {
  220. if (*incx == 1) {
  221. for (j = *n; j >= 1; --j) {
  222. i__1 = j;
  223. if (x[i__1].r != 0.f || x[i__1].i != 0.f) {
  224. i__1 = j;
  225. temp.r = x[i__1].r, temp.i = x[i__1].i;
  226. i__1 = j + 1;
  227. for (i__ = *n; i__ >= i__1; --i__) {
  228. i__2 = i__;
  229. i__3 = i__;
  230. i__4 = i__ + j * a_dim1;
  231. q__2.r = temp.r * a[i__4].r - temp.i * a[i__4].i,
  232. q__2.i = temp.r * a[i__4].i + temp.i * a[
  233. i__4].r;
  234. q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i +
  235. q__2.i;
  236. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  237. /* L50: */
  238. }
  239. if (nounit) {
  240. i__1 = j;
  241. i__2 = j;
  242. i__3 = j + j * a_dim1;
  243. q__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
  244. i__3].i, q__1.i = x[i__2].r * a[i__3].i +
  245. x[i__2].i * a[i__3].r;
  246. x[i__1].r = q__1.r, x[i__1].i = q__1.i;
  247. }
  248. }
  249. /* L60: */
  250. }
  251. } else {
  252. kx += (*n - 1) * *incx;
  253. jx = kx;
  254. for (j = *n; j >= 1; --j) {
  255. i__1 = jx;
  256. if (x[i__1].r != 0.f || x[i__1].i != 0.f) {
  257. i__1 = jx;
  258. temp.r = x[i__1].r, temp.i = x[i__1].i;
  259. ix = kx;
  260. i__1 = j + 1;
  261. for (i__ = *n; i__ >= i__1; --i__) {
  262. i__2 = ix;
  263. i__3 = ix;
  264. i__4 = i__ + j * a_dim1;
  265. q__2.r = temp.r * a[i__4].r - temp.i * a[i__4].i,
  266. q__2.i = temp.r * a[i__4].i + temp.i * a[
  267. i__4].r;
  268. q__1.r = x[i__3].r + q__2.r, q__1.i = x[i__3].i +
  269. q__2.i;
  270. x[i__2].r = q__1.r, x[i__2].i = q__1.i;
  271. ix -= *incx;
  272. /* L70: */
  273. }
  274. if (nounit) {
  275. i__1 = jx;
  276. i__2 = jx;
  277. i__3 = j + j * a_dim1;
  278. q__1.r = x[i__2].r * a[i__3].r - x[i__2].i * a[
  279. i__3].i, q__1.i = x[i__2].r * a[i__3].i +
  280. x[i__2].i * a[i__3].r;
  281. x[i__1].r = q__1.r, x[i__1].i = q__1.i;
  282. }
  283. }
  284. jx -= *incx;
  285. /* L80: */
  286. }
  287. }
  288. }
  289. } else {
  290. /* Form x := A'*x or x := conjg( A' )*x. */
  291. if (lsame_(uplo, "U")) {
  292. if (*incx == 1) {
  293. for (j = *n; j >= 1; --j) {
  294. i__1 = j;
  295. temp.r = x[i__1].r, temp.i = x[i__1].i;
  296. if (noconj) {
  297. if (nounit) {
  298. i__1 = j + j * a_dim1;
  299. q__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
  300. q__1.i = temp.r * a[i__1].i + temp.i * a[
  301. i__1].r;
  302. temp.r = q__1.r, temp.i = q__1.i;
  303. }
  304. for (i__ = j - 1; i__ >= 1; --i__) {
  305. i__1 = i__ + j * a_dim1;
  306. i__2 = i__;
  307. q__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
  308. i__2].i, q__2.i = a[i__1].r * x[i__2].i +
  309. a[i__1].i * x[i__2].r;
  310. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  311. q__2.i;
  312. temp.r = q__1.r, temp.i = q__1.i;
  313. /* L90: */
  314. }
  315. } else {
  316. if (nounit) {
  317. r_cnjg(&q__2, &a[j + j * a_dim1]);
  318. q__1.r = temp.r * q__2.r - temp.i * q__2.i,
  319. q__1.i = temp.r * q__2.i + temp.i *
  320. q__2.r;
  321. temp.r = q__1.r, temp.i = q__1.i;
  322. }
  323. for (i__ = j - 1; i__ >= 1; --i__) {
  324. r_cnjg(&q__3, &a[i__ + j * a_dim1]);
  325. i__1 = i__;
  326. q__2.r = q__3.r * x[i__1].r - q__3.i * x[i__1].i,
  327. q__2.i = q__3.r * x[i__1].i + q__3.i * x[
  328. i__1].r;
  329. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  330. q__2.i;
  331. temp.r = q__1.r, temp.i = q__1.i;
  332. /* L100: */
  333. }
  334. }
  335. i__1 = j;
  336. x[i__1].r = temp.r, x[i__1].i = temp.i;
  337. /* L110: */
  338. }
  339. } else {
  340. jx = kx + (*n - 1) * *incx;
  341. for (j = *n; j >= 1; --j) {
  342. i__1 = jx;
  343. temp.r = x[i__1].r, temp.i = x[i__1].i;
  344. ix = jx;
  345. if (noconj) {
  346. if (nounit) {
  347. i__1 = j + j * a_dim1;
  348. q__1.r = temp.r * a[i__1].r - temp.i * a[i__1].i,
  349. q__1.i = temp.r * a[i__1].i + temp.i * a[
  350. i__1].r;
  351. temp.r = q__1.r, temp.i = q__1.i;
  352. }
  353. for (i__ = j - 1; i__ >= 1; --i__) {
  354. ix -= *incx;
  355. i__1 = i__ + j * a_dim1;
  356. i__2 = ix;
  357. q__2.r = a[i__1].r * x[i__2].r - a[i__1].i * x[
  358. i__2].i, q__2.i = a[i__1].r * x[i__2].i +
  359. a[i__1].i * x[i__2].r;
  360. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  361. q__2.i;
  362. temp.r = q__1.r, temp.i = q__1.i;
  363. /* L120: */
  364. }
  365. } else {
  366. if (nounit) {
  367. r_cnjg(&q__2, &a[j + j * a_dim1]);
  368. q__1.r = temp.r * q__2.r - temp.i * q__2.i,
  369. q__1.i = temp.r * q__2.i + temp.i *
  370. q__2.r;
  371. temp.r = q__1.r, temp.i = q__1.i;
  372. }
  373. for (i__ = j - 1; i__ >= 1; --i__) {
  374. ix -= *incx;
  375. r_cnjg(&q__3, &a[i__ + j * a_dim1]);
  376. i__1 = ix;
  377. q__2.r = q__3.r * x[i__1].r - q__3.i * x[i__1].i,
  378. q__2.i = q__3.r * x[i__1].i + q__3.i * x[
  379. i__1].r;
  380. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  381. q__2.i;
  382. temp.r = q__1.r, temp.i = q__1.i;
  383. /* L130: */
  384. }
  385. }
  386. i__1 = jx;
  387. x[i__1].r = temp.r, x[i__1].i = temp.i;
  388. jx -= *incx;
  389. /* L140: */
  390. }
  391. }
  392. } else {
  393. if (*incx == 1) {
  394. i__1 = *n;
  395. for (j = 1; j <= i__1; ++j) {
  396. i__2 = j;
  397. temp.r = x[i__2].r, temp.i = x[i__2].i;
  398. if (noconj) {
  399. if (nounit) {
  400. i__2 = j + j * a_dim1;
  401. q__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
  402. q__1.i = temp.r * a[i__2].i + temp.i * a[
  403. i__2].r;
  404. temp.r = q__1.r, temp.i = q__1.i;
  405. }
  406. i__2 = *n;
  407. for (i__ = j + 1; i__ <= i__2; ++i__) {
  408. i__3 = i__ + j * a_dim1;
  409. i__4 = i__;
  410. q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
  411. i__4].i, q__2.i = a[i__3].r * x[i__4].i +
  412. a[i__3].i * x[i__4].r;
  413. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  414. q__2.i;
  415. temp.r = q__1.r, temp.i = q__1.i;
  416. /* L150: */
  417. }
  418. } else {
  419. if (nounit) {
  420. r_cnjg(&q__2, &a[j + j * a_dim1]);
  421. q__1.r = temp.r * q__2.r - temp.i * q__2.i,
  422. q__1.i = temp.r * q__2.i + temp.i *
  423. q__2.r;
  424. temp.r = q__1.r, temp.i = q__1.i;
  425. }
  426. i__2 = *n;
  427. for (i__ = j + 1; i__ <= i__2; ++i__) {
  428. r_cnjg(&q__3, &a[i__ + j * a_dim1]);
  429. i__3 = i__;
  430. q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i,
  431. q__2.i = q__3.r * x[i__3].i + q__3.i * x[
  432. i__3].r;
  433. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  434. q__2.i;
  435. temp.r = q__1.r, temp.i = q__1.i;
  436. /* L160: */
  437. }
  438. }
  439. i__2 = j;
  440. x[i__2].r = temp.r, x[i__2].i = temp.i;
  441. /* L170: */
  442. }
  443. } else {
  444. jx = kx;
  445. i__1 = *n;
  446. for (j = 1; j <= i__1; ++j) {
  447. i__2 = jx;
  448. temp.r = x[i__2].r, temp.i = x[i__2].i;
  449. ix = jx;
  450. if (noconj) {
  451. if (nounit) {
  452. i__2 = j + j * a_dim1;
  453. q__1.r = temp.r * a[i__2].r - temp.i * a[i__2].i,
  454. q__1.i = temp.r * a[i__2].i + temp.i * a[
  455. i__2].r;
  456. temp.r = q__1.r, temp.i = q__1.i;
  457. }
  458. i__2 = *n;
  459. for (i__ = j + 1; i__ <= i__2; ++i__) {
  460. ix += *incx;
  461. i__3 = i__ + j * a_dim1;
  462. i__4 = ix;
  463. q__2.r = a[i__3].r * x[i__4].r - a[i__3].i * x[
  464. i__4].i, q__2.i = a[i__3].r * x[i__4].i +
  465. a[i__3].i * x[i__4].r;
  466. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  467. q__2.i;
  468. temp.r = q__1.r, temp.i = q__1.i;
  469. /* L180: */
  470. }
  471. } else {
  472. if (nounit) {
  473. r_cnjg(&q__2, &a[j + j * a_dim1]);
  474. q__1.r = temp.r * q__2.r - temp.i * q__2.i,
  475. q__1.i = temp.r * q__2.i + temp.i *
  476. q__2.r;
  477. temp.r = q__1.r, temp.i = q__1.i;
  478. }
  479. i__2 = *n;
  480. for (i__ = j + 1; i__ <= i__2; ++i__) {
  481. ix += *incx;
  482. r_cnjg(&q__3, &a[i__ + j * a_dim1]);
  483. i__3 = ix;
  484. q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i,
  485. q__2.i = q__3.r * x[i__3].i + q__3.i * x[
  486. i__3].r;
  487. q__1.r = temp.r + q__2.r, q__1.i = temp.i +
  488. q__2.i;
  489. temp.r = q__1.r, temp.i = q__1.i;
  490. /* L190: */
  491. }
  492. }
  493. i__2 = jx;
  494. x[i__2].r = temp.r, x[i__2].i = temp.i;
  495. jx += *incx;
  496. /* L200: */
  497. }
  498. }
  499. }
  500. }
  501. return 0;
  502. /* End of CTRMV . */
  503. } /* ctrmv_ */