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
|
/*
* 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>
#include "protos.h"
int gels(A, R, M, EPS, AUX)
double A[], R[];
int M;
double EPS;
double AUX[];
{
int I = 0, J = 0, 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);
}
|