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