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.

cppequ.f 6.5 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240
  1. *> \brief \b CPPEQU
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CPPEQU + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cppequ.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cppequ.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cppequ.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFO )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER UPLO
  25. * INTEGER INFO, N
  26. * REAL AMAX, SCOND
  27. * ..
  28. * .. Array Arguments ..
  29. * REAL S( * )
  30. * COMPLEX AP( * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> CPPEQU computes row and column scalings intended to equilibrate a
  40. *> Hermitian positive definite matrix A in packed storage and reduce
  41. *> its condition number (with respect to the two-norm). S contains the
  42. *> scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix
  43. *> B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.
  44. *> This choice of S puts the condition number of B within a factor N of
  45. *> the smallest possible condition number over all possible diagonal
  46. *> scalings.
  47. *> \endverbatim
  48. *
  49. * Arguments:
  50. * ==========
  51. *
  52. *> \param[in] UPLO
  53. *> \verbatim
  54. *> UPLO is CHARACTER*1
  55. *> = 'U': Upper triangle of A is stored;
  56. *> = 'L': Lower triangle of A is stored.
  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 array, dimension (N*(N+1)/2)
  68. *> The upper or lower triangle of the Hermitian matrix A, packed
  69. *> columnwise in a linear array. The j-th column of A is stored
  70. *> in the array AP as follows:
  71. *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  72. *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  73. *> \endverbatim
  74. *>
  75. *> \param[out] S
  76. *> \verbatim
  77. *> S is REAL array, dimension (N)
  78. *> If INFO = 0, S contains the scale factors for A.
  79. *> \endverbatim
  80. *>
  81. *> \param[out] SCOND
  82. *> \verbatim
  83. *> SCOND is REAL
  84. *> If INFO = 0, S contains the ratio of the smallest S(i) to
  85. *> the largest S(i). If SCOND >= 0.1 and AMAX is neither too
  86. *> large nor too small, it is not worth scaling by S.
  87. *> \endverbatim
  88. *>
  89. *> \param[out] AMAX
  90. *> \verbatim
  91. *> AMAX is REAL
  92. *> Absolute value of largest matrix element. If AMAX is very
  93. *> close to overflow or very close to underflow, the matrix
  94. *> should be scaled.
  95. *> \endverbatim
  96. *>
  97. *> \param[out] INFO
  98. *> \verbatim
  99. *> INFO is INTEGER
  100. *> = 0: successful exit
  101. *> < 0: if INFO = -i, the i-th argument had an illegal value
  102. *> > 0: if INFO = i, the i-th diagonal element is nonpositive.
  103. *> \endverbatim
  104. *
  105. * Authors:
  106. * ========
  107. *
  108. *> \author Univ. of Tennessee
  109. *> \author Univ. of California Berkeley
  110. *> \author Univ. of Colorado Denver
  111. *> \author NAG Ltd.
  112. *
  113. *> \date November 2011
  114. *
  115. *> \ingroup complexOTHERcomputational
  116. *
  117. * =====================================================================
  118. SUBROUTINE CPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFO )
  119. *
  120. * -- LAPACK computational routine (version 3.4.0) --
  121. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  122. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  123. * November 2011
  124. *
  125. * .. Scalar Arguments ..
  126. CHARACTER UPLO
  127. INTEGER INFO, N
  128. REAL AMAX, SCOND
  129. * ..
  130. * .. Array Arguments ..
  131. REAL S( * )
  132. COMPLEX AP( * )
  133. * ..
  134. *
  135. * =====================================================================
  136. *
  137. * .. Parameters ..
  138. REAL ONE, ZERO
  139. PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  140. * ..
  141. * .. Local Scalars ..
  142. LOGICAL UPPER
  143. INTEGER I, JJ
  144. REAL SMIN
  145. * ..
  146. * .. External Functions ..
  147. LOGICAL LSAME
  148. EXTERNAL LSAME
  149. * ..
  150. * .. External Subroutines ..
  151. EXTERNAL XERBLA
  152. * ..
  153. * .. Intrinsic Functions ..
  154. INTRINSIC MAX, MIN, REAL, SQRT
  155. * ..
  156. * .. Executable Statements ..
  157. *
  158. * Test the input parameters.
  159. *
  160. INFO = 0
  161. UPPER = LSAME( UPLO, 'U' )
  162. IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  163. INFO = -1
  164. ELSE IF( N.LT.0 ) THEN
  165. INFO = -2
  166. END IF
  167. IF( INFO.NE.0 ) THEN
  168. CALL XERBLA( 'CPPEQU', -INFO )
  169. RETURN
  170. END IF
  171. *
  172. * Quick return if possible
  173. *
  174. IF( N.EQ.0 ) THEN
  175. SCOND = ONE
  176. AMAX = ZERO
  177. RETURN
  178. END IF
  179. *
  180. * Initialize SMIN and AMAX.
  181. *
  182. S( 1 ) = REAL( AP( 1 ) )
  183. SMIN = S( 1 )
  184. AMAX = S( 1 )
  185. *
  186. IF( UPPER ) THEN
  187. *
  188. * UPLO = 'U': Upper triangle of A is stored.
  189. * Find the minimum and maximum diagonal elements.
  190. *
  191. JJ = 1
  192. DO 10 I = 2, N
  193. JJ = JJ + I
  194. S( I ) = REAL( AP( JJ ) )
  195. SMIN = MIN( SMIN, S( I ) )
  196. AMAX = MAX( AMAX, S( I ) )
  197. 10 CONTINUE
  198. *
  199. ELSE
  200. *
  201. * UPLO = 'L': Lower triangle of A is stored.
  202. * Find the minimum and maximum diagonal elements.
  203. *
  204. JJ = 1
  205. DO 20 I = 2, N
  206. JJ = JJ + N - I + 2
  207. S( I ) = REAL( AP( JJ ) )
  208. SMIN = MIN( SMIN, S( I ) )
  209. AMAX = MAX( AMAX, S( I ) )
  210. 20 CONTINUE
  211. END IF
  212. *
  213. IF( SMIN.LE.ZERO ) THEN
  214. *
  215. * Find the first non-positive diagonal element and return.
  216. *
  217. DO 30 I = 1, N
  218. IF( S( I ).LE.ZERO ) THEN
  219. INFO = I
  220. RETURN
  221. END IF
  222. 30 CONTINUE
  223. ELSE
  224. *
  225. * Set the scale factors to the reciprocals
  226. * of the diagonal elements.
  227. *
  228. DO 40 I = 1, N
  229. S( I ) = ONE / SQRT( S( I ) )
  230. 40 CONTINUE
  231. *
  232. * Compute SCOND = min(S(I)) / max(S(I))
  233. *
  234. SCOND = SQRT( SMIN ) / SQRT( AMAX )
  235. END IF
  236. RETURN
  237. *
  238. * End of CPPEQU
  239. *
  240. END