4
* -- SuperLU routine (version 2.0) --
5
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
6
* and Lawrence Berkeley National Lab.
14
dlacon_(int *n, double *v, double *x, int *isgn, double *est, int *kase)
21
DLACON estimates the 1-norm of a square matrix A.
22
Reverse communication is used for evaluating matrix-vector products.
29
The order of the matrix. N >= 1.
31
V (workspace) DOUBLE PRECISION array, dimension (N)
32
On the final return, V = A*W, where EST = norm(V)/norm(W)
35
X (input/output) DOUBLE PRECISION array, dimension (N)
36
On an intermediate return, X should be overwritten by
39
and DLACON must be re-called with all the other parameters
42
ISGN (workspace) INT array, dimension (N)
44
EST (output) DOUBLE PRECISION
45
An estimate (a lower bound) for norm(A).
47
KASE (input/output) INT
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' * X.
51
On the final return from DLACON, KASE will again be 0.
56
Contributed by Nick Higham, University of Manchester.
57
Originally named CONEST, dated March 16, 1988.
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
=====================================================================
65
/* Table of constant values */
72
static int jump, jlast;
73
static double altsgn, estold;
77
extern int ISAMAX(int *, double *, int *);
78
extern double SASUM(int *, double *, int *);
79
extern int SCOPY(int *, double *, int *, double *, int *);
81
extern int idamax_(int *, double *, int *);
82
extern double dasum_(int *, double *, int *);
83
extern int dcopy_(int *, double *, int *, double *, int *);
85
#define d_sign(a, b) (b >= 0 ? fabs(a) : -fabs(a)) /* Copy sign */
87
( a>=0 ? floor(a+.5) : -floor(.5-a) ) /* Round to nearest integer */
90
for (i = 0; i < *n; ++i) {
91
x[i] = 1. / (double) (*n);
106
/* ................ ENTRY (JUMP = 1)
107
FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. */
116
*est = SASUM(n, x, &c__1);
118
*est = dasum_(n, x, &c__1);
121
for (i = 0; i < *n; ++i) {
122
x[i] = d_sign(one, x[i]);
123
isgn[i] = i_dnnt(x[i]);
129
/* ................ ENTRY (JUMP = 2)
130
FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
133
j = ISAMAX(n, &x[0], &c__1);
135
j = idamax_(n, &x[0], &c__1);
140
/* MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
142
for (i = 0; i < *n; ++i) x[i] = zero;
148
/* ................ ENTRY (JUMP = 3)
149
X HAS BEEN OVERWRITTEN BY A*X. */
152
SCOPY(n, x, &c__1, v, &c__1);
154
dcopy_(n, x, &c__1, v, &c__1);
158
*est = SASUM(n, v, &c__1);
160
*est = dasum_(n, v, &c__1);
163
for (i = 0; i < *n; ++i)
164
if (i_dnnt(d_sign(one, x[i])) != isgn[i])
167
/* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
171
/* TEST FOR CYCLING. */
172
if (*est <= estold) goto L120;
174
for (i = 0; i < *n; ++i) {
175
x[i] = d_sign(one, x[i]);
176
isgn[i] = i_dnnt(x[i]);
182
/* ................ ENTRY (JUMP = 4)
183
X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
187
j = ISAMAX(n, &x[0], &c__1);
189
j = idamax_(n, &x[0], &c__1);
192
if (x[jlast] != fabs(x[j]) && iter < 5) {
197
/* ITERATION COMPLETE. FINAL STAGE. */
200
for (i = 1; i <= *n; ++i) {
201
x[i-1] = altsgn * ((double)(i - 1) / (double)(*n - 1) + 1.);
208
/* ................ ENTRY (JUMP = 5)
209
X HAS BEEN OVERWRITTEN BY A*X. */
212
temp = SASUM(n, x, &c__1) / (double)(*n * 3) * 2.;
214
temp = dasum_(n, x, &c__1) / (double)(*n * 3) * 2.;
218
SCOPY(n, &x[0], &c__1, &v[0], &c__1);
220
dcopy_(n, &x[0], &c__1, &v[0], &c__1);