diff options
Diffstat (limited to 'libm/double/gels.c')
| -rw-r--r-- | libm/double/gels.c | 232 |
1 files changed, 0 insertions, 232 deletions
diff --git a/libm/double/gels.c b/libm/double/gels.c deleted file mode 100644 index 4d548d050..000000000 --- a/libm/double/gels.c +++ /dev/null @@ -1,232 +0,0 @@ -/* -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 ); -} |