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.

zqrt14.f 7.1 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256
  1. *> \brief \b ZQRT14
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. * Definition:
  9. * ===========
  10. *
  11. * DOUBLE PRECISION FUNCTION ZQRT14( TRANS, M, N, NRHS, A, LDA, X,
  12. * LDX, WORK, LWORK )
  13. *
  14. * .. Scalar Arguments ..
  15. * CHARACTER TRANS
  16. * INTEGER LDA, LDX, LWORK, M, N, NRHS
  17. * ..
  18. * .. Array Arguments ..
  19. * COMPLEX*16 A( LDA, * ), WORK( LWORK ), X( LDX, * )
  20. * ..
  21. *
  22. *
  23. *> \par Purpose:
  24. * =============
  25. *>
  26. *> \verbatim
  27. *>
  28. *> ZQRT14 checks whether X is in the row space of A or A'. It does so
  29. *> by scaling both X and A such that their norms are in the range
  30. *> [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X]
  31. *> (if TRANS = 'C') or an LQ factorization of [A',X]' (if TRANS = 'N'),
  32. *> and returning the norm of the trailing triangle, scaled by
  33. *> MAX(M,N,NRHS)*eps.
  34. *> \endverbatim
  35. *
  36. * Arguments:
  37. * ==========
  38. *
  39. *> \param[in] TRANS
  40. *> \verbatim
  41. *> TRANS is CHARACTER*1
  42. *> = 'N': No transpose, check for X in the row space of A
  43. *> = 'C': Conjugate transpose, check for X in row space of A'.
  44. *> \endverbatim
  45. *>
  46. *> \param[in] M
  47. *> \verbatim
  48. *> M is INTEGER
  49. *> The number of rows of the matrix A.
  50. *> \endverbatim
  51. *>
  52. *> \param[in] N
  53. *> \verbatim
  54. *> N is INTEGER
  55. *> The number of columns of the matrix A.
  56. *> \endverbatim
  57. *>
  58. *> \param[in] NRHS
  59. *> \verbatim
  60. *> NRHS is INTEGER
  61. *> The number of right hand sides, i.e., the number of columns
  62. *> of X.
  63. *> \endverbatim
  64. *>
  65. *> \param[in] A
  66. *> \verbatim
  67. *> A is COMPLEX*16 array, dimension (LDA,N)
  68. *> The M-by-N matrix A.
  69. *> \endverbatim
  70. *>
  71. *> \param[in] LDA
  72. *> \verbatim
  73. *> LDA is INTEGER
  74. *> The leading dimension of the array A.
  75. *> \endverbatim
  76. *>
  77. *> \param[in] X
  78. *> \verbatim
  79. *> X is COMPLEX*16 array, dimension (LDX,NRHS)
  80. *> If TRANS = 'N', the N-by-NRHS matrix X.
  81. *> IF TRANS = 'C', the M-by-NRHS matrix X.
  82. *> \endverbatim
  83. *>
  84. *> \param[in] LDX
  85. *> \verbatim
  86. *> LDX is INTEGER
  87. *> The leading dimension of the array X.
  88. *> \endverbatim
  89. *>
  90. *> \param[out] WORK
  91. *> \verbatim
  92. *> WORK is COMPLEX*16 array dimension (LWORK)
  93. *> \endverbatim
  94. *>
  95. *> \param[in] LWORK
  96. *> \verbatim
  97. *> LWORK is INTEGER
  98. *> length of workspace array required
  99. *> If TRANS = 'N', LWORK >= (M+NRHS)*(N+2);
  100. *> if TRANS = 'C', LWORK >= (N+NRHS)*(M+2).
  101. *> \endverbatim
  102. *
  103. * Authors:
  104. * ========
  105. *
  106. *> \author Univ. of Tennessee
  107. *> \author Univ. of California Berkeley
  108. *> \author Univ. of Colorado Denver
  109. *> \author NAG Ltd.
  110. *
  111. *> \ingroup complex16_lin
  112. *
  113. * =====================================================================
  114. DOUBLE PRECISION FUNCTION ZQRT14( TRANS, M, N, NRHS, A, LDA, X,
  115. $ LDX, WORK, LWORK )
  116. *
  117. * -- LAPACK test routine --
  118. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  119. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  120. *
  121. * .. Scalar Arguments ..
  122. CHARACTER TRANS
  123. INTEGER LDA, LDX, LWORK, M, N, NRHS
  124. * ..
  125. * .. Array Arguments ..
  126. COMPLEX*16 A( LDA, * ), WORK( LWORK ), X( LDX, * )
  127. * ..
  128. *
  129. * =====================================================================
  130. *
  131. * .. Parameters ..
  132. DOUBLE PRECISION ZERO, ONE
  133. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
  134. * ..
  135. * .. Local Scalars ..
  136. LOGICAL TPSD
  137. INTEGER I, INFO, J, LDWORK
  138. DOUBLE PRECISION ANRM, ERR, XNRM
  139. * ..
  140. * .. Local Arrays ..
  141. DOUBLE PRECISION RWORK( 1 )
  142. * ..
  143. * .. External Functions ..
  144. LOGICAL LSAME
  145. DOUBLE PRECISION DLAMCH, ZLANGE
  146. EXTERNAL LSAME, DLAMCH, ZLANGE
  147. * ..
  148. * .. External Subroutines ..
  149. EXTERNAL XERBLA, ZGELQ2, ZGEQR2, ZLACPY, ZLASCL
  150. * ..
  151. * .. Intrinsic Functions ..
  152. INTRINSIC ABS, DBLE, DCONJG, MAX, MIN
  153. * ..
  154. * .. Executable Statements ..
  155. *
  156. ZQRT14 = ZERO
  157. IF( LSAME( TRANS, 'N' ) ) THEN
  158. LDWORK = M + NRHS
  159. TPSD = .FALSE.
  160. IF( LWORK.LT.( M+NRHS )*( N+2 ) ) THEN
  161. CALL XERBLA( 'ZQRT14', 10 )
  162. RETURN
  163. ELSE IF( N.LE.0 .OR. NRHS.LE.0 ) THEN
  164. RETURN
  165. END IF
  166. ELSE IF( LSAME( TRANS, 'C' ) ) THEN
  167. LDWORK = M
  168. TPSD = .TRUE.
  169. IF( LWORK.LT.( N+NRHS )*( M+2 ) ) THEN
  170. CALL XERBLA( 'ZQRT14', 10 )
  171. RETURN
  172. ELSE IF( M.LE.0 .OR. NRHS.LE.0 ) THEN
  173. RETURN
  174. END IF
  175. ELSE
  176. CALL XERBLA( 'ZQRT14', 1 )
  177. RETURN
  178. END IF
  179. *
  180. * Copy and scale A
  181. *
  182. CALL ZLACPY( 'All', M, N, A, LDA, WORK, LDWORK )
  183. ANRM = ZLANGE( 'M', M, N, WORK, LDWORK, RWORK )
  184. IF( ANRM.NE.ZERO )
  185. $ CALL ZLASCL( 'G', 0, 0, ANRM, ONE, M, N, WORK, LDWORK, INFO )
  186. *
  187. * Copy X or X' into the right place and scale it
  188. *
  189. IF( TPSD ) THEN
  190. *
  191. * Copy X into columns n+1:n+nrhs of work
  192. *
  193. CALL ZLACPY( 'All', M, NRHS, X, LDX, WORK( N*LDWORK+1 ),
  194. $ LDWORK )
  195. XNRM = ZLANGE( 'M', M, NRHS, WORK( N*LDWORK+1 ), LDWORK,
  196. $ RWORK )
  197. IF( XNRM.NE.ZERO )
  198. $ CALL ZLASCL( 'G', 0, 0, XNRM, ONE, M, NRHS,
  199. $ WORK( N*LDWORK+1 ), LDWORK, INFO )
  200. *
  201. * Compute QR factorization of X
  202. *
  203. CALL ZGEQR2( M, N+NRHS, WORK, LDWORK,
  204. $ WORK( LDWORK*( N+NRHS )+1 ),
  205. $ WORK( LDWORK*( N+NRHS )+MIN( M, N+NRHS )+1 ),
  206. $ INFO )
  207. *
  208. * Compute largest entry in upper triangle of
  209. * work(n+1:m,n+1:n+nrhs)
  210. *
  211. ERR = ZERO
  212. DO 20 J = N + 1, N + NRHS
  213. DO 10 I = N + 1, MIN( M, J )
  214. ERR = MAX( ERR, ABS( WORK( I+( J-1 )*M ) ) )
  215. 10 CONTINUE
  216. 20 CONTINUE
  217. *
  218. ELSE
  219. *
  220. * Copy X' into rows m+1:m+nrhs of work
  221. *
  222. DO 40 I = 1, N
  223. DO 30 J = 1, NRHS
  224. WORK( M+J+( I-1 )*LDWORK ) = DCONJG( X( I, J ) )
  225. 30 CONTINUE
  226. 40 CONTINUE
  227. *
  228. XNRM = ZLANGE( 'M', NRHS, N, WORK( M+1 ), LDWORK, RWORK )
  229. IF( XNRM.NE.ZERO )
  230. $ CALL ZLASCL( 'G', 0, 0, XNRM, ONE, NRHS, N, WORK( M+1 ),
  231. $ LDWORK, INFO )
  232. *
  233. * Compute LQ factorization of work
  234. *
  235. CALL ZGELQ2( LDWORK, N, WORK, LDWORK, WORK( LDWORK*N+1 ),
  236. $ WORK( LDWORK*( N+1 )+1 ), INFO )
  237. *
  238. * Compute largest entry in lower triangle in
  239. * work(m+1:m+nrhs,m+1:n)
  240. *
  241. ERR = ZERO
  242. DO 60 J = M + 1, N
  243. DO 50 I = J, LDWORK
  244. ERR = MAX( ERR, ABS( WORK( I+( J-1 )*LDWORK ) ) )
  245. 50 CONTINUE
  246. 60 CONTINUE
  247. *
  248. END IF
  249. *
  250. ZQRT14 = ERR / ( DBLE( MAX( M, N, NRHS ) )*DLAMCH( 'Epsilon' ) )
  251. *
  252. RETURN
  253. *
  254. * End of ZQRT14
  255. *
  256. END