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.

sqrt15.f 8.3 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315
  1. *> \brief \b SQRT15
  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 SQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S,
  12. * RANK, NORMA, NORMB, ISEED, WORK, LWORK )
  13. *
  14. * .. Scalar Arguments ..
  15. * INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
  16. * REAL NORMA, NORMB
  17. * ..
  18. * .. Array Arguments ..
  19. * INTEGER ISEED( 4 )
  20. * REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
  21. * ..
  22. *
  23. *
  24. *> \par Purpose:
  25. * =============
  26. *>
  27. *> \verbatim
  28. *>
  29. *> SQRT15 generates a matrix with full or deficient rank and of various
  30. *> norms.
  31. *> \endverbatim
  32. *
  33. * Arguments:
  34. * ==========
  35. *
  36. *> \param[in] SCALE
  37. *> \verbatim
  38. *> SCALE is INTEGER
  39. *> SCALE = 1: normally scaled matrix
  40. *> SCALE = 2: matrix scaled up
  41. *> SCALE = 3: matrix scaled down
  42. *> \endverbatim
  43. *>
  44. *> \param[in] RKSEL
  45. *> \verbatim
  46. *> RKSEL is INTEGER
  47. *> RKSEL = 1: full rank matrix
  48. *> RKSEL = 2: rank-deficient matrix
  49. *> \endverbatim
  50. *>
  51. *> \param[in] M
  52. *> \verbatim
  53. *> M is INTEGER
  54. *> The number of rows of the matrix A.
  55. *> \endverbatim
  56. *>
  57. *> \param[in] N
  58. *> \verbatim
  59. *> N is INTEGER
  60. *> The number of columns of A.
  61. *> \endverbatim
  62. *>
  63. *> \param[in] NRHS
  64. *> \verbatim
  65. *> NRHS is INTEGER
  66. *> The number of columns of B.
  67. *> \endverbatim
  68. *>
  69. *> \param[out] A
  70. *> \verbatim
  71. *> A is REAL array, dimension (LDA,N)
  72. *> The M-by-N matrix A.
  73. *> \endverbatim
  74. *>
  75. *> \param[in] LDA
  76. *> \verbatim
  77. *> LDA is INTEGER
  78. *> The leading dimension of the array A.
  79. *> \endverbatim
  80. *>
  81. *> \param[out] B
  82. *> \verbatim
  83. *> B is REAL array, dimension (LDB, NRHS)
  84. *> A matrix that is in the range space of matrix A.
  85. *> \endverbatim
  86. *>
  87. *> \param[in] LDB
  88. *> \verbatim
  89. *> LDB is INTEGER
  90. *> The leading dimension of the array B.
  91. *> \endverbatim
  92. *>
  93. *> \param[out] S
  94. *> \verbatim
  95. *> S is REAL array, dimension MIN(M,N)
  96. *> Singular values of A.
  97. *> \endverbatim
  98. *>
  99. *> \param[out] RANK
  100. *> \verbatim
  101. *> RANK is INTEGER
  102. *> number of nonzero singular values of A.
  103. *> \endverbatim
  104. *>
  105. *> \param[out] NORMA
  106. *> \verbatim
  107. *> NORMA is REAL
  108. *> one-norm of A.
  109. *> \endverbatim
  110. *>
  111. *> \param[out] NORMB
  112. *> \verbatim
  113. *> NORMB is REAL
  114. *> one-norm of B.
  115. *> \endverbatim
  116. *>
  117. *> \param[in,out] ISEED
  118. *> \verbatim
  119. *> ISEED is integer array, dimension (4)
  120. *> seed for random number generator.
  121. *> \endverbatim
  122. *>
  123. *> \param[out] WORK
  124. *> \verbatim
  125. *> WORK is REAL array, dimension (LWORK)
  126. *> \endverbatim
  127. *>
  128. *> \param[in] LWORK
  129. *> \verbatim
  130. *> LWORK is INTEGER
  131. *> length of work space required.
  132. *> LWORK >= MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M)
  133. *> \endverbatim
  134. *
  135. * Authors:
  136. * ========
  137. *
  138. *> \author Univ. of Tennessee
  139. *> \author Univ. of California Berkeley
  140. *> \author Univ. of Colorado Denver
  141. *> \author NAG Ltd.
  142. *
  143. *> \date December 2016
  144. *
  145. *> \ingroup single_lin
  146. *
  147. * =====================================================================
  148. SUBROUTINE SQRT15( SCALE, RKSEL, M, N, NRHS, A, LDA, B, LDB, S,
  149. $ RANK, NORMA, NORMB, ISEED, WORK, LWORK )
  150. *
  151. * -- LAPACK test routine (version 3.7.0) --
  152. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  153. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  154. * December 2016
  155. *
  156. * .. Scalar Arguments ..
  157. INTEGER LDA, LDB, LWORK, M, N, NRHS, RANK, RKSEL, SCALE
  158. REAL NORMA, NORMB
  159. * ..
  160. * .. Array Arguments ..
  161. INTEGER ISEED( 4 )
  162. REAL A( LDA, * ), B( LDB, * ), S( * ), WORK( LWORK )
  163. * ..
  164. *
  165. * =====================================================================
  166. *
  167. * .. Parameters ..
  168. REAL ZERO, ONE, TWO, SVMIN
  169. PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
  170. $ SVMIN = 0.1E0 )
  171. * ..
  172. * .. Local Scalars ..
  173. INTEGER INFO, J, MN
  174. REAL BIGNUM, EPS, SMLNUM, TEMP
  175. * ..
  176. * .. Local Arrays ..
  177. REAL DUMMY( 1 )
  178. * ..
  179. * .. External Functions ..
  180. REAL SASUM, SLAMCH, SLANGE, SLARND, SNRM2
  181. EXTERNAL SASUM, SLAMCH, SLANGE, SLARND, SNRM2
  182. * ..
  183. * .. External Subroutines ..
  184. EXTERNAL SGEMM, SLAORD, SLARF, SLARNV, SLAROR, SLASCL,
  185. $ SLASET, SSCAL, XERBLA
  186. * ..
  187. * .. Intrinsic Functions ..
  188. INTRINSIC ABS, MAX, MIN
  189. * ..
  190. * .. Executable Statements ..
  191. *
  192. MN = MIN( M, N )
  193. IF( LWORK.LT.MAX( M+MN, MN*NRHS, 2*N+M ) ) THEN
  194. CALL XERBLA( 'SQRT15', 16 )
  195. RETURN
  196. END IF
  197. *
  198. SMLNUM = SLAMCH( 'Safe minimum' )
  199. BIGNUM = ONE / SMLNUM
  200. EPS = SLAMCH( 'Epsilon' )
  201. SMLNUM = ( SMLNUM / EPS ) / EPS
  202. BIGNUM = ONE / SMLNUM
  203. *
  204. * Determine rank and (unscaled) singular values
  205. *
  206. IF( RKSEL.EQ.1 ) THEN
  207. RANK = MN
  208. ELSE IF( RKSEL.EQ.2 ) THEN
  209. RANK = ( 3*MN ) / 4
  210. DO 10 J = RANK + 1, MN
  211. S( J ) = ZERO
  212. 10 CONTINUE
  213. ELSE
  214. CALL XERBLA( 'SQRT15', 2 )
  215. END IF
  216. *
  217. IF( RANK.GT.0 ) THEN
  218. *
  219. * Nontrivial case
  220. *
  221. S( 1 ) = ONE
  222. DO 30 J = 2, RANK
  223. 20 CONTINUE
  224. TEMP = SLARND( 1, ISEED )
  225. IF( TEMP.GT.SVMIN ) THEN
  226. S( J ) = ABS( TEMP )
  227. ELSE
  228. GO TO 20
  229. END IF
  230. 30 CONTINUE
  231. CALL SLAORD( 'Decreasing', RANK, S, 1 )
  232. *
  233. * Generate 'rank' columns of a random orthogonal matrix in A
  234. *
  235. CALL SLARNV( 2, ISEED, M, WORK )
  236. CALL SSCAL( M, ONE / SNRM2( M, WORK, 1 ), WORK, 1 )
  237. CALL SLASET( 'Full', M, RANK, ZERO, ONE, A, LDA )
  238. CALL SLARF( 'Left', M, RANK, WORK, 1, TWO, A, LDA,
  239. $ WORK( M+1 ) )
  240. *
  241. * workspace used: m+mn
  242. *
  243. * Generate consistent rhs in the range space of A
  244. *
  245. CALL SLARNV( 2, ISEED, RANK*NRHS, WORK )
  246. CALL SGEMM( 'No transpose', 'No transpose', M, NRHS, RANK, ONE,
  247. $ A, LDA, WORK, RANK, ZERO, B, LDB )
  248. *
  249. * work space used: <= mn *nrhs
  250. *
  251. * generate (unscaled) matrix A
  252. *
  253. DO 40 J = 1, RANK
  254. CALL SSCAL( M, S( J ), A( 1, J ), 1 )
  255. 40 CONTINUE
  256. IF( RANK.LT.N )
  257. $ CALL SLASET( 'Full', M, N-RANK, ZERO, ZERO, A( 1, RANK+1 ),
  258. $ LDA )
  259. CALL SLAROR( 'Right', 'No initialization', M, N, A, LDA, ISEED,
  260. $ WORK, INFO )
  261. *
  262. ELSE
  263. *
  264. * work space used 2*n+m
  265. *
  266. * Generate null matrix and rhs
  267. *
  268. DO 50 J = 1, MN
  269. S( J ) = ZERO
  270. 50 CONTINUE
  271. CALL SLASET( 'Full', M, N, ZERO, ZERO, A, LDA )
  272. CALL SLASET( 'Full', M, NRHS, ZERO, ZERO, B, LDB )
  273. *
  274. END IF
  275. *
  276. * Scale the matrix
  277. *
  278. IF( SCALE.NE.1 ) THEN
  279. NORMA = SLANGE( 'Max', M, N, A, LDA, DUMMY )
  280. IF( NORMA.NE.ZERO ) THEN
  281. IF( SCALE.EQ.2 ) THEN
  282. *
  283. * matrix scaled up
  284. *
  285. CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, M, N, A,
  286. $ LDA, INFO )
  287. CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, MN, 1, S,
  288. $ MN, INFO )
  289. CALL SLASCL( 'General', 0, 0, NORMA, BIGNUM, M, NRHS, B,
  290. $ LDB, INFO )
  291. ELSE IF( SCALE.EQ.3 ) THEN
  292. *
  293. * matrix scaled down
  294. *
  295. CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, M, N, A,
  296. $ LDA, INFO )
  297. CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, MN, 1, S,
  298. $ MN, INFO )
  299. CALL SLASCL( 'General', 0, 0, NORMA, SMLNUM, M, NRHS, B,
  300. $ LDB, INFO )
  301. ELSE
  302. CALL XERBLA( 'SQRT15', 1 )
  303. RETURN
  304. END IF
  305. END IF
  306. END IF
  307. *
  308. NORMA = SASUM( MN, S, 1 )
  309. NORMB = SLANGE( 'One-norm', M, NRHS, B, LDB, DUMMY )
  310. *
  311. RETURN
  312. *
  313. * End of SQRT15
  314. *
  315. END