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