1
2 /*! @file zlacon.c
3 * \brief Estimates the 1-norm
4 *
5 * <pre>
6 * -- SuperLU routine (version 2.0) --
7 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
8 * and Lawrence Berkeley National Lab.
9 * November 15, 1997
10 * </pre>
11 */
12 #include <math.h>
13 #include "slu_Cnames.h"
14 #include "slu_dcomplex.h"
15
16 /*! \brief
17 *
18 * <pre>
19 * Purpose
20 * =======
21 *
22 * ZLACON estimates the 1-norm of a square matrix A.
23 * Reverse communication is used for evaluating matrix-vector products.
24 *
25 *
26 * Arguments
27 * =========
28 *
29 * N (input) INT
30 * The order of the matrix. N >= 1.
31 *
32 * V (workspace) DOUBLE COMPLEX PRECISION array, dimension (N)
33 * On the final return, V = A*W, where EST = norm(V)/norm(W)
34 * (W is not returned).
35 *
36 * X (input/output) DOUBLE COMPLEX PRECISION array, dimension (N)
37 * On an intermediate return, X should be overwritten by
38 * A * X, if KASE=1,
39 * A' * X, if KASE=2,
40 * where A' is the conjugate transpose of A,
41 * and ZLACON must be re-called with all the other parameters
42 * unchanged.
43 *
44 *
45 * EST (output) DOUBLE PRECISION
46 * An estimate (a lower bound) for norm(A).
47 *
48 * KASE (input/output) INT
49 * On the initial call to ZLACON, KASE should be 0.
50 * On an intermediate return, KASE will be 1 or 2, indicating
51 * whether X should be overwritten by A * X or A' * X.
52 * On the final return from ZLACON, KASE will again be 0.
53 *
54 * Further Details
55 * ======= =======
56 *
57 * Contributed by Nick Higham, University of Manchester.
58 * Originally named CONEST, dated March 16, 1988.
59 *
60 * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
61 * a real or complex matrix, with applications to condition estimation",
62 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
63 * =====================================================================
64 * </pre>
65 */
66
67 int
zlacon_(int * n,doublecomplex * v,doublecomplex * x,double * est,int * kase)68 zlacon_(int *n, doublecomplex *v, doublecomplex *x, double *est, int *kase)
69
70 {
71
72
73 /* Table of constant values */
74 int c__1 = 1;
75 doublecomplex zero = {0.0, 0.0};
76 doublecomplex one = {1.0, 0.0};
77
78 /* System generated locals */
79 double d__1;
80
81 /* Local variables */
82 static int iter;
83 static int jump, jlast;
84 static double altsgn, estold;
85 static int i, j;
86 double temp;
87 double safmin;
88 extern double dlamch_(char *);
89 extern int izmax1_(int *, doublecomplex *, int *);
90 extern double dzsum1_(int *, doublecomplex *, int *);
91
92 safmin = dlamch_("Safe minimum");
93 if ( *kase == 0 ) {
94 for (i = 0; i < *n; ++i) {
95 x[i].r = 1. / (double) (*n);
96 x[i].i = 0.;
97 }
98 *kase = 1;
99 jump = 1;
100 return 0;
101 }
102
103 switch (jump) {
104 case 1: goto L20;
105 case 2: goto L40;
106 case 3: goto L70;
107 case 4: goto L110;
108 case 5: goto L140;
109 }
110
111 /* ................ ENTRY (JUMP = 1)
112 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
113 L20:
114 if (*n == 1) {
115 v[0] = x[0];
116 *est = z_abs(&v[0]);
117 /* ... QUIT */
118 goto L150;
119 }
120 *est = dzsum1_(n, x, &c__1);
121
122 for (i = 0; i < *n; ++i) {
123 d__1 = z_abs(&x[i]);
124 if (d__1 > safmin) {
125 d__1 = 1 / d__1;
126 x[i].r *= d__1;
127 x[i].i *= d__1;
128 } else {
129 x[i] = one;
130 }
131 }
132 *kase = 2;
133 jump = 2;
134 return 0;
135
136 /* ................ ENTRY (JUMP = 2)
137 FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
138 L40:
139 j = izmax1_(n, &x[0], &c__1);
140 --j;
141 iter = 2;
142
143 /* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
144 L50:
145 for (i = 0; i < *n; ++i) x[i] = zero;
146 x[j] = one;
147 *kase = 1;
148 jump = 3;
149 return 0;
150
151 /* ................ ENTRY (JUMP = 3)
152 X HAS BEEN OVERWRITTEN BY A*X. */
153 L70:
154 #ifdef _CRAY
155 CCOPY(n, x, &c__1, v, &c__1);
156 #else
157 zcopy_(n, x, &c__1, v, &c__1);
158 #endif
159 estold = *est;
160 *est = dzsum1_(n, v, &c__1);
161
162
163 L90:
164 /* TEST FOR CYCLING. */
165 if (*est <= estold) goto L120;
166
167 for (i = 0; i < *n; ++i) {
168 d__1 = z_abs(&x[i]);
169 if (d__1 > safmin) {
170 d__1 = 1 / d__1;
171 x[i].r *= d__1;
172 x[i].i *= d__1;
173 } else {
174 x[i] = one;
175 }
176 }
177 *kase = 2;
178 jump = 4;
179 return 0;
180
181 /* ................ ENTRY (JUMP = 4)
182 X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
183 L110:
184 jlast = j;
185 j = izmax1_(n, &x[0], &c__1);
186 --j;
187 if (x[jlast].r != (d__1 = x[j].r, fabs(d__1)) && iter < 5) {
188 ++iter;
189 goto L50;
190 }
191
192 /* ITERATION COMPLETE. FINAL STAGE. */
193 L120:
194 altsgn = 1.;
195 for (i = 1; i <= *n; ++i) {
196 x[i-1].r = altsgn * ((double)(i - 1) / (double)(*n - 1) + 1.);
197 x[i-1].i = 0.;
198 altsgn = -altsgn;
199 }
200 *kase = 1;
201 jump = 5;
202 return 0;
203
204 /* ................ ENTRY (JUMP = 5)
205 X HAS BEEN OVERWRITTEN BY A*X. */
206 L140:
207 temp = dzsum1_(n, x, &c__1) / (double)(*n * 3) * 2.;
208 if (temp > *est) {
209 #ifdef _CRAY
210 CCOPY(n, &x[0], &c__1, &v[0], &c__1);
211 #else
212 zcopy_(n, &x[0], &c__1, &v[0], &c__1);
213 #endif
214 *est = temp;
215 }
216
217 L150:
218 *kase = 0;
219 return 0;
220
221 } /* zlacon_ */
222