← Back to branch summary

~ubuntu-branches/ubuntu/raring/python-scipy/raring-proposed

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/slacon.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
 
 
13
 
int
14
 
slacon_(int *n, float *v, float *x, int *isgn, float *est, int *kase)
15
 
 
16
 
{
17
 
/*
18
 
    Purpose   
19
 
    =======   
20
 
 
21
 
    SLACON estimates the 1-norm of a square matrix A.   
22
 
    Reverse communication is used for evaluating matrix-vector products. 
23
 
  
24
 
 
25
 
    Arguments   
26
 
    =========   
27
 
 
28
 
    N      (input) INT
29
 
           The order of the matrix.  N >= 1.   
30
 
 
31
 
    V      (workspace) FLOAT PRECISION array, dimension (N)   
32
 
           On the final return, V = A*W,  where  EST = norm(V)/norm(W)   
33
 
           (W is not returned).   
34
 
 
35
 
    X      (input/output) FLOAT PRECISION array, dimension (N)   
36
 
           On an intermediate return, X should be overwritten by   
37
 
                 A * X,   if KASE=1,   
38
 
                 A' * X,  if KASE=2,
39
 
           and SLACON must be re-called with all the other parameters   
40
 
           unchanged.   
41
 
 
42
 
    ISGN   (workspace) INT array, dimension (N)
43
 
 
44
 
    EST    (output) FLOAT PRECISION   
45
 
           An estimate (a lower bound) for norm(A).   
46
 
 
47
 
    KASE   (input/output) INT
48
 
           On the initial call to SLACON, 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 SLACON, KASE will again be 0.   
52
 
 
53
 
    Further Details   
54
 
    ======= =======   
55
 
 
56
 
    Contributed by Nick Higham, University of Manchester.   
57
 
    Originally named CONEST, dated March 16, 1988.   
58
 
 
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
 
    ===================================================================== 
63
 
*/
64
 
 
65
 
    /* Table of constant values */
66
 
    int c__1 = 1;
67
 
    float      zero = 0.0;
68
 
    float      one = 1.0;
69
 
    
70
 
    /* Local variables */
71
 
    static int iter;
72
 
    static int jump, jlast;
73
 
    static float altsgn, estold;
74
 
    static int i, j;
75
 
    float temp;
76
 
#ifdef _CRAY
77
 
    extern int ISAMAX(int *, float *, int *);
78
 
    extern float SASUM(int *, float *, int *);
79
 
    extern int SCOPY(int *, float *, int *, float *, int *);
80
 
#else
81
 
    extern int isamax_(int *, float *, int *);
82
 
    extern float sasum_(int *, float *, int *);
83
 
    extern int scopy_(int *, float *, int *, float *, int *);
84
 
#endif
85
 
#define d_sign(a, b) (b >= 0 ? fabs(a) : -fabs(a))    /* Copy sign */
86
 
#define i_dnnt(a) \
87
 
        ( a>=0 ? floor(a+.5) : -floor(.5-a) ) /* Round to nearest integer */
88
 
 
89
 
    if ( *kase == 0 ) {
90
 
        for (i = 0; i < *n; ++i) {
91
 
            x[i] = 1. / (float) (*n);
92
 
        }
93
 
        *kase = 1;
94
 
        jump = 1;
95
 
        return 0;
96
 
    }
97
 
 
98
 
    switch (jump) {
99
 
        case 1:  goto L20;
100
 
        case 2:  goto L40;
101
 
        case 3:  goto L70;
102
 
        case 4:  goto L110;
103
 
        case 5:  goto L140;
104
 
    }
105
 
 
106
 
    /*     ................ ENTRY   (JUMP = 1)   
107
 
           FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X. */
108
 
  L20:
109
 
    if (*n == 1) {
110
 
        v[0] = x[0];
111
 
        *est = fabs(v[0]);
112
 
        /*        ... QUIT */
113
 
        goto L150;
114
 
    }
115
 
#ifdef _CRAY
116
 
    *est = SASUM(n, x, &c__1);
117
 
#else
118
 
    *est = sasum_(n, x, &c__1);
119
 
#endif
120
 
 
121
 
    for (i = 0; i < *n; ++i) {
122
 
        x[i] = d_sign(one, x[i]);
123
 
        isgn[i] = i_dnnt(x[i]);
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
 
#ifdef _CRAY
133
 
    j = ISAMAX(n, &x[0], &c__1);
134
 
#else
135
 
    j = isamax_(n, &x[0], &c__1);
136
 
#endif
137
 
    --j;
138
 
    iter = 2;
139
 
 
140
 
    /*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX. */
141
 
L50:
142
 
    for (i = 0; i < *n; ++i) x[i] = zero;
143
 
    x[j] = one;
144
 
    *kase = 1;
145
 
    jump = 3;
146
 
    return 0;
147
 
 
148
 
    /*     ................ ENTRY   (JUMP = 3)   
149
 
           X HAS BEEN OVERWRITTEN BY A*X. */
150
 
L70:
151
 
#ifdef _CRAY
152
 
    SCOPY(n, x, &c__1, v, &c__1);
153
 
#else
154
 
    scopy_(n, x, &c__1, v, &c__1);
155
 
#endif
156
 
    estold = *est;
157
 
#ifdef _CRAY
158
 
    *est = SASUM(n, v, &c__1);
159
 
#else
160
 
    *est = sasum_(n, v, &c__1);
161
 
#endif
162
 
 
163
 
    for (i = 0; i < *n; ++i)
164
 
        if (i_dnnt(d_sign(one, x[i])) != isgn[i])
165
 
            goto L90;
166
 
 
167
 
    /*     REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. */
168
 
    goto L120;
169
 
 
170
 
L90:
171
 
    /*     TEST FOR CYCLING. */
172
 
    if (*est <= estold) goto L120;
173
 
 
174
 
    for (i = 0; i < *n; ++i) {
175
 
        x[i] = d_sign(one, x[i]);
176
 
        isgn[i] = i_dnnt(x[i]);
177
 
    }
178
 
    *kase = 2;
179
 
    jump = 4;
180
 
    return 0;
181
 
 
182
 
    /*     ................ ENTRY   (JUMP = 4)   
183
 
           X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X. */
184
 
L110:
185
 
    jlast = j;
186
 
#ifdef _CRAY
187
 
    j = ISAMAX(n, &x[0], &c__1);
188
 
#else
189
 
    j = isamax_(n, &x[0], &c__1);
190
 
#endif
191
 
    --j;
192
 
    if (x[jlast] != fabs(x[j]) && iter < 5) {
193
 
        ++iter;
194
 
        goto L50;
195
 
    }
196
 
 
197
 
    /*     ITERATION COMPLETE.  FINAL STAGE. */
198
 
L120:
199
 
    altsgn = 1.;
200
 
    for (i = 1; i <= *n; ++i) {
201
 
        x[i-1] = altsgn * ((float)(i - 1) / (float)(*n - 1) + 1.);
202
 
        altsgn = -altsgn;
203
 
    }
204
 
    *kase = 1;
205
 
    jump = 5;
206
 
    return 0;
207
 
    
208
 
    /*     ................ ENTRY   (JUMP = 5)   
209
 
           X HAS BEEN OVERWRITTEN BY A*X. */
210
 
L140:
211
 
#ifdef _CRAY
212
 
    temp = SASUM(n, x, &c__1) / (float)(*n * 3) * 2.;
213
 
#else
214
 
    temp = sasum_(n, x, &c__1) / (float)(*n * 3) * 2.;
215
 
#endif
216
 
    if (temp > *est) {
217
 
#ifdef _CRAY
218
 
        SCOPY(n, &x[0], &c__1, &v[0], &c__1);
219
 
#else
220
 
        scopy_(n, &x[0], &c__1, &v[0], &c__1);
221
 
#endif
222
 
        *est = temp;
223
 
    }
224
 
 
225
 
L150:
226
 
    *kase = 0;
227
 
    return 0;
228
 
 
229
 
} /* slacon_ */