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.

sspcon.f 6.3 kB

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