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.

cla_syrcond_x.f 7.8 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298
  1. *> \brief \b CLA_SYRCOND_X computes the infinity norm condition number of op(A)*diag(x) for symmetric indefinite matrices.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CLA_SYRCOND_X + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cla_syrcond_x.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cla_syrcond_x.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cla_syrcond_x.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * REAL FUNCTION CLA_SYRCOND_X( UPLO, N, A, LDA, AF, LDAF, IPIV, X,
  22. * INFO, WORK, RWORK )
  23. *
  24. * .. Scalar Arguments ..
  25. * CHARACTER UPLO
  26. * INTEGER N, LDA, LDAF, INFO
  27. * ..
  28. * .. Array Arguments ..
  29. * INTEGER IPIV( * )
  30. * COMPLEX A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
  31. * REAL RWORK( * )
  32. * ..
  33. *
  34. *
  35. *> \par Purpose:
  36. * =============
  37. *>
  38. *> \verbatim
  39. *>
  40. *> CLA_SYRCOND_X Computes the infinity norm condition number of
  41. *> op(A) * diag(X) where X is a COMPLEX vector.
  42. *> \endverbatim
  43. *
  44. * Arguments:
  45. * ==========
  46. *
  47. *> \param[in] UPLO
  48. *> \verbatim
  49. *> UPLO is CHARACTER*1
  50. *> = 'U': Upper triangle of A is stored;
  51. *> = 'L': Lower triangle of A is stored.
  52. *> \endverbatim
  53. *>
  54. *> \param[in] N
  55. *> \verbatim
  56. *> N is INTEGER
  57. *> The number of linear equations, i.e., the order of the
  58. *> matrix A. N >= 0.
  59. *> \endverbatim
  60. *>
  61. *> \param[in] A
  62. *> \verbatim
  63. *> A is COMPLEX array, dimension (LDA,N)
  64. *> On entry, the N-by-N matrix A.
  65. *> \endverbatim
  66. *>
  67. *> \param[in] LDA
  68. *> \verbatim
  69. *> LDA is INTEGER
  70. *> The leading dimension of the array A. LDA >= max(1,N).
  71. *> \endverbatim
  72. *>
  73. *> \param[in] AF
  74. *> \verbatim
  75. *> AF is COMPLEX array, dimension (LDAF,N)
  76. *> The block diagonal matrix D and the multipliers used to
  77. *> obtain the factor U or L as computed by CSYTRF.
  78. *> \endverbatim
  79. *>
  80. *> \param[in] LDAF
  81. *> \verbatim
  82. *> LDAF is INTEGER
  83. *> The leading dimension of the array AF. LDAF >= max(1,N).
  84. *> \endverbatim
  85. *>
  86. *> \param[in] IPIV
  87. *> \verbatim
  88. *> IPIV is INTEGER array, dimension (N)
  89. *> Details of the interchanges and the block structure of D
  90. *> as determined by CSYTRF.
  91. *> \endverbatim
  92. *>
  93. *> \param[in] X
  94. *> \verbatim
  95. *> X is COMPLEX array, dimension (N)
  96. *> The vector X in the formula op(A) * diag(X).
  97. *> \endverbatim
  98. *>
  99. *> \param[out] INFO
  100. *> \verbatim
  101. *> INFO is INTEGER
  102. *> = 0: Successful exit.
  103. *> i > 0: The ith argument is invalid.
  104. *> \endverbatim
  105. *>
  106. *> \param[in] WORK
  107. *> \verbatim
  108. *> WORK is COMPLEX array, dimension (2*N).
  109. *> Workspace.
  110. *> \endverbatim
  111. *>
  112. *> \param[in] RWORK
  113. *> \verbatim
  114. *> RWORK is REAL array, dimension (N).
  115. *> Workspace.
  116. *> \endverbatim
  117. *
  118. * Authors:
  119. * ========
  120. *
  121. *> \author Univ. of Tennessee
  122. *> \author Univ. of California Berkeley
  123. *> \author Univ. of Colorado Denver
  124. *> \author NAG Ltd.
  125. *
  126. *> \date December 2016
  127. *
  128. *> \ingroup complexSYcomputational
  129. *
  130. * =====================================================================
  131. REAL FUNCTION CLA_SYRCOND_X( UPLO, N, A, LDA, AF, LDAF, IPIV, X,
  132. $ INFO, WORK, RWORK )
  133. *
  134. * -- LAPACK computational routine (version 3.7.0) --
  135. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  136. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  137. * December 2016
  138. *
  139. * .. Scalar Arguments ..
  140. CHARACTER UPLO
  141. INTEGER N, LDA, LDAF, INFO
  142. * ..
  143. * .. Array Arguments ..
  144. INTEGER IPIV( * )
  145. COMPLEX A( LDA, * ), AF( LDAF, * ), WORK( * ), X( * )
  146. REAL RWORK( * )
  147. * ..
  148. *
  149. * =====================================================================
  150. *
  151. * .. Local Scalars ..
  152. INTEGER KASE
  153. REAL AINVNM, ANORM, TMP
  154. INTEGER I, J
  155. LOGICAL UP, UPPER
  156. COMPLEX ZDUM
  157. * ..
  158. * .. Local Arrays ..
  159. INTEGER ISAVE( 3 )
  160. * ..
  161. * .. External Functions ..
  162. LOGICAL LSAME
  163. EXTERNAL LSAME
  164. * ..
  165. * .. External Subroutines ..
  166. EXTERNAL CLACN2, CSYTRS, XERBLA
  167. * ..
  168. * .. Intrinsic Functions ..
  169. INTRINSIC ABS, MAX
  170. * ..
  171. * .. Statement Functions ..
  172. REAL CABS1
  173. * ..
  174. * .. Statement Function Definitions ..
  175. CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )
  176. * ..
  177. * .. Executable Statements ..
  178. *
  179. CLA_SYRCOND_X = 0.0E+0
  180. *
  181. INFO = 0
  182. UPPER = LSAME( UPLO, 'U' )
  183. IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  184. INFO = -1
  185. ELSE IF ( N.LT.0 ) THEN
  186. INFO = -2
  187. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  188. INFO = -4
  189. ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
  190. INFO = -6
  191. END IF
  192. IF( INFO.NE.0 ) THEN
  193. CALL XERBLA( 'CLA_SYRCOND_X', -INFO )
  194. RETURN
  195. END IF
  196. UP = .FALSE.
  197. IF ( LSAME( UPLO, 'U' ) ) UP = .TRUE.
  198. *
  199. * Compute norm of op(A)*op2(C).
  200. *
  201. ANORM = 0.0
  202. IF ( UP ) THEN
  203. DO I = 1, N
  204. TMP = 0.0E+0
  205. DO J = 1, I
  206. TMP = TMP + CABS1( A( J, I ) * X( J ) )
  207. END DO
  208. DO J = I+1, N
  209. TMP = TMP + CABS1( A( I, J ) * X( J ) )
  210. END DO
  211. RWORK( I ) = TMP
  212. ANORM = MAX( ANORM, TMP )
  213. END DO
  214. ELSE
  215. DO I = 1, N
  216. TMP = 0.0E+0
  217. DO J = 1, I
  218. TMP = TMP + CABS1( A( I, J ) * X( J ) )
  219. END DO
  220. DO J = I+1, N
  221. TMP = TMP + CABS1( A( J, I ) * X( J ) )
  222. END DO
  223. RWORK( I ) = TMP
  224. ANORM = MAX( ANORM, TMP )
  225. END DO
  226. END IF
  227. *
  228. * Quick return if possible.
  229. *
  230. IF( N.EQ.0 ) THEN
  231. CLA_SYRCOND_X = 1.0E+0
  232. RETURN
  233. ELSE IF( ANORM .EQ. 0.0E+0 ) THEN
  234. RETURN
  235. END IF
  236. *
  237. * Estimate the norm of inv(op(A)).
  238. *
  239. AINVNM = 0.0E+0
  240. *
  241. KASE = 0
  242. 10 CONTINUE
  243. CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  244. IF( KASE.NE.0 ) THEN
  245. IF( KASE.EQ.2 ) THEN
  246. *
  247. * Multiply by R.
  248. *
  249. DO I = 1, N
  250. WORK( I ) = WORK( I ) * RWORK( I )
  251. END DO
  252. *
  253. IF ( UP ) THEN
  254. CALL CSYTRS( 'U', N, 1, AF, LDAF, IPIV,
  255. $ WORK, N, INFO )
  256. ELSE
  257. CALL CSYTRS( 'L', N, 1, AF, LDAF, IPIV,
  258. $ WORK, N, INFO )
  259. ENDIF
  260. *
  261. * Multiply by inv(X).
  262. *
  263. DO I = 1, N
  264. WORK( I ) = WORK( I ) / X( I )
  265. END DO
  266. ELSE
  267. *
  268. * Multiply by inv(X**T).
  269. *
  270. DO I = 1, N
  271. WORK( I ) = WORK( I ) / X( I )
  272. END DO
  273. *
  274. IF ( UP ) THEN
  275. CALL CSYTRS( 'U', N, 1, AF, LDAF, IPIV,
  276. $ WORK, N, INFO )
  277. ELSE
  278. CALL CSYTRS( 'L', N, 1, AF, LDAF, IPIV,
  279. $ WORK, N, INFO )
  280. END IF
  281. *
  282. * Multiply by R.
  283. *
  284. DO I = 1, N
  285. WORK( I ) = WORK( I ) * RWORK( I )
  286. END DO
  287. END IF
  288. GO TO 10
  289. END IF
  290. *
  291. * Compute the estimate of the reciprocal condition number.
  292. *
  293. IF( AINVNM .NE. 0.0E+0 )
  294. $ CLA_SYRCOND_X = 1.0E+0 / AINVNM
  295. *
  296. RETURN
  297. *
  298. END