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.

zspcon.f 6.2 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228
  1. *> \brief \b ZSPCON
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download ZSPCON + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zspcon.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zspcon.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zspcon.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZSPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER UPLO
  25. * INTEGER INFO, N
  26. * DOUBLE PRECISION ANORM, RCOND
  27. * ..
  28. * .. Array Arguments ..
  29. * INTEGER IPIV( * )
  30. * COMPLEX*16 AP( * ), WORK( * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> ZSPCON estimates the reciprocal of the condition number (in the
  40. *> 1-norm) of a complex symmetric packed matrix A using the
  41. *> factorization A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.
  42. *>
  43. *> An estimate is obtained for norm(inv(A)), and the reciprocal of the
  44. *> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
  45. *> \endverbatim
  46. *
  47. * Arguments:
  48. * ==========
  49. *
  50. *> \param[in] UPLO
  51. *> \verbatim
  52. *> UPLO is CHARACTER*1
  53. *> Specifies whether the details of the factorization are stored
  54. *> as an upper or lower triangular matrix.
  55. *> = 'U': Upper triangular, form is A = U*D*U**T;
  56. *> = 'L': Lower triangular, form is A = L*D*L**T.
  57. *> \endverbatim
  58. *>
  59. *> \param[in] N
  60. *> \verbatim
  61. *> N is INTEGER
  62. *> The order of the matrix A. N >= 0.
  63. *> \endverbatim
  64. *>
  65. *> \param[in] AP
  66. *> \verbatim
  67. *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
  68. *> The block diagonal matrix D and the multipliers used to
  69. *> obtain the factor U or L as computed by ZSPTRF, stored as a
  70. *> packed triangular matrix.
  71. *> \endverbatim
  72. *>
  73. *> \param[in] IPIV
  74. *> \verbatim
  75. *> IPIV is INTEGER array, dimension (N)
  76. *> Details of the interchanges and the block structure of D
  77. *> as determined by ZSPTRF.
  78. *> \endverbatim
  79. *>
  80. *> \param[in] ANORM
  81. *> \verbatim
  82. *> ANORM is DOUBLE PRECISION
  83. *> The 1-norm of the original matrix A.
  84. *> \endverbatim
  85. *>
  86. *> \param[out] RCOND
  87. *> \verbatim
  88. *> RCOND is DOUBLE PRECISION
  89. *> The reciprocal of the condition number of the matrix A,
  90. *> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
  91. *> estimate of the 1-norm of inv(A) computed in this routine.
  92. *> \endverbatim
  93. *>
  94. *> \param[out] WORK
  95. *> \verbatim
  96. *> WORK is COMPLEX*16 array, dimension (2*N)
  97. *> \endverbatim
  98. *>
  99. *> \param[out] INFO
  100. *> \verbatim
  101. *> INFO is INTEGER
  102. *> = 0: successful exit
  103. *> < 0: if INFO = -i, the i-th argument had an illegal value
  104. *> \endverbatim
  105. *
  106. * Authors:
  107. * ========
  108. *
  109. *> \author Univ. of Tennessee
  110. *> \author Univ. of California Berkeley
  111. *> \author Univ. of Colorado Denver
  112. *> \author NAG Ltd.
  113. *
  114. *> \ingroup complex16OTHERcomputational
  115. *
  116. * =====================================================================
  117. SUBROUTINE ZSPCON( UPLO, N, AP, IPIV, ANORM, RCOND, WORK, INFO )
  118. *
  119. * -- LAPACK computational routine --
  120. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  121. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  122. *
  123. * .. Scalar Arguments ..
  124. CHARACTER UPLO
  125. INTEGER INFO, N
  126. DOUBLE PRECISION ANORM, RCOND
  127. * ..
  128. * .. Array Arguments ..
  129. INTEGER IPIV( * )
  130. COMPLEX*16 AP( * ), WORK( * )
  131. * ..
  132. *
  133. * =====================================================================
  134. *
  135. * .. Parameters ..
  136. DOUBLE PRECISION ONE, ZERO
  137. PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  138. * ..
  139. * .. Local Scalars ..
  140. LOGICAL UPPER
  141. INTEGER I, IP, KASE
  142. DOUBLE PRECISION AINVNM
  143. * ..
  144. * .. Local Arrays ..
  145. INTEGER ISAVE( 3 )
  146. * ..
  147. * .. External Functions ..
  148. LOGICAL LSAME
  149. EXTERNAL LSAME
  150. * ..
  151. * .. External Subroutines ..
  152. EXTERNAL XERBLA, ZLACN2, ZSPTRS
  153. * ..
  154. * .. Executable Statements ..
  155. *
  156. * Test the input parameters.
  157. *
  158. INFO = 0
  159. UPPER = LSAME( UPLO, 'U' )
  160. IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  161. INFO = -1
  162. ELSE IF( N.LT.0 ) THEN
  163. INFO = -2
  164. ELSE IF( ANORM.LT.ZERO ) THEN
  165. INFO = -5
  166. END IF
  167. IF( INFO.NE.0 ) THEN
  168. CALL XERBLA( 'ZSPCON', -INFO )
  169. RETURN
  170. END IF
  171. *
  172. * Quick return if possible
  173. *
  174. RCOND = ZERO
  175. IF( N.EQ.0 ) THEN
  176. RCOND = ONE
  177. RETURN
  178. ELSE IF( ANORM.LE.ZERO ) THEN
  179. RETURN
  180. END IF
  181. *
  182. * Check that the diagonal matrix D is nonsingular.
  183. *
  184. IF( UPPER ) THEN
  185. *
  186. * Upper triangular storage: examine D from bottom to top
  187. *
  188. IP = N*( N+1 ) / 2
  189. DO 10 I = N, 1, -1
  190. IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
  191. $ RETURN
  192. IP = IP - I
  193. 10 CONTINUE
  194. ELSE
  195. *
  196. * Lower triangular storage: examine D from top to bottom.
  197. *
  198. IP = 1
  199. DO 20 I = 1, N
  200. IF( IPIV( I ).GT.0 .AND. AP( IP ).EQ.ZERO )
  201. $ RETURN
  202. IP = IP + N - I + 1
  203. 20 CONTINUE
  204. END IF
  205. *
  206. * Estimate the 1-norm of the inverse.
  207. *
  208. KASE = 0
  209. 30 CONTINUE
  210. CALL ZLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE )
  211. IF( KASE.NE.0 ) THEN
  212. *
  213. * Multiply by inv(L*D*L**T) or inv(U*D*U**T).
  214. *
  215. CALL ZSPTRS( UPLO, N, 1, AP, IPIV, WORK, N, INFO )
  216. GO TO 30
  217. END IF
  218. *
  219. * Compute the estimate of the reciprocal condition number.
  220. *
  221. IF( AINVNM.NE.ZERO )
  222. $ RCOND = ( ONE / AINVNM ) / ANORM
  223. *
  224. RETURN
  225. *
  226. * End of ZSPCON
  227. *
  228. END