diff options
Diffstat (limited to 'libm/ldouble/gelsl.c')
| -rw-r--r-- | libm/ldouble/gelsl.c | 240 |
1 files changed, 240 insertions, 0 deletions
diff --git a/libm/ldouble/gelsl.c b/libm/ldouble/gelsl.c new file mode 100644 index 000000000..d66ad55e9 --- /dev/null +++ b/libm/ldouble/gelsl.c @@ -0,0 +1,240 @@ +/* +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 <stdio.h> +#define fabsl(x) ( (x) < 0.0L ? -(x) : (x) ) + +int gels( A, R, M, EPS, AUX ) +long double A[],R[]; +int M; +long double EPS; +long double AUX[]; +{ +int I, J, K, L, IER; +int II, LL, LLD, LR, LT, LST, LLST, LEND; +long double tb, piv, tol, pivi; + +IER = 0; +if( M <= 0 ) + { +fatal: + IER = -1; + goto done; + } +/* SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT */ + +/* Diagonal elements are at A(i,i) = 0, 2, 5, 9, 14, ... + * A(i,j) = A( i(i-1)/2 + j - 1 ) + */ +piv = 0.0L; +I = 0; +J = 0; +L = 0; +for( K=1; K<=M; K++ ) + { + L += K; + tb = fabsl( 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.0L ) + { + printf( "gels: piv <= 0 at K = %d\n", K ); + goto fatal; + } + if( IER == 0 ) + { + if( piv <= tol ) + { + IER = K; +/* + goto done; +*/ + } + } + LT = J - K; + LST += K; + +/* PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R */ + pivi = 1.0L / 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.0L; + 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 =fabsl( 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 ) + { + printf( "gels: LEND = %d\n", LEND ); + if( LEND < 0 ) + goto fatal; + goto done; + } +II = M; +for( I=2; I<=M; I++ ) + { + LST -= II; + II -= 1; + L = A[LST-1] + 0.5L; + 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: +if( IER ) + printf( "gels error %d!\n", IER ); +return( IER ); +} |