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.

sqrt12.f 5.8 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  1. *> \brief \b SQRT12
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. * Definition:
  9. * ===========
  10. *
  11. * REAL FUNCTION SQRT12( M, N, A, LDA, S, WORK, LWORK )
  12. *
  13. * .. Scalar Arguments ..
  14. * INTEGER LDA, LWORK, M, N
  15. * ..
  16. * .. Array Arguments ..
  17. * REAL A( LDA, * ), S( * ), WORK( LWORK )
  18. * ..
  19. *
  20. *
  21. *> \par Purpose:
  22. * =============
  23. *>
  24. *> \verbatim
  25. *>
  26. *> SQRT12 computes the singular values `svlues' of the upper trapezoid
  27. *> of A(1:M,1:N) and returns the ratio
  28. *>
  29. *> || s - svlues||/(||svlues||*eps*max(M,N))
  30. *> \endverbatim
  31. *
  32. * Arguments:
  33. * ==========
  34. *
  35. *> \param[in] M
  36. *> \verbatim
  37. *> M is INTEGER
  38. *> The number of rows of the matrix A.
  39. *> \endverbatim
  40. *>
  41. *> \param[in] N
  42. *> \verbatim
  43. *> N is INTEGER
  44. *> The number of columns of the matrix A.
  45. *> \endverbatim
  46. *>
  47. *> \param[in] A
  48. *> \verbatim
  49. *> A is REAL array, dimension (LDA,N)
  50. *> The M-by-N matrix A. Only the upper trapezoid is referenced.
  51. *> \endverbatim
  52. *>
  53. *> \param[in] LDA
  54. *> \verbatim
  55. *> LDA is INTEGER
  56. *> The leading dimension of the array A.
  57. *> \endverbatim
  58. *>
  59. *> \param[in] S
  60. *> \verbatim
  61. *> S is REAL array, dimension (min(M,N))
  62. *> The singular values of the matrix A.
  63. *> \endverbatim
  64. *>
  65. *> \param[out] WORK
  66. *> \verbatim
  67. *> WORK is REAL array, dimension (LWORK)
  68. *> \endverbatim
  69. *>
  70. *> \param[in] LWORK
  71. *> \verbatim
  72. *> LWORK is INTEGER
  73. *> The length of the array WORK. LWORK >= max(M*N + 4*min(M,N) +
  74. *> max(M,N), M*N+2*MIN( M, N )+4*N).
  75. *> \endverbatim
  76. *
  77. * Authors:
  78. * ========
  79. *
  80. *> \author Univ. of Tennessee
  81. *> \author Univ. of California Berkeley
  82. *> \author Univ. of Colorado Denver
  83. *> \author NAG Ltd.
  84. *
  85. *> \date December 2016
  86. *
  87. *> \ingroup single_lin
  88. *
  89. * =====================================================================
  90. REAL FUNCTION SQRT12( M, N, A, LDA, S, WORK, LWORK )
  91. *
  92. * -- LAPACK test routine (version 3.7.0) --
  93. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  94. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  95. * December 2016
  96. *
  97. * .. Scalar Arguments ..
  98. INTEGER LDA, LWORK, M, N
  99. * ..
  100. * .. Array Arguments ..
  101. REAL A( LDA, * ), S( * ), WORK( LWORK )
  102. * ..
  103. *
  104. * =====================================================================
  105. *
  106. * .. Parameters ..
  107. REAL ZERO, ONE
  108. PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
  109. * ..
  110. * .. Local Scalars ..
  111. INTEGER I, INFO, ISCL, J, MN
  112. REAL ANRM, BIGNUM, NRMSVL, SMLNUM
  113. * ..
  114. * .. External Functions ..
  115. REAL SASUM, SLAMCH, SLANGE, SNRM2
  116. EXTERNAL SASUM, SLAMCH, SLANGE, SNRM2
  117. * ..
  118. * .. External Subroutines ..
  119. EXTERNAL SAXPY, SBDSQR, SGEBD2, SLABAD, SLASCL, SLASET,
  120. $ XERBLA
  121. * ..
  122. * .. Intrinsic Functions ..
  123. INTRINSIC MAX, MIN, REAL
  124. * ..
  125. * .. Local Arrays ..
  126. REAL DUMMY( 1 )
  127. * ..
  128. * .. Executable Statements ..
  129. *
  130. SQRT12 = ZERO
  131. *
  132. * Test that enough workspace is supplied
  133. *
  134. IF( LWORK.LT.MAX( M*N+4*MIN( M, N )+MAX( M, N ),
  135. $ M*N+2*MIN( M, N )+4*N) ) THEN
  136. CALL XERBLA( 'SQRT12', 7 )
  137. RETURN
  138. END IF
  139. *
  140. * Quick return if possible
  141. *
  142. MN = MIN( M, N )
  143. IF( MN.LE.ZERO )
  144. $ RETURN
  145. *
  146. NRMSVL = SNRM2( MN, S, 1 )
  147. *
  148. * Copy upper triangle of A into work
  149. *
  150. CALL SLASET( 'Full', M, N, ZERO, ZERO, WORK, M )
  151. DO 20 J = 1, N
  152. DO 10 I = 1, MIN( J, M )
  153. WORK( ( J-1 )*M+I ) = A( I, J )
  154. 10 CONTINUE
  155. 20 CONTINUE
  156. *
  157. * Get machine parameters
  158. *
  159. SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
  160. BIGNUM = ONE / SMLNUM
  161. CALL SLABAD( SMLNUM, BIGNUM )
  162. *
  163. * Scale work if max entry outside range [SMLNUM,BIGNUM]
  164. *
  165. ANRM = SLANGE( 'M', M, N, WORK, M, DUMMY )
  166. ISCL = 0
  167. IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
  168. *
  169. * Scale matrix norm up to SMLNUM
  170. *
  171. CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, WORK, M, INFO )
  172. ISCL = 1
  173. ELSE IF( ANRM.GT.BIGNUM ) THEN
  174. *
  175. * Scale matrix norm down to BIGNUM
  176. *
  177. CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, WORK, M, INFO )
  178. ISCL = 1
  179. END IF
  180. *
  181. IF( ANRM.NE.ZERO ) THEN
  182. *
  183. * Compute SVD of work
  184. *
  185. CALL SGEBD2( M, N, WORK, M, WORK( M*N+1 ), WORK( M*N+MN+1 ),
  186. $ WORK( M*N+2*MN+1 ), WORK( M*N+3*MN+1 ),
  187. $ WORK( M*N+4*MN+1 ), INFO )
  188. CALL SBDSQR( 'Upper', MN, 0, 0, 0, WORK( M*N+1 ),
  189. $ WORK( M*N+MN+1 ), DUMMY, MN, DUMMY, 1, DUMMY, MN,
  190. $ WORK( M*N+2*MN+1 ), INFO )
  191. *
  192. IF( ISCL.EQ.1 ) THEN
  193. IF( ANRM.GT.BIGNUM ) THEN
  194. CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MN, 1,
  195. $ WORK( M*N+1 ), MN, INFO )
  196. END IF
  197. IF( ANRM.LT.SMLNUM ) THEN
  198. CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MN, 1,
  199. $ WORK( M*N+1 ), MN, INFO )
  200. END IF
  201. END IF
  202. *
  203. ELSE
  204. *
  205. DO 30 I = 1, MN
  206. WORK( M*N+I ) = ZERO
  207. 30 CONTINUE
  208. END IF
  209. *
  210. * Compare s and singular values of work
  211. *
  212. CALL SAXPY( MN, -ONE, S, 1, WORK( M*N+1 ), 1 )
  213. SQRT12 = SASUM( MN, WORK( M*N+1 ), 1 ) /
  214. $ ( SLAMCH( 'Epsilon' )*REAL( MAX( M, N ) ) )
  215. IF( NRMSVL.NE.ZERO )
  216. $ SQRT12 = SQRT12 / NRMSVL
  217. *
  218. RETURN
  219. *
  220. * End of SQRT12
  221. *
  222. END