← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/cephes/simq.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/special/cephes/simq.c
1
 
/*                                                      simq.c
2
 
 *
3
 
 *      Solution of simultaneous linear equations AX = B
4
 
 *      by Gaussian elimination with partial pivoting
5
 
 *
6
 
 *
7
 
 *
8
 
 * SYNOPSIS:
9
 
 *
10
 
 * double A[n*n], B[n], X[n];
11
 
 * int n, flag;
12
 
 * int IPS[];
13
 
 * int simq();
14
 
 *
15
 
 * ercode = simq( A, B, X, n, flag, IPS );
16
 
 *
17
 
 *
18
 
 *
19
 
 * DESCRIPTION:
20
 
 *
21
 
 * B, X, IPS are vectors of length n.
22
 
 * A is an n x n matrix (i.e., a vector of length n*n),
23
 
 * stored row-wise: that is, A(i,j) = A[ij],
24
 
 * where ij = i*n + j, which is the transpose of the normal
25
 
 * column-wise storage.
26
 
 *
27
 
 * The contents of matrix A are destroyed.
28
 
 *
29
 
 * Set flag=0 to solve.
30
 
 * Set flag=-1 to do a new back substitution for different B vector
31
 
 * using the same A matrix previously reduced when flag=0.
32
 
 *
33
 
 * The routine returns nonzero on error; messages are printed.
34
 
 *
35
 
 *
36
 
 * ACCURACY:
37
 
 *
38
 
 * Depends on the conditioning (range of eigenvalues) of matrix A.
39
 
 *
40
 
 *
41
 
 * REFERENCE:
42
 
 *
43
 
 * Computer Solution of Linear Algebraic Systems,
44
 
 * by George E. Forsythe and Cleve B. Moler; Prentice-Hall, 1967.
45
 
 *
46
 
 */
47
 
 
48
 
/*                                                      simq    2 */
49
 
 
50
 
#include <stdio.h>
51
 
#define ANSIPROT
52
 
#ifdef ANSIPROT
53
 
int simq(double [], double [], double [], int, int, int [] );
54
 
#endif
55
 
 
56
 
#define fabs(x) ((x) < 0 ? -(x) : (x))
57
 
 
58
 
int simq( A, B, X, n, flag, IPS )
59
 
double A[], B[], X[];
60
 
int n, flag;
61
 
int IPS[];
62
 
{
63
 
int i, j, ij, ip, ipj, ipk, ipn;
64
 
int idxpiv, iback;
65
 
int k, kp, kp1, kpk, kpn;
66
 
int nip, nkp, nm1;
67
 
double em, q, rownrm, big, size, pivot, sum;
68
 
 
69
 
nm1 = n-1;
70
 
if( flag < 0 )
71
 
        goto solve;
72
 
 
73
 
/*      Initialize IPS and X    */
74
 
 
75
 
ij=0;
76
 
for( i=0; i<n; i++ )
77
 
        {
78
 
        IPS[i] = i;
79
 
        rownrm = 0.0;
80
 
        for( j=0; j<n; j++ )
81
 
                {
82
 
                q = fabs( A[ij] );
83
 
                if( rownrm < q )
84
 
                        rownrm = q;
85
 
                ++ij;
86
 
                }
87
 
        if( rownrm == 0.0 )
88
 
                {
89
 
                puts("SIMQ ROWNRM=0");
90
 
                return(1);
91
 
                }
92
 
        X[i] = 1.0/rownrm;
93
 
        }
94
 
 
95
 
/*                                                      simq    3 */
96
 
/*      Gaussian elimination with partial pivoting      */
97
 
 
98
 
for( k=0; k<nm1; k++ )
99
 
        {
100
 
        big= 0.0;
101
 
        idxpiv = 0;
102
 
        for( i=k; i<n; i++ )
103
 
                {
104
 
                ip = IPS[i];
105
 
                ipk = n*ip + k;
106
 
                size = fabs( A[ipk] ) * X[ip];
107
 
                if( size > big )
108
 
                        {
109
 
                        big = size;
110
 
                        idxpiv = i;
111
 
                        }
112
 
                }
113
 
 
114
 
        if( big == 0.0 )
115
 
                {
116
 
                puts( "SIMQ BIG=0" );
117
 
                return(2);
118
 
                }
119
 
        if( idxpiv != k )
120
 
                {
121
 
                j = IPS[k];
122
 
                IPS[k] = IPS[idxpiv];
123
 
                IPS[idxpiv] = j;
124
 
                }
125
 
        kp = IPS[k];
126
 
        kpk = n*kp + k;
127
 
        pivot = A[kpk];
128
 
        kp1 = k+1;
129
 
        for( i=kp1; i<n; i++ )
130
 
                {
131
 
                ip = IPS[i];
132
 
                ipk = n*ip + k;
133
 
                em = -A[ipk]/pivot;
134
 
                A[ipk] = -em;
135
 
                nip = n*ip;
136
 
                nkp = n*kp;
137
 
                for( j=kp1; j<n; j++ )
138
 
                        {
139
 
                        ipj = nip + j;
140
 
                        A[ipj] = A[ipj] + em * A[nkp + j];
141
 
                        }
142
 
                }
143
 
        }
144
 
kpn = n * IPS[n-1] + n - 1;     /* last element of IPS[n] th row */
145
 
if( A[kpn] == 0.0 )
146
 
        {
147
 
        puts( "SIMQ A[kpn]=0");
148
 
        return(3);
149
 
        }
150
 
 
151
 
/*                                                      simq 4 */
152
 
/*      back substitution       */
153
 
 
154
 
solve:
155
 
ip = IPS[0];
156
 
X[0] = B[ip];
157
 
for( i=1; i<n; i++ )
158
 
        {
159
 
        ip = IPS[i];
160
 
        ipj = n * ip;
161
 
        sum = 0.0;
162
 
        for( j=0; j<i; j++ )
163
 
                {
164
 
                sum += A[ipj] * X[j];
165
 
                ++ipj;
166
 
                }
167
 
        X[i] = B[ip] - sum;
168
 
        }
169
 
 
170
 
ipn = n * IPS[n-1] + n - 1;
171
 
X[n-1] = X[n-1]/A[ipn];
172
 
 
173
 
for( iback=1; iback<n; iback++ )
174
 
        {
175
 
/* i goes (n-1),...,1   */
176
 
        i = nm1 - iback;
177
 
        ip = IPS[i];
178
 
        nip = n*ip;
179
 
        sum = 0.0;
180
 
        for( j=i+1; j<n; j++ )
181
 
                sum += A[nip+j] * X[j];
182
 
        X[i] = (X[i] - sum)/A[nip+i];
183
 
        }
184
 
return(0);
185
 
}