|
7034
|
1 SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE ) |
|
|
2 * |
|
|
3 * -- LAPACK auxiliary routine (version 3.1) -- |
|
|
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
5 * November 2006 |
|
|
6 * |
|
|
7 * .. Scalar Arguments .. |
|
|
8 INTEGER KASE, N |
|
|
9 DOUBLE PRECISION EST |
|
|
10 * .. |
|
|
11 * .. Array Arguments .. |
|
|
12 INTEGER ISGN( * ), ISAVE( 3 ) |
|
|
13 DOUBLE PRECISION V( * ), X( * ) |
|
|
14 * .. |
|
|
15 * |
|
|
16 * Purpose |
|
|
17 * ======= |
|
|
18 * |
|
|
19 * DLACN2 estimates the 1-norm of a square, real 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 * and DLACN2 must be re-called with all the other parameters |
|
|
37 * unchanged. |
|
|
38 * |
|
|
39 * ISGN (workspace) INTEGER array, dimension (N) |
|
|
40 * |
|
|
41 * EST (input/output) DOUBLE PRECISION |
|
|
42 * On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be |
|
|
43 * unchanged from the previous call to DLACN2. |
|
|
44 * On exit, EST is an estimate (a lower bound) for norm(A). |
|
|
45 * |
|
|
46 * KASE (input/output) INTEGER |
|
|
47 * On the initial call to DLACN2, KASE should be 0. |
|
|
48 * On an intermediate return, KASE will be 1 or 2, indicating |
|
|
49 * whether X should be overwritten by A * X or A' * X. |
|
|
50 * On the final return from DLACN2, KASE will again be 0. |
|
|
51 * |
|
|
52 * ISAVE (input/output) INTEGER array, dimension (3) |
|
|
53 * ISAVE is used to save variables between calls to DLACN2 |
|
|
54 * |
|
|
55 * Further Details |
|
|
56 * ======= ======= |
|
|
57 * |
|
|
58 * Contributed by Nick Higham, University of Manchester. |
|
|
59 * Originally named SONEST, dated March 16, 1988. |
|
|
60 * |
|
|
61 * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of |
|
|
62 * a real or complex matrix, with applications to condition estimation", |
|
|
63 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. |
|
|
64 * |
|
|
65 * This is a thread safe version of DLACON, which uses the array ISAVE |
|
|
66 * in place of a SAVE statement, as follows: |
|
|
67 * |
|
|
68 * DLACON DLACN2 |
|
|
69 * JUMP ISAVE(1) |
|
|
70 * J ISAVE(2) |
|
|
71 * ITER ISAVE(3) |
|
|
72 * |
|
|
73 * ===================================================================== |
|
|
74 * |
|
|
75 * .. Parameters .. |
|
|
76 INTEGER ITMAX |
|
|
77 PARAMETER ( ITMAX = 5 ) |
|
|
78 DOUBLE PRECISION ZERO, ONE, TWO |
|
|
79 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) |
|
|
80 * .. |
|
|
81 * .. Local Scalars .. |
|
|
82 INTEGER I, JLAST |
|
|
83 DOUBLE PRECISION ALTSGN, ESTOLD, TEMP |
|
|
84 * .. |
|
|
85 * .. External Functions .. |
|
|
86 INTEGER IDAMAX |
|
|
87 DOUBLE PRECISION DASUM |
|
|
88 EXTERNAL IDAMAX, DASUM |
|
|
89 * .. |
|
|
90 * .. External Subroutines .. |
|
|
91 EXTERNAL DCOPY |
|
|
92 * .. |
|
|
93 * .. Intrinsic Functions .. |
|
|
94 INTRINSIC ABS, DBLE, NINT, SIGN |
|
|
95 * .. |
|
|
96 * .. Executable Statements .. |
|
|
97 * |
|
|
98 IF( KASE.EQ.0 ) THEN |
|
|
99 DO 10 I = 1, N |
|
|
100 X( I ) = ONE / DBLE( N ) |
|
|
101 10 CONTINUE |
|
|
102 KASE = 1 |
|
|
103 ISAVE( 1 ) = 1 |
|
|
104 RETURN |
|
|
105 END IF |
|
|
106 * |
|
|
107 GO TO ( 20, 40, 70, 110, 140 )ISAVE( 1 ) |
|
|
108 * |
|
|
109 * ................ ENTRY (ISAVE( 1 ) = 1) |
|
|
110 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. |
|
|
111 * |
|
|
112 20 CONTINUE |
|
|
113 IF( N.EQ.1 ) THEN |
|
|
114 V( 1 ) = X( 1 ) |
|
|
115 EST = ABS( V( 1 ) ) |
|
|
116 * ... QUIT |
|
|
117 GO TO 150 |
|
|
118 END IF |
|
|
119 EST = DASUM( N, X, 1 ) |
|
|
120 * |
|
|
121 DO 30 I = 1, N |
|
|
122 X( I ) = SIGN( ONE, X( I ) ) |
|
|
123 ISGN( I ) = NINT( X( I ) ) |
|
|
124 30 CONTINUE |
|
|
125 KASE = 2 |
|
|
126 ISAVE( 1 ) = 2 |
|
|
127 RETURN |
|
|
128 * |
|
|
129 * ................ ENTRY (ISAVE( 1 ) = 2) |
|
|
130 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. |
|
|
131 * |
|
|
132 40 CONTINUE |
|
|
133 ISAVE( 2 ) = IDAMAX( N, X, 1 ) |
|
|
134 ISAVE( 3 ) = 2 |
|
|
135 * |
|
|
136 * MAIN LOOP - ITERATIONS 2,3,...,ITMAX. |
|
|
137 * |
|
|
138 50 CONTINUE |
|
|
139 DO 60 I = 1, N |
|
|
140 X( I ) = ZERO |
|
|
141 60 CONTINUE |
|
|
142 X( ISAVE( 2 ) ) = ONE |
|
|
143 KASE = 1 |
|
|
144 ISAVE( 1 ) = 3 |
|
|
145 RETURN |
|
|
146 * |
|
|
147 * ................ ENTRY (ISAVE( 1 ) = 3) |
|
|
148 * X HAS BEEN OVERWRITTEN BY A*X. |
|
|
149 * |
|
|
150 70 CONTINUE |
|
|
151 CALL DCOPY( N, X, 1, V, 1 ) |
|
|
152 ESTOLD = EST |
|
|
153 EST = DASUM( N, V, 1 ) |
|
|
154 DO 80 I = 1, N |
|
|
155 IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) ) |
|
|
156 $ GO TO 90 |
|
|
157 80 CONTINUE |
|
|
158 * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. |
|
|
159 GO TO 120 |
|
|
160 * |
|
|
161 90 CONTINUE |
|
|
162 * TEST FOR CYCLING. |
|
|
163 IF( EST.LE.ESTOLD ) |
|
|
164 $ GO TO 120 |
|
|
165 * |
|
|
166 DO 100 I = 1, N |
|
|
167 X( I ) = SIGN( ONE, X( I ) ) |
|
|
168 ISGN( I ) = NINT( X( I ) ) |
|
|
169 100 CONTINUE |
|
|
170 KASE = 2 |
|
|
171 ISAVE( 1 ) = 4 |
|
|
172 RETURN |
|
|
173 * |
|
|
174 * ................ ENTRY (ISAVE( 1 ) = 4) |
|
|
175 * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. |
|
|
176 * |
|
|
177 110 CONTINUE |
|
|
178 JLAST = ISAVE( 2 ) |
|
|
179 ISAVE( 2 ) = IDAMAX( N, X, 1 ) |
|
|
180 IF( ( X( JLAST ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND. |
|
|
181 $ ( ISAVE( 3 ).LT.ITMAX ) ) THEN |
|
|
182 ISAVE( 3 ) = ISAVE( 3 ) + 1 |
|
|
183 GO TO 50 |
|
|
184 END IF |
|
|
185 * |
|
|
186 * ITERATION COMPLETE. FINAL STAGE. |
|
|
187 * |
|
|
188 120 CONTINUE |
|
|
189 ALTSGN = ONE |
|
|
190 DO 130 I = 1, N |
|
|
191 X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) |
|
|
192 ALTSGN = -ALTSGN |
|
|
193 130 CONTINUE |
|
|
194 KASE = 1 |
|
|
195 ISAVE( 1 ) = 5 |
|
|
196 RETURN |
|
|
197 * |
|
|
198 * ................ ENTRY (ISAVE( 1 ) = 5) |
|
|
199 * X HAS BEEN OVERWRITTEN BY A*X. |
|
|
200 * |
|
|
201 140 CONTINUE |
|
|
202 TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) ) |
|
|
203 IF( TEMP.GT.EST ) THEN |
|
|
204 CALL DCOPY( N, X, 1, V, 1 ) |
|
|
205 EST = TEMP |
|
|
206 END IF |
|
|
207 * |
|
|
208 150 CONTINUE |
|
|
209 KASE = 0 |
|
|
210 RETURN |
|
|
211 * |
|
|
212 * End of DLACN2 |
|
|
213 * |
|
|
214 END |
