ssyrk.f 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299
  1. *DECK SSYRK
  2. SUBROUTINE SSYRK (UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
  3. C***BEGIN PROLOGUE SSYRK
  4. C***PURPOSE Perform symmetric rank k update of a real symmetric matrix.
  5. C***LIBRARY SLATEC (BLAS)
  6. C***CATEGORY D1B6
  7. C***TYPE SINGLE PRECISION (SSYRK-S, DSYRK-D, CSYRK-C)
  8. C***KEYWORDS LEVEL 3 BLAS, LINEAR ALGEBRA
  9. C***AUTHOR Dongarra, J., (ANL)
  10. C Duff, I., (AERE)
  11. C Du Croz, J., (NAG)
  12. C Hammarling, S. (NAG)
  13. C***DESCRIPTION
  14. C
  15. C SSYRK performs one of the symmetric rank k operations
  16. C
  17. C C := alpha*A*A' + beta*C,
  18. C
  19. C or
  20. C
  21. C C := alpha*A'*A + beta*C,
  22. C
  23. C where alpha and beta are scalars, C is an n by n symmetric matrix
  24. C and A is an n by k matrix in the first case and a k by n matrix
  25. C in the second case.
  26. C
  27. C Parameters
  28. C ==========
  29. C
  30. C UPLO - CHARACTER*1.
  31. C On entry, UPLO specifies whether the upper or lower
  32. C triangular part of the array C is to be referenced as
  33. C follows:
  34. C
  35. C UPLO = 'U' or 'u' Only the upper triangular part of C
  36. C is to be referenced.
  37. C
  38. C UPLO = 'L' or 'l' Only the lower triangular part of C
  39. C is to be referenced.
  40. C
  41. C Unchanged on exit.
  42. C
  43. C TRANS - CHARACTER*1.
  44. C On entry, TRANS specifies the operation to be performed as
  45. C follows:
  46. C
  47. C TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
  48. C
  49. C TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
  50. C
  51. C TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
  52. C
  53. C Unchanged on exit.
  54. C
  55. C N - INTEGER.
  56. C On entry, N specifies the order of the matrix C. N must be
  57. C at least zero.
  58. C Unchanged on exit.
  59. C
  60. C K - INTEGER.
  61. C On entry with TRANS = 'N' or 'n', K specifies the number
  62. C of columns of the matrix A, and on entry with
  63. C TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
  64. C of rows of the matrix A. K must be at least zero.
  65. C Unchanged on exit.
  66. C
  67. C ALPHA - REAL .
  68. C On entry, ALPHA specifies the scalar alpha.
  69. C Unchanged on exit.
  70. C
  71. C A - REAL array of DIMENSION ( LDA, ka ), where ka is
  72. C k when TRANS = 'N' or 'n', and is n otherwise.
  73. C Before entry with TRANS = 'N' or 'n', the leading n by k
  74. C part of the array A must contain the matrix A, otherwise
  75. C the leading k by n part of the array A must contain the
  76. C matrix A.
  77. C Unchanged on exit.
  78. C
  79. C LDA - INTEGER.
  80. C On entry, LDA specifies the first dimension of A as declared
  81. C in the calling (sub) program. When TRANS = 'N' or 'n'
  82. C then LDA must be at least max( 1, n ), otherwise LDA must
  83. C be at least max( 1, k ).
  84. C Unchanged on exit.
  85. C
  86. C BETA - REAL .
  87. C On entry, BETA specifies the scalar beta.
  88. C Unchanged on exit.
  89. C
  90. C C - REAL array of DIMENSION ( LDC, n ).
  91. C Before entry with UPLO = 'U' or 'u', the leading n by n
  92. C upper triangular part of the array C must contain the upper
  93. C triangular part of the symmetric matrix and the strictly
  94. C lower triangular part of C is not referenced. On exit, the
  95. C upper triangular part of the array C is overwritten by the
  96. C upper triangular part of the updated matrix.
  97. C Before entry with UPLO = 'L' or 'l', the leading n by n
  98. C lower triangular part of the array C must contain the lower
  99. C triangular part of the symmetric matrix and the strictly
  100. C upper triangular part of C is not referenced. On exit, the
  101. C lower triangular part of the array C is overwritten by the
  102. C lower triangular part of the updated matrix.
  103. C
  104. C LDC - INTEGER.
  105. C On entry, LDC specifies the first dimension of C as declared
  106. C in the calling (sub) program. LDC must be at least
  107. C max( 1, n ).
  108. C Unchanged on exit.
  109. C
  110. C***REFERENCES Dongarra, J., Du Croz, J., Duff, I., and Hammarling, S.
  111. C A set of level 3 basic linear algebra subprograms.
  112. C ACM TOMS, Vol. 16, No. 1, pp. 1-17, March 1990.
  113. C***ROUTINES CALLED LSAME, XERBLA
  114. C***REVISION HISTORY (YYMMDD)
  115. C 890208 DATE WRITTEN
  116. C 910605 Modified to meet SLATEC prologue standards. Only comment
  117. C lines were modified. (BKS)
  118. C***END PROLOGUE SSYRK
  119. C .. Scalar Arguments ..
  120. CHARACTER*1 UPLO, TRANS
  121. INTEGER N, K, LDA, LDC
  122. REAL ALPHA, BETA
  123. C .. Array Arguments ..
  124. REAL A( LDA, * ), C( LDC, * )
  125. C .. External Functions ..
  126. LOGICAL LSAME
  127. EXTERNAL LSAME
  128. C .. External Subroutines ..
  129. EXTERNAL XERBLA
  130. C .. Intrinsic Functions ..
  131. INTRINSIC MAX
  132. C .. Local Scalars ..
  133. LOGICAL UPPER
  134. INTEGER I, INFO, J, L, NROWA
  135. REAL TEMP
  136. C .. Parameters ..
  137. REAL ONE , ZERO
  138. PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  139. C***FIRST EXECUTABLE STATEMENT SSYRK
  140. C
  141. C Test the input parameters.
  142. C
  143. IF( LSAME( TRANS, 'N' ) )THEN
  144. NROWA = N
  145. ELSE
  146. NROWA = K
  147. END IF
  148. UPPER = LSAME( UPLO, 'U' )
  149. C
  150. INFO = 0
  151. IF( ( .NOT.UPPER ).AND.
  152. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
  153. INFO = 1
  154. ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
  155. $ ( .NOT.LSAME( TRANS, 'T' ) ).AND.
  156. $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN
  157. INFO = 2
  158. ELSE IF( N .LT.0 )THEN
  159. INFO = 3
  160. ELSE IF( K .LT.0 )THEN
  161. INFO = 4
  162. ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  163. INFO = 7
  164. ELSE IF( LDC.LT.MAX( 1, N ) )THEN
  165. INFO = 10
  166. END IF
  167. IF( INFO.NE.0 )THEN
  168. CALL XERBLA( 'SSYRK ', INFO )
  169. RETURN
  170. END IF
  171. C
  172. C Quick return if possible.
  173. C
  174. IF( ( N.EQ.0 ).OR.
  175. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
  176. $ RETURN
  177. C
  178. C And when alpha.eq.zero.
  179. C
  180. IF( ALPHA.EQ.ZERO )THEN
  181. IF( UPPER )THEN
  182. IF( BETA.EQ.ZERO )THEN
  183. DO 20, J = 1, N
  184. DO 10, I = 1, J
  185. C( I, J ) = ZERO
  186. 10 CONTINUE
  187. 20 CONTINUE
  188. ELSE
  189. DO 40, J = 1, N
  190. DO 30, I = 1, J
  191. C( I, J ) = BETA*C( I, J )
  192. 30 CONTINUE
  193. 40 CONTINUE
  194. END IF
  195. ELSE
  196. IF( BETA.EQ.ZERO )THEN
  197. DO 60, J = 1, N
  198. DO 50, I = J, N
  199. C( I, J ) = ZERO
  200. 50 CONTINUE
  201. 60 CONTINUE
  202. ELSE
  203. DO 80, J = 1, N
  204. DO 70, I = J, N
  205. C( I, J ) = BETA*C( I, J )
  206. 70 CONTINUE
  207. 80 CONTINUE
  208. END IF
  209. END IF
  210. RETURN
  211. END IF
  212. C
  213. C Start the operations.
  214. C
  215. IF( LSAME( TRANS, 'N' ) )THEN
  216. C
  217. C Form C := alpha*A*A' + beta*C.
  218. C
  219. IF( UPPER )THEN
  220. DO 130, J = 1, N
  221. IF( BETA.EQ.ZERO )THEN
  222. DO 90, I = 1, J
  223. C( I, J ) = ZERO
  224. 90 CONTINUE
  225. ELSE IF( BETA.NE.ONE )THEN
  226. DO 100, I = 1, J
  227. C( I, J ) = BETA*C( I, J )
  228. 100 CONTINUE
  229. END IF
  230. DO 120, L = 1, K
  231. IF( A( J, L ).NE.ZERO )THEN
  232. TEMP = ALPHA*A( J, L )
  233. DO 110, I = 1, J
  234. C( I, J ) = C( I, J ) + TEMP*A( I, L )
  235. 110 CONTINUE
  236. END IF
  237. 120 CONTINUE
  238. 130 CONTINUE
  239. ELSE
  240. DO 180, J = 1, N
  241. IF( BETA.EQ.ZERO )THEN
  242. DO 140, I = J, N
  243. C( I, J ) = ZERO
  244. 140 CONTINUE
  245. ELSE IF( BETA.NE.ONE )THEN
  246. DO 150, I = J, N
  247. C( I, J ) = BETA*C( I, J )
  248. 150 CONTINUE
  249. END IF
  250. DO 170, L = 1, K
  251. IF( A( J, L ).NE.ZERO )THEN
  252. TEMP = ALPHA*A( J, L )
  253. DO 160, I = J, N
  254. C( I, J ) = C( I, J ) + TEMP*A( I, L )
  255. 160 CONTINUE
  256. END IF
  257. 170 CONTINUE
  258. 180 CONTINUE
  259. END IF
  260. ELSE
  261. C
  262. C Form C := alpha*A'*A + beta*C.
  263. C
  264. IF( UPPER )THEN
  265. DO 210, J = 1, N
  266. DO 200, I = 1, J
  267. TEMP = ZERO
  268. DO 190, L = 1, K
  269. TEMP = TEMP + A( L, I )*A( L, J )
  270. 190 CONTINUE
  271. IF( BETA.EQ.ZERO )THEN
  272. C( I, J ) = ALPHA*TEMP
  273. ELSE
  274. C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  275. END IF
  276. 200 CONTINUE
  277. 210 CONTINUE
  278. ELSE
  279. DO 240, J = 1, N
  280. DO 230, I = J, N
  281. TEMP = ZERO
  282. DO 220, L = 1, K
  283. TEMP = TEMP + A( L, I )*A( L, J )
  284. 220 CONTINUE
  285. IF( BETA.EQ.ZERO )THEN
  286. C( I, J ) = ALPHA*TEMP
  287. ELSE
  288. C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
  289. END IF
  290. 230 CONTINUE
  291. 240 CONTINUE
  292. END IF
  293. END IF
  294. C
  295. RETURN
  296. C
  297. C End of SSYRK .
  298. C
  299. END