|
/*
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 );
}
¤ Dauer der Verarbeitung: 0.19 Sekunden
(vorverarbeitet)
¤
|
Laden
Fehler beim Verzeichnis:
in der Quellcodebibliothek suchen
Die farbliche Syntaxdarstellung ist noch experimentell.
|
