← Back to branch summary

~ubuntu-branches/ubuntu/edgy/k3d/edgy-proposed

« back to all changes in this revision

Viewing changes to modules/mesh/superlu/src/dlacon.c

  • Committer: Bazaar Package Importer
  • Author(s): David Martínez Moreno
  • Date: 2005-04-28 18:38:10 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050428183810-kvc6nz46z6gheo0k
Tags: 0.4.5.0-5
* AAAAAAAARGH, *patching* configure.ac broke again the build process! Fixed
  (I hope).
* Removed more cruft of shaderpreprocessor/ from configure.

Show diffs side-by-side

added added

removed removed

Lines of Context:
modules/mesh/superlu/src/dlacon.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
dlacon_(int *n, double *v, double *x, int *isgn, double *est, int *kase)
 
15
 
 
16
{
 
17
/*
 
18
    Purpose   
 
19
    =======   
 
20
 
 
21
    DLACON 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) DOUBLE 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) DOUBLE 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 DLACON must be re-called with all the other parameters   
 
40
           unchanged.   
 
41
 
 
42
    ISGN   (workspace) INT array, dimension (N)
 
43
 
 
44
    EST    (output) DOUBLE PRECISION   
 
45
           An estimate (a lower bound) for norm(A).   
 
46
 
 
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.   
 
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
    double      zero = 0.0;
 
68
    double      one = 1.0;
 
69
    
 
70
    /* Local variables */
 
71
    static int iter;
 
72
    static int jump, jlast;
 
73
    static double altsgn, estold;
 
74
    static int i, j;
 
75
    double temp;
 
76
#ifdef _CRAY
 
77
    extern int ISAMAX(int *, double *, int *);
 
78
    extern double SASUM(int *, double *, int *);
 
79
    extern int SCOPY(int *, double *, int *, double *, int *);
 
80
#else
 
81
    extern int idamax_(int *, double *, int *);
 
82
    extern double dasum_(int *, double *, int *);
 
83
    extern int dcopy_(int *, double *, int *, double *, 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. / (double) (*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 = dasum_(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 = idamax_(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
    dcopy_(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 = dasum_(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 = idamax_(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 * ((double)(i - 1) / (double)(*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) / (double)(*n * 3) * 2.;
 
213
#else
 
214
    temp = dasum_(n, x, &c__1) / (double)(*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
        dcopy_(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
} /* dlacon_ */