You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

306 lines
9.4 KiB

  1. SUBROUTINE SSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
  2. * .. Scalar Arguments ..
  3. REAL ALPHA,BETA
  4. INTEGER INCX,INCY,K,LDA,N
  5. CHARACTER UPLO
  6. * ..
  7. * .. Array Arguments ..
  8. REAL A(LDA,*),X(*),Y(*)
  9. * ..
  10. *
  11. * Purpose
  12. * =======
  13. *
  14. * SSBMV performs the matrix-vector operation
  15. *
  16. * y := alpha*A*x + beta*y,
  17. *
  18. * where alpha and beta are scalars, x and y are n element vectors and
  19. * A is an n by n symmetric band matrix, with k super-diagonals.
  20. *
  21. * Arguments
  22. * ==========
  23. *
  24. * UPLO - CHARACTER*1.
  25. * On entry, UPLO specifies whether the upper or lower
  26. * triangular part of the band matrix A is being supplied as
  27. * follows:
  28. *
  29. * UPLO = 'U' or 'u' The upper triangular part of A is
  30. * being supplied.
  31. *
  32. * UPLO = 'L' or 'l' The lower triangular part of A is
  33. * being supplied.
  34. *
  35. * Unchanged on exit.
  36. *
  37. * N - INTEGER.
  38. * On entry, N specifies the order of the matrix A.
  39. * N must be at least zero.
  40. * Unchanged on exit.
  41. *
  42. * K - INTEGER.
  43. * On entry, K specifies the number of super-diagonals of the
  44. * matrix A. K must satisfy 0 .le. K.
  45. * Unchanged on exit.
  46. *
  47. * ALPHA - REAL .
  48. * On entry, ALPHA specifies the scalar alpha.
  49. * Unchanged on exit.
  50. *
  51. * A - REAL array of DIMENSION ( LDA, n ).
  52. * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
  53. * by n part of the array A must contain the upper triangular
  54. * band part of the symmetric matrix, supplied column by
  55. * column, with the leading diagonal of the matrix in row
  56. * ( k + 1 ) of the array, the first super-diagonal starting at
  57. * position 2 in row k, and so on. The top left k by k triangle
  58. * of the array A is not referenced.
  59. * The following program segment will transfer the upper
  60. * triangular part of a symmetric band matrix from conventional
  61. * full matrix storage to band storage:
  62. *
  63. * DO 20, J = 1, N
  64. * M = K + 1 - J
  65. * DO 10, I = MAX( 1, J - K ), J
  66. * A( M + I, J ) = matrix( I, J )
  67. * 10 CONTINUE
  68. * 20 CONTINUE
  69. *
  70. * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
  71. * by n part of the array A must contain the lower triangular
  72. * band part of the symmetric matrix, supplied column by
  73. * column, with the leading diagonal of the matrix in row 1 of
  74. * the array, the first sub-diagonal starting at position 1 in
  75. * row 2, and so on. The bottom right k by k triangle of the
  76. * array A is not referenced.
  77. * The following program segment will transfer the lower
  78. * triangular part of a symmetric band matrix from conventional
  79. * full matrix storage to band storage:
  80. *
  81. * DO 20, J = 1, N
  82. * M = 1 - J
  83. * DO 10, I = J, MIN( N, J + K )
  84. * A( M + I, J ) = matrix( I, J )
  85. * 10 CONTINUE
  86. * 20 CONTINUE
  87. *
  88. * Unchanged on exit.
  89. *
  90. * LDA - INTEGER.
  91. * On entry, LDA specifies the first dimension of A as declared
  92. * in the calling (sub) program. LDA must be at least
  93. * ( k + 1 ).
  94. * Unchanged on exit.
  95. *
  96. * X - REAL array of DIMENSION at least
  97. * ( 1 + ( n - 1 )*abs( INCX ) ).
  98. * Before entry, the incremented array X must contain the
  99. * vector x.
  100. * Unchanged on exit.
  101. *
  102. * INCX - INTEGER.
  103. * On entry, INCX specifies the increment for the elements of
  104. * X. INCX must not be zero.
  105. * Unchanged on exit.
  106. *
  107. * BETA - REAL .
  108. * On entry, BETA specifies the scalar beta.
  109. * Unchanged on exit.
  110. *
  111. * Y - REAL array of DIMENSION at least
  112. * ( 1 + ( n - 1 )*abs( INCY ) ).
  113. * Before entry, the incremented array Y must contain the
  114. * vector y. On exit, Y is overwritten by the updated vector y.
  115. *
  116. * INCY - INTEGER.
  117. * On entry, INCY specifies the increment for the elements of
  118. * Y. INCY must not be zero.
  119. * Unchanged on exit.
  120. *
  121. * Further Details
  122. * ===============
  123. *
  124. * Level 2 Blas routine.
  125. *
  126. * -- Written on 22-October-1986.
  127. * Jack Dongarra, Argonne National Lab.
  128. * Jeremy Du Croz, Nag Central Office.
  129. * Sven Hammarling, Nag Central Office.
  130. * Richard Hanson, Sandia National Labs.
  131. *
  132. * =====================================================================
  133. *
  134. * .. Parameters ..
  135. REAL ONE,ZERO
  136. PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
  137. * ..
  138. * .. Local Scalars ..
  139. REAL TEMP1,TEMP2
  140. INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
  141. * ..
  142. * .. External Functions ..
  143. LOGICAL LSAME
  144. EXTERNAL LSAME
  145. * ..
  146. * .. External Subroutines ..
  147. EXTERNAL XERBLA
  148. * ..
  149. * .. Intrinsic Functions ..
  150. INTRINSIC MAX,MIN
  151. * ..
  152. *
  153. * Test the input parameters.
  154. *
  155. INFO = 0
  156. IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  157. INFO = 1
  158. ELSE IF (N.LT.0) THEN
  159. INFO = 2
  160. ELSE IF (K.LT.0) THEN
  161. INFO = 3
  162. ELSE IF (LDA.LT. (K+1)) THEN
  163. INFO = 6
  164. ELSE IF (INCX.EQ.0) THEN
  165. INFO = 8
  166. ELSE IF (INCY.EQ.0) THEN
  167. INFO = 11
  168. END IF
  169. IF (INFO.NE.0) THEN
  170. CALL XERBLA('SSBMV ',INFO)
  171. RETURN
  172. END IF
  173. *
  174. * Quick return if possible.
  175. *
  176. IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
  177. *
  178. * Set up the start points in X and Y.
  179. *
  180. IF (INCX.GT.0) THEN
  181. KX = 1
  182. ELSE
  183. KX = 1 - (N-1)*INCX
  184. END IF
  185. IF (INCY.GT.0) THEN
  186. KY = 1
  187. ELSE
  188. KY = 1 - (N-1)*INCY
  189. END IF
  190. *
  191. * Start the operations. In this version the elements of the array A
  192. * are accessed sequentially with one pass through A.
  193. *
  194. * First form y := beta*y.
  195. *
  196. IF (BETA.NE.ONE) THEN
  197. IF (INCY.EQ.1) THEN
  198. IF (BETA.EQ.ZERO) THEN
  199. DO 10 I = 1,N
  200. Y(I) = ZERO
  201. 10 CONTINUE
  202. ELSE
  203. DO 20 I = 1,N
  204. Y(I) = BETA*Y(I)
  205. 20 CONTINUE
  206. END IF
  207. ELSE
  208. IY = KY
  209. IF (BETA.EQ.ZERO) THEN
  210. DO 30 I = 1,N
  211. Y(IY) = ZERO
  212. IY = IY + INCY
  213. 30 CONTINUE
  214. ELSE
  215. DO 40 I = 1,N
  216. Y(IY) = BETA*Y(IY)
  217. IY = IY + INCY
  218. 40 CONTINUE
  219. END IF
  220. END IF
  221. END IF
  222. IF (ALPHA.EQ.ZERO) RETURN
  223. IF (LSAME(UPLO,'U')) THEN
  224. *
  225. * Form y when upper triangle of A is stored.
  226. *
  227. KPLUS1 = K + 1
  228. IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  229. DO 60 J = 1,N
  230. TEMP1 = ALPHA*X(J)
  231. TEMP2 = ZERO
  232. L = KPLUS1 - J
  233. DO 50 I = MAX(1,J-K),J - 1
  234. Y(I) = Y(I) + TEMP1*A(L+I,J)
  235. TEMP2 = TEMP2 + A(L+I,J)*X(I)
  236. 50 CONTINUE
  237. Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
  238. 60 CONTINUE
  239. ELSE
  240. JX = KX
  241. JY = KY
  242. DO 80 J = 1,N
  243. TEMP1 = ALPHA*X(JX)
  244. TEMP2 = ZERO
  245. IX = KX
  246. IY = KY
  247. L = KPLUS1 - J
  248. DO 70 I = MAX(1,J-K),J - 1
  249. Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  250. TEMP2 = TEMP2 + A(L+I,J)*X(IX)
  251. IX = IX + INCX
  252. IY = IY + INCY
  253. 70 CONTINUE
  254. Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
  255. JX = JX + INCX
  256. JY = JY + INCY
  257. IF (J.GT.K) THEN
  258. KX = KX + INCX
  259. KY = KY + INCY
  260. END IF
  261. 80 CONTINUE
  262. END IF
  263. ELSE
  264. *
  265. * Form y when lower triangle of A is stored.
  266. *
  267. IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  268. DO 100 J = 1,N
  269. TEMP1 = ALPHA*X(J)
  270. TEMP2 = ZERO
  271. Y(J) = Y(J) + TEMP1*A(1,J)
  272. L = 1 - J
  273. DO 90 I = J + 1,MIN(N,J+K)
  274. Y(I) = Y(I) + TEMP1*A(L+I,J)
  275. TEMP2 = TEMP2 + A(L+I,J)*X(I)
  276. 90 CONTINUE
  277. Y(J) = Y(J) + ALPHA*TEMP2
  278. 100 CONTINUE
  279. ELSE
  280. JX = KX
  281. JY = KY
  282. DO 120 J = 1,N
  283. TEMP1 = ALPHA*X(JX)
  284. TEMP2 = ZERO
  285. Y(JY) = Y(JY) + TEMP1*A(1,J)
  286. L = 1 - J
  287. IX = JX
  288. IY = JY
  289. DO 110 I = J + 1,MIN(N,J+K)
  290. IX = IX + INCX
  291. IY = IY + INCY
  292. Y(IY) = Y(IY) + TEMP1*A(L+I,J)
  293. TEMP2 = TEMP2 + A(L+I,J)*X(IX)
  294. 110 CONTINUE
  295. Y(JY) = Y(JY) + ALPHA*TEMP2
  296. JX = JX + INCX
  297. JY = JY + INCY
  298. 120 CONTINUE
  299. END IF
  300. END IF
  301. *
  302. RETURN
  303. *
  304. * End of SSBMV .
  305. *
  306. END