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

ssyr2f.f 7.3 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230
  1. SUBROUTINE SSYR2F ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  2. * .. Scalar Arguments ..
  3. REAL ALPHA
  4. INTEGER INCX, INCY, LDA, N
  5. CHARACTER*1 UPLO
  6. * .. Array Arguments ..
  7. REAL A( LDA, * ), X( * ), Y( * )
  8. * ..
  9. *
  10. * Purpose
  11. * =======
  12. *
  13. * SSYR2 performs the symmetric rank 2 operation
  14. *
  15. * A := alpha*x*y' + alpha*y*x' + A,
  16. *
  17. * where alpha is a scalar, x and y are n element vectors and A is an n
  18. * by n symmetric matrix.
  19. *
  20. * Parameters
  21. * ==========
  22. *
  23. * UPLO - CHARACTER*1.
  24. * On entry, UPLO specifies whether the upper or lower
  25. * triangular part of the array A is to be referenced as
  26. * follows:
  27. *
  28. * UPLO = 'U' or 'u' Only the upper triangular part of A
  29. * is to be referenced.
  30. *
  31. * UPLO = 'L' or 'l' Only the lower triangular part of A
  32. * is to be referenced.
  33. *
  34. * Unchanged on exit.
  35. *
  36. * N - INTEGER.
  37. * On entry, N specifies the order of the matrix A.
  38. * N must be at least zero.
  39. * Unchanged on exit.
  40. *
  41. * ALPHA - REAL .
  42. * On entry, ALPHA specifies the scalar alpha.
  43. * Unchanged on exit.
  44. *
  45. * X - REAL array of dimension at least
  46. * ( 1 + ( n - 1 )*abs( INCX ) ).
  47. * Before entry, the incremented array X must contain the n
  48. * element vector x.
  49. * Unchanged on exit.
  50. *
  51. * INCX - INTEGER.
  52. * On entry, INCX specifies the increment for the elements of
  53. * X. INCX must not be zero.
  54. * Unchanged on exit.
  55. *
  56. * Y - REAL array of dimension at least
  57. * ( 1 + ( n - 1 )*abs( INCY ) ).
  58. * Before entry, the incremented array Y must contain the n
  59. * element vector y.
  60. * Unchanged on exit.
  61. *
  62. * INCY - INTEGER.
  63. * On entry, INCY specifies the increment for the elements of
  64. * Y. INCY must not be zero.
  65. * Unchanged on exit.
  66. *
  67. * A - REAL array of DIMENSION ( LDA, n ).
  68. * Before entry with UPLO = 'U' or 'u', the leading n by n
  69. * upper triangular part of the array A must contain the upper
  70. * triangular part of the symmetric matrix and the strictly
  71. * lower triangular part of A is not referenced. On exit, the
  72. * upper triangular part of the array A is overwritten by the
  73. * upper triangular part of the updated matrix.
  74. * Before entry with UPLO = 'L' or 'l', the leading n by n
  75. * lower triangular part of the array A must contain the lower
  76. * triangular part of the symmetric matrix and the strictly
  77. * upper triangular part of A is not referenced. On exit, the
  78. * lower triangular part of the array A is overwritten by the
  79. * lower triangular part of the updated matrix.
  80. *
  81. * LDA - INTEGER.
  82. * On entry, LDA specifies the first dimension of A as declared
  83. * in the calling (sub) program. LDA must be at least
  84. * max( 1, n ).
  85. * Unchanged on exit.
  86. *
  87. *
  88. * Level 2 Blas routine.
  89. *
  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. *
  96. *
  97. * .. Parameters ..
  98. REAL ZERO
  99. PARAMETER ( ZERO = 0.0E+0 )
  100. * .. Local Scalars ..
  101. REAL TEMP1, TEMP2
  102. INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY
  103. * .. External Functions ..
  104. LOGICAL LSAME
  105. EXTERNAL LSAME
  106. * .. External Subroutines ..
  107. EXTERNAL XERBLA
  108. * .. Intrinsic Functions ..
  109. INTRINSIC MAX
  110. * ..
  111. * .. Executable Statements ..
  112. *
  113. * Test the input parameters.
  114. *
  115. INFO = 0
  116. IF ( .NOT.LSAME( UPLO, 'U' ).AND.
  117. $ .NOT.LSAME( UPLO, 'L' ) )THEN
  118. INFO = 1
  119. ELSE IF( N.LT.0 )THEN
  120. INFO = 2
  121. ELSE IF( INCX.EQ.0 )THEN
  122. INFO = 5
  123. ELSE IF( INCY.EQ.0 )THEN
  124. INFO = 7
  125. ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  126. INFO = 9
  127. END IF
  128. IF( INFO.NE.0 )THEN
  129. CALL XERBLA( 'SSYR2 ', INFO )
  130. RETURN
  131. END IF
  132. *
  133. * Quick return if possible.
  134. *
  135. IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  136. $ RETURN
  137. *
  138. * Set up the start points in X and Y if the increments are not both
  139. * unity.
  140. *
  141. IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
  142. IF( INCX.GT.0 )THEN
  143. KX = 1
  144. ELSE
  145. KX = 1 - ( N - 1 )*INCX
  146. END IF
  147. IF( INCY.GT.0 )THEN
  148. KY = 1
  149. ELSE
  150. KY = 1 - ( N - 1 )*INCY
  151. END IF
  152. JX = KX
  153. JY = KY
  154. END IF
  155. *
  156. * Start the operations. In this version the elements of A are
  157. * accessed sequentially with one pass through the triangular part
  158. * of A.
  159. *
  160. IF( LSAME( UPLO, 'U' ) )THEN
  161. *
  162. * Form A when A is stored in the upper triangle.
  163. *
  164. IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  165. DO 20, J = 1, N
  166. IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  167. TEMP1 = ALPHA*Y( J )
  168. TEMP2 = ALPHA*X( J )
  169. DO 10, I = 1, J
  170. A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  171. 10 CONTINUE
  172. END IF
  173. 20 CONTINUE
  174. ELSE
  175. DO 40, J = 1, N
  176. IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  177. TEMP1 = ALPHA*Y( JY )
  178. TEMP2 = ALPHA*X( JX )
  179. IX = KX
  180. IY = KY
  181. DO 30, I = 1, J
  182. A( I, J ) = A( I, J ) + X( IX )*TEMP1
  183. $ + Y( IY )*TEMP2
  184. IX = IX + INCX
  185. IY = IY + INCY
  186. 30 CONTINUE
  187. END IF
  188. JX = JX + INCX
  189. JY = JY + INCY
  190. 40 CONTINUE
  191. END IF
  192. ELSE
  193. *
  194. * Form A when A is stored in the lower triangle.
  195. *
  196. IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  197. DO 60, J = 1, N
  198. IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  199. TEMP1 = ALPHA*Y( J )
  200. TEMP2 = ALPHA*X( J )
  201. DO 50, I = J, N
  202. A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  203. 50 CONTINUE
  204. END IF
  205. 60 CONTINUE
  206. ELSE
  207. DO 80, J = 1, N
  208. IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  209. TEMP1 = ALPHA*Y( JY )
  210. TEMP2 = ALPHA*X( JX )
  211. IX = JX
  212. IY = JY
  213. DO 70, I = J, N
  214. A( I, J ) = A( I, J ) + X( IX )*TEMP1
  215. $ + Y( IY )*TEMP2
  216. IX = IX + INCX
  217. IY = IY + INCY
  218. 70 CONTINUE
  219. END IF
  220. JX = JX + INCX
  221. JY = JY + INCY
  222. 80 CONTINUE
  223. END IF
  224. END IF
  225. *
  226. RETURN
  227. *
  228. * End of SSYR2 .
  229. *
  230. END