summaryrefslogtreecommitdiff
path: root/libm/ldouble/gelsl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/gelsl.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/ldouble/gelsl.c')
-rw-r--r--libm/ldouble/gelsl.c240
1 files changed, 0 insertions, 240 deletions
diff --git a/libm/ldouble/gelsl.c b/libm/ldouble/gelsl.c
deleted file mode 100644
index d66ad55e9..000000000
--- a/libm/ldouble/gelsl.c
+++ /dev/null
@@ -1,240 +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 <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 );
-}