C SUBROUTINE 'GELS' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE. C C NOTE: TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS C (1) HAVE BEEN CHANGED TO (*). 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 SUBROUTINE GELS(R,A,M,N,EPS,IER,AUX) C C DIMENSION A(*),R(*),AUX(*) IF(M)24,24,1 C C SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT 1 IER=0 PIV=0. L=0 DO 3 K=1,M L=L+K TB=ABS(A(L)) IF(TB-PIV)3,3,2 2 PIV=TB I=L J=K 3 CONTINUE 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). C C C START ELIMINATION LOOP LST=0 NM=N*M LEND=M-1 DO 18 K=1,M C C TEST ON USEFULNESS OF SYMMETRIC ALGORITHM IF(PIV)24,24,4 4 IF(IER)7,5,7 5 IF(PIV-TOL)6,6,7 6 IER=K-1 7 LT=J-K LST=LST+K C C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R PIVI=1./A(I) DO 8 L=K,NM,M LL=L+LT TB=PIVI*R(LL) R(LL)=R(L) 8 R(L)=TB C C IS ELIMINATION TERMINATED IF(K-M)9,19,19 C C ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A. C ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX. 9 LR=LST+(LT*(K+J-1))/2 LL=LR L=LST DO 14 II=K,LEND L=L+II LL=LL+1 IF(L-LR)12,10,11 10 A(LL)=A(LST) TB=A(L) GO TO 13 11 LL=L+LT 12 TB=A(LL) A(LL)=A(L) 13 AUX(II)=TB 14 A(L)=PIVI*TB C C SAVE COLUMN INTERCHANGE INFORMATION A(LST)=LT C C ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT PIV=0. LLST=LST LT=0 DO 18 II=K,LEND PIVI=-AUX(II) LL=LLST LT=LT+1 DO 15 LLD=II,LEND LL=LL+LLD L=LL+LT 15 A(L)=A(L)+PIVI*A(LL) LLST=LLST+II LR=LLST+LT TB=ABS(A(LR)) IF(TB-PIV)17,17,16 16 PIV=TB I=LR J=II+1 17 DO 18 LR=K,NM,M LL=LR+LT 18 R(LL)=R(LL)+PIVI*R(LR) C END OF ELIMINATION LOOP C C C BACK SUBSTITUTION AND BACK INTERCHANGE 19 IF(LEND)24,23,20 20 II=M DO 22 I=2,M LST=LST-II II=II-1 L=A(LST)+.5 DO 22 J=II,NM,M TB=R(J) LL=J K=LST DO 21 LT=II,LEND LL=LL+1 K=K+LT 21 TB=TB-A(K)*R(LL) K=J+L R(J)=R(K) 22 R(K)=TB 23 RETURN C C C ERROR RETURN 24 IER=-1 RETURN END C C======================================================================= C