summaryrefslogtreecommitdiff
path: root/libm/double/gels.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
committerEric Andersen <andersen@codepoet.org>2001-05-10 00:40:28 +0000
commit1077fa4d772832f77a677ce7fb7c2d513b959e3f (patch)
tree579bee13fb0b58d2800206366ec2caecbb15f3fc /libm/double/gels.c
parent22358dd7ce7bb49792204b698f01a6f69b9c8e08 (diff)
uClibc now has a math library. muahahahaha!
-Erik
Diffstat (limited to 'libm/double/gels.c')
-rw-r--r--libm/double/gels.c232
1 files changed, 232 insertions, 0 deletions
diff --git a/libm/double/gels.c b/libm/double/gels.c
new file mode 100644
index 000000000..4d548d050
--- /dev/null
+++ b/libm/double/gels.c
@@ -0,0 +1,232 @@
+/*
+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 );
+}