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