Mercurial > octave-nkf

comparison libcruft/lapack/zlacn2.f @ 7034:68db500cb558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-16 18:54:19 by jwe]
author jwe
date Tue, 16 Oct 2007 18:54:23 +0000
parents
children
comparison
equal deleted inserted replaced
7033:f0142f2afdc6 7034:68db500cb558
1 SUBROUTINE ZLACN2( N, V, X, 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 ISAVE( 3 )
13 COMPLEX*16 V( * ), X( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZLACN2 estimates the 1-norm of a square, complex 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) COMPLEX*16 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) COMPLEX*16 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 * where A' is the conjugate transpose of A, and ZLACN2 must be
37 * re-called with all the other parameters unchanged.
38 *
39 * EST (input/output) DOUBLE PRECISION
40 * On entry with KASE = 1 or 2 and ISAVE(1) = 3, EST should be
41 * unchanged from the previous call to ZLACN2.
42 * On exit, EST is an estimate (a lower bound) for norm(A).
43 *
44 * KASE (input/output) INTEGER
45 * On the initial call to ZLACN2, KASE should be 0.
46 * On an intermediate return, KASE will be 1 or 2, indicating
47 * whether X should be overwritten by A * X or A' * X.
48 * On the final return from ZLACN2, KASE will again be 0.
49 *
50 * ISAVE (input/output) INTEGER array, dimension (3)
51 * ISAVE is used to save variables between calls to ZLACN2
52 *
53 * Further Details
54 * ======= =======
55 *
56 * Contributed by Nick Higham, University of Manchester.
57 * Originally named CONEST, 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 * Last modified: April, 1999
64 *
65 * This is a thread safe version of ZLACON, which uses the array ISAVE
66 * in place of a SAVE statement, as follows:
67 *
68 * ZLACON ZLACN2
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 ONE, TWO
79 PARAMETER ( ONE = 1.0D0, TWO = 2.0D0 )
80 COMPLEX*16 CZERO, CONE
81 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
82 $ CONE = ( 1.0D0, 0.0D0 ) )
83 * ..
84 * .. Local Scalars ..
85 INTEGER I, JLAST
86 DOUBLE PRECISION ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
87 * ..
88 * .. External Functions ..
89 INTEGER IZMAX1
90 DOUBLE PRECISION DLAMCH, DZSUM1
91 EXTERNAL IZMAX1, DLAMCH, DZSUM1
92 * ..
93 * .. External Subroutines ..
94 EXTERNAL ZCOPY
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC ABS, DBLE, DCMPLX, DIMAG
98 * ..
99 * .. Executable Statements ..
100 *
101 SAFMIN = DLAMCH( 'Safe minimum' )
102 IF( KASE.EQ.0 ) THEN
103 DO 10 I = 1, N
104 X( I ) = DCMPLX( ONE / DBLE( N ) )
105 10 CONTINUE
106 KASE = 1
107 ISAVE( 1 ) = 1
108 RETURN
109 END IF
110 *
111 GO TO ( 20, 40, 70, 90, 120 )ISAVE( 1 )
112 *
113 * ................ ENTRY (ISAVE( 1 ) = 1)
114 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
115 *
116 20 CONTINUE
117 IF( N.EQ.1 ) THEN
118 V( 1 ) = X( 1 )
119 EST = ABS( V( 1 ) )
120 * ... QUIT
121 GO TO 130
122 END IF
123 EST = DZSUM1( N, X, 1 )
124 *
125 DO 30 I = 1, N
126 ABSXI = ABS( X( I ) )
127 IF( ABSXI.GT.SAFMIN ) THEN
128 X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
129 $ DIMAG( X( I ) ) / ABSXI )
130 ELSE
131 X( I ) = CONE
132 END IF
133 30 CONTINUE
134 KASE = 2
135 ISAVE( 1 ) = 2
136 RETURN
137 *
138 * ................ ENTRY (ISAVE( 1 ) = 2)
139 * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
140 *
141 40 CONTINUE
142 ISAVE( 2 ) = IZMAX1( N, X, 1 )
143 ISAVE( 3 ) = 2
144 *
145 * MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
146 *
147 50 CONTINUE
148 DO 60 I = 1, N
149 X( I ) = CZERO
150 60 CONTINUE
151 X( ISAVE( 2 ) ) = CONE
152 KASE = 1
153 ISAVE( 1 ) = 3
154 RETURN
155 *
156 * ................ ENTRY (ISAVE( 1 ) = 3)
157 * X HAS BEEN OVERWRITTEN BY A*X.
158 *
159 70 CONTINUE
160 CALL ZCOPY( N, X, 1, V, 1 )
161 ESTOLD = EST
162 EST = DZSUM1( N, V, 1 )
163 *
164 * TEST FOR CYCLING.
165 IF( EST.LE.ESTOLD )
166 $ GO TO 100
167 *
168 DO 80 I = 1, N
169 ABSXI = ABS( X( I ) )
170 IF( ABSXI.GT.SAFMIN ) THEN
171 X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
172 $ DIMAG( X( I ) ) / ABSXI )
173 ELSE
174 X( I ) = CONE
175 END IF
176 80 CONTINUE
177 KASE = 2
178 ISAVE( 1 ) = 4
179 RETURN
180 *
181 * ................ ENTRY (ISAVE( 1 ) = 4)
182 * X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
183 *
184 90 CONTINUE
185 JLAST = ISAVE( 2 )
186 ISAVE( 2 ) = IZMAX1( N, X, 1 )
187 IF( ( ABS( X( JLAST ) ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND.
188 $ ( ISAVE( 3 ).LT.ITMAX ) ) THEN
189 ISAVE( 3 ) = ISAVE( 3 ) + 1
190 GO TO 50
191 END IF
192 *
193 * ITERATION COMPLETE. FINAL STAGE.
194 *
195 100 CONTINUE
196 ALTSGN = ONE
197 DO 110 I = 1, N
198 X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
199 ALTSGN = -ALTSGN
200 110 CONTINUE
201 KASE = 1
202 ISAVE( 1 ) = 5
203 RETURN
204 *
205 * ................ ENTRY (ISAVE( 1 ) = 5)
206 * X HAS BEEN OVERWRITTEN BY A*X.
207 *
208 120 CONTINUE
209 TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
210 IF( TEMP.GT.EST ) THEN
211 CALL ZCOPY( N, X, 1, V, 1 )
212 EST = TEMP
213 END IF
214 *
215 130 CONTINUE
216 KASE = 0
217 RETURN
218 *
219 * End of ZLACN2
220 *
221 END