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.

clansp.f 9.0 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294
  1. *> \brief \b CLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CLANSP + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clansp.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clansp.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clansp.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK )
  22. *
  23. * .. Scalar Arguments ..
  24. * CHARACTER NORM, UPLO
  25. * INTEGER N
  26. * ..
  27. * .. Array Arguments ..
  28. * REAL WORK( * )
  29. * COMPLEX AP( * )
  30. * ..
  31. *
  32. *
  33. *> \par Purpose:
  34. * =============
  35. *>
  36. *> \verbatim
  37. *>
  38. *> CLANSP returns the value of the one norm, or the Frobenius norm, or
  39. *> the infinity norm, or the element of largest absolute value of a
  40. *> complex symmetric matrix A, supplied in packed form.
  41. *> \endverbatim
  42. *>
  43. *> \return CLANSP
  44. *> \verbatim
  45. *>
  46. *> CLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
  47. *> (
  48. *> ( norm1(A), NORM = '1', 'O' or 'o'
  49. *> (
  50. *> ( normI(A), NORM = 'I' or 'i'
  51. *> (
  52. *> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
  53. *>
  54. *> where norm1 denotes the one norm of a matrix (maximum column sum),
  55. *> normI denotes the infinity norm of a matrix (maximum row sum) and
  56. *> normF denotes the Frobenius norm of a matrix (square root of sum of
  57. *> squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.
  58. *> \endverbatim
  59. *
  60. * Arguments:
  61. * ==========
  62. *
  63. *> \param[in] NORM
  64. *> \verbatim
  65. *> NORM is CHARACTER*1
  66. *> Specifies the value to be returned in CLANSP as described
  67. *> above.
  68. *> \endverbatim
  69. *>
  70. *> \param[in] UPLO
  71. *> \verbatim
  72. *> UPLO is CHARACTER*1
  73. *> Specifies whether the upper or lower triangular part of the
  74. *> symmetric matrix A is supplied.
  75. *> = 'U': Upper triangular part of A is supplied
  76. *> = 'L': Lower triangular part of A is supplied
  77. *> \endverbatim
  78. *>
  79. *> \param[in] N
  80. *> \verbatim
  81. *> N is INTEGER
  82. *> The order of the matrix A. N >= 0. When N = 0, CLANSP is
  83. *> set to zero.
  84. *> \endverbatim
  85. *>
  86. *> \param[in] AP
  87. *> \verbatim
  88. *> AP is COMPLEX array, dimension (N*(N+1)/2)
  89. *> The upper or lower triangle of the symmetric matrix A, packed
  90. *> columnwise in a linear array. The j-th column of A is stored
  91. *> in the array AP as follows:
  92. *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  93. *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  94. *> \endverbatim
  95. *>
  96. *> \param[out] WORK
  97. *> \verbatim
  98. *> WORK is REAL array, dimension (MAX(1,LWORK)),
  99. *> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
  100. *> WORK is not referenced.
  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. *> \date December 2016
  112. *
  113. *> \ingroup complexOTHERauxiliary
  114. *
  115. * =====================================================================
  116. REAL FUNCTION CLANSP( NORM, UPLO, N, AP, WORK )
  117. *
  118. * -- LAPACK auxiliary routine (version 3.7.0) --
  119. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  120. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  121. * December 2016
  122. *
  123. IMPLICIT NONE
  124. * .. Scalar Arguments ..
  125. CHARACTER NORM, UPLO
  126. INTEGER N
  127. * ..
  128. * .. Array Arguments ..
  129. REAL WORK( * )
  130. COMPLEX AP( * )
  131. * ..
  132. *
  133. * =====================================================================
  134. *
  135. * .. Parameters ..
  136. REAL ONE, ZERO
  137. PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  138. * ..
  139. * .. Local Scalars ..
  140. INTEGER I, J, K
  141. REAL ABSA, SUM, VALUE
  142. * ..
  143. * .. Local Arrays ..
  144. REAL SSQ( 2 ), COLSSQ( 2 )
  145. * ..
  146. * .. External Functions ..
  147. LOGICAL LSAME, SISNAN
  148. EXTERNAL LSAME, SISNAN
  149. * ..
  150. * .. External Subroutines ..
  151. EXTERNAL CLASSQ, SCOMBSSQ
  152. * ..
  153. * .. Intrinsic Functions ..
  154. INTRINSIC ABS, AIMAG, REAL, SQRT
  155. * ..
  156. * .. Executable Statements ..
  157. *
  158. IF( N.EQ.0 ) THEN
  159. VALUE = ZERO
  160. ELSE IF( LSAME( NORM, 'M' ) ) THEN
  161. *
  162. * Find max(abs(A(i,j))).
  163. *
  164. VALUE = ZERO
  165. IF( LSAME( UPLO, 'U' ) ) THEN
  166. K = 1
  167. DO 20 J = 1, N
  168. DO 10 I = K, K + J - 1
  169. SUM = ABS( AP( I ) )
  170. IF( VALUE .LT. SUM .OR. SISNAN( SUM ) ) VALUE = SUM
  171. 10 CONTINUE
  172. K = K + J
  173. 20 CONTINUE
  174. ELSE
  175. K = 1
  176. DO 40 J = 1, N
  177. DO 30 I = K, K + N - J
  178. SUM = ABS( AP( I ) )
  179. IF( VALUE .LT. SUM .OR. SISNAN( SUM ) ) VALUE = SUM
  180. 30 CONTINUE
  181. K = K + N - J + 1
  182. 40 CONTINUE
  183. END IF
  184. ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
  185. $ ( NORM.EQ.'1' ) ) THEN
  186. *
  187. * Find normI(A) ( = norm1(A), since A is symmetric).
  188. *
  189. VALUE = ZERO
  190. K = 1
  191. IF( LSAME( UPLO, 'U' ) ) THEN
  192. DO 60 J = 1, N
  193. SUM = ZERO
  194. DO 50 I = 1, J - 1
  195. ABSA = ABS( AP( K ) )
  196. SUM = SUM + ABSA
  197. WORK( I ) = WORK( I ) + ABSA
  198. K = K + 1
  199. 50 CONTINUE
  200. WORK( J ) = SUM + ABS( AP( K ) )
  201. K = K + 1
  202. 60 CONTINUE
  203. DO 70 I = 1, N
  204. SUM = WORK( I )
  205. IF( VALUE .LT. SUM .OR. SISNAN( SUM ) ) VALUE = SUM
  206. 70 CONTINUE
  207. ELSE
  208. DO 80 I = 1, N
  209. WORK( I ) = ZERO
  210. 80 CONTINUE
  211. DO 100 J = 1, N
  212. SUM = WORK( J ) + ABS( AP( K ) )
  213. K = K + 1
  214. DO 90 I = J + 1, N
  215. ABSA = ABS( AP( K ) )
  216. SUM = SUM + ABSA
  217. WORK( I ) = WORK( I ) + ABSA
  218. K = K + 1
  219. 90 CONTINUE
  220. IF( VALUE .LT. SUM .OR. SISNAN( SUM ) ) VALUE = SUM
  221. 100 CONTINUE
  222. END IF
  223. ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
  224. *
  225. * Find normF(A).
  226. * SSQ(1) is scale
  227. * SSQ(2) is sum-of-squares
  228. * For better accuracy, sum each column separately.
  229. *
  230. SSQ( 1 ) = ZERO
  231. SSQ( 2 ) = ONE
  232. *
  233. * Sum off-diagonals
  234. *
  235. K = 2
  236. IF( LSAME( UPLO, 'U' ) ) THEN
  237. DO 110 J = 2, N
  238. COLSSQ( 1 ) = ZERO
  239. COLSSQ( 2 ) = ONE
  240. CALL CLASSQ( J-1, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
  241. CALL SCOMBSSQ( SSQ, COLSSQ )
  242. K = K + J
  243. 110 CONTINUE
  244. ELSE
  245. DO 120 J = 1, N - 1
  246. COLSSQ( 1 ) = ZERO
  247. COLSSQ( 2 ) = ONE
  248. CALL CLASSQ( N-J, AP( K ), 1, COLSSQ( 1 ), COLSSQ( 2 ) )
  249. CALL SCOMBSSQ( SSQ, COLSSQ )
  250. K = K + N - J + 1
  251. 120 CONTINUE
  252. END IF
  253. SSQ( 2 ) = 2*SSQ( 2 )
  254. *
  255. * Sum diagonal
  256. *
  257. K = 1
  258. COLSSQ( 1 ) = ZERO
  259. COLSSQ( 2 ) = ONE
  260. DO 130 I = 1, N
  261. IF( REAL( AP( K ) ).NE.ZERO ) THEN
  262. ABSA = ABS( REAL( AP( K ) ) )
  263. IF( COLSSQ( 1 ).LT.ABSA ) THEN
  264. COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
  265. COLSSQ( 1 ) = ABSA
  266. ELSE
  267. COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
  268. END IF
  269. END IF
  270. IF( AIMAG( AP( K ) ).NE.ZERO ) THEN
  271. ABSA = ABS( AIMAG( AP( K ) ) )
  272. IF( COLSSQ( 1 ).LT.ABSA ) THEN
  273. COLSSQ( 2 ) = ONE + COLSSQ(2)*( COLSSQ(1) / ABSA )**2
  274. COLSSQ( 1 ) = ABSA
  275. ELSE
  276. COLSSQ( 2 ) = COLSSQ( 2 ) + ( ABSA / COLSSQ( 1 ) )**2
  277. END IF
  278. END IF
  279. IF( LSAME( UPLO, 'U' ) ) THEN
  280. K = K + I + 1
  281. ELSE
  282. K = K + N - I + 1
  283. END IF
  284. 130 CONTINUE
  285. CALL SCOMBSSQ( SSQ, COLSSQ )
  286. VALUE = SSQ( 1 )*SQRT( SSQ( 2 ) )
  287. END IF
  288. *
  289. CLANSP = VALUE
  290. RETURN
  291. *
  292. * End of CLANSP
  293. *
  294. END