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.

clagsy.f 7.9 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286
  1. *> \brief \b CLAGSY
  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 CLAGSY( N, K, D, A, LDA, ISEED, WORK, INFO )
  12. *
  13. * .. Scalar Arguments ..
  14. * INTEGER INFO, K, LDA, N
  15. * ..
  16. * .. Array Arguments ..
  17. * INTEGER ISEED( 4 )
  18. * REAL D( * )
  19. * COMPLEX A( LDA, * ), WORK( * )
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> CLAGSY generates a complex symmetric matrix A, by pre- and post-
  29. *> multiplying a real diagonal matrix D with a random unitary matrix:
  30. *> A = U*D*U**T. The semi-bandwidth may then be reduced to k by
  31. *> additional unitary transformations.
  32. *> \endverbatim
  33. *
  34. * Arguments:
  35. * ==========
  36. *
  37. *> \param[in] N
  38. *> \verbatim
  39. *> N is INTEGER
  40. *> The order of the matrix A. N >= 0.
  41. *> \endverbatim
  42. *>
  43. *> \param[in] K
  44. *> \verbatim
  45. *> K is INTEGER
  46. *> The number of nonzero subdiagonals within the band of A.
  47. *> 0 <= K <= N-1.
  48. *> \endverbatim
  49. *>
  50. *> \param[in] D
  51. *> \verbatim
  52. *> D is REAL array, dimension (N)
  53. *> The diagonal elements of the diagonal matrix D.
  54. *> \endverbatim
  55. *>
  56. *> \param[out] A
  57. *> \verbatim
  58. *> A is COMPLEX array, dimension (LDA,N)
  59. *> The generated n by n symmetric matrix A (the full matrix is
  60. *> stored).
  61. *> \endverbatim
  62. *>
  63. *> \param[in] LDA
  64. *> \verbatim
  65. *> LDA is INTEGER
  66. *> The leading dimension of the array A. LDA >= N.
  67. *> \endverbatim
  68. *>
  69. *> \param[in,out] ISEED
  70. *> \verbatim
  71. *> ISEED is INTEGER array, dimension (4)
  72. *> On entry, the seed of the random number generator; the array
  73. *> elements must be between 0 and 4095, and ISEED(4) must be
  74. *> odd.
  75. *> On exit, the seed is updated.
  76. *> \endverbatim
  77. *>
  78. *> \param[out] WORK
  79. *> \verbatim
  80. *> WORK is COMPLEX array, dimension (2*N)
  81. *> \endverbatim
  82. *>
  83. *> \param[out] INFO
  84. *> \verbatim
  85. *> INFO is INTEGER
  86. *> = 0: successful exit
  87. *> < 0: if INFO = -i, the i-th argument had an illegal value
  88. *> \endverbatim
  89. *
  90. * Authors:
  91. * ========
  92. *
  93. *> \author Univ. of Tennessee
  94. *> \author Univ. of California Berkeley
  95. *> \author Univ. of Colorado Denver
  96. *> \author NAG Ltd.
  97. *
  98. *> \date December 2016
  99. *
  100. *> \ingroup complex_matgen
  101. *
  102. * =====================================================================
  103. SUBROUTINE CLAGSY( N, K, D, A, LDA, ISEED, WORK, INFO )
  104. *
  105. * -- LAPACK auxiliary routine (version 3.7.0) --
  106. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  107. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  108. * December 2016
  109. *
  110. * .. Scalar Arguments ..
  111. INTEGER INFO, K, LDA, N
  112. * ..
  113. * .. Array Arguments ..
  114. INTEGER ISEED( 4 )
  115. REAL D( * )
  116. COMPLEX A( LDA, * ), WORK( * )
  117. * ..
  118. *
  119. * =====================================================================
  120. *
  121. * .. Parameters ..
  122. COMPLEX ZERO, ONE, HALF
  123. PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ),
  124. $ ONE = ( 1.0E+0, 0.0E+0 ),
  125. $ HALF = ( 0.5E+0, 0.0E+0 ) )
  126. * ..
  127. * .. Local Scalars ..
  128. INTEGER I, II, J, JJ
  129. REAL WN
  130. COMPLEX ALPHA, TAU, WA, WB
  131. * ..
  132. * .. External Subroutines ..
  133. EXTERNAL CAXPY, CGEMV, CGERC, CLACGV, CLARNV, CSCAL,
  134. $ CSYMV, XERBLA
  135. * ..
  136. * .. External Functions ..
  137. REAL SCNRM2
  138. COMPLEX CDOTC
  139. EXTERNAL SCNRM2, CDOTC
  140. * ..
  141. * .. Intrinsic Functions ..
  142. INTRINSIC ABS, MAX, REAL
  143. * ..
  144. * .. Executable Statements ..
  145. *
  146. * Test the input arguments
  147. *
  148. INFO = 0
  149. IF( N.LT.0 ) THEN
  150. INFO = -1
  151. ELSE IF( K.LT.0 .OR. K.GT.N-1 ) THEN
  152. INFO = -2
  153. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  154. INFO = -5
  155. END IF
  156. IF( INFO.LT.0 ) THEN
  157. CALL XERBLA( 'CLAGSY', -INFO )
  158. RETURN
  159. END IF
  160. *
  161. * initialize lower triangle of A to diagonal matrix
  162. *
  163. DO 20 J = 1, N
  164. DO 10 I = J + 1, N
  165. A( I, J ) = ZERO
  166. 10 CONTINUE
  167. 20 CONTINUE
  168. DO 30 I = 1, N
  169. A( I, I ) = D( I )
  170. 30 CONTINUE
  171. *
  172. * Generate lower triangle of symmetric matrix
  173. *
  174. DO 60 I = N - 1, 1, -1
  175. *
  176. * generate random reflection
  177. *
  178. CALL CLARNV( 3, ISEED, N-I+1, WORK )
  179. WN = SCNRM2( N-I+1, WORK, 1 )
  180. WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 )
  181. IF( WN.EQ.ZERO ) THEN
  182. TAU = ZERO
  183. ELSE
  184. WB = WORK( 1 ) + WA
  185. CALL CSCAL( N-I, ONE / WB, WORK( 2 ), 1 )
  186. WORK( 1 ) = ONE
  187. TAU = REAL( WB / WA )
  188. END IF
  189. *
  190. * apply random reflection to A(i:n,i:n) from the left
  191. * and the right
  192. *
  193. * compute y := tau * A * conjg(u)
  194. *
  195. CALL CLACGV( N-I+1, WORK, 1 )
  196. CALL CSYMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO,
  197. $ WORK( N+1 ), 1 )
  198. CALL CLACGV( N-I+1, WORK, 1 )
  199. *
  200. * compute v := y - 1/2 * tau * ( u, y ) * u
  201. *
  202. ALPHA = -HALF*TAU*CDOTC( N-I+1, WORK, 1, WORK( N+1 ), 1 )
  203. CALL CAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 )
  204. *
  205. * apply the transformation as a rank-2 update to A(i:n,i:n)
  206. *
  207. * CALL CSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1,
  208. * $ A( I, I ), LDA )
  209. *
  210. DO 50 JJ = I, N
  211. DO 40 II = JJ, N
  212. A( II, JJ ) = A( II, JJ ) -
  213. $ WORK( II-I+1 )*WORK( N+JJ-I+1 ) -
  214. $ WORK( N+II-I+1 )*WORK( JJ-I+1 )
  215. 40 CONTINUE
  216. 50 CONTINUE
  217. 60 CONTINUE
  218. *
  219. * Reduce number of subdiagonals to K
  220. *
  221. DO 100 I = 1, N - 1 - K
  222. *
  223. * generate reflection to annihilate A(k+i+1:n,i)
  224. *
  225. WN = SCNRM2( N-K-I+1, A( K+I, I ), 1 )
  226. WA = ( WN / ABS( A( K+I, I ) ) )*A( K+I, I )
  227. IF( WN.EQ.ZERO ) THEN
  228. TAU = ZERO
  229. ELSE
  230. WB = A( K+I, I ) + WA
  231. CALL CSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 )
  232. A( K+I, I ) = ONE
  233. TAU = REAL( WB / WA )
  234. END IF
  235. *
  236. * apply reflection to A(k+i:n,i+1:k+i-1) from the left
  237. *
  238. CALL CGEMV( 'Conjugate transpose', N-K-I+1, K-1, ONE,
  239. $ A( K+I, I+1 ), LDA, A( K+I, I ), 1, ZERO, WORK, 1 )
  240. CALL CGERC( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1,
  241. $ A( K+I, I+1 ), LDA )
  242. *
  243. * apply reflection to A(k+i:n,k+i:n) from the left and the right
  244. *
  245. * compute y := tau * A * conjg(u)
  246. *
  247. CALL CLACGV( N-K-I+1, A( K+I, I ), 1 )
  248. CALL CSYMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA,
  249. $ A( K+I, I ), 1, ZERO, WORK, 1 )
  250. CALL CLACGV( N-K-I+1, A( K+I, I ), 1 )
  251. *
  252. * compute v := y - 1/2 * tau * ( u, y ) * u
  253. *
  254. ALPHA = -HALF*TAU*CDOTC( N-K-I+1, A( K+I, I ), 1, WORK, 1 )
  255. CALL CAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 )
  256. *
  257. * apply symmetric rank-2 update to A(k+i:n,k+i:n)
  258. *
  259. * CALL CSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
  260. * $ A( K+I, K+I ), LDA )
  261. *
  262. DO 80 JJ = K + I, N
  263. DO 70 II = JJ, N
  264. A( II, JJ ) = A( II, JJ ) - A( II, I )*WORK( JJ-K-I+1 ) -
  265. $ WORK( II-K-I+1 )*A( JJ, I )
  266. 70 CONTINUE
  267. 80 CONTINUE
  268. *
  269. A( K+I, I ) = -WA
  270. DO 90 J = K + I + 1, N
  271. A( J, I ) = ZERO
  272. 90 CONTINUE
  273. 100 CONTINUE
  274. *
  275. * Store full symmetric matrix
  276. *
  277. DO 120 J = 1, N
  278. DO 110 I = J + 1, N
  279. A( J, I ) = A( I, J )
  280. 110 CONTINUE
  281. 120 CONTINUE
  282. RETURN
  283. *
  284. * End of CLAGSY
  285. *
  286. END