← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/superlu/clacon.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/pysparse/superlu/clacon.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 "scomplex.h"
13
 
 
14
 
int
15
 
clacon_(int *n, complex *v, complex *x, float *est, int *kase)
16
 
 
17
 
{
18
 
/*
19
 
    Purpose   
20
 
    =======   
21
 
 
22
 
    CLACON 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) 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) 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 CLACON must be re-called with all the other parameters   
42
 
           unchanged.   
43
 
 
44
 
 
45
 
    EST    (output) FLOAT PRECISION   
46
 
           An estimate (a lower bound) for norm(A).   
47
 
 
48
 
    KASE   (input/output) INT
49
 
           On the initial call to CLACON, 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 CLACON, 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
 
    complex      zero = {0.0, 0.0};
69
 
    complex      one = {1.0, 0.0};
70
 
 
71
 
    /* System generated locals */
72
 
    float d__1;
73
 
    
74
 
    /* Local variables */
75
 
    static int iter;
76
 
    static int jump, jlast;
77
 
    static float altsgn, estold;
78
 
    static int i, j;
79
 
    float temp;
80
 
    float safmin;
81
 
    extern double slamch_(char *);
82
 
    extern int icmax1_(int *, complex *, int *);
83
 
    extern double scsum1_(int *, complex *, int *);
84
 
 
85
 
    safmin = slamch_("Safe minimum");
86
 
    if ( *kase == 0 ) {
87
 
        for (i = 0; i < *n; ++i) {
88
 
            x[i].r = 1. / (float) (*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 = c_abs(&v[0]);
110
 
        /*        ... QUIT */
111
 
        goto L150;
112
 
    }
113
 
    *est = scsum1_(n, x, &c__1);
114
 
 
115
 
    for (i = 0; i < *n; ++i) {
116
 
        d__1 = c_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 = icmax1_(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
 
    ccopy_(n, x, &c__1, v, &c__1);
151
 
#endif
152
 
    estold = *est;
153
 
    *est = scsum1_(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 = c_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 = icmax1_(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 * ((float)(i - 1) / (float)(*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 = scsum1_(n, x, &c__1) / (float)(*n * 3) * 2.;
201
 
    if (temp > *est) {
202
 
#ifdef _CRAY
203
 
        CCOPY(n, &x[0], &c__1, &v[0], &c__1);
204
 
#else
205
 
        ccopy_(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
 
} /* clacon_ */