ssymv.c 7.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313
  1. /* ssymv.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 ssymv_(char *uplo, integer *n, real *alpha, real *a,
  14. integer *lda, real *x, integer *incx, real *beta, real *y, integer *
  15. 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 temp1, temp2;
  22. extern logical lsame_(char *, char *);
  23. extern /* Subroutine */ int xerbla_(char *, integer *);
  24. /* .. Scalar Arguments .. */
  25. /* .. */
  26. /* .. Array Arguments .. */
  27. /* .. */
  28. /* Purpose */
  29. /* ======= */
  30. /* SSYMV performs the matrix-vector operation */
  31. /* y := alpha*A*x + beta*y, */
  32. /* where alpha and beta are scalars, x and y are n element vectors and */
  33. /* A is an n by n symmetric matrix. */
  34. /* Arguments */
  35. /* ========== */
  36. /* UPLO - CHARACTER*1. */
  37. /* On entry, UPLO specifies whether the upper or lower */
  38. /* triangular part of the array A is to be referenced as */
  39. /* follows: */
  40. /* UPLO = 'U' or 'u' Only the upper triangular part of A */
  41. /* is to be referenced. */
  42. /* UPLO = 'L' or 'l' Only the lower triangular part of A */
  43. /* is to be referenced. */
  44. /* Unchanged on exit. */
  45. /* N - INTEGER. */
  46. /* On entry, N specifies the order of the matrix A. */
  47. /* N must be at least zero. */
  48. /* Unchanged on exit. */
  49. /* ALPHA - REAL . */
  50. /* On entry, ALPHA specifies the scalar alpha. */
  51. /* Unchanged on exit. */
  52. /* A - REAL array of DIMENSION ( LDA, n ). */
  53. /* Before entry with UPLO = 'U' or 'u', the leading n by n */
  54. /* upper triangular part of the array A must contain the upper */
  55. /* triangular part of the symmetric matrix and the strictly */
  56. /* lower triangular part of A is not referenced. */
  57. /* Before entry with UPLO = 'L' or 'l', the leading n by n */
  58. /* lower triangular part of the array A must contain the lower */
  59. /* triangular part of the symmetric matrix and the strictly */
  60. /* upper triangular part of A is not referenced. */
  61. /* Unchanged on exit. */
  62. /* LDA - INTEGER. */
  63. /* On entry, LDA specifies the first dimension of A as declared */
  64. /* in the calling (sub) program. LDA must be at least */
  65. /* max( 1, n ). */
  66. /* Unchanged on exit. */
  67. /* X - REAL array of dimension at least */
  68. /* ( 1 + ( n - 1 )*abs( INCX ) ). */
  69. /* Before entry, the incremented array X must contain the n */
  70. /* element vector x. */
  71. /* Unchanged on exit. */
  72. /* INCX - INTEGER. */
  73. /* On entry, INCX specifies the increment for the elements of */
  74. /* X. INCX must not be zero. */
  75. /* Unchanged on exit. */
  76. /* BETA - REAL . */
  77. /* On entry, BETA specifies the scalar beta. When BETA is */
  78. /* supplied as zero then Y need not be set on input. */
  79. /* Unchanged on exit. */
  80. /* Y - REAL array of dimension at least */
  81. /* ( 1 + ( n - 1 )*abs( INCY ) ). */
  82. /* Before entry, the incremented array Y must contain the n */
  83. /* element vector y. On exit, Y is overwritten by the updated */
  84. /* vector y. */
  85. /* INCY - INTEGER. */
  86. /* On entry, INCY specifies the increment for the elements of */
  87. /* Y. INCY 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. --y;
  112. /* Function Body */
  113. info = 0;
  114. if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
  115. info = 1;
  116. } else if (*n < 0) {
  117. info = 2;
  118. } else if (*lda < max(1,*n)) {
  119. info = 5;
  120. } else if (*incx == 0) {
  121. info = 7;
  122. } else if (*incy == 0) {
  123. info = 10;
  124. }
  125. if (info != 0) {
  126. xerbla_("SSYMV ", &info);
  127. return 0;
  128. }
  129. /* Quick return if possible. */
  130. if (*n == 0 || *alpha == 0.f && *beta == 1.f) {
  131. return 0;
  132. }
  133. /* Set up the start points in X and Y. */
  134. if (*incx > 0) {
  135. kx = 1;
  136. } else {
  137. kx = 1 - (*n - 1) * *incx;
  138. }
  139. if (*incy > 0) {
  140. ky = 1;
  141. } else {
  142. ky = 1 - (*n - 1) * *incy;
  143. }
  144. /* Start the operations. In this version the elements of A are */
  145. /* accessed sequentially with one pass through the triangular part */
  146. /* of A. */
  147. /* First form y := beta*y. */
  148. if (*beta != 1.f) {
  149. if (*incy == 1) {
  150. if (*beta == 0.f) {
  151. i__1 = *n;
  152. for (i__ = 1; i__ <= i__1; ++i__) {
  153. y[i__] = 0.f;
  154. /* L10: */
  155. }
  156. } else {
  157. i__1 = *n;
  158. for (i__ = 1; i__ <= i__1; ++i__) {
  159. y[i__] = *beta * y[i__];
  160. /* L20: */
  161. }
  162. }
  163. } else {
  164. iy = ky;
  165. if (*beta == 0.f) {
  166. i__1 = *n;
  167. for (i__ = 1; i__ <= i__1; ++i__) {
  168. y[iy] = 0.f;
  169. iy += *incy;
  170. /* L30: */
  171. }
  172. } else {
  173. i__1 = *n;
  174. for (i__ = 1; i__ <= i__1; ++i__) {
  175. y[iy] = *beta * y[iy];
  176. iy += *incy;
  177. /* L40: */
  178. }
  179. }
  180. }
  181. }
  182. if (*alpha == 0.f) {
  183. return 0;
  184. }
  185. if (lsame_(uplo, "U")) {
  186. /* Form y when A is stored in upper triangle. */
  187. if (*incx == 1 && *incy == 1) {
  188. i__1 = *n;
  189. for (j = 1; j <= i__1; ++j) {
  190. temp1 = *alpha * x[j];
  191. temp2 = 0.f;
  192. i__2 = j - 1;
  193. for (i__ = 1; i__ <= i__2; ++i__) {
  194. y[i__] += temp1 * a[i__ + j * a_dim1];
  195. temp2 += a[i__ + j * a_dim1] * x[i__];
  196. /* L50: */
  197. }
  198. y[j] = y[j] + temp1 * a[j + j * a_dim1] + *alpha * temp2;
  199. /* L60: */
  200. }
  201. } else {
  202. jx = kx;
  203. jy = ky;
  204. i__1 = *n;
  205. for (j = 1; j <= i__1; ++j) {
  206. temp1 = *alpha * x[jx];
  207. temp2 = 0.f;
  208. ix = kx;
  209. iy = ky;
  210. i__2 = j - 1;
  211. for (i__ = 1; i__ <= i__2; ++i__) {
  212. y[iy] += temp1 * a[i__ + j * a_dim1];
  213. temp2 += a[i__ + j * a_dim1] * x[ix];
  214. ix += *incx;
  215. iy += *incy;
  216. /* L70: */
  217. }
  218. y[jy] = y[jy] + temp1 * a[j + j * a_dim1] + *alpha * temp2;
  219. jx += *incx;
  220. jy += *incy;
  221. /* L80: */
  222. }
  223. }
  224. } else {
  225. /* Form y when A is stored in lower triangle. */
  226. if (*incx == 1 && *incy == 1) {
  227. i__1 = *n;
  228. for (j = 1; j <= i__1; ++j) {
  229. temp1 = *alpha * x[j];
  230. temp2 = 0.f;
  231. y[j] += temp1 * a[j + j * a_dim1];
  232. i__2 = *n;
  233. for (i__ = j + 1; i__ <= i__2; ++i__) {
  234. y[i__] += temp1 * a[i__ + j * a_dim1];
  235. temp2 += a[i__ + j * a_dim1] * x[i__];
  236. /* L90: */
  237. }
  238. y[j] += *alpha * temp2;
  239. /* L100: */
  240. }
  241. } else {
  242. jx = kx;
  243. jy = ky;
  244. i__1 = *n;
  245. for (j = 1; j <= i__1; ++j) {
  246. temp1 = *alpha * x[jx];
  247. temp2 = 0.f;
  248. y[jy] += temp1 * a[j + j * a_dim1];
  249. ix = jx;
  250. iy = jy;
  251. i__2 = *n;
  252. for (i__ = j + 1; i__ <= i__2; ++i__) {
  253. ix += *incx;
  254. iy += *incy;
  255. y[iy] += temp1 * a[i__ + j * a_dim1];
  256. temp2 += a[i__ + j * a_dim1] * x[ix];
  257. /* L110: */
  258. }
  259. y[jy] += *alpha * temp2;
  260. jx += *incx;
  261. jy += *incy;
  262. /* L120: */
  263. }
  264. }
  265. }
  266. return 0;
  267. /* End of SSYMV . */
  268. } /* ssymv_ */