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.

sspr.f 7.2 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
  1. *> \brief \b SSPR
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. * Definition:
  9. * ===========
  10. *
  11. * SUBROUTINE SSPR(UPLO,N,ALPHA,X,INCX,AP)
  12. *
  13. * .. Scalar Arguments ..
  14. * REAL ALPHA
  15. * INTEGER INCX,N
  16. * CHARACTER UPLO
  17. * ..
  18. * .. Array Arguments ..
  19. * REAL AP(*),X(*)
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> SSPR performs the symmetric rank 1 operation
  29. *>
  30. *> A := alpha*x*x**T + A,
  31. *>
  32. *> where alpha is a real scalar, x is an n element vector and A is an
  33. *> n by n symmetric matrix, supplied in packed form.
  34. *> \endverbatim
  35. *
  36. * Arguments:
  37. * ==========
  38. *
  39. *> \param[in] UPLO
  40. *> \verbatim
  41. *> UPLO is CHARACTER*1
  42. *> On entry, UPLO specifies whether the upper or lower
  43. *> triangular part of the matrix A is supplied in the packed
  44. *> array AP as follows:
  45. *>
  46. *> UPLO = 'U' or 'u' The upper triangular part of A is
  47. *> supplied in AP.
  48. *>
  49. *> UPLO = 'L' or 'l' The lower triangular part of A is
  50. *> supplied in AP.
  51. *> \endverbatim
  52. *>
  53. *> \param[in] N
  54. *> \verbatim
  55. *> N is INTEGER
  56. *> On entry, N specifies the order of the matrix A.
  57. *> N must be at least zero.
  58. *> \endverbatim
  59. *>
  60. *> \param[in] ALPHA
  61. *> \verbatim
  62. *> ALPHA is REAL
  63. *> On entry, ALPHA specifies the scalar alpha.
  64. *> \endverbatim
  65. *>
  66. *> \param[in] X
  67. *> \verbatim
  68. *> X is REAL array of dimension at least
  69. *> ( 1 + ( n - 1 )*abs( INCX ) ).
  70. *> Before entry, the incremented array X must contain the n
  71. *> element vector x.
  72. *> \endverbatim
  73. *>
  74. *> \param[in] INCX
  75. *> \verbatim
  76. *> INCX is INTEGER
  77. *> On entry, INCX specifies the increment for the elements of
  78. *> X. INCX must not be zero.
  79. *> \endverbatim
  80. *>
  81. *> \param[in,out] AP
  82. *> \verbatim
  83. *> AP is REAL array of DIMENSION at least
  84. *> ( ( n*( n + 1 ) )/2 ).
  85. *> Before entry with UPLO = 'U' or 'u', the array AP must
  86. *> contain the upper triangular part of the symmetric matrix
  87. *> packed sequentially, column by column, so that AP( 1 )
  88. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  89. *> and a( 2, 2 ) respectively, and so on. On exit, the array
  90. *> AP is overwritten by the upper triangular part of the
  91. *> updated matrix.
  92. *> Before entry with UPLO = 'L' or 'l', the array AP must
  93. *> contain the lower triangular part of the symmetric matrix
  94. *> packed sequentially, column by column, so that AP( 1 )
  95. *> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  96. *> and a( 3, 1 ) respectively, and so on. On exit, the array
  97. *> AP is overwritten by the lower triangular part of the
  98. *> updated matrix.
  99. *> \endverbatim
  100. *
  101. * Authors:
  102. * ========
  103. *
  104. *> \author Univ. of Tennessee
  105. *> \author Univ. of California Berkeley
  106. *> \author Univ. of Colorado Denver
  107. *> \author NAG Ltd.
  108. *
  109. *> \date November 2011
  110. *
  111. *> \ingroup single_blas_level2
  112. *
  113. *> \par Further Details:
  114. * =====================
  115. *>
  116. *> \verbatim
  117. *>
  118. *> Level 2 Blas routine.
  119. *>
  120. *> -- Written on 22-October-1986.
  121. *> Jack Dongarra, Argonne National Lab.
  122. *> Jeremy Du Croz, Nag Central Office.
  123. *> Sven Hammarling, Nag Central Office.
  124. *> Richard Hanson, Sandia National Labs.
  125. *> \endverbatim
  126. *>
  127. * =====================================================================
  128. SUBROUTINE SSPR(UPLO,N,ALPHA,X,INCX,AP)
  129. *
  130. * -- Reference BLAS level2 routine (version 3.4.0) --
  131. * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
  132. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  133. * November 2011
  134. *
  135. * .. Scalar Arguments ..
  136. REAL ALPHA
  137. INTEGER INCX,N
  138. CHARACTER UPLO
  139. * ..
  140. * .. Array Arguments ..
  141. REAL AP(*),X(*)
  142. * ..
  143. *
  144. * =====================================================================
  145. *
  146. * .. Parameters ..
  147. REAL ZERO
  148. PARAMETER (ZERO=0.0E+0)
  149. * ..
  150. * .. Local Scalars ..
  151. REAL TEMP
  152. INTEGER I,INFO,IX,J,JX,K,KK,KX
  153. * ..
  154. * .. External Functions ..
  155. LOGICAL LSAME
  156. EXTERNAL LSAME
  157. * ..
  158. * .. External Subroutines ..
  159. EXTERNAL XERBLA
  160. * ..
  161. *
  162. * Test the input parameters.
  163. *
  164. INFO = 0
  165. IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  166. INFO = 1
  167. ELSE IF (N.LT.0) THEN
  168. INFO = 2
  169. ELSE IF (INCX.EQ.0) THEN
  170. INFO = 5
  171. END IF
  172. IF (INFO.NE.0) THEN
  173. CALL XERBLA('SSPR ',INFO)
  174. RETURN
  175. END IF
  176. *
  177. * Quick return if possible.
  178. *
  179. IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  180. *
  181. * Set the start point in X if the increment is not unity.
  182. *
  183. IF (INCX.LE.0) THEN
  184. KX = 1 - (N-1)*INCX
  185. ELSE IF (INCX.NE.1) THEN
  186. KX = 1
  187. END IF
  188. *
  189. * Start the operations. In this version the elements of the array AP
  190. * are accessed sequentially with one pass through AP.
  191. *
  192. KK = 1
  193. IF (LSAME(UPLO,'U')) THEN
  194. *
  195. * Form A when upper triangle is stored in AP.
  196. *
  197. IF (INCX.EQ.1) THEN
  198. DO 20 J = 1,N
  199. IF (X(J).NE.ZERO) THEN
  200. TEMP = ALPHA*X(J)
  201. K = KK
  202. DO 10 I = 1,J
  203. AP(K) = AP(K) + X(I)*TEMP
  204. K = K + 1
  205. 10 CONTINUE
  206. END IF
  207. KK = KK + J
  208. 20 CONTINUE
  209. ELSE
  210. JX = KX
  211. DO 40 J = 1,N
  212. IF (X(JX).NE.ZERO) THEN
  213. TEMP = ALPHA*X(JX)
  214. IX = KX
  215. DO 30 K = KK,KK + J - 1
  216. AP(K) = AP(K) + X(IX)*TEMP
  217. IX = IX + INCX
  218. 30 CONTINUE
  219. END IF
  220. JX = JX + INCX
  221. KK = KK + J
  222. 40 CONTINUE
  223. END IF
  224. ELSE
  225. *
  226. * Form A when lower triangle is stored in AP.
  227. *
  228. IF (INCX.EQ.1) THEN
  229. DO 60 J = 1,N
  230. IF (X(J).NE.ZERO) THEN
  231. TEMP = ALPHA*X(J)
  232. K = KK
  233. DO 50 I = J,N
  234. AP(K) = AP(K) + X(I)*TEMP
  235. K = K + 1
  236. 50 CONTINUE
  237. END IF
  238. KK = KK + N - J + 1
  239. 60 CONTINUE
  240. ELSE
  241. JX = KX
  242. DO 80 J = 1,N
  243. IF (X(JX).NE.ZERO) THEN
  244. TEMP = ALPHA*X(JX)
  245. IX = JX
  246. DO 70 K = KK,KK + N - J
  247. AP(K) = AP(K) + X(IX)*TEMP
  248. IX = IX + INCX
  249. 70 CONTINUE
  250. END IF
  251. JX = JX + INCX
  252. KK = KK + N - J + 1
  253. 80 CONTINUE
  254. END IF
  255. END IF
  256. *
  257. RETURN
  258. *
  259. * End of SSPR .
  260. *
  261. END