1: SUBROUTINE ZLACON( N, V, X, EST, KASE )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: INTEGER KASE, N
10: DOUBLE PRECISION EST
11: * ..
12: * .. Array Arguments ..
13: COMPLEX*16 V( N ), X( N )
14: * ..
15: *
16: * Purpose
17: * =======
18: *
19: * ZLACON estimates the 1-norm of a square, complex matrix A.
20: * Reverse communication is used for evaluating matrix-vector products.
21: *
22: * Arguments
23: * =========
24: *
25: * N (input) INTEGER
26: * The order of the matrix. N >= 1.
27: *
28: * V (workspace) COMPLEX*16 array, dimension (N)
29: * On the final return, V = A*W, where EST = norm(V)/norm(W)
30: * (W is not returned).
31: *
32: * X (input/output) COMPLEX*16 array, dimension (N)
33: * On an intermediate return, X should be overwritten by
34: * A * X, if KASE=1,
35: * A' * X, if KASE=2,
36: * where A' is the conjugate transpose of A, and ZLACON must be
37: * re-called with all the other parameters unchanged.
38: *
39: * EST (input/output) DOUBLE PRECISION
40: * On entry with KASE = 1 or 2 and JUMP = 3, EST should be
41: * unchanged from the previous call to ZLACON.
42: * On exit, EST is an estimate (a lower bound) for norm(A).
43: *
44: * KASE (input/output) INTEGER
45: * On the initial call to ZLACON, KASE should be 0.
46: * On an intermediate return, KASE will be 1 or 2, indicating
47: * whether X should be overwritten by A * X or A' * X.
48: * On the final return from ZLACON, KASE will again be 0.
49: *
50: * Further Details
51: * ======= =======
52: *
53: * Contributed by Nick Higham, University of Manchester.
54: * Originally named CONEST, dated March 16, 1988.
55: *
56: * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
57: * a real or complex matrix, with applications to condition estimation",
58: * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
59: *
60: * Last modified: April, 1999
61: *
62: * =====================================================================
63: *
64: * .. Parameters ..
65: INTEGER ITMAX
66: PARAMETER ( ITMAX = 5 )
67: DOUBLE PRECISION ONE, TWO
68: PARAMETER ( ONE = 1.0D0, TWO = 2.0D0 )
69: COMPLEX*16 CZERO, CONE
70: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
71: $ CONE = ( 1.0D0, 0.0D0 ) )
72: * ..
73: * .. Local Scalars ..
74: INTEGER I, ITER, J, JLAST, JUMP
75: DOUBLE PRECISION ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
76: * ..
77: * .. External Functions ..
78: INTEGER IZMAX1
79: DOUBLE PRECISION DLAMCH, DZSUM1
80: EXTERNAL IZMAX1, DLAMCH, DZSUM1
81: * ..
82: * .. External Subroutines ..
83: EXTERNAL ZCOPY
84: * ..
85: * .. Intrinsic Functions ..
86: INTRINSIC ABS, DBLE, DCMPLX, DIMAG
87: * ..
88: * .. Save statement ..
89: SAVE
90: * ..
91: * .. Executable Statements ..
92: *
93: SAFMIN = DLAMCH( 'Safe minimum' )
94: IF( KASE.EQ.0 ) THEN
95: DO 10 I = 1, N
96: X( I ) = DCMPLX( ONE / DBLE( N ) )
97: 10 CONTINUE
98: KASE = 1
99: JUMP = 1
100: RETURN
101: END IF
102: *
103: GO TO ( 20, 40, 70, 90, 120 )JUMP
104: *
105: * ................ ENTRY (JUMP = 1)
106: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
107: *
108: 20 CONTINUE
109: IF( N.EQ.1 ) THEN
110: V( 1 ) = X( 1 )
111: EST = ABS( V( 1 ) )
112: * ... QUIT
113: GO TO 130
114: END IF
115: EST = DZSUM1( N, X, 1 )
116: *
117: DO 30 I = 1, N
118: ABSXI = ABS( X( I ) )
119: IF( ABSXI.GT.SAFMIN ) THEN
120: X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
121: $ DIMAG( X( I ) ) / ABSXI )
122: ELSE
123: X( I ) = CONE
124: END IF
125: 30 CONTINUE
126: KASE = 2
127: JUMP = 2
128: RETURN
129: *
130: * ................ ENTRY (JUMP = 2)
131: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
132: *
133: 40 CONTINUE
134: J = IZMAX1( N, X, 1 )
135: ITER = 2
136: *
137: * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
138: *
139: 50 CONTINUE
140: DO 60 I = 1, N
141: X( I ) = CZERO
142: 60 CONTINUE
143: X( J ) = CONE
144: KASE = 1
145: JUMP = 3
146: RETURN
147: *
148: * ................ ENTRY (JUMP = 3)
149: * X HAS BEEN OVERWRITTEN BY A*X.
150: *
151: 70 CONTINUE
152: CALL ZCOPY( N, X, 1, V, 1 )
153: ESTOLD = EST
154: EST = DZSUM1( N, V, 1 )
155: *
156: * TEST FOR CYCLING.
157: IF( EST.LE.ESTOLD )
158: $ GO TO 100
159: *
160: DO 80 I = 1, N
161: ABSXI = ABS( X( I ) )
162: IF( ABSXI.GT.SAFMIN ) THEN
163: X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
164: $ DIMAG( X( I ) ) / ABSXI )
165: ELSE
166: X( I ) = CONE
167: END IF
168: 80 CONTINUE
169: KASE = 2
170: JUMP = 4
171: RETURN
172: *
173: * ................ ENTRY (JUMP = 4)
174: * X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
175: *
176: 90 CONTINUE
177: JLAST = J
178: J = IZMAX1( N, X, 1 )
179: IF( ( ABS( X( JLAST ) ).NE.ABS( X( J ) ) ) .AND.
180: $ ( ITER.LT.ITMAX ) ) THEN
181: ITER = ITER + 1
182: GO TO 50
183: END IF
184: *
185: * ITERATION COMPLETE. FINAL STAGE.
186: *
187: 100 CONTINUE
188: ALTSGN = ONE
189: DO 110 I = 1, N
190: X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
191: ALTSGN = -ALTSGN
192: 110 CONTINUE
193: KASE = 1
194: JUMP = 5
195: RETURN
196: *
197: * ................ ENTRY (JUMP = 5)
198: * X HAS BEEN OVERWRITTEN BY A*X.
199: *
200: 120 CONTINUE
201: TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
202: IF( TEMP.GT.EST ) THEN
203: CALL ZCOPY( N, X, 1, V, 1 )
204: EST = TEMP
205: END IF
206: *
207: 130 CONTINUE
208: KASE = 0
209: RETURN
210: *
211: * End of ZLACON
212: *
213: END
CVSweb interface <joel.bertrand@systella.fr>