uClibc now has a math library. muahahahaha!

 -Erik

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