sgbmv.c 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368
  1. /* sgbmv.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 sgbmv_(char *trans, integer *m, integer *n, integer *kl,
  14. integer *ku, real *alpha, real *a, integer *lda, real *x, integer *
  15. incx, real *beta, real *y, integer *incy)
  16. {
  17. /* System generated locals */
  18. integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
  19. /* Local variables */
  20. integer i__, j, k, ix, iy, jx, jy, kx, ky, kup1, 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. /* SGBMV 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 band matrix, with kl sub-diagonals and ku super-diagonals. */
  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. /* KL - INTEGER. */
  53. /* On entry, KL specifies the number of sub-diagonals of the */
  54. /* matrix A. KL must satisfy 0 .le. KL. */
  55. /* Unchanged on exit. */
  56. /* KU - INTEGER. */
  57. /* On entry, KU specifies the number of super-diagonals of the */
  58. /* matrix A. KU must satisfy 0 .le. KU. */
  59. /* Unchanged on exit. */
  60. /* ALPHA - REAL . */
  61. /* On entry, ALPHA specifies the scalar alpha. */
  62. /* Unchanged on exit. */
  63. /* A - REAL array of DIMENSION ( LDA, n ). */
  64. /* Before entry, the leading ( kl + ku + 1 ) by n part of the */
  65. /* array A must contain the matrix of coefficients, supplied */
  66. /* column by column, with the leading diagonal of the matrix in */
  67. /* row ( ku + 1 ) of the array, the first super-diagonal */
  68. /* starting at position 2 in row ku, the first sub-diagonal */
  69. /* starting at position 1 in row ( ku + 2 ), and so on. */
  70. /* Elements in the array A that do not correspond to elements */
  71. /* in the band matrix (such as the top left ku by ku triangle) */
  72. /* are not referenced. */
  73. /* The following program segment will transfer a band matrix */
  74. /* from conventional full matrix storage to band storage: */
  75. /* DO 20, J = 1, N */
  76. /* K = KU + 1 - J */
  77. /* DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) */
  78. /* A( K + I, J ) = matrix( I, J ) */
  79. /* 10 CONTINUE */
  80. /* 20 CONTINUE */
  81. /* Unchanged on exit. */
  82. /* LDA - INTEGER. */
  83. /* On entry, LDA specifies the first dimension of A as declared */
  84. /* in the calling (sub) program. LDA must be at least */
  85. /* ( kl + ku + 1 ). */
  86. /* Unchanged on exit. */
  87. /* X - REAL array of DIMENSION at least */
  88. /* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' */
  89. /* and at least */
  90. /* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. */
  91. /* Before entry, the incremented array X must contain the */
  92. /* vector x. */
  93. /* Unchanged on exit. */
  94. /* INCX - INTEGER. */
  95. /* On entry, INCX specifies the increment for the elements of */
  96. /* X. INCX must not be zero. */
  97. /* Unchanged on exit. */
  98. /* BETA - REAL . */
  99. /* On entry, BETA specifies the scalar beta. When BETA is */
  100. /* supplied as zero then Y need not be set on input. */
  101. /* Unchanged on exit. */
  102. /* Y - REAL array of DIMENSION at least */
  103. /* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' */
  104. /* and at least */
  105. /* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. */
  106. /* Before entry, the incremented array Y must contain the */
  107. /* vector y. On exit, Y is overwritten by the updated vector y. */
  108. /* INCY - INTEGER. */
  109. /* On entry, INCY specifies the increment for the elements of */
  110. /* Y. INCY must not be zero. */
  111. /* Unchanged on exit. */
  112. /* Level 2 Blas routine. */
  113. /* -- Written on 22-October-1986. */
  114. /* Jack Dongarra, Argonne National Lab. */
  115. /* Jeremy Du Croz, Nag Central Office. */
  116. /* Sven Hammarling, Nag Central Office. */
  117. /* Richard Hanson, Sandia National Labs. */
  118. /* .. Parameters .. */
  119. /* .. */
  120. /* .. Local Scalars .. */
  121. /* .. */
  122. /* .. External Functions .. */
  123. /* .. */
  124. /* .. External Subroutines .. */
  125. /* .. */
  126. /* .. Intrinsic Functions .. */
  127. /* .. */
  128. /* Test the input parameters. */
  129. /* Parameter adjustments */
  130. a_dim1 = *lda;
  131. a_offset = 1 + a_dim1;
  132. a -= a_offset;
  133. --x;
  134. --y;
  135. /* Function Body */
  136. info = 0;
  137. if (! lsame_(trans, "N") && ! lsame_(trans, "T") && ! lsame_(trans, "C")
  138. ) {
  139. info = 1;
  140. } else if (*m < 0) {
  141. info = 2;
  142. } else if (*n < 0) {
  143. info = 3;
  144. } else if (*kl < 0) {
  145. info = 4;
  146. } else if (*ku < 0) {
  147. info = 5;
  148. } else if (*lda < *kl + *ku + 1) {
  149. info = 8;
  150. } else if (*incx == 0) {
  151. info = 10;
  152. } else if (*incy == 0) {
  153. info = 13;
  154. }
  155. if (info != 0) {
  156. xerbla_("SGBMV ", &info);
  157. return 0;
  158. }
  159. /* Quick return if possible. */
  160. if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
  161. return 0;
  162. }
  163. /* Set LENX and LENY, the lengths of the vectors x and y, and set */
  164. /* up the start points in X and Y. */
  165. if (lsame_(trans, "N")) {
  166. lenx = *n;
  167. leny = *m;
  168. } else {
  169. lenx = *m;
  170. leny = *n;
  171. }
  172. if (*incx > 0) {
  173. kx = 1;
  174. } else {
  175. kx = 1 - (lenx - 1) * *incx;
  176. }
  177. if (*incy > 0) {
  178. ky = 1;
  179. } else {
  180. ky = 1 - (leny - 1) * *incy;
  181. }
  182. /* Start the operations. In this version the elements of A are */
  183. /* accessed sequentially with one pass through the band part of A. */
  184. /* First form y := beta*y. */
  185. if (*beta != 1.f) {
  186. if (*incy == 1) {
  187. if (*beta == 0.f) {
  188. i__1 = leny;
  189. for (i__ = 1; i__ <= i__1; ++i__) {
  190. y[i__] = 0.f;
  191. /* L10: */
  192. }
  193. } else {
  194. i__1 = leny;
  195. for (i__ = 1; i__ <= i__1; ++i__) {
  196. y[i__] = *beta * y[i__];
  197. /* L20: */
  198. }
  199. }
  200. } else {
  201. iy = ky;
  202. if (*beta == 0.f) {
  203. i__1 = leny;
  204. for (i__ = 1; i__ <= i__1; ++i__) {
  205. y[iy] = 0.f;
  206. iy += *incy;
  207. /* L30: */
  208. }
  209. } else {
  210. i__1 = leny;
  211. for (i__ = 1; i__ <= i__1; ++i__) {
  212. y[iy] = *beta * y[iy];
  213. iy += *incy;
  214. /* L40: */
  215. }
  216. }
  217. }
  218. }
  219. if (*alpha == 0.f) {
  220. return 0;
  221. }
  222. kup1 = *ku + 1;
  223. if (lsame_(trans, "N")) {
  224. /* Form y := alpha*A*x + y. */
  225. jx = kx;
  226. if (*incy == 1) {
  227. i__1 = *n;
  228. for (j = 1; j <= i__1; ++j) {
  229. if (x[jx] != 0.f) {
  230. temp = *alpha * x[jx];
  231. k = kup1 - j;
  232. /* Computing MAX */
  233. i__2 = 1, i__3 = j - *ku;
  234. /* Computing MIN */
  235. i__5 = *m, i__6 = j + *kl;
  236. i__4 = min(i__5,i__6);
  237. for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
  238. y[i__] += temp * a[k + i__ + j * a_dim1];
  239. /* L50: */
  240. }
  241. }
  242. jx += *incx;
  243. /* L60: */
  244. }
  245. } else {
  246. i__1 = *n;
  247. for (j = 1; j <= i__1; ++j) {
  248. if (x[jx] != 0.f) {
  249. temp = *alpha * x[jx];
  250. iy = ky;
  251. k = kup1 - j;
  252. /* Computing MAX */
  253. i__4 = 1, i__2 = j - *ku;
  254. /* Computing MIN */
  255. i__5 = *m, i__6 = j + *kl;
  256. i__3 = min(i__5,i__6);
  257. for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
  258. y[iy] += temp * a[k + i__ + j * a_dim1];
  259. iy += *incy;
  260. /* L70: */
  261. }
  262. }
  263. jx += *incx;
  264. if (j > *ku) {
  265. ky += *incy;
  266. }
  267. /* L80: */
  268. }
  269. }
  270. } else {
  271. /* Form y := alpha*A'*x + y. */
  272. jy = ky;
  273. if (*incx == 1) {
  274. i__1 = *n;
  275. for (j = 1; j <= i__1; ++j) {
  276. temp = 0.f;
  277. k = kup1 - j;
  278. /* Computing MAX */
  279. i__3 = 1, i__4 = j - *ku;
  280. /* Computing MIN */
  281. i__5 = *m, i__6 = j + *kl;
  282. i__2 = min(i__5,i__6);
  283. for (i__ = max(i__3,i__4); i__ <= i__2; ++i__) {
  284. temp += a[k + i__ + j * a_dim1] * x[i__];
  285. /* L90: */
  286. }
  287. y[jy] += *alpha * temp;
  288. jy += *incy;
  289. /* L100: */
  290. }
  291. } else {
  292. i__1 = *n;
  293. for (j = 1; j <= i__1; ++j) {
  294. temp = 0.f;
  295. ix = kx;
  296. k = kup1 - j;
  297. /* Computing MAX */
  298. i__2 = 1, i__3 = j - *ku;
  299. /* Computing MIN */
  300. i__5 = *m, i__6 = j + *kl;
  301. i__4 = min(i__5,i__6);
  302. for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
  303. temp += a[k + i__ + j * a_dim1] * x[ix];
  304. ix += *incx;
  305. /* L110: */
  306. }
  307. y[jy] += *alpha * temp;
  308. jy += *incy;
  309. if (j > *ku) {
  310. kx += *incx;
  311. }
  312. /* L120: */
  313. }
  314. }
  315. }
  316. return 0;
  317. /* End of SGBMV . */
  318. } /* sgbmv_ */