KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlacon.f
Go to the documentation of this file.
1  SUBROUTINE dlacon( N, V, X, ISGN, EST, KASE )
2 *
3 * -- LAPACK auxiliary routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * February 29, 1992
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' * 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 (output) DOUBLE PRECISION
43 * An estimate (a lower bound) for norm(A).
44 *
45 * KASE (input/output) INTEGER
46 * On the initial call to DLACON, KASE should be 0.
47 * On an intermediate return, KASE will be 1 or 2, indicating
48 * whether X should be overwritten by A * X or A' * X.
49 * On the final return from DLACON, KASE will again be 0.
50 *
51 * Further Details
52 * ======= =======
53 *
54 * Contributed by Nick Higham, University of Manchester.
55 * Originally named SONEST, dated March 16, 1988.
56 *
57 * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
58 * a real or complex matrix, with applications to condition estimation",
59 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
60 *
61 * =====================================================================
62 *
63 * .. Parameters ..
64  INTEGER ITMAX
65  parameter( itmax = 5 )
66  DOUBLE PRECISION ZERO, ONE, TWO
67  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
68 * ..
69 * .. Local Scalars ..
70  INTEGER I, ITER, J, JLAST, JUMP
71  DOUBLE PRECISION ALTSGN, ESTOLD, TEMP
72 * ..
73 * .. External Functions ..
74  INTEGER IDAMAX
75  DOUBLE PRECISION DASUM
76  EXTERNAL idamax, dasum
77 * ..
78 * .. External Subroutines ..
79  EXTERNAL dcopy
80 * ..
81 * .. Intrinsic Functions ..
82  INTRINSIC abs, dble, nint, sign
83 * ..
84 * .. Save statement ..
85  SAVE
86 * ..
87 * .. Executable Statements ..
88 *
89  IF( kase.EQ.0 ) THEN
90  DO 10 i = 1, n
91  x( i ) = one / dble( n )
92  10 CONTINUE
93  kase = 1
94  jump = 1
95  RETURN
96  END IF
97 *
98  GO TO ( 20, 40, 70, 110, 140 )jump
99 *
100 * ................ ENTRY (JUMP = 1)
101 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
102 *
103  20 CONTINUE
104  IF( n.EQ.1 ) THEN
105  v( 1 ) = x( 1 )
106  est = abs( v( 1 ) )
107 * ... QUIT
108  GO TO 150
109  END IF
110  est = dasum( n, x, 1 )
111 *
112  DO 30 i = 1, n
113  x( i ) = sign( one, x( i ) )
114  isgn( i ) = nint( x( i ) )
115  30 CONTINUE
116  kase = 2
117  jump = 2
118  RETURN
119 *
120 * ................ ENTRY (JUMP = 2)
121 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
122 *
123  40 CONTINUE
124  j = idamax( n, x, 1 )
125  iter = 2
126 *
127 * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
128 *
129  50 CONTINUE
130  DO 60 i = 1, n
131  x( i ) = zero
132  60 CONTINUE
133  x( j ) = one
134  kase = 1
135  jump = 3
136  RETURN
137 *
138 * ................ ENTRY (JUMP = 3)
139 * X HAS BEEN OVERWRITTEN BY A*X.
140 *
141  70 CONTINUE
142  CALL dcopy( n, x, 1, v, 1 )
143  estold = est
144  est = dasum( n, v, 1 )
145  DO 80 i = 1, n
146  IF( nint( sign( one, x( i ) ) ).NE.isgn( i ) )
147  $ GO TO 90
148  80 CONTINUE
149 * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
150  GO TO 120
151 *
152  90 CONTINUE
153 * TEST FOR CYCLING.
154  IF( est.LE.estold )
155  $ GO TO 120
156 *
157  DO 100 i = 1, n
158  x( i ) = sign( one, x( i ) )
159  isgn( i ) = nint( x( i ) )
160  100 CONTINUE
161  kase = 2
162  jump = 4
163  RETURN
164 *
165 * ................ ENTRY (JUMP = 4)
166 * X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
167 *
168  110 CONTINUE
169  jlast = j
170  j = idamax( n, x, 1 )
171  IF( ( x( jlast ).NE.abs( x( j ) ) ) .AND. ( iter.LT.itmax ) ) THEN
172  iter = iter + 1
173  GO TO 50
174  END IF
175 *
176 * ITERATION COMPLETE. FINAL STAGE.
177 *
178  120 CONTINUE
179  altsgn = one
180  DO 130 i = 1, n
181  x( i ) = altsgn*( one+dble( i-1 ) / dble( n-1 ) )
182  altsgn = -altsgn
183  130 CONTINUE
184  kase = 1
185  jump = 5
186  RETURN
187 *
188 * ................ ENTRY (JUMP = 5)
189 * X HAS BEEN OVERWRITTEN BY A*X.
190 *
191  140 CONTINUE
192  temp = two*( dasum( n, x, 1 ) / dble( 3*n ) )
193  IF( temp.GT.est ) THEN
194  CALL dcopy( n, x, 1, v, 1 )
195  est = temp
196  END IF
197 *
198  150 CONTINUE
199  kase = 0
200  RETURN
201 *
202 * End of DLACON
203 *
204  END
subroutine dcopy(n, dx, incx, dy, incy)
Definition: dcopy.f:2
subroutine dlacon(N, V, X, ISGN, EST, KASE)
Definition: dlacon.f:2