C C .................................................................. C C SUBROUTINE GELB C C PURPOSE C TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH A C COEFFICIENT MATRIX OF BAND STRUCTURE. C C USAGE C CALL GELB(R,A,M,N,MUD,MLD,EPS,IER) 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 - M BY M COEFFICIENT MATRIX WITH BAND STRUCTURE C (DESTROYED). C M - THE NUMBER OF EQUATIONS IN THE SYSTEM. C N - THE NUMBER OF RIGHT HAND SIDE VECTORS. C MUD - THE NUMBER OF UPPER CODIAGONALS (THAT MEANS C CODIAGONALS ABOVE MAIN DIAGONAL). C MLD - THE NUMBER OF LOWER CODIAGONALS (THAT MEANS C CODIAGONALS BELOW MAIN DIAGONAL). 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 WRONG INPUT PARAME- C TERS M,MUD,MLD OR BECAUSE OF PIVOT ELEMENT C AT ANY ELIMINATION STEP 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 ELEMENT OF MATRIX A. C C REMARKS C BAND MATRIX A IS ASSUMED TO BE STORED ROWWISE IN THE FIRST C ME SUCCESSIVE STORAGE LOCATIONS OF TOTALLY NEEDED MA C STORAGE LOCATIONS, WHERE C MA=M*MC-ML*(ML+1)/2 AND ME=MA-MU*(MU+1)/2 WITH C MC=MIN(M,1+MUD+MLD), ML=MC-1-MLD, MU=MC-1-MUD. C RIGHT HAND SIDE MATRIX R IS ASSUMED TO BE STORED COLUMNWISE C IN N*M SUCCESSIVE STORAGE LOCATIONS. ON RETURN SOLUTION C MATRIX R IS STORED COLUMNWISE TOO. C INPUT PARAMETERS M, MUD, MLD SHOULD SATISFY THE FOLLOWING C RESTRICTIONS MUD NOT LESS THAN ZERO C MLD NOT LESS THAN ZERO C MUD+MLD NOT GREATER THAN 2*M-2. C NO ACTION BESIDES ERROR MESSAGE IER=-1 TAKES PLACE IF THESE C RESTRICTIONS ARE NOT SATISFIED. C THE PROCEDURE GIVES RESULTS IF THE RESTRICTIONS ON INPUT C PARAMETERS ARE SATISFIED AND IF PIVOT ELEMENTS AT ALL C ELIMINATION STEPS ARE DIFFERENT FROM 0. HOWEVER WARNING C IER=K - IF GIVEN - INDICATES POSSIBLE LOSS OF SIGNIFICANCE. C IN CASE OF A WELL SCALED MATRIX A AND APPROPRIATE TOLERANCE C EPS, IER=K MAY BE INTERPRETED THAT MATRIX A HAS THE RANK K. C NO WARNING IS GIVEN IF MATRIX A HAS NO LOWER CODIAGONAL. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SOLUTION IS DONE BY MEANS OF GAUSS ELIMINATION WITH C COLUMN PIVOTING ONLY, IN ORDER TO PRESERVE BAND STRUCTURE C IN REMAINING COEFFICIENT MATRICES. C C .................................................................. C SUBROUTINE GELB(R,A,M,N,MUD,MLD,EPS,IER) C C DIMENSION R(1),A(1) C C TEST ON WRONG INPUT PARAMETERS IF(MLD)47,1,1 1 IF(MUD)47,2,2 2 MC=1+MLD+MUD IF(MC+1-M-M)3,3,47 C C PREPARE INTEGER PARAMETERS C MC=NUMBER OF COLUMNS IN MATRIX A C MU=NUMBER OF ZEROS TO BE INSERTED IN FIRST ROW OF MATRIX A C ML=NUMBER OF MISSING ELEMENTS IN LAST ROW OF MATRIX A C MR=INDEX OF LAST ROW IN MATRIX A WITH MC ELEMENTS C MZ=TOTAL NUMBER OF ZEROS TO BE INSERTED IN MATRIX A C MA=TOTAL NUMBER OF STORAGE LOCATIONS NECESSARY FOR MATRIX A C NM=NUMBER OF ELEMENTS IN MATRIX R 3 IF(MC-M)5,5,4 4 MC=M 5 MU=MC-MUD-1 ML=MC-MLD-1 MR=M-ML MZ=(MU*(MU+1))/2 MA=M*MC-(ML*(ML+1))/2 NM=N*M C C MOVE ELEMENTS BACKWARD AND SEARCH FOR ABSOLUTELY GREATEST ELEMENT C (NOT NECESSARY IN CASE OF A MATRIX WITHOUT LOWER CODIAGONALS) IER=0 PIV=0. IF(MLD)14,14,6 6 JJ=MA J=MA-MZ KST=J DO 9 K=1,KST TB=A(J) A(JJ)=TB TB=ABS(TB) IF(TB-PIV)8,8,7 7 PIV=TB 8 J=J-1 9 JJ=JJ-1 C C INSERT ZEROS IN FIRST MU ROWS (NOT NECESSARY IN CASE MZ=0) IF(MZ)14,14,10 10 JJ=1 J=1+MZ IC=1+MUD DO 13 I=1,MU DO 12 K=1,MC A(JJ)=0. IF(K-IC)11,11,12 11 A(JJ)=A(J) J=J+1 12 JJ=JJ+1 13 IC=IC+1 C C GENERATE TEST VALUE FOR SINGULARITY 14 TOL=EPS*PIV C C C START DECOMPOSITION LOOP KST=1 IDST=MC IC=MC-1 DO 38 K=1,M IF(K-MR-1)16,16,15 15 IDST=IDST-1 16 ID=IDST ILR=K+MLD IF(ILR-M)18,18,17 17 ILR=M 18 II=KST C C PIVOT SEARCH IN FIRST COLUMN (ROW INDEXES FROM I=K UP TO I=ILR) PIV=0. DO 22 I=K,ILR TB=ABS(A(II)) IF(TB-PIV)20,20,19 19 PIV=TB J=I JJ=II 20 IF(I-MR)22,22,21 21 ID=ID-1 22 II=II+ID C C TEST ON SINGULARITY IF(PIV)47,47,23 23 IF(IER)26,24,26 24 IF(PIV-TOL)25,25,26 25 IER=K-1 26 PIV=1./A(JJ) C C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R ID=J-K DO 27 I=K,NM,M II=I+ID TB=PIV*R(II) R(II)=R(I) 27 R(I)=TB C C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN COEFFICIENT MATRIX A II=KST J=JJ+IC DO 28 I=JJ,J TB=PIV*A(I) A(I)=A(II) A(II)=TB 28 II=II+1 C C ELEMENT REDUCTION IF(K-ILR)29,34,34 29 ID=KST II=K+1 MU=KST+1 MZ=KST+IC DO 33 I=II,ILR C C IN MATRIX A ID=ID+MC JJ=I-MR-1 IF(JJ)31,31,30 30 ID=ID-JJ 31 PIV=-A(ID) J=ID+1 DO 32 JJ=MU,MZ A(J-1)=A(J)+PIV*A(JJ) 32 J=J+1 A(J-1)=0. C C IN MATRIX R J=K DO 33 JJ=I,NM,M R(JJ)=R(JJ)+PIV*R(J) 33 J=J+M 34 KST=KST+MC IF(ILR-MR)36,35,35 35 IC=IC-1 36 ID=K-MR IF(ID)38,38,37 37 KST=KST-ID 38 CONTINUE C END OF DECOMPOSITION LOOP C C C BACK SUBSTITUTION IF(MC-1)46,46,39 39 IC=2 KST=MA+ML-MC+2 II=M DO 45 I=2,M KST=KST-MC II=II-1 J=II-MR IF(J)41,41,40 40 KST=KST+J 41 DO 43 J=II,NM,M TB=R(J) MZ=KST+IC-2 ID=J DO 42 JJ=KST,MZ ID=ID+1 42 TB=TB-A(JJ)*R(ID) 43 R(J)=TB IF(IC-MC)44,45,45 44 IC=IC+1 45 CONTINUE 46 RETURN C C C ERROR RETURN 47 IER=-1 RETURN END