1: *> \brief \b DLACN2 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 DLACN2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlacn2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlacn2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlacn2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER KASE, N
25: * DOUBLE PRECISION EST
26: * ..
27: * .. Array Arguments ..
28: * INTEGER ISGN( * ), ISAVE( 3 )
29: * DOUBLE PRECISION V( * ), X( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLACN2 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 DLACN2 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 ISAVE(1) = 3, EST should be
77: *> unchanged from the previous call to DLACN2.
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 DLACN2, 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 DLACN2, KASE will again be 0.
88: *> \endverbatim
89: *>
90: *> \param[in,out] ISAVE
91: *> \verbatim
92: *> ISAVE is INTEGER array, dimension (3)
93: *> ISAVE is used to save variables between calls to DLACN2
94: *> \endverbatim
95: *
96: * Authors:
97: * ========
98: *
99: *> \author Univ. of Tennessee
100: *> \author Univ. of California Berkeley
101: *> \author Univ. of Colorado Denver
102: *> \author NAG Ltd.
103: *
104: *> \ingroup doubleOTHERauxiliary
105: *
106: *> \par Further Details:
107: * =====================
108: *>
109: *> \verbatim
110: *>
111: *> Originally named SONEST, dated March 16, 1988.
112: *>
113: *> This is a thread safe version of DLACON, which uses the array ISAVE
114: *> in place of a SAVE statement, as follows:
115: *>
116: *> DLACON DLACN2
117: *> JUMP ISAVE(1)
118: *> J ISAVE(2)
119: *> ITER ISAVE(3)
120: *> \endverbatim
121: *
122: *> \par Contributors:
123: * ==================
124: *>
125: *> Nick Higham, University of Manchester
126: *
127: *> \par References:
128: * ================
129: *>
130: *> N.J. Higham, "FORTRAN codes for estimating the one-norm of
131: *> a real or complex matrix, with applications to condition estimation",
132: *> ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
133: *>
134: * =====================================================================
135: SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE )
136: *
137: * -- LAPACK auxiliary routine --
138: * -- LAPACK is a software package provided by Univ. of Tennessee, --
139: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
140: *
141: * .. Scalar Arguments ..
142: INTEGER KASE, N
143: DOUBLE PRECISION EST
144: * ..
145: * .. Array Arguments ..
146: INTEGER ISGN( * ), ISAVE( 3 )
147: DOUBLE PRECISION V( * ), X( * )
148: * ..
149: *
150: * =====================================================================
151: *
152: * .. Parameters ..
153: INTEGER ITMAX
154: PARAMETER ( ITMAX = 5 )
155: DOUBLE PRECISION ZERO, ONE, TWO
156: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
157: * ..
158: * .. Local Scalars ..
159: INTEGER I, JLAST
160: DOUBLE PRECISION ALTSGN, ESTOLD, TEMP, XS
161: * ..
162: * .. External Functions ..
163: INTEGER IDAMAX
164: DOUBLE PRECISION DASUM
165: EXTERNAL IDAMAX, DASUM
166: * ..
167: * .. External Subroutines ..
168: EXTERNAL DCOPY
169: * ..
170: * .. Intrinsic Functions ..
171: INTRINSIC ABS, DBLE, NINT
172: * ..
173: * .. Executable Statements ..
174: *
175: IF( KASE.EQ.0 ) THEN
176: DO 10 I = 1, N
177: X( I ) = ONE / DBLE( N )
178: 10 CONTINUE
179: KASE = 1
180: ISAVE( 1 ) = 1
181: RETURN
182: END IF
183: *
184: GO TO ( 20, 40, 70, 110, 140 )ISAVE( 1 )
185: *
186: * ................ ENTRY (ISAVE( 1 ) = 1)
187: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
188: *
189: 20 CONTINUE
190: IF( N.EQ.1 ) THEN
191: V( 1 ) = X( 1 )
192: EST = ABS( V( 1 ) )
193: * ... QUIT
194: GO TO 150
195: END IF
196: EST = DASUM( N, X, 1 )
197: *
198: DO 30 I = 1, N
199: IF( X(I).GE.ZERO ) THEN
200: X(I) = ONE
201: ELSE
202: X(I) = -ONE
203: END IF
204: ISGN( I ) = NINT( X( I ) )
205: 30 CONTINUE
206: KASE = 2
207: ISAVE( 1 ) = 2
208: RETURN
209: *
210: * ................ ENTRY (ISAVE( 1 ) = 2)
211: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
212: *
213: 40 CONTINUE
214: ISAVE( 2 ) = IDAMAX( N, X, 1 )
215: ISAVE( 3 ) = 2
216: *
217: * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
218: *
219: 50 CONTINUE
220: DO 60 I = 1, N
221: X( I ) = ZERO
222: 60 CONTINUE
223: X( ISAVE( 2 ) ) = ONE
224: KASE = 1
225: ISAVE( 1 ) = 3
226: RETURN
227: *
228: * ................ ENTRY (ISAVE( 1 ) = 3)
229: * X HAS BEEN OVERWRITTEN BY A*X.
230: *
231: 70 CONTINUE
232: CALL DCOPY( N, X, 1, V, 1 )
233: ESTOLD = EST
234: EST = DASUM( N, V, 1 )
235: DO 80 I = 1, N
236: IF( X(I).GE.ZERO ) THEN
237: XS = ONE
238: ELSE
239: XS = -ONE
240: END IF
241: IF( NINT( XS ).NE.ISGN( I ) )
242: $ GO TO 90
243: 80 CONTINUE
244: * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
245: GO TO 120
246: *
247: 90 CONTINUE
248: * TEST FOR CYCLING.
249: IF( EST.LE.ESTOLD )
250: $ GO TO 120
251: *
252: DO 100 I = 1, N
253: IF( X(I).GE.ZERO ) THEN
254: X(I) = ONE
255: ELSE
256: X(I) = -ONE
257: END IF
258: ISGN( I ) = NINT( X( I ) )
259: 100 CONTINUE
260: KASE = 2
261: ISAVE( 1 ) = 4
262: RETURN
263: *
264: * ................ ENTRY (ISAVE( 1 ) = 4)
265: * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
266: *
267: 110 CONTINUE
268: JLAST = ISAVE( 2 )
269: ISAVE( 2 ) = IDAMAX( N, X, 1 )
270: IF( ( X( JLAST ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND.
271: $ ( ISAVE( 3 ).LT.ITMAX ) ) THEN
272: ISAVE( 3 ) = ISAVE( 3 ) + 1
273: GO TO 50
274: END IF
275: *
276: * ITERATION COMPLETE. FINAL STAGE.
277: *
278: 120 CONTINUE
279: ALTSGN = ONE
280: DO 130 I = 1, N
281: X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) )
282: ALTSGN = -ALTSGN
283: 130 CONTINUE
284: KASE = 1
285: ISAVE( 1 ) = 5
286: RETURN
287: *
288: * ................ ENTRY (ISAVE( 1 ) = 5)
289: * X HAS BEEN OVERWRITTEN BY A*X.
290: *
291: 140 CONTINUE
292: TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
293: IF( TEMP.GT.EST ) THEN
294: CALL DCOPY( N, X, 1, V, 1 )
295: EST = TEMP
296: END IF
297: *
298: 150 CONTINUE
299: KASE = 0
300: RETURN
301: *
302: * End of DLACN2
303: *
304: END
CVSweb interface <joel.bertrand@systella.fr>