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.

slaror.f 9.1 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304
  1. *> \brief \b SLAROR
  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 SLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
  12. *
  13. * .. Scalar Arguments ..
  14. * CHARACTER INIT, SIDE
  15. * INTEGER INFO, LDA, M, N
  16. * ..
  17. * .. Array Arguments ..
  18. * INTEGER ISEED( 4 )
  19. * REAL A( LDA, * ), X( * )
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> SLAROR pre- or post-multiplies an M by N matrix A by a random
  29. *> orthogonal matrix U, overwriting A. A may optionally be initialized
  30. *> to the identity matrix before multiplying by U. U is generated using
  31. *> the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
  32. *> \endverbatim
  33. *
  34. * Arguments:
  35. * ==========
  36. *
  37. *> \param[in] SIDE
  38. *> \verbatim
  39. *> SIDE is CHARACTER*1
  40. *> Specifies whether A is multiplied on the left or right by U.
  41. *> = 'L': Multiply A on the left (premultiply) by U
  42. *> = 'R': Multiply A on the right (postmultiply) by U'
  43. *> = 'C' or 'T': Multiply A on the left by U and the right
  44. *> by U' (Here, U' means U-transpose.)
  45. *> \endverbatim
  46. *>
  47. *> \param[in] INIT
  48. *> \verbatim
  49. *> INIT is CHARACTER*1
  50. *> Specifies whether or not A should be initialized to the
  51. *> identity matrix.
  52. *> = 'I': Initialize A to (a section of) the identity matrix
  53. *> before applying U.
  54. *> = 'N': No initialization. Apply U to the input matrix A.
  55. *>
  56. *> INIT = 'I' may be used to generate square or rectangular
  57. *> orthogonal matrices:
  58. *>
  59. *> For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
  60. *> to each other, as will the columns.
  61. *>
  62. *> If M < N, SIDE = 'R' produces a dense matrix whose rows are
  63. *> orthogonal and whose columns are not, while SIDE = 'L'
  64. *> produces a matrix whose rows are orthogonal, and whose first
  65. *> M columns are orthogonal, and whose remaining columns are
  66. *> zero.
  67. *>
  68. *> If M > N, SIDE = 'L' produces a dense matrix whose columns
  69. *> are orthogonal and whose rows are not, while SIDE = 'R'
  70. *> produces a matrix whose columns are orthogonal, and whose
  71. *> first M rows are orthogonal, and whose remaining rows are
  72. *> zero.
  73. *> \endverbatim
  74. *>
  75. *> \param[in] M
  76. *> \verbatim
  77. *> M is INTEGER
  78. *> The number of rows of A.
  79. *> \endverbatim
  80. *>
  81. *> \param[in] N
  82. *> \verbatim
  83. *> N is INTEGER
  84. *> The number of columns of A.
  85. *> \endverbatim
  86. *>
  87. *> \param[in,out] A
  88. *> \verbatim
  89. *> A is REAL array, dimension (LDA, N)
  90. *> On entry, the array A.
  91. *> On exit, overwritten by U A ( if SIDE = 'L' ),
  92. *> or by A U ( if SIDE = 'R' ),
  93. *> or by U A U' ( if SIDE = 'C' or 'T').
  94. *> \endverbatim
  95. *>
  96. *> \param[in] LDA
  97. *> \verbatim
  98. *> LDA is INTEGER
  99. *> The leading dimension of the array A. LDA >= max(1,M).
  100. *> \endverbatim
  101. *>
  102. *> \param[in,out] ISEED
  103. *> \verbatim
  104. *> ISEED is INTEGER array, dimension (4)
  105. *> On entry ISEED specifies the seed of the random number
  106. *> generator. The array elements should be between 0 and 4095;
  107. *> if not they will be reduced mod 4096. Also, ISEED(4) must
  108. *> be odd. The random number generator uses a linear
  109. *> congruential sequence limited to small integers, and so
  110. *> should produce machine independent random numbers. The
  111. *> values of ISEED are changed on exit, and can be used in the
  112. *> next call to SLAROR to continue the same random number
  113. *> sequence.
  114. *> \endverbatim
  115. *>
  116. *> \param[out] X
  117. *> \verbatim
  118. *> X is REAL array, dimension (3*MAX( M, N ))
  119. *> Workspace of length
  120. *> 2*M + N if SIDE = 'L',
  121. *> 2*N + M if SIDE = 'R',
  122. *> 3*N if SIDE = 'C' or 'T'.
  123. *> \endverbatim
  124. *>
  125. *> \param[out] INFO
  126. *> \verbatim
  127. *> INFO is INTEGER
  128. *> An error flag. It is set to:
  129. *> = 0: normal return
  130. *> < 0: if INFO = -k, the k-th argument had an illegal value
  131. *> = 1: if the random numbers generated by SLARND are bad.
  132. *> \endverbatim
  133. *
  134. * Authors:
  135. * ========
  136. *
  137. *> \author Univ. of Tennessee
  138. *> \author Univ. of California Berkeley
  139. *> \author Univ. of Colorado Denver
  140. *> \author NAG Ltd.
  141. *
  142. *> \date November 2011
  143. *
  144. *> \ingroup real_matgen
  145. *
  146. * =====================================================================
  147. SUBROUTINE SLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
  148. *
  149. * -- LAPACK auxiliary routine (version 3.4.0) --
  150. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  151. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  152. * November 2011
  153. *
  154. * .. Scalar Arguments ..
  155. CHARACTER INIT, SIDE
  156. INTEGER INFO, LDA, M, N
  157. * ..
  158. * .. Array Arguments ..
  159. INTEGER ISEED( 4 )
  160. REAL A( LDA, * ), X( * )
  161. * ..
  162. *
  163. * =====================================================================
  164. *
  165. * .. Parameters ..
  166. REAL ZERO, ONE, TOOSML
  167. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0,
  168. $ TOOSML = 1.0E-20 )
  169. * ..
  170. * .. Local Scalars ..
  171. INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
  172. REAL FACTOR, XNORM, XNORMS
  173. * ..
  174. * .. External Functions ..
  175. LOGICAL LSAME
  176. REAL SLARND, SNRM2
  177. EXTERNAL LSAME, SLARND, SNRM2
  178. * ..
  179. * .. External Subroutines ..
  180. EXTERNAL SGEMV, SGER, SLASET, SSCAL, XERBLA
  181. * ..
  182. * .. Intrinsic Functions ..
  183. INTRINSIC ABS, SIGN
  184. * ..
  185. * .. Executable Statements ..
  186. *
  187. INFO = 0
  188. IF( N.EQ.0 .OR. M.EQ.0 )
  189. $ RETURN
  190. *
  191. ITYPE = 0
  192. IF( LSAME( SIDE, 'L' ) ) THEN
  193. ITYPE = 1
  194. ELSE IF( LSAME( SIDE, 'R' ) ) THEN
  195. ITYPE = 2
  196. ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN
  197. ITYPE = 3
  198. END IF
  199. *
  200. * Check for argument errors.
  201. *
  202. IF( ITYPE.EQ.0 ) THEN
  203. INFO = -1
  204. ELSE IF( M.LT.0 ) THEN
  205. INFO = -3
  206. ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
  207. INFO = -4
  208. ELSE IF( LDA.LT.M ) THEN
  209. INFO = -6
  210. END IF
  211. IF( INFO.NE.0 ) THEN
  212. CALL XERBLA( 'SLAROR', -INFO )
  213. RETURN
  214. END IF
  215. *
  216. IF( ITYPE.EQ.1 ) THEN
  217. NXFRM = M
  218. ELSE
  219. NXFRM = N
  220. END IF
  221. *
  222. * Initialize A to the identity matrix if desired
  223. *
  224. IF( LSAME( INIT, 'I' ) )
  225. $ CALL SLASET( 'Full', M, N, ZERO, ONE, A, LDA )
  226. *
  227. * If no rotation possible, multiply by random +/-1
  228. *
  229. * Compute rotation by computing Householder transformations
  230. * H(2), H(3), ..., H(nhouse)
  231. *
  232. DO 10 J = 1, NXFRM
  233. X( J ) = ZERO
  234. 10 CONTINUE
  235. *
  236. DO 30 IXFRM = 2, NXFRM
  237. KBEG = NXFRM - IXFRM + 1
  238. *
  239. * Generate independent normal( 0, 1 ) random numbers
  240. *
  241. DO 20 J = KBEG, NXFRM
  242. X( J ) = SLARND( 3, ISEED )
  243. 20 CONTINUE
  244. *
  245. * Generate a Householder transformation from the random vector X
  246. *
  247. XNORM = SNRM2( IXFRM, X( KBEG ), 1 )
  248. XNORMS = SIGN( XNORM, X( KBEG ) )
  249. X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) )
  250. FACTOR = XNORMS*( XNORMS+X( KBEG ) )
  251. IF( ABS( FACTOR ).LT.TOOSML ) THEN
  252. INFO = 1
  253. CALL XERBLA( 'SLAROR', INFO )
  254. RETURN
  255. ELSE
  256. FACTOR = ONE / FACTOR
  257. END IF
  258. X( KBEG ) = X( KBEG ) + XNORMS
  259. *
  260. * Apply Householder transformation to A
  261. *
  262. IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
  263. *
  264. * Apply H(k) from the left.
  265. *
  266. CALL SGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA,
  267. $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
  268. CALL SGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ),
  269. $ 1, A( KBEG, 1 ), LDA )
  270. *
  271. END IF
  272. *
  273. IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
  274. *
  275. * Apply H(k) from the right.
  276. *
  277. CALL SGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA,
  278. $ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
  279. CALL SGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ),
  280. $ 1, A( 1, KBEG ), LDA )
  281. *
  282. END IF
  283. 30 CONTINUE
  284. *
  285. X( 2*NXFRM ) = SIGN( ONE, SLARND( 3, ISEED ) )
  286. *
  287. * Scale the matrix A by D.
  288. *
  289. IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
  290. DO 40 IROW = 1, M
  291. CALL SSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA )
  292. 40 CONTINUE
  293. END IF
  294. *
  295. IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
  296. DO 50 JCOL = 1, N
  297. CALL SSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
  298. 50 CONTINUE
  299. END IF
  300. RETURN
  301. *
  302. * End of SLAROR
  303. *
  304. END