1 /* ./src_f77/clacon.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11
clacon_(integer * n,complex * v,complex * x,real * est,integer * kase)12 /* Subroutine */ int clacon_(integer *n, complex *v, complex *x, real *est,
13 integer *kase)
14 {
15 /* System generated locals */
16 integer i__1, i__2, i__3;
17 real r__1, r__2;
18 complex q__1;
19
20 /* Builtin functions */
21 double c_abs(complex *), r_imag(complex *);
22
23 /* Local variables */
24 static integer i__, j, iter;
25 static real temp;
26 static integer jump;
27 static real absxi;
28 static integer jlast;
29 extern /* Subroutine */ int ccopy_(integer *, complex *, integer *,
30 complex *, integer *);
31 extern integer icmax1_(integer *, complex *, integer *);
32 extern doublereal scsum1_(integer *, complex *, integer *), slamch_(char *
33 , ftnlen);
34 static real safmin, altsgn, estold;
35
36
37 /* -- LAPACK auxiliary routine (version 3.0) -- */
38 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
39 /* Courant Institute, Argonne National Lab, and Rice University */
40 /* June 30, 1999 */
41
42 /* .. Scalar Arguments .. */
43 /* .. */
44 /* .. Array Arguments .. */
45 /* .. */
46
47 /* Purpose */
48 /* ======= */
49
50 /* CLACON estimates the 1-norm of a square, complex matrix A. */
51 /* Reverse communication is used for evaluating matrix-vector products. */
52
53 /* Arguments */
54 /* ========= */
55
56 /* N (input) INTEGER */
57 /* The order of the matrix. N >= 1. */
58
59 /* V (workspace) COMPLEX array, dimension (N) */
60 /* On the final return, V = A*W, where EST = norm(V)/norm(W) */
61 /* (W is not returned). */
62
63 /* X (input/output) COMPLEX array, dimension (N) */
64 /* On an intermediate return, X should be overwritten by */
65 /* A * X, if KASE=1, */
66 /* A' * X, if KASE=2, */
67 /* where A' is the conjugate transpose of A, and CLACON must be */
68 /* re-called with all the other parameters unchanged. */
69
70 /* EST (output) REAL */
71 /* An estimate (a lower bound) for norm(A). */
72
73 /* KASE (input/output) INTEGER */
74 /* On the initial call to CLACON, KASE should be 0. */
75 /* On an intermediate return, KASE will be 1 or 2, indicating */
76 /* whether X should be overwritten by A * X or A' * X. */
77 /* On the final return from CLACON, KASE will again be 0. */
78
79 /* Further Details */
80 /* ======= ======= */
81
82 /* Contributed by Nick Higham, University of Manchester. */
83 /* Originally named CONEST, dated March 16, 1988. */
84
85 /* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of */
86 /* a real or complex matrix, with applications to condition estimation", */
87 /* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. */
88
89 /* Last modified: April, 1999 */
90
91 /* ===================================================================== */
92
93 /* .. Parameters .. */
94 /* .. */
95 /* .. Local Scalars .. */
96 /* .. */
97 /* .. External Functions .. */
98 /* .. */
99 /* .. External Subroutines .. */
100 /* .. */
101 /* .. Intrinsic Functions .. */
102 /* .. */
103 /* .. Save statement .. */
104 /* .. */
105 /* .. Executable Statements .. */
106
107 /* Parameter adjustments */
108 --x;
109 --v;
110
111 /* Function Body */
112 safmin = slamch_("Safe minimum", (ftnlen)12);
113 if (*kase == 0) {
114 i__1 = *n;
115 for (i__ = 1; i__ <= i__1; ++i__) {
116 i__2 = i__;
117 r__1 = 1.f / (real) (*n);
118 q__1.r = r__1, q__1.i = 0.f;
119 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
120 /* L10: */
121 }
122 *kase = 1;
123 jump = 1;
124 return 0;
125 }
126
127 switch (jump) {
128 case 1: goto L20;
129 case 2: goto L40;
130 case 3: goto L70;
131 case 4: goto L90;
132 case 5: goto L120;
133 }
134
135 /* ................ ENTRY (JUMP = 1) */
136 /* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
137
138 L20:
139 if (*n == 1) {
140 v[1].r = x[1].r, v[1].i = x[1].i;
141 *est = c_abs(&v[1]);
142 /* ... QUIT */
143 goto L130;
144 }
145 *est = scsum1_(n, &x[1], &c__1);
146
147 i__1 = *n;
148 for (i__ = 1; i__ <= i__1; ++i__) {
149 absxi = c_abs(&x[i__]);
150 if (absxi > safmin) {
151 i__2 = i__;
152 i__3 = i__;
153 r__1 = x[i__3].r / absxi;
154 r__2 = r_imag(&x[i__]) / absxi;
155 q__1.r = r__1, q__1.i = r__2;
156 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
157 } else {
158 i__2 = i__;
159 x[i__2].r = 1.f, x[i__2].i = 0.f;
160 }
161 /* L30: */
162 }
163 *kase = 2;
164 jump = 2;
165 return 0;
166
167 /* ................ ENTRY (JUMP = 2) */
168 /* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. */
169
170 L40:
171 j = icmax1_(n, &x[1], &c__1);
172 iter = 2;
173
174 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
175
176 L50:
177 i__1 = *n;
178 for (i__ = 1; i__ <= i__1; ++i__) {
179 i__2 = i__;
180 x[i__2].r = 0.f, x[i__2].i = 0.f;
181 /* L60: */
182 }
183 i__1 = j;
184 x[i__1].r = 1.f, x[i__1].i = 0.f;
185 *kase = 1;
186 jump = 3;
187 return 0;
188
189 /* ................ ENTRY (JUMP = 3) */
190 /* X HAS BEEN OVERWRITTEN BY A*X. */
191
192 L70:
193 ccopy_(n, &x[1], &c__1, &v[1], &c__1);
194 estold = *est;
195 *est = scsum1_(n, &v[1], &c__1);
196
197 /* TEST FOR CYCLING. */
198 if (*est <= estold) {
199 goto L100;
200 }
201
202 i__1 = *n;
203 for (i__ = 1; i__ <= i__1; ++i__) {
204 absxi = c_abs(&x[i__]);
205 if (absxi > safmin) {
206 i__2 = i__;
207 i__3 = i__;
208 r__1 = x[i__3].r / absxi;
209 r__2 = r_imag(&x[i__]) / absxi;
210 q__1.r = r__1, q__1.i = r__2;
211 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
212 } else {
213 i__2 = i__;
214 x[i__2].r = 1.f, x[i__2].i = 0.f;
215 }
216 /* L80: */
217 }
218 *kase = 2;
219 jump = 4;
220 return 0;
221
222 /* ................ ENTRY (JUMP = 4) */
223 /* X HAS BEEN OVERWRITTEN BY CTRANS(A)*X. */
224
225 L90:
226 jlast = j;
227 j = icmax1_(n, &x[1], &c__1);
228 if (c_abs(&x[jlast]) != c_abs(&x[j]) && iter < 5) {
229 ++iter;
230 goto L50;
231 }
232
233 /* ITERATION COMPLETE. FINAL STAGE. */
234
235 L100:
236 altsgn = 1.f;
237 i__1 = *n;
238 for (i__ = 1; i__ <= i__1; ++i__) {
239 i__2 = i__;
240 r__1 = altsgn * ((real) (i__ - 1) / (real) (*n - 1) + 1.f);
241 q__1.r = r__1, q__1.i = 0.f;
242 x[i__2].r = q__1.r, x[i__2].i = q__1.i;
243 altsgn = -altsgn;
244 /* L110: */
245 }
246 *kase = 1;
247 jump = 5;
248 return 0;
249
250 /* ................ ENTRY (JUMP = 5) */
251 /* X HAS BEEN OVERWRITTEN BY A*X. */
252
253 L120:
254 temp = scsum1_(n, &x[1], &c__1) / (real) (*n * 3) * 2.f;
255 if (temp > *est) {
256 ccopy_(n, &x[1], &c__1, &v[1], &c__1);
257 *est = temp;
258 }
259
260 L130:
261 *kase = 0;
262 return 0;
263
264 /* End of CLACON */
265
266 } /* clacon_ */
267
268