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

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 );
-}