sgemv.c 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312
  1. /* sgemv.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 sgemv_(char *trans, integer *m, integer *n, real *alpha,
  14. real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
  15. integer *incy)
  16. {
  17. /* System generated locals */
  18. integer a_dim1, a_offset, i__1, i__2;
  19. /* Local variables */
  20. integer i__, j, ix, iy, jx, jy, kx, ky, info;
  21. real temp;
  22. integer lenx, leny;
  23. extern logical lsame_(char *, char *);
  24. extern /* Subroutine */ int xerbla_(char *, integer *);
  25. /* .. Scalar Arguments .. */
  26. /* .. */
  27. /* .. Array Arguments .. */
  28. /* .. */
  29. /* Purpose */
  30. /* ======= */
  31. /* SGEMV performs one of the matrix-vector operations */
  32. /* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, */
  33. /* where alpha and beta are scalars, x and y are vectors and A is an */
  34. /* m by n matrix. */
  35. /* Arguments */
  36. /* ========== */
  37. /* TRANS - CHARACTER*1. */
  38. /* On entry, TRANS specifies the operation to be performed as */
  39. /* follows: */
  40. /* TRANS = 'N' or 'n' y := alpha*A*x + beta*y. */
  41. /* TRANS = 'T' or 't' y := alpha*A'*x + beta*y. */
  42. /* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. */
  43. /* Unchanged on exit. */
  44. /* M - INTEGER. */
  45. /* On entry, M specifies the number of rows of the matrix A. */
  46. /* M must be at least zero. */
  47. /* Unchanged on exit. */
  48. /* N - INTEGER. */
  49. /* On entry, N specifies the number of columns of the matrix A. */
  50. /* N must be at least zero. */
  51. /* Unchanged on exit. */
  52. /* ALPHA - REAL . */
  53. /* On entry, ALPHA specifies the scalar alpha. */
  54. /* Unchanged on exit. */
  55. /* A - REAL array of DIMENSION ( LDA, n ). */
  56. /* Before entry, the leading m by n part of the array A must */
  57. /* contain the matrix of coefficients. */
  58. /* Unchanged on exit. */
  59. /* LDA - INTEGER. */
  60. /* On entry, LDA specifies the first dimension of A as declared */
  61. /* in the calling (sub) program. LDA must be at least */
  62. /* max( 1, m ). */
  63. /* Unchanged on exit. */
  64. /* X - REAL array of DIMENSION at least */
  65. /* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
  66. /* and at least */
  67. /* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
  68. /* Before entry, the incremented array X must contain the */
  69. /* vector x. */
  70. /* Unchanged on exit. */
  71. /* INCX - INTEGER. */
  72. /* On entry, INCX specifies the increment for the elements of */
  73. /* X. INCX must not be zero. */
  74. /* Unchanged on exit. */
  75. /* BETA - REAL . */
  76. /* On entry, BETA specifies the scalar beta. When BETA is */
  77. /* supplied as zero then Y need not be set on input. */
  78. /* Unchanged on exit. */
  79. /* Y - REAL array of DIMENSION at least */
  80. /* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
  81. /* and at least */
  82. /* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
  83. /* Before entry with BETA non-zero, the incremented array Y */
  84. /* must contain the vector y. On exit, Y is overwritten by the */
  85. /* updated vector y. */
  86. /* INCY - INTEGER. */
  87. /* On entry, INCY specifies the increment for the elements of */
  88. /* Y. INCY must not be zero. */
  89. /* Unchanged on exit. */
  90. /* Level 2 Blas routine. */
  91. /* -- Written on 22-October-1986. */
  92. /* Jack Dongarra, Argonne National Lab. */
  93. /* Jeremy Du Croz, Nag Central Office. */
  94. /* Sven Hammarling, Nag Central Office. */
  95. /* Richard Hanson, Sandia National Labs. */
  96. /* .. Parameters .. */
  97. /* .. */
  98. /* .. Local Scalars .. */
  99. /* .. */
  100. /* .. External Functions .. */
  101. /* .. */
  102. /* .. External Subroutines .. */
  103. /* .. */
  104. /* .. Intrinsic Functions .. */
  105. /* .. */
  106. /* Test the input parameters. */
  107. /* Parameter adjustments */
  108. a_dim1 = *lda;
  109. a_offset = 1 + a_dim1;
  110. a -= a_offset;
  111. --x;
  112. --y;
  113. /* Function Body */
  114. info = 0;
  115. if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
  116. ) {
  117. info = 1;
  118. } else if (*m < 0) {
  119. info = 2;
  120. } else if (*n < 0) {
  121. info = 3;
  122. } else if (*lda < max(1,*m)) {
  123. info = 6;
  124. } else if (*incx == 0) {
  125. info = 8;
  126. } else if (*incy == 0) {
  127. info = 11;
  128. }
  129. if (info != 0) {
  130. xerbla_("SGEMV ", &info);
  131. return 0;
  132. }
  133. /* Quick return if possible. */
  134. if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
  135. return 0;
  136. }
  137. /* Set LENX and LENY, the lengths of the vectors x and y, and set */
  138. /* up the start points in X and Y. */
  139. if (lsame_(trans, "N")) {
  140. lenx = *n;
  141. leny = *m;
  142. } else {
  143. lenx = *m;
  144. leny = *n;
  145. }
  146. if (*incx > 0) {
  147. kx = 1;
  148. } else {
  149. kx = 1 - (lenx - 1) * *incx;
  150. }
  151. if (*incy > 0) {
  152. ky = 1;
  153. } else {
  154. ky = 1 - (leny - 1) * *incy;
  155. }
  156. /* Start the operations. In this version the elements of A are */
  157. /* accessed sequentially with one pass through A. */
  158. /* First form y := beta*y. */
  159. if (*beta != 1.f) {
  160. if (*incy == 1) {
  161. if (*beta == 0.f) {
  162. i__1 = leny;
  163. for (i__ = 1; i__ <= i__1; ++i__) {
  164. y[i__] = 0.f;
  165. /* L10: */
  166. }
  167. } else {
  168. i__1 = leny;
  169. for (i__ = 1; i__ <= i__1; ++i__) {
  170. y[i__] = *beta * y[i__];
  171. /* L20: */
  172. }
  173. }
  174. } else {
  175. iy = ky;
  176. if (*beta == 0.f) {
  177. i__1 = leny;
  178. for (i__ = 1; i__ <= i__1; ++i__) {
  179. y[iy] = 0.f;
  180. iy += *incy;
  181. /* L30: */
  182. }
  183. } else {
  184. i__1 = leny;
  185. for (i__ = 1; i__ <= i__1; ++i__) {
  186. y[iy] = *beta * y[iy];
  187. iy += *incy;
  188. /* L40: */
  189. }
  190. }
  191. }
  192. }
  193. if (*alpha == 0.f) {
  194. return 0;
  195. }
  196. if (lsame_(trans, "N")) {
  197. /* Form y := alpha*A*x + y. */
  198. jx = kx;
  199. if (*incy == 1) {
  200. i__1 = *n;
  201. for (j = 1; j <= i__1; ++j) {
  202. if (x[jx] != 0.f) {
  203. temp = *alpha * x[jx];
  204. i__2 = *m;
  205. for (i__ = 1; i__ <= i__2; ++i__) {
  206. y[i__] += temp * a[i__ + j * a_dim1];
  207. /* L50: */
  208. }
  209. }
  210. jx += *incx;
  211. /* L60: */
  212. }
  213. } else {
  214. i__1 = *n;
  215. for (j = 1; j <= i__1; ++j) {
  216. if (x[jx] != 0.f) {
  217. temp = *alpha * x[jx];
  218. iy = ky;
  219. i__2 = *m;
  220. for (i__ = 1; i__ <= i__2; ++i__) {
  221. y[iy] += temp * a[i__ + j * a_dim1];
  222. iy += *incy;
  223. /* L70: */
  224. }
  225. }
  226. jx += *incx;
  227. /* L80: */
  228. }
  229. }
  230. } else {
  231. /* Form y := alpha*A'*x + y. */
  232. jy = ky;
  233. if (*incx == 1) {
  234. i__1 = *n;
  235. for (j = 1; j <= i__1; ++j) {
  236. temp = 0.f;
  237. i__2 = *m;
  238. for (i__ = 1; i__ <= i__2; ++i__) {
  239. temp += a[i__ + j * a_dim1] * x[i__];
  240. /* L90: */
  241. }
  242. y[jy] += *alpha * temp;
  243. jy += *incy;
  244. /* L100: */
  245. }
  246. } else {
  247. i__1 = *n;
  248. for (j = 1; j <= i__1; ++j) {
  249. temp = 0.f;
  250. ix = kx;
  251. i__2 = *m;
  252. for (i__ = 1; i__ <= i__2; ++i__) {
  253. temp += a[i__ + j * a_dim1] * x[ix];
  254. ix += *incx;
  255. /* L110: */
  256. }
  257. y[jy] += *alpha * temp;
  258. jy += *incy;
  259. /* L120: */
  260. }
  261. }
  262. }
  263. return 0;
  264. /* End of SGEMV . */
  265. } /* sgemv_ */