← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/zlacon.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/zlacon.c
 
1
 
 
2
 
 
3
/*
 
4
 * -- SuperLU routine (version 2.0) --
 
5
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
6
 * and Lawrence Berkeley National Lab.
 
7
 * November 15, 1997
 
8
 *
 
9
 */
 
10
#include <math.h>
 
11
#include "Cnames.h"
 
12
#include "dcomplex.h"
 
13
 
 
14
int
 
15
zlacon_(int *n, doublecomplex *v, doublecomplex *x, double *est, int *kase)
 
16
 
 
17
{
 
18
/*
 
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
*/
 
65
 
 
66
    /* Table of constant values */
 
67
    int c__1 = 1;
 
68
    doublecomplex      zero = {0.0, 0.0};
 
69
    doublecomplex      one = {1.0, 0.0};
 
70
 
 
71
    /* System generated locals */
 
72
    double d__1;
 
73
    
 
74
    /* Local variables */
 
75
    static int iter;
 
76
    static int jump, jlast;
 
77
    static double altsgn, estold;
 
78
    static int i, j;
 
79
    double temp;
 
80
    double safmin;
 
81
    extern double dlamch_(char *);
 
82
    extern int izmax1_(int *, doublecomplex *, int *);
 
83
    extern double dzsum1_(int *, doublecomplex *, int *);
 
84
 
 
85
    safmin = dlamch_("Safe minimum");
 
86
    if ( *kase == 0 ) {
 
87
        for (i = 0; i < *n; ++i) {
 
88
            x[i].r = 1. / (double) (*n);
 
89
            x[i].i = 0.;
 
90
        }
 
91
        *kase = 1;
 
92
        jump = 1;
 
93
        return 0;
 
94
    }
 
95
 
 
96
    switch (jump) {
 
97
        case 1:  goto L20;
 
98
        case 2:  goto L40;
 
99
        case 3:  goto L70;
 
100
        case 4:  goto L110;
 
101
        case 5:  goto L140;
 
102
    }
 
103
 
 
104
    /*     ................ ENTRY   (JUMP = 1)   
 
105
           FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
 
106
  L20:
 
107
    if (*n == 1) {
 
108
        v[0] = x[0];
 
109
        *est = z_abs(&v[0]);
 
110
        /*        ... QUIT */
 
111
        goto L150;
 
112
    }
 
113
    *est = dzsum1_(n, x, &c__1);
 
114
 
 
115
    for (i = 0; i < *n; ++i) {
 
116
        d__1 = z_abs(&x[i]);
 
117
        if (d__1 > safmin) {
 
118
            d__1 = 1 / d__1;
 
119
            x[i].r *= d__1;
 
120
            x[i].i *= d__1;
 
121
        } else {
 
122
            x[i] = one;
 
123
        }
 
124
    }
 
125
    *kase = 2;
 
126
    jump = 2;
 
127
    return 0;
 
128
 
 
129
    /*     ................ ENTRY   (JUMP = 2)   
 
130
           FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. */
 
131
L40:
 
132
    j = izmax1_(n, &x[0], &c__1);
 
133
    --j;
 
134
    iter = 2;
 
135
 
 
136
    /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
 
137
L50:
 
138
    for (i = 0; i < *n; ++i) x[i] = zero;
 
139
    x[j] = one;
 
140
    *kase = 1;
 
141
    jump = 3;
 
142
    return 0;
 
143
 
 
144
    /*     ................ ENTRY   (JUMP = 3)   
 
145
           X HAS BEEN OVERWRITTEN BY A*X. */
 
146
L70:
 
147
#ifdef _CRAY
 
148
    CCOPY(n, x, &c__1, v, &c__1);
 
149
#else
 
150
    zcopy_(n, x, &c__1, v, &c__1);
 
151
#endif
 
152
    estold = *est;
 
153
    *est = dzsum1_(n, v, &c__1);
 
154
 
 
155
 
 
156
L90:
 
157
    /*     TEST FOR CYCLING. */
 
158
    if (*est <= estold) goto L120;
 
159
 
 
160
    for (i = 0; i < *n; ++i) {
 
161
        d__1 = z_abs(&x[i]);
 
162
        if (d__1 > safmin) {
 
163
            d__1 = 1 / d__1;
 
164
            x[i].r *= d__1;
 
165
            x[i].i *= d__1;
 
166
        } else {
 
167
            x[i] = one;
 
168
        }
 
169
    }
 
170
    *kase = 2;
 
171
    jump = 4;
 
172
    return 0;
 
173
 
 
174
    /*     ................ ENTRY   (JUMP = 4)   
 
175
           X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
 
176
L110:
 
177
    jlast = j;
 
178
    j = izmax1_(n, &x[0], &c__1);
 
179
    --j;
 
180
    if (x[jlast].r != (d__1 = x[j].r, fabs(d__1)) && iter < 5) {
 
181
        ++iter;
 
182
        goto L50;
 
183
    }
 
184
 
 
185
    /*     ITERATION COMPLETE.  FINAL STAGE. */
 
186
L120:
 
187
    altsgn = 1.;
 
188
    for (i = 1; i <= *n; ++i) {
 
189
        x[i-1].r = altsgn * ((double)(i - 1) / (double)(*n - 1) + 1.);
 
190
        x[i-1].i = 0.;
 
191
        altsgn = -altsgn;
 
192
    }
 
193
    *kase = 1;
 
194
    jump = 5;
 
195
    return 0;
 
196
    
 
197
    /*     ................ ENTRY   (JUMP = 5)   
 
198
           X HAS BEEN OVERWRITTEN BY A*X. */
 
199
L140:
 
200
    temp = dzsum1_(n, x, &c__1) / (double)(*n * 3) * 2.;
 
201
    if (temp > *est) {
 
202
#ifdef _CRAY
 
203
        CCOPY(n, &x[0], &c__1, &v[0], &c__1);
 
204
#else
 
205
        zcopy_(n, &x[0], &c__1, &v[0], &c__1);
 
206
#endif
 
207
        *est = temp;
 
208
    }
 
209
 
 
210
L150:
 
211
    *kase = 0;
 
212
    return 0;
 
213
 
 
214
} /* zlacon_ */