dtrmv.c 8.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345
  1. /* dtrmv.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 dtrmv_(char *uplo, char *trans, char *diag, integer *n,
  14. doublereal *a, integer *lda, doublereal *x, integer *incx)
  15. {
  16. /* System generated locals */
  17. integer a_dim1, a_offset, i__1, i__2;
  18. /* Local variables */
  19. integer i__, j, ix, jx, kx, info;
  20. doublereal temp;
  21. extern logical lsame_(char *, char *);
  22. extern /* Subroutine */ int xerbla_(char *, integer *);
  23. logical nounit;
  24. /* .. Scalar Arguments .. */
  25. /* .. */
  26. /* .. Array Arguments .. */
  27. /* .. */
  28. /* Purpose */
  29. /* ======= */
  30. /* DTRMV performs one of the matrix-vector operations */
  31. /* x := A*x, or x := A'*x, */
  32. /* where x is an n element vector and A is an n by n unit, or non-unit, */
  33. /* upper or lower triangular matrix. */
  34. /* Arguments */
  35. /* ========== */
  36. /* UPLO - CHARACTER*1. */
  37. /* On entry, UPLO specifies whether the matrix is an upper or */
  38. /* lower triangular matrix as follows: */
  39. /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
  40. /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
  41. /* Unchanged on exit. */
  42. /* TRANS - CHARACTER*1. */
  43. /* On entry, TRANS specifies the operation to be performed as */
  44. /* follows: */
  45. /* TRANS = 'N' or 'n' x := A*x. */
  46. /* TRANS = 'T' or 't' x := A'*x. */
  47. /* TRANS = 'C' or 'c' x := A'*x. */
  48. /* Unchanged on exit. */
  49. /* DIAG - CHARACTER*1. */
  50. /* On entry, DIAG specifies whether or not A is unit */
  51. /* triangular as follows: */
  52. /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
  53. /* DIAG = 'N' or 'n' A is not assumed to be unit */
  54. /* triangular. */
  55. /* Unchanged on exit. */
  56. /* N - INTEGER. */
  57. /* On entry, N specifies the order of the matrix A. */
  58. /* N must be at least zero. */
  59. /* Unchanged on exit. */
  60. /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
  61. /* Before entry with UPLO = 'U' or 'u', the leading n by n */
  62. /* upper triangular part of the array A must contain the upper */
  63. /* triangular matrix and the strictly lower triangular part of */
  64. /* A is not referenced. */
  65. /* Before entry with UPLO = 'L' or 'l', the leading n by n */
  66. /* lower triangular part of the array A must contain the lower */
  67. /* triangular matrix and the strictly upper triangular part of */
  68. /* A is not referenced. */
  69. /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
  70. /* A are not referenced either, but are assumed to be unity. */
  71. /* Unchanged on exit. */
  72. /* LDA - INTEGER. */
  73. /* On entry, LDA specifies the first dimension of A as declared */
  74. /* in the calling (sub) program. LDA must be at least */
  75. /* max( 1, n ). */
  76. /* Unchanged on exit. */
  77. /* X - DOUBLE PRECISION array of dimension at least */
  78. /* ( 1 + ( n - 1 )*abs( INCX ) ). */
  79. /* Before entry, the incremented array X must contain the n */
  80. /* element vector x. On exit, X is overwritten with the */
  81. /* tranformed vector x. */
  82. /* INCX - INTEGER. */
  83. /* On entry, INCX specifies the increment for the elements of */
  84. /* X. INCX must not be zero. */
  85. /* Unchanged on exit. */
  86. /* Level 2 Blas routine. */
  87. /* -- Written on 22-October-1986. */
  88. /* Jack Dongarra, Argonne National Lab. */
  89. /* Jeremy Du Croz, Nag Central Office. */
  90. /* Sven Hammarling, Nag Central Office. */
  91. /* Richard Hanson, Sandia National Labs. */
  92. /* .. Parameters .. */
  93. /* .. */
  94. /* .. Local Scalars .. */
  95. /* .. */
  96. /* .. External Functions .. */
  97. /* .. */
  98. /* .. External Subroutines .. */
  99. /* .. */
  100. /* .. Intrinsic Functions .. */
  101. /* .. */
  102. /* Test the input parameters. */
  103. /* Parameter adjustments */
  104. a_dim1 = *lda;
  105. a_offset = 1 + a_dim1;
  106. a -= a_offset;
  107. --x;
  108. /* Function Body */
  109. info = 0;
  110. if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
  111. info = 1;
  112. } else if (! lsame_(trans, "N") && ! lsame_(trans,
  113. "T") && ! lsame_(trans, "C")) {
  114. info = 2;
  115. } else if (! lsame_(diag, "U") && ! lsame_(diag,
  116. "N")) {
  117. info = 3;
  118. } else if (*n < 0) {
  119. info = 4;
  120. } else if (*lda < max(1,*n)) {
  121. info = 6;
  122. } else if (*incx == 0) {
  123. info = 8;
  124. }
  125. if (info != 0) {
  126. xerbla_("DTRMV ", &info);
  127. return 0;
  128. }
  129. /* Quick return if possible. */
  130. if (*n == 0) {
  131. return 0;
  132. }
  133. nounit = lsame_(diag, "N");
  134. /* Set up the start point in X if the increment is not unity. This */
  135. /* will be ( N - 1 )*INCX too small for descending loops. */
  136. if (*incx <= 0) {
  137. kx = 1 - (*n - 1) * *incx;
  138. } else if (*incx != 1) {
  139. kx = 1;
  140. }
  141. /* Start the operations. In this version the elements of A are */
  142. /* accessed sequentially with one pass through A. */
  143. if (lsame_(trans, "N")) {
  144. /* Form x := A*x. */
  145. if (lsame_(uplo, "U")) {
  146. if (*incx == 1) {
  147. i__1 = *n;
  148. for (j = 1; j <= i__1; ++j) {
  149. if (x[j] != 0.) {
  150. temp = x[j];
  151. i__2 = j - 1;
  152. for (i__ = 1; i__ <= i__2; ++i__) {
  153. x[i__] += temp * a[i__ + j * a_dim1];
  154. /* L10: */
  155. }
  156. if (nounit) {
  157. x[j] *= a[j + j * a_dim1];
  158. }
  159. }
  160. /* L20: */
  161. }
  162. } else {
  163. jx = kx;
  164. i__1 = *n;
  165. for (j = 1; j <= i__1; ++j) {
  166. if (x[jx] != 0.) {
  167. temp = x[jx];
  168. ix = kx;
  169. i__2 = j - 1;
  170. for (i__ = 1; i__ <= i__2; ++i__) {
  171. x[ix] += temp * a[i__ + j * a_dim1];
  172. ix += *incx;
  173. /* L30: */
  174. }
  175. if (nounit) {
  176. x[jx] *= a[j + j * a_dim1];
  177. }
  178. }
  179. jx += *incx;
  180. /* L40: */
  181. }
  182. }
  183. } else {
  184. if (*incx == 1) {
  185. for (j = *n; j >= 1; --j) {
  186. if (x[j] != 0.) {
  187. temp = x[j];
  188. i__1 = j + 1;
  189. for (i__ = *n; i__ >= i__1; --i__) {
  190. x[i__] += temp * a[i__ + j * a_dim1];
  191. /* L50: */
  192. }
  193. if (nounit) {
  194. x[j] *= a[j + j * a_dim1];
  195. }
  196. }
  197. /* L60: */
  198. }
  199. } else {
  200. kx += (*n - 1) * *incx;
  201. jx = kx;
  202. for (j = *n; j >= 1; --j) {
  203. if (x[jx] != 0.) {
  204. temp = x[jx];
  205. ix = kx;
  206. i__1 = j + 1;
  207. for (i__ = *n; i__ >= i__1; --i__) {
  208. x[ix] += temp * a[i__ + j * a_dim1];
  209. ix -= *incx;
  210. /* L70: */
  211. }
  212. if (nounit) {
  213. x[jx] *= a[j + j * a_dim1];
  214. }
  215. }
  216. jx -= *incx;
  217. /* L80: */
  218. }
  219. }
  220. }
  221. } else {
  222. /* Form x := A'*x. */
  223. if (lsame_(uplo, "U")) {
  224. if (*incx == 1) {
  225. for (j = *n; j >= 1; --j) {
  226. temp = x[j];
  227. if (nounit) {
  228. temp *= a[j + j * a_dim1];
  229. }
  230. for (i__ = j - 1; i__ >= 1; --i__) {
  231. temp += a[i__ + j * a_dim1] * x[i__];
  232. /* L90: */
  233. }
  234. x[j] = temp;
  235. /* L100: */
  236. }
  237. } else {
  238. jx = kx + (*n - 1) * *incx;
  239. for (j = *n; j >= 1; --j) {
  240. temp = x[jx];
  241. ix = jx;
  242. if (nounit) {
  243. temp *= a[j + j * a_dim1];
  244. }
  245. for (i__ = j - 1; i__ >= 1; --i__) {
  246. ix -= *incx;
  247. temp += a[i__ + j * a_dim1] * x[ix];
  248. /* L110: */
  249. }
  250. x[jx] = temp;
  251. jx -= *incx;
  252. /* L120: */
  253. }
  254. }
  255. } else {
  256. if (*incx == 1) {
  257. i__1 = *n;
  258. for (j = 1; j <= i__1; ++j) {
  259. temp = x[j];
  260. if (nounit) {
  261. temp *= a[j + j * a_dim1];
  262. }
  263. i__2 = *n;
  264. for (i__ = j + 1; i__ <= i__2; ++i__) {
  265. temp += a[i__ + j * a_dim1] * x[i__];
  266. /* L130: */
  267. }
  268. x[j] = temp;
  269. /* L140: */
  270. }
  271. } else {
  272. jx = kx;
  273. i__1 = *n;
  274. for (j = 1; j <= i__1; ++j) {
  275. temp = x[jx];
  276. ix = jx;
  277. if (nounit) {
  278. temp *= a[j + j * a_dim1];
  279. }
  280. i__2 = *n;
  281. for (i__ = j + 1; i__ <= i__2; ++i__) {
  282. ix += *incx;
  283. temp += a[i__ + j * a_dim1] * x[ix];
  284. /* L150: */
  285. }
  286. x[jx] = temp;
  287. jx += *incx;
  288. /* L160: */
  289. }
  290. }
  291. }
  292. }
  293. return 0;
  294. /* End of DTRMV . */
  295. } /* dtrmv_ */