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.

clacon.f 7.3 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281
  1. *> \brief \b CLACON estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vector products.
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CLACON + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clacon.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clacon.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clacon.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CLACON( N, V, X, EST, KASE )
  22. *
  23. * .. Scalar Arguments ..
  24. * INTEGER KASE, N
  25. * REAL EST
  26. * ..
  27. * .. Array Arguments ..
  28. * COMPLEX V( N ), X( N )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> CLACON estimates the 1-norm of a square, complex matrix A.
  38. *> Reverse communication is used for evaluating matrix-vector products.
  39. *> \endverbatim
  40. *
  41. * Arguments:
  42. * ==========
  43. *
  44. *> \param[in] N
  45. *> \verbatim
  46. *> N is INTEGER
  47. *> The order of the matrix. N >= 1.
  48. *> \endverbatim
  49. *>
  50. *> \param[out] V
  51. *> \verbatim
  52. *> V is COMPLEX array, dimension (N)
  53. *> On the final return, V = A*W, where EST = norm(V)/norm(W)
  54. *> (W is not returned).
  55. *> \endverbatim
  56. *>
  57. *> \param[in,out] X
  58. *> \verbatim
  59. *> X is COMPLEX array, dimension (N)
  60. *> On an intermediate return, X should be overwritten by
  61. *> A * X, if KASE=1,
  62. *> A**H * X, if KASE=2,
  63. *> where A**H is the conjugate transpose of A, and CLACON must be
  64. *> re-called with all the other parameters unchanged.
  65. *> \endverbatim
  66. *>
  67. *> \param[in,out] EST
  68. *> \verbatim
  69. *> EST is REAL
  70. *> On entry with KASE = 1 or 2 and JUMP = 3, EST should be
  71. *> unchanged from the previous call to CLACON.
  72. *> On exit, EST is an estimate (a lower bound) for norm(A).
  73. *> \endverbatim
  74. *>
  75. *> \param[in,out] KASE
  76. *> \verbatim
  77. *> KASE is INTEGER
  78. *> On the initial call to CLACON, KASE should be 0.
  79. *> On an intermediate return, KASE will be 1 or 2, indicating
  80. *> whether X should be overwritten by A * X or A**H * X.
  81. *> On the final return from CLACON, KASE will again be 0.
  82. *> \endverbatim
  83. *
  84. * Authors:
  85. * ========
  86. *
  87. *> \author Univ. of Tennessee
  88. *> \author Univ. of California Berkeley
  89. *> \author Univ. of Colorado Denver
  90. *> \author NAG Ltd.
  91. *
  92. *> \date December 2016
  93. *
  94. *> \ingroup complexOTHERauxiliary
  95. *
  96. *> \par Further Details:
  97. * =====================
  98. *>
  99. *> Originally named CONEST, dated March 16, 1988. \n
  100. *> Last modified: April, 1999
  101. *
  102. *> \par Contributors:
  103. * ==================
  104. *>
  105. *> Nick Higham, University of Manchester
  106. *
  107. *> \par References:
  108. * ================
  109. *>
  110. *> N.J. Higham, "FORTRAN codes for estimating the one-norm of
  111. *> a real or complex matrix, with applications to condition estimation",
  112. *> ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
  113. *>
  114. * =====================================================================
  115. SUBROUTINE CLACON( N, V, X, EST, KASE )
  116. *
  117. * -- LAPACK auxiliary routine (version 3.7.0) --
  118. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  119. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  120. * December 2016
  121. *
  122. * .. Scalar Arguments ..
  123. INTEGER KASE, N
  124. REAL EST
  125. * ..
  126. * .. Array Arguments ..
  127. COMPLEX V( N ), X( N )
  128. * ..
  129. *
  130. * =====================================================================
  131. *
  132. * .. Parameters ..
  133. INTEGER ITMAX
  134. PARAMETER ( ITMAX = 5 )
  135. REAL ONE, TWO
  136. PARAMETER ( ONE = 1.0E0, TWO = 2.0E0 )
  137. COMPLEX CZERO, CONE
  138. PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ),
  139. $ CONE = ( 1.0E0, 0.0E0 ) )
  140. * ..
  141. * .. Local Scalars ..
  142. INTEGER I, ITER, J, JLAST, JUMP
  143. REAL ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
  144. * ..
  145. * .. External Functions ..
  146. INTEGER ICMAX1
  147. REAL SCSUM1, SLAMCH
  148. EXTERNAL ICMAX1, SCSUM1, SLAMCH
  149. * ..
  150. * .. External Subroutines ..
  151. EXTERNAL CCOPY
  152. * ..
  153. * .. Intrinsic Functions ..
  154. INTRINSIC ABS, AIMAG, CMPLX, REAL
  155. * ..
  156. * .. Save statement ..
  157. SAVE
  158. * ..
  159. * .. Executable Statements ..
  160. *
  161. SAFMIN = SLAMCH( 'Safe minimum' )
  162. IF( KASE.EQ.0 ) THEN
  163. DO 10 I = 1, N
  164. X( I ) = CMPLX( ONE / REAL( N ) )
  165. 10 CONTINUE
  166. KASE = 1
  167. JUMP = 1
  168. RETURN
  169. END IF
  170. *
  171. GO TO ( 20, 40, 70, 90, 120 )JUMP
  172. *
  173. * ................ ENTRY (JUMP = 1)
  174. * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
  175. *
  176. 20 CONTINUE
  177. IF( N.EQ.1 ) THEN
  178. V( 1 ) = X( 1 )
  179. EST = ABS( V( 1 ) )
  180. * ... QUIT
  181. GO TO 130
  182. END IF
  183. EST = SCSUM1( N, X, 1 )
  184. *
  185. DO 30 I = 1, N
  186. ABSXI = ABS( X( I ) )
  187. IF( ABSXI.GT.SAFMIN ) THEN
  188. X( I ) = CMPLX( REAL( X( I ) ) / ABSXI,
  189. $ AIMAG( X( I ) ) / ABSXI )
  190. ELSE
  191. X( I ) = CONE
  192. END IF
  193. 30 CONTINUE
  194. KASE = 2
  195. JUMP = 2
  196. RETURN
  197. *
  198. * ................ ENTRY (JUMP = 2)
  199. * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
  200. *
  201. 40 CONTINUE
  202. J = ICMAX1( N, X, 1 )
  203. ITER = 2
  204. *
  205. * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
  206. *
  207. 50 CONTINUE
  208. DO 60 I = 1, N
  209. X( I ) = CZERO
  210. 60 CONTINUE
  211. X( J ) = CONE
  212. KASE = 1
  213. JUMP = 3
  214. RETURN
  215. *
  216. * ................ ENTRY (JUMP = 3)
  217. * X HAS BEEN OVERWRITTEN BY A*X.
  218. *
  219. 70 CONTINUE
  220. CALL CCOPY( N, X, 1, V, 1 )
  221. ESTOLD = EST
  222. EST = SCSUM1( N, V, 1 )
  223. *
  224. * TEST FOR CYCLING.
  225. IF( EST.LE.ESTOLD )
  226. $ GO TO 100
  227. *
  228. DO 80 I = 1, N
  229. ABSXI = ABS( X( I ) )
  230. IF( ABSXI.GT.SAFMIN ) THEN
  231. X( I ) = CMPLX( REAL( X( I ) ) / ABSXI,
  232. $ AIMAG( X( I ) ) / ABSXI )
  233. ELSE
  234. X( I ) = CONE
  235. END IF
  236. 80 CONTINUE
  237. KASE = 2
  238. JUMP = 4
  239. RETURN
  240. *
  241. * ................ ENTRY (JUMP = 4)
  242. * X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
  243. *
  244. 90 CONTINUE
  245. JLAST = J
  246. J = ICMAX1( N, X, 1 )
  247. IF( ( ABS( X( JLAST ) ).NE.ABS( X( J ) ) ) .AND.
  248. $ ( ITER.LT.ITMAX ) ) THEN
  249. ITER = ITER + 1
  250. GO TO 50
  251. END IF
  252. *
  253. * ITERATION COMPLETE. FINAL STAGE.
  254. *
  255. 100 CONTINUE
  256. ALTSGN = ONE
  257. DO 110 I = 1, N
  258. X( I ) = CMPLX( ALTSGN*( ONE+REAL( I-1 ) / REAL( N-1 ) ) )
  259. ALTSGN = -ALTSGN
  260. 110 CONTINUE
  261. KASE = 1
  262. JUMP = 5
  263. RETURN
  264. *
  265. * ................ ENTRY (JUMP = 5)
  266. * X HAS BEEN OVERWRITTEN BY A*X.
  267. *
  268. 120 CONTINUE
  269. TEMP = TWO*( SCSUM1( N, X, 1 ) / REAL( 3*N ) )
  270. IF( TEMP.GT.EST ) THEN
  271. CALL CCOPY( N, X, 1, V, 1 )
  272. EST = TEMP
  273. END IF
  274. *
  275. 130 CONTINUE
  276. KASE = 0
  277. RETURN
  278. *
  279. * End of CLACON
  280. *
  281. END