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.

zlacon.f 7.3 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278
  1. *> \brief \b ZLACON 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 ZLACON + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlacon.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlacon.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlacon.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE ZLACON( N, V, X, EST, KASE )
  22. *
  23. * .. Scalar Arguments ..
  24. * INTEGER KASE, N
  25. * DOUBLE PRECISION EST
  26. * ..
  27. * .. Array Arguments ..
  28. * COMPLEX*16 V( N ), X( N )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> ZLACON 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*16 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*16 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 ZLACON must be
  64. *> re-called with all the other parameters unchanged.
  65. *> \endverbatim
  66. *>
  67. *> \param[in,out] EST
  68. *> \verbatim
  69. *> EST is DOUBLE PRECISION
  70. *> On entry with KASE = 1 or 2 and JUMP = 3, EST should be
  71. *> unchanged from the previous call to ZLACON.
  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 ZLACON, 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 ZLACON, 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. *> \ingroup complex16OTHERauxiliary
  93. *
  94. *> \par Further Details:
  95. * =====================
  96. *>
  97. *> Originally named CONEST, dated March 16, 1988. \n
  98. *> Last modified: April, 1999
  99. *
  100. *> \par Contributors:
  101. * ==================
  102. *>
  103. *> Nick Higham, University of Manchester
  104. *
  105. *> \par References:
  106. * ================
  107. *>
  108. *> N.J. Higham, "FORTRAN codes for estimating the one-norm of
  109. *> a real or complex matrix, with applications to condition estimation",
  110. *> ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
  111. *>
  112. * =====================================================================
  113. SUBROUTINE ZLACON( N, V, X, EST, KASE )
  114. *
  115. * -- LAPACK auxiliary routine --
  116. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  117. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  118. *
  119. * .. Scalar Arguments ..
  120. INTEGER KASE, N
  121. DOUBLE PRECISION EST
  122. * ..
  123. * .. Array Arguments ..
  124. COMPLEX*16 V( N ), X( N )
  125. * ..
  126. *
  127. * =====================================================================
  128. *
  129. * .. Parameters ..
  130. INTEGER ITMAX
  131. PARAMETER ( ITMAX = 5 )
  132. DOUBLE PRECISION ONE, TWO
  133. PARAMETER ( ONE = 1.0D0, TWO = 2.0D0 )
  134. COMPLEX*16 CZERO, CONE
  135. PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
  136. $ CONE = ( 1.0D0, 0.0D0 ) )
  137. * ..
  138. * .. Local Scalars ..
  139. INTEGER I, ITER, J, JLAST, JUMP
  140. DOUBLE PRECISION ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
  141. * ..
  142. * .. External Functions ..
  143. INTEGER IZMAX1
  144. DOUBLE PRECISION DLAMCH, DZSUM1
  145. EXTERNAL IZMAX1, DLAMCH, DZSUM1
  146. * ..
  147. * .. External Subroutines ..
  148. EXTERNAL ZCOPY
  149. * ..
  150. * .. Intrinsic Functions ..
  151. INTRINSIC ABS, DBLE, DCMPLX, DIMAG
  152. * ..
  153. * .. Save statement ..
  154. SAVE
  155. * ..
  156. * .. Executable Statements ..
  157. *
  158. SAFMIN = DLAMCH( 'Safe minimum' )
  159. IF( KASE.EQ.0 ) THEN
  160. DO 10 I = 1, N
  161. X( I ) = DCMPLX( ONE / DBLE( N ) )
  162. 10 CONTINUE
  163. KASE = 1
  164. JUMP = 1
  165. RETURN
  166. END IF
  167. *
  168. GO TO ( 20, 40, 70, 90, 120 )JUMP
  169. *
  170. * ................ ENTRY (JUMP = 1)
  171. * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
  172. *
  173. 20 CONTINUE
  174. IF( N.EQ.1 ) THEN
  175. V( 1 ) = X( 1 )
  176. EST = ABS( V( 1 ) )
  177. * ... QUIT
  178. GO TO 130
  179. END IF
  180. EST = DZSUM1( N, X, 1 )
  181. *
  182. DO 30 I = 1, N
  183. ABSXI = ABS( X( I ) )
  184. IF( ABSXI.GT.SAFMIN ) THEN
  185. X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
  186. $ DIMAG( X( I ) ) / ABSXI )
  187. ELSE
  188. X( I ) = CONE
  189. END IF
  190. 30 CONTINUE
  191. KASE = 2
  192. JUMP = 2
  193. RETURN
  194. *
  195. * ................ ENTRY (JUMP = 2)
  196. * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
  197. *
  198. 40 CONTINUE
  199. J = IZMAX1( N, X, 1 )
  200. ITER = 2
  201. *
  202. * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
  203. *
  204. 50 CONTINUE
  205. DO 60 I = 1, N
  206. X( I ) = CZERO
  207. 60 CONTINUE
  208. X( J ) = CONE
  209. KASE = 1
  210. JUMP = 3
  211. RETURN
  212. *
  213. * ................ ENTRY (JUMP = 3)
  214. * X HAS BEEN OVERWRITTEN BY A*X.
  215. *
  216. 70 CONTINUE
  217. CALL ZCOPY( N, X, 1, V, 1 )
  218. ESTOLD = EST
  219. EST = DZSUM1( N, V, 1 )
  220. *
  221. * TEST FOR CYCLING.
  222. IF( EST.LE.ESTOLD )
  223. $ GO TO 100
  224. *
  225. DO 80 I = 1, N
  226. ABSXI = ABS( X( I ) )
  227. IF( ABSXI.GT.SAFMIN ) THEN
  228. X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
  229. $ DIMAG( X( I ) ) / ABSXI )
  230. ELSE
  231. X( I ) = CONE
  232. END IF
  233. 80 CONTINUE
  234. KASE = 2
  235. JUMP = 4
  236. RETURN
  237. *
  238. * ................ ENTRY (JUMP = 4)
  239. * X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
  240. *
  241. 90 CONTINUE
  242. JLAST = J
  243. J = IZMAX1( N, X, 1 )
  244. IF( ( ABS( X( JLAST ) ).NE.ABS( X( J ) ) ) .AND.
  245. $ ( ITER.LT.ITMAX ) ) THEN
  246. ITER = ITER + 1
  247. GO TO 50
  248. END IF
  249. *
  250. * ITERATION COMPLETE. FINAL STAGE.
  251. *
  252. 100 CONTINUE
  253. ALTSGN = ONE
  254. DO 110 I = 1, N
  255. X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
  256. ALTSGN = -ALTSGN
  257. 110 CONTINUE
  258. KASE = 1
  259. JUMP = 5
  260. RETURN
  261. *
  262. * ................ ENTRY (JUMP = 5)
  263. * X HAS BEEN OVERWRITTEN BY A*X.
  264. *
  265. 120 CONTINUE
  266. TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
  267. IF( TEMP.GT.EST ) THEN
  268. CALL ZCOPY( N, X, 1, V, 1 )
  269. EST = TEMP
  270. END IF
  271. *
  272. 130 CONTINUE
  273. KASE = 0
  274. RETURN
  275. *
  276. * End of ZLACON
  277. *
  278. END