summaryrefslogtreecommitdiff
path: root/libm/double/gels.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/double/gels.c')
-rw-r--r--libm/double/gels.c232
1 files changed, 0 insertions, 232 deletions
diff --git a/libm/double/gels.c b/libm/double/gels.c
deleted file mode 100644
index 4d548d050..000000000
--- a/libm/double/gels.c
+++ /dev/null
@@ -1,232 +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 <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 );
-}