← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/special/cephes/gels.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:
scipy/special/cephes/gels.c
 
1
/*
 
2
C
 
3
C     ..................................................................
 
4
C
 
5
C        SUBROUTINE GELS
 
6
C
 
7
C        PURPOSE
 
8
C           TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
 
9
C           SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
 
10
C           IS ASSUMED TO BE STORED COLUMNWISE.
 
11
C
 
12
C        USAGE
 
13
C           CALL GELS(R,A,M,N,EPS,IER,AUX)
 
14
C
 
15
C        DESCRIPTION OF PARAMETERS
 
16
C           R      - M BY N RIGHT HAND SIDE MATRIX.  (DESTROYED)
 
17
C                    ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
 
18
C           A      - UPPER TRIANGULAR PART OF THE SYMMETRIC
 
19
C                    M BY M COEFFICIENT MATRIX.  (DESTROYED)
 
20
C           M      - THE NUMBER OF EQUATIONS IN THE SYSTEM.
 
21
C           N      - THE NUMBER OF RIGHT HAND SIDE VECTORS.
 
22
C           EPS    - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
 
23
C                    TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
 
24
C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
 
25
C                    IER=0  - NO ERROR,
 
26
C                    IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
 
27
C                             PIVOT ELEMENT AT ANY ELIMINATION STEP
 
28
C                             EQUAL TO 0,
 
29
C                    IER=K  - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
 
30
C                             CANCE INDICATED AT ELIMINATION STEP K+1,
 
31
C                             WHERE PIVOT ELEMENT WAS LESS THAN OR
 
32
C                             EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
 
33
C                             ABSOLUTELY GREATEST MAIN DIAGONAL
 
34
C                             ELEMENT OF MATRIX A.
 
35
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH DIMENSION M-1.
 
36
C
 
37
C        REMARKS
 
38
C           UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
 
39
C           COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
 
40
C           HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
 
41
C           LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
 
42
C           TOO.
 
43
C           THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
 
44
C           GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
 
45
C           ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
 
46
C           INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
 
47
C           SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
 
48
C           INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
 
49
C           GIVEN IN CASE M=1.
 
50
C           ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
 
51
C           MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
 
52
C           ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE GELG (WHICH
 
53
C           WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
 
54
C
 
55
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
 
56
C           NONE
 
57
C
 
58
C        METHOD
 
59
C           SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
 
60
C           PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
 
61
C           SYMMETRY IN REMAINING COEFFICIENT MATRICES.
 
62
C
 
63
C     ..................................................................
 
64
C
 
65
*/
 
66
#define ANSIPROT
 
67
#ifndef ANSIPROT
 
68
double fabs();
 
69
#else
 
70
extern double fabs(double);
 
71
int gels( double [], double [], int, double, double [] );
 
72
#endif
 
73
 
 
74
int
 
75
gels( A, R, M, EPS, AUX )
 
76
double A[],R[];
 
77
int M;
 
78
double EPS;
 
79
double AUX[];
 
80
{
 
81
int I = 0, J = 0, K, L, IER;
 
82
int II, LL, LLD, LR, LT, LST, LLST, LEND;
 
83
double tb, piv, tol, pivi;
 
84
 
 
85
if( M <= 0 )
 
86
        {
 
87
fatal:
 
88
        IER = -1;
 
89
        goto done;
 
90
        }
 
91
/* SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT */
 
92
 
 
93
/*  Diagonal elements are at A(i,i) = 1, 3, 6, 10, ...
 
94
 *  A(i,j) = A( i(i-1)/2 + j )
 
95
 */
 
96
IER = 0;
 
97
piv = 0.0;
 
98
L = 0;
 
99
for( K=1; K<=M; K++ )
 
100
        {
 
101
        L += K;
 
102
        tb = fabs( A[L-1] );
 
103
        if( tb > piv )
 
104
                {
 
105
                piv = tb;
 
106
                I = L;
 
107
                J = K;
 
108
                }
 
109
        }
 
110
tol = EPS * piv;
 
111
 
 
112
/*
 
113
C     MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
 
114
C     PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
 
115
*/
 
116
 
 
117
/*     START ELIMINATION LOOP */
 
118
LST = 0;
 
119
LEND = M - 1;
 
120
for( K=1; K<=M; K++ )
 
121
        {
 
122
/*     TEST ON USEFULNESS OF SYMMETRIC ALGORITHM */
 
123
        if( piv <= 0.0 )
 
124
                goto fatal;
 
125
        if( IER == 0 )
 
126
                {
 
127
                if( piv <= tol )
 
128
                        {
 
129
                        IER = K - 1;
 
130
                        }
 
131
                }
 
132
        LT = J - K;
 
133
        LST += K;
 
134
 
 
135
/*  PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R */
 
136
        pivi = 1.0 / A[I-1];
 
137
        L = K;
 
138
        LL = L + LT;
 
139
        tb = pivi * R[LL-1];
 
140
        R[LL-1] = R[L-1];
 
141
        R[L-1] = tb;
 
142
/* IS ELIMINATION TERMINATED */
 
143
        if( K >= M )
 
144
                break;
 
145
/*
 
146
C     ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
 
147
C     ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
 
148
*/
 
149
        LR = LST + (LT*(K+J-1))/2;
 
150
        LL = LR;
 
151
        L=LST;
 
152
        for( II=K; II<=LEND; II++ )
 
153
                {
 
154
                L += II;
 
155
                LL += 1;
 
156
                if( L == LR )
 
157
                        {
 
158
                        A[LL-1] = A[LST-1];
 
159
                        tb = A[L-1];
 
160
                        goto lab13;
 
161
                        }
 
162
                if( L > LR )
 
163
                        LL = L + LT;
 
164
 
 
165
                tb = A[LL-1];
 
166
                A[LL-1] = A[L-1];
 
167
lab13:
 
168
                AUX[II-1] = tb;
 
169
                A[L-1] = pivi * tb;
 
170
                }
 
171
/* SAVE COLUMN INTERCHANGE INFORMATION */
 
172
        A[LST-1] = LT;
 
173
/* ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT */
 
174
        piv = 0.0;
 
175
        LLST = LST;
 
176
        LT = 0;
 
177
        for( II=K; II<=LEND; II++ )
 
178
                {
 
179
                pivi = -AUX[II-1];
 
180
                LL = LLST;
 
181
                LT += 1;
 
182
                for( LLD=II; LLD<=LEND; LLD++ )
 
183
                        {
 
184
                        LL += LLD;
 
185
                        L = LL + LT;
 
186
                        A[L-1] += pivi * A[LL-1];
 
187
                        }
 
188
                LLST += II;
 
189
                LR = LLST + LT;
 
190
                tb =fabs( A[LR-1] );
 
191
                if( tb > piv )
 
192
                        {
 
193
                        piv = tb;
 
194
                        I = LR;
 
195
                        J = II + 1;
 
196
                        }
 
197
                LR = K;
 
198
                LL = LR + LT;
 
199
                R[LL-1] += pivi * R[LR-1];
 
200
                }
 
201
        }
 
202
/* END OF ELIMINATION LOOP */
 
203
 
 
204
/* BACK SUBSTITUTION AND BACK INTERCHANGE */
 
205
 
 
206
if( LEND <= 0 )
 
207
        {
 
208
        if( LEND < 0 )
 
209
                goto fatal;
 
210
        goto done;
 
211
        }
 
212
II = M;
 
213
for( I=2; I<=M; I++ )
 
214
        {
 
215
        LST -= II;
 
216
        II -= 1;
 
217
        L = A[LST-1] + 0.5;
 
218
        J = II;
 
219
        tb = R[J-1];
 
220
        LL = J;
 
221
        K = LST;
 
222
        for( LT=II; LT<=LEND; LT++ )
 
223
                {
 
224
                LL += 1;
 
225
                K += LT;
 
226
                tb -= A[K-1] * R[LL-1];
 
227
                }
 
228
        K = J + L;
 
229
        R[J-1] = R[K-1];
 
230
        R[K-1] = tb;
 
231
        }
 
232
done:
 
233
return( IER );
 
234
}