1 SUBROUTINE DLACON( N, V, X, ISGN, 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 INTEGER ISGN( * )
14 DOUBLE PRECISION V( * ), X( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLACON 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 DLACON 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 JUMP = 3, EST should be
44 * unchanged from the previous call to DLACON.
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 DLACON, 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 DLACON, KASE will again be 0.
52 *
53 * Further Details
54 * ======= =======
55 *
56 * Contributed by Nick Higham, University of Manchester.
57 * Originally named SONEST, dated March 16, 1988.
58 *
59 * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
60 * a real or complex matrix, with applications to condition estimation",
61 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
62 *
63 * =====================================================================
64 *
65 * .. Parameters ..
66 INTEGER ITMAX
67 PARAMETER ( ITMAX = 5 )
68 DOUBLE PRECISION ZERO, ONE, TWO
69 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
70 * ..
71 * .. Local Scalars ..
72 INTEGER I, ITER, J, JLAST, JUMP
73 DOUBLE PRECISION ALTSGN, ESTOLD, TEMP
74 * ..
75 * .. External Functions ..
76 INTEGER IDAMAX
77 DOUBLE PRECISION DASUM
78 EXTERNAL IDAMAX, DASUM
79 * ..
80 * .. External Subroutines ..
81 EXTERNAL DCOPY
82 * ..
83 * .. Intrinsic Functions ..
84 INTRINSIC ABS, DBLE, NINT, SIGN
85 * ..
86 * .. Save statement ..
87 SAVE
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( KASE.EQ.0 ) THEN
92 DO 10 I = 1, N
93 X( I ) = ONE / DBLE( N )
94 10 CONTINUE
95 KASE = 1
96 JUMP = 1
97 RETURN
98 END IF
99 *
100 GO TO ( 20, 40, 70, 110, 140 )JUMP
101 *
102 * ................ ENTRY (JUMP = 1)
103 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
104 *
105 20 CONTINUE
106 IF( N.EQ.1 ) THEN
107 V( 1 ) = X( 1 )
108 EST = ABS( V( 1 ) )
109 * ... QUIT
110 GO TO 150
111 END IF
112 EST = DASUM( N, X, 1 )
113 *
114 DO 30 I = 1, N
115 X( I ) = SIGN( ONE, X( I ) )
116 ISGN( I ) = NINT( X( I ) )
117 30 CONTINUE
118 KASE = 2
119 JUMP = 2
120 RETURN
121 *
122 * ................ ENTRY (JUMP = 2)
123 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
124 *
125 40 CONTINUE
126 J = IDAMAX( N, X, 1 )
127 ITER = 2
128 *
129 * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
130 *
131 50 CONTINUE
132 DO 60 I = 1, N
133 X( I ) = ZERO
134 60 CONTINUE
135 X( J ) = ONE
136 KASE = 1
137 JUMP = 3
138 RETURN
139 *
140 * ................ ENTRY (JUMP = 3)
141 * X HAS BEEN OVERWRITTEN BY A*X.
142 *
143 70 CONTINUE
144 CALL DCOPY( N, X, 1, V, 1 )
145 ESTOLD = EST
146 EST = DASUM( N, V, 1 )
147 DO 80 I = 1, N
148 IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) )
149 $ GO TO 90
150 80 CONTINUE
151 * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
152 GO TO 120
153 *
154 90 CONTINUE
155 * TEST FOR CYCLING.
156 IF( EST.LE.ESTOLD )
157 $ GO TO 120
158 *
159 DO 100 I = 1, N
160 X( I ) = SIGN( ONE, X( I ) )
161 ISGN( I ) = NINT( X( I ) )
162 100 CONTINUE
163 KASE = 2
164 JUMP = 4
165 RETURN
166 *
167 * ................ ENTRY (JUMP = 4)
168 * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
169 *
170 110 CONTINUE
171 JLAST = J
172 J = IDAMAX( N, X, 1 )
173 IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN
174 ITER = ITER + 1
175 GO TO 50
176 END IF
177 *
178 * ITERATION COMPLETE. FINAL STAGE.
179 *
180 120 CONTINUE
181 ALTSGN = ONE
182 DO 130 I = 1, N
183 X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) )
184 ALTSGN = -ALTSGN
185 130 CONTINUE
186 KASE = 1
187 JUMP = 5
188 RETURN
189 *
190 * ................ ENTRY (JUMP = 5)
191 * X HAS BEEN OVERWRITTEN BY A*X.
192 *
193 140 CONTINUE
194 TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
195 IF( TEMP.GT.EST ) THEN
196 CALL DCOPY( N, X, 1, V, 1 )
197 EST = TEMP
198 END IF
199 *
200 150 CONTINUE
201 KASE = 0
202 RETURN
203 *
204 * End of DLACON
205 *
206 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( * )
14 DOUBLE PRECISION V( * ), X( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DLACON 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 DLACON 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 JUMP = 3, EST should be
44 * unchanged from the previous call to DLACON.
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 DLACON, 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 DLACON, KASE will again be 0.
52 *
53 * Further Details
54 * ======= =======
55 *
56 * Contributed by Nick Higham, University of Manchester.
57 * Originally named SONEST, dated March 16, 1988.
58 *
59 * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
60 * a real or complex matrix, with applications to condition estimation",
61 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
62 *
63 * =====================================================================
64 *
65 * .. Parameters ..
66 INTEGER ITMAX
67 PARAMETER ( ITMAX = 5 )
68 DOUBLE PRECISION ZERO, ONE, TWO
69 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
70 * ..
71 * .. Local Scalars ..
72 INTEGER I, ITER, J, JLAST, JUMP
73 DOUBLE PRECISION ALTSGN, ESTOLD, TEMP
74 * ..
75 * .. External Functions ..
76 INTEGER IDAMAX
77 DOUBLE PRECISION DASUM
78 EXTERNAL IDAMAX, DASUM
79 * ..
80 * .. External Subroutines ..
81 EXTERNAL DCOPY
82 * ..
83 * .. Intrinsic Functions ..
84 INTRINSIC ABS, DBLE, NINT, SIGN
85 * ..
86 * .. Save statement ..
87 SAVE
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( KASE.EQ.0 ) THEN
92 DO 10 I = 1, N
93 X( I ) = ONE / DBLE( N )
94 10 CONTINUE
95 KASE = 1
96 JUMP = 1
97 RETURN
98 END IF
99 *
100 GO TO ( 20, 40, 70, 110, 140 )JUMP
101 *
102 * ................ ENTRY (JUMP = 1)
103 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
104 *
105 20 CONTINUE
106 IF( N.EQ.1 ) THEN
107 V( 1 ) = X( 1 )
108 EST = ABS( V( 1 ) )
109 * ... QUIT
110 GO TO 150
111 END IF
112 EST = DASUM( N, X, 1 )
113 *
114 DO 30 I = 1, N
115 X( I ) = SIGN( ONE, X( I ) )
116 ISGN( I ) = NINT( X( I ) )
117 30 CONTINUE
118 KASE = 2
119 JUMP = 2
120 RETURN
121 *
122 * ................ ENTRY (JUMP = 2)
123 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
124 *
125 40 CONTINUE
126 J = IDAMAX( N, X, 1 )
127 ITER = 2
128 *
129 * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
130 *
131 50 CONTINUE
132 DO 60 I = 1, N
133 X( I ) = ZERO
134 60 CONTINUE
135 X( J ) = ONE
136 KASE = 1
137 JUMP = 3
138 RETURN
139 *
140 * ................ ENTRY (JUMP = 3)
141 * X HAS BEEN OVERWRITTEN BY A*X.
142 *
143 70 CONTINUE
144 CALL DCOPY( N, X, 1, V, 1 )
145 ESTOLD = EST
146 EST = DASUM( N, V, 1 )
147 DO 80 I = 1, N
148 IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) )
149 $ GO TO 90
150 80 CONTINUE
151 * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
152 GO TO 120
153 *
154 90 CONTINUE
155 * TEST FOR CYCLING.
156 IF( EST.LE.ESTOLD )
157 $ GO TO 120
158 *
159 DO 100 I = 1, N
160 X( I ) = SIGN( ONE, X( I ) )
161 ISGN( I ) = NINT( X( I ) )
162 100 CONTINUE
163 KASE = 2
164 JUMP = 4
165 RETURN
166 *
167 * ................ ENTRY (JUMP = 4)
168 * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X.
169 *
170 110 CONTINUE
171 JLAST = J
172 J = IDAMAX( N, X, 1 )
173 IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN
174 ITER = ITER + 1
175 GO TO 50
176 END IF
177 *
178 * ITERATION COMPLETE. FINAL STAGE.
179 *
180 120 CONTINUE
181 ALTSGN = ONE
182 DO 130 I = 1, N
183 X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) )
184 ALTSGN = -ALTSGN
185 130 CONTINUE
186 KASE = 1
187 JUMP = 5
188 RETURN
189 *
190 * ................ ENTRY (JUMP = 5)
191 * X HAS BEEN OVERWRITTEN BY A*X.
192 *
193 140 CONTINUE
194 TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
195 IF( TEMP.GT.EST ) THEN
196 CALL DCOPY( N, X, 1, V, 1 )
197 EST = TEMP
198 END IF
199 *
200 150 CONTINUE
201 KASE = 0
202 RETURN
203 *
204 * End of DLACON
205 *
206 END