C C Программа Lin3dsysccmSolver - обращения трехдиагональных матриц C_3 и C решения систем линейных уравнений C_3X=Y C (метод критических компонент [88]) C C Общие блоки: COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN C Обpащение: CALL Lin2dsysccmSolver(M,R,N,INF,IM,JM,B,A,DET) C Входные данные: Блок F499_RCONST содержит вещественные константы C e_8,e_0,e_1,e_2=e_1/beta, определенные в INIT_CONST C (или заранее известные); C M - (integer) поpядок квадратной матpицы C_3; C N - (integer) первая размерность массива B (N=M); C R - (real*8) одномерный (или двумеpный) массив C pазмеpности 3M. На входе в массиве R размещается C исходная матрица C_3, если она несимметричная, в виде: C C R=[*,q_1,p_2,r_2,q_2,...,p_M,r_M,q_M]. (***) C C Если она симметричная, то в виде (*). Здесь символом * C отмечены ячейки, которые при вводе не заполняются; C INF - (integer): INF=1, если задача решается с симметричной C трехдиагональной матрицей; INF=-1, если задача решается C с несимметричной трехдиагональной матpицей; C IM,JM - (integer). Если JM=0, то решается система C_3X=Y. При C этом IM - число правых частей. Если же JM=k (k<=IM), C то вычисляются элементы обратной к C_3 матрицы, C находящиеся в столбцах с номерами k,k+1,...,IM. При C этом k - номер последнего вычисляемого столбца C обратной матрицы; C B - (real*8) двумеpный массив (если отведенная память C компьютера не достаточна для хранения массива B C размерности [N,IM], то обращение к подпрограмме C Lin3dsysccmSolver организуется самим пользователем в C цикле с соответствующим выделением из B блоков C нужных размеров. В частности, B может быть одномерным, C если IM=1). Если решается система C_3X=Y, т.е. JM=0, C то массив B имеет размерность [N,IM]. При этом в B C хранятся векторы правых частей Y_j в виде (**). C Если вычисляются элементы обратной к C_3 матрицы, C т.е. JM!=0, то массив B имеет размерность [N,IM-JM+1]; C A - (real*8) рабочий массив размерности [3,M]; C DET - (real*8) если на входе DET=1, то на выходе C DET=det(C_3) - определитель матрицы C_3, если на входе C DET=0, то вычисления определителя не производится. C Выходные данные: INF - (integer): C INF=0 - ноpмальное завеpшение pаботы подпpогpамм; C INF=-1 - исходная матpица выpождена; C B - (real*8) массив содеpжит матpицу pешений и элементов C обратной матрицы в виде (**); C DET - (real*8) содержит вычисленный определитель матрицы C_3, C если на входе DET=1 C SUBROUTINE Lin3dsysccmSolver(M,R,N,INF,IM,JM,B,A,DET) IMPLICIT REAL*8(A-H, O-Z) DIMENSION R(*),B(N,*),A(N,*) COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN * call Init_Const ! Sap: to hide this actin from User IF(M.LE.0.OR.N.LE.0) write(*,100)M,N IS=1 DET1=DET Q1=2D0*EPS IN1=IABS(INF)+1 IQ=(5-ISIGN(1,INF))/2 JQ=2-IQ INF0=INF INF=-1 71 IS=-IS M2=M+1 I=M2*(1-IS)/2+IS M2=M2-I IN=IN1-IS NQ=JQ+IS MQ=JQ*(1-IS)+IS IL=0 I1=I A(I,IN)=R(IQ*I+JQ) 1 KQ=IQ*I IF(I.EQ.M2) GOTO 4 I=I+IS IF(DABS(A(I-IS,IN)).GE.UNDFLO) GOTO 2 A(I,IN)=-R(KQ+MQ)*R(KQ+NQ) DET=A(I,IN) IF(DABS(DET).LT.UNDFLO) RETURN IL=1 IF(I.EQ.M2) GOTO 69 I=I+IS A(I,IN)=R(IQ*I+JQ) GOTO 1 2 A(I-IS,IN1)=R(KQ+MQ)/A(I-IS,IN) IF(DABS(A(I-IS,IN1)).GT.1D0) IL=1 IF(DABS(R(KQ+MQ)).GT.DABS(R(KQ+NQ))) GOTO 3 A(I,IN)=R(IQ*I+JQ)-A(I-IS,IN1)*R(KQ+NQ) GOTO 1 3 A(I,IN)=R(IQ*I+JQ)-R(KQ+NQ)/A(I-IS,IN)*R(KQ+MQ) GOTO 1 4 DET=A(I,IN) IF(DABS(DET).LT.UNDFLO) A(I,IN)=EPS IF(IL.NE.0) GOTO 69 IF(JM.NE.0) GOTO 8 C SOLUTION WITH THE SWEEP METHOD DO 7 J=1,IM B(I1,J)=B(I1,J)/A(I1,IN) DO 5 L=I1+IS,I,IS 5 B(L,J)=(B(L,J)-R(IQ*(L-IS)+NQ)*B(L-IS,J))/A(L,IN) DO 6 L=I-IS,I1,-IS 6 B(L,J)=B(L,J)-A(L,IN1)*B(L+IS,J) 7 CONTINUE GOTO 72 C INVERSION WITH THE SWEEP METHOD 8 J=0 DO 12 J1=JM,IM J=J+1 DO 9 L=I1,J1-IS,IS 9 B(L,J)=0D0 B(J1,J)=1D0/A(J1,IN) DO 10 L=J1+IS,I,IS 10 B(L,J)=-R(IQ*(L-IS)+NQ)/A(L,IN)*B(L-IS,J) DO 11 L=I-IS,J1,-IS 11 B(L,J)=B(L,J)-A(L,IN1)*B(L+IS,J) DO 12 L=J1-IS,I1,-IS B(L,J)=-A(L,IN1)*B(L+IS,J) 12 CONTINUE GOTO 72 69 IF(IS.EQ.-1) GOTO 71 C M1=M2-IS NI=IN1+IS NQ1=NQ NQ=JQ-IS MQ=JQ*(1+IS)-IS IF(JM.NE.0) GOTO 38 DO 37 J=1,IM D=B(I1,J) D0=B(M2,J) DO 13 L=I1+IS,M2,IS A(L,IN1)=-R(IQ*(L-IS)+NQ1)*D D=A(L,IN1) IF(DABS(A(L-IS,IN)).LT.UNDFLO) GOTO 13 A(L,IN1)=A(L,IN1)/A(L-IS,IN) D=B(L,J)+A(L,IN1) 13 CONTINUE B(M2,J)=D/A(M2,IN) D2=DABS(D0) L=M2 LK=M2 14 KQ=L*IQ L=L-IS D0=-R(KQ+NQ)*D0 IF(DABS(A(L+IS,NI)).LT.UNDFLO) GOTO 20 D0=B(L,J)+D0/A(L+IS,NI) IF(L.EQ.I1) GOTO 36 IF(DABS(A(L-IS,IN)).LT.DABS(A(L+IS,NI))) GOTO 16 D1=A(L,IN) IF(L.EQ.M1) GOTO 15 IF(DABS(A(L+2*IS,NI)).LT.UNDFLO) GOTO 19 15 D1=D1-R(KQ+NQ)*R(KQ+MQ)/A(L+IS,NI) GOTO 19 16 IF(DABS(A(L-IS,IN)).GE.UNDFLO) GOTO 17 IF(DABS(A(L,IN)).LT.UNDFLO) RETURN D=A(L,IN1)/A(L,IN) GOTO 21 17 D1=A(L,NI) IF(L.EQ.I1+IS) GOTO 18 IF(DABS(A(L-2*IS,IN)).LT.UNDFLO) GOTO 19 18 D1=D1-R(IQ*L+NQ)*R(IQ*L+MQ)/A(L-IS,IN) 19 IF(DABS(D1).LT.UNDFLO) RETURN D=(D0+A(L,IN1))/D1 GOTO 21 20 IF(L.EQ.I1) GOTO 36 IF(DABS(A(L,NI)).LT.UNDFLO) RETURN D=D0/A(L,NI) 21 D3=0D0 IF(L.NE.M1) D3=R(KQ+NQ+IQ)*B(L+2*IS,J) IF(RL(R(KQ+MQ)*D,R(KQ+JQ)*B(L+IS,J),D3,D2).LE.Q1) GOTO 35 22 IF(DABS(A(L,IN)).GE.UNDFLO) GOTO 23 B(L,J)=D IF(LK.NE.M2) B(L,J)=DD+D1 IF(L.EQ.I1) GOTO 37 L=L-IS 23 D2=R(IQ*L+JQ) LK=L-IS D=-R(IQ*(L+IS)+NQ)/A(L,IN)*B(L+IS,J) D0=B(L,J) DD=(A(L,IN1)+D0)/A(L,IN) D1=D 24 KQ=L*IQ A(L,IN1)=DD DD2=DABS(B(L,J)) B(L,J)=DD+D1 IF(L.EQ.I1) GOTO 37 L=L-IS D4=D2 D=-R(KQ+NQ)*D IF(DABS(A(L,IN)).GE.UNDFLO) D=D/A(L,IN) IF(DABS(D).GT.1D0/EPS) GOTO 22 IF(DABS(D2).GE.UNDFLO) GOTO 25 D2=R(KQ+MQ) IF(L.EQ.I1) GOTO 34 IF(DABS(D2).LT.UNDFLO) RETURN DD=D0/D2 D1=D DD3=0D0 IF(L.NE.LK) DD3=R(KQ+NQ+IQ)*A(L+2*IS,IN1) IF(RL(R(KQ+MQ)*DD,R(KQ+JQ)*A(L+IS,IN1),DD3,DD2).GT.Q1) GOTO 22 DD2=DABS(B(L,J)) A(L,IN1)=DD B(L,J)=DD+D KQ=L*IQ L=L-IS D=-R(KQ+NQ)/A(L,IN)*D D0=B(L,J)-R(KQ+NQ)*DD D2=R(IQ*L+JQ) IF(L.EQ.I1) GOTO 34 IF(DABS(D).GT.1D0/EPS) GOTO 22 IF(DABS(A(L-IS,IN)).LT.UNDFLO) GOTO 29 D1=A(L,IN) GOTO 32 25 D0=B(L,J)-R(KQ+NQ)*D0/D2 IF(DABS(R(KQ+NQ)).LT.DABS(R(KQ+MQ))) GOTO 26 D3=R(KQ+MQ)/D2*R(KQ+NQ) GOTO 27 26 D3=R(KQ+NQ)/D2*R(KQ+MQ) 27 D2=R(IQ*L+JQ)-D3 IF(L.EQ.I1) GOTO 34 IF(DABS(A(L-IS,IN)).LT.DABS(D4)) GOTO 28 D1=A(L,IN)-D3 GOTO 32 28 IF(DABS(A(L-IS,IN)).GE.UNDFLO) GOTO 30 29 IF(DABS(A(L,IN)).LT.UNDFLO) RETURN DD=A(L,IN1)/A(L,IN) D1=0D0 GOTO 33 30 D1=D2 IF(L.EQ.I1+IS) GOTO 31 IF(DABS(A(L-2*IS,IN)).LT.UNDFLO) GOTO 32 31 D1=D1-R(IQ*L+NQ)*R(IQ*L+MQ)/A(L-IS,IN) 32 IF(DABS(D1).LT.UNDFLO) RETURN DD=(A(L,IN1)+D0)/D1 D1=D 33 DD3=0D0 IF(L.NE.LK) DD3=R(KQ+NQ+IQ)*A(L+2*IS,IN1) IF(RL(R(KQ+MQ)*DD,R(KQ+JQ)*A(L+IS,IN1),DD3,DD2).GT.Q1) GOTO 22 GOTO 24 34 IF(DABS(D2).LT.UNDFLO) RETURN B(L,J)=D0/D2+D GOTO 37 35 D2=DABS(B(L,J)) B(L,J)=D IF (L.NE.I1) GOTO 14 36 B(I1,J)=D0/A(I1,NI) 37 CONTINUE GOTO 72 C 38 J=0 DO 68 J1=JM,IM J=J+1 L=M2 B(J1,J)=1D0 IF(J1.EQ.I1) GOTO 40 IF(DABS(A(J1-IS,IN)).GE.UNDFLO) GOTO 40 DO 39 I=J1,M2,IS 39 B(I,J)=0D0 L=J1 GOTO 42 40 DO 41 I=J1+IS,M2,IS B(I,J)=-R(IQ*(I-IS)+NQ1)*B(I-IS,J) IF(DABS(A(I-IS,IN)).LT.UNDFLO) GOTO 41 B(I,J)=B(I,J)/A(I-IS,IN) 41 CONTINUE B(M2,J)=B(M2,J)/A(M2,IN) 42 LK=M2 D0=1D0 LM=I1 IF(J1.EQ.M2) GOTO 44 IF(DABS(A(J1+IS,NI)).GE.UNDFLO) GOTO 44 DO 43 I=I1,J1,IS 43 B(I,J)=0D0 LM=J1 IF(L.EQ.LM) GOTO 68 44 KQ=IQ*L L=L-IS D=0D0 IF(L.GE.J1) GOTO 45 D1=A(L,NI) B(L,J)=-R(KQ+NQ)*D0 D0=B(L,J) IF(DABS(A(L+IS,NI)).LT.UNDFLO) GOTO 50 B(L,J)=B(L,J)/A(L+IS,NI) D0=B(L,J) 45 IF(L.EQ.I1) GOTO 67 IF(DABS(A(L-IS,IN)).LT.DABS(A(L+IS,NI))) GOTO 47 IF(DABS(A(L+IS,NI)).LT.UNDFLO) GOTO 66 D1=A(L,IN) IF(L.EQ.M1) GOTO 46 IF(DABS(A(L+2*IS,NI)).LT.UNDFLO) GOTO 50 46 D1=D1-R(KQ+NQ)*R(KQ+MQ)/A(L+IS,NI) GOTO 50 47 D1=A(L,IN) IF(DABS(A(L-IS,IN)).GE.UNDFLO) GOTO 51 IF(L.GE.J1) GOTO 50 GOTO 66 51 D1=A(L,NI) IF(L.EQ.I1+IS) GOTO 49 IF(DABS(A(L-2*IS,IN)).LT.UNDFLO) GOTO 50 49 D1=D1-R(IQ*L+NQ)*R(IQ*L+MQ)/A(L-IS,IN) 50 IF(DABS(D1).LT.UNDFLO) RETURN D=B(L,J)/D1 D3=0D0 IF(L.NE.M1) D3=R(KQ+NQ+IQ)*B(L+2*IS,J) D2=0D0 IF(L+IS.EQ.J1) D2=1D0 If(RL(R(KQ+MQ)*D,R(KQ+JQ)*B(L+IS,J),D3,D2).LE.Q1) GOTO 66 52 IF(DABS(A(L,IN)).GE.UNDFLO) GOTO 53 B(L,J)=D IF(LK.NE.M2) B(L,J)=DD+D1 IF(L.EQ.LM) GOTO 68 L=L-IS 53 D2=R(IQ*L+JQ) DD1=0D0 LK=L D=-R(KQ+NQ)/A(L,IN)*B(L+IS,J) DD=0D0 IF(LK.GE.J1) DD=B(L,J)/A(L,IN) 54 D1=D 55 KQ=L*IQ D0=B(L,J) B(L,J)=DD+D1 IF(L.EQ.LM) GOTO 68 L=L-IS D4=D2 D=-R(KQ+NQ)*D IF(DABS(A(L,IN)).GE.UNDFLO) D=D/A(L,IN) IF(LK.LT.J1) GOTO 54 IF(DABS(D2).GE.UNDFLO) GOTO 56 D2=R(KQ+MQ) IF(L.EQ.I1) GOTO 65 IF(DABS(D2).LT.UNDFLO) RETURN DD2=DD1 DD1=DD DD=0D0 IF(L.LT.J1) DD=D0/D2 DD3=DABS(R(KQ+MQ)*DD+R(KQ+JQ)*DD1+R(KQ+NQ+IQ)*DD2) IF(L+IS.EQ.J1) DD3=DABS(1D0-D3) IF(DD3.GT.Q1) GOTO 52 D0=B(L,J) B(L,J)=DD+D KQ=L*IQ L=L-IS D=-R(KQ+NQ)/A(L,IN)*D IF(DABS(D).GT.1D0/EPS) GOTO 52 IF(L.LT.J1) B(L,J)=-R(KQ+NQ)*DD D2=R(IQ*L+JQ) IF(L.EQ.I1) GOTO 65 IF(DABS(A(L-IS,IN)).LT.UNDFLO) GOTO 60 D1=A(L,IN) GOTO 63 56 IF(L.LT.J1) B(L,J)=-R(KQ+NQ)*D0/D2 IF(DABS(R(KQ+NQ)).LT.DABS(R(KQ+MQ))) GOTO 57 D3=R(KQ+MQ)/D2*R(KQ+NQ) GOTO 58 57 D3=R(KQ+NQ)/D2*R(KQ+MQ) 58 D2=R(IQ*L+JQ)-D3 IF(L.EQ.I1) GOTO 65 IF(DABS(D).GT.1D0/EPS) GOTO 52 IF(DABS(A(L-IS,IN)).LT.DABS(D4)) GOTO 59 D1=A(L,IN)-D3 GOTO 63 59 IF(DABS(A(L-IS,IN)).GE.UNDFLO) GOTO 61 60 DD2=DD1 DD1=DD IF(DABS(A(L,IN)).LT.UNDFLO) RETURN DD=B(L,J)/A(L,IN) D1=0D0 GOTO 64 61 D1=D2 IF(L.EQ.I1+IS) GOTO 62 IF(DABS(A(L-2*IS,IN)).LT.UNDFLO) GOTO 63 62 D1=D1-R(IQ*L+NQ)*R(IQ*L+MQ)/A(L-IS,IN) 63 IF(DABS(D1).LT.UNDFLO) RETURN DD2=DD1 DD1=DD DD=B(L,J)/D1 D1=D 64 DD3=DABS(R(KQ+MQ)*DD+R(KQ+JQ)*DD1+R(KQ+NQ+IQ)*DD2) IF(L+IS.EQ.J1) DD3=DABS(1D0-DD3) IF(DD3.LE.Q1) GOTO 55 GOTO 52 65 IF(DABS(D2).LT.UNDFLO) RETURN B(L,J)=B(L,J)/D2+D GOTO 68 66 B(L,J)=D IF (L.GT.LM) GOTO 44 GOTO 68 67 B(L,J)=B(L,J)/A(L,NI) 68 CONTINUE 72 INF=0 IF(DET1.EQ.0D0) DET=0D0 IF(DABS(DET).LT.UNDFLO) RETURN DET=1D0 IE=0 DO 73 I=M2,I1,-IS IF(DABS(A(I,IN)).LT.UNDFLO) GOTO 73 DET1=GTEXP(A(I,IN),JE) DET=DET*DET1 DET=GTEXP(DET,IJE) IE=IE+JE+IJE 73 CONTINUE DET=CHEXP(DET,IE) RETURN 100 FORMAT(' M OR N LESS THEN 0 ON PROGRAM Lin3dsysccmSolver: M=', * I2,' N=',I2) END C C Функция RL - вычисления величину D=|D1+D2+D3|-|Y| C C Обpащение: RL(D1,D2,D3,Y) C Входные данные: C D1,D2,D3,Y - (real*8) вещественные данные. C Выходные данные: C RL - (real*8) содержить величину D. C REAL*8 FUNCTION RL(D1,D2,D3,Y) REAL*8 D1,D2,D3,D,Y IF(DABS(D1).LT.DABS(D2)) GOTO 3 IF(DABS(D2).LT.DABS(D3)) GOTO 2 1 D=DABS(D1+D2+D3) GOTO 4 2 D=DABS(D1+D3+D2) GOTO 4 3 IF(DABS(D3).LT.DABS(D1)) GOTO 1 D=DABS(D2+D3+D1) 4 IF(Y.GT.1D0) GOTO 5 RL=DABS(Y-D) RETURN 5 RL=DABS(1D0-D/Y) RETURN END