C SUBROUTINE Invert3dMatrix(M,A,INF,R1,R2,R,IZ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(M,*),R1(*),R2(*),R(*),IZ(*) C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C* * C* AUTHORS :Emel'yanenko G.A., Rakhmonov T.T., Dushanov E.B. * C* REPRESENTATIVE: A.P.Sapoznikov (LCTA JINR) * C* * C* * C* SUBROUTINE Invert3dMatrix(M,A,INF,R1,R2,R,IZ) * C* INVERTS THE TRIDIAGONAL MATRICES [1] OF GENERAL * C* FORM (INCLUDING ILL-POSED ONES) * C* BY USING THE SUBPROGRAM Invert3dwpMatrix. * C* * C* M (INTEGER) DIMENSION OF QUADRATIC MATRIX C; * C* A (REAL*8) A TWO-DIMENSIONAL ARRAY OF DIMENSION M*M,THAT AT * C* THE INPUT CONTAINS THE MATRIX C, AND AT THE OUTPUT B=C^{-1} C* AT THE INPUT, THE MATRIX C LIES IN THE ARRAY A IN THE FORM * C* i i * C* ³ * q{1} r{2} * ... * ³ * C* ³ p{2} q{2} r{3} * ... * ³ * C* A = ³ : : : * ... * ³, * C* ³ p{M-1} q{M-1} r{M} * ... * ³ * C* ³ p{M} q{M} * * ... * ³ * C* i i * C* WHERE SYMBOL * MARKS CELLS OF ARRAY A, WHOSE CONTENTS IS * C* NOT IMPORTANT AT THE INPUT, AND AT THE OUTPUT B=C^{-1} IS * C* LOCATED BY COLUMNS IN THE FIRST M COLUMNS AND M ROWS OF * C* THE ARRAY A; * C* INF (INTEGER) AN OUTPUT PARAMETER: * C* INF= 0 NORMAL COMPLETION OF THE WORK; * C* INF=-1 THE INITIAL MATRIX IS SINGULAR; * C* R1,R2,R (REAL*8) ONE-DIMENSIONAL WORKING ARRAYS OF DIMENSION M. * C* IZ (INTEGER) ONE-DIMENSIONAL WORKING ARRAY OF DIMENSION M. * C* * C* * C* REFERENCE: * C* 1. Emel'yanenko G.A., Rakhmonov T.T., Dushanov E.B. Algorithms and * C* programs of the critical-component method of inversion tridiagonal* C* matrices and solution of systems of linear equations. JINR * C* preprint, E11-96-106, Dubna, 1996. * C* * C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DATA ZERO/0.D0/ C IF(M.LE.1) GO TO 38 M1=M-1 I=0 K=0 K1=4 K2=0 2 I=I+1 K=K+1 IZ(K)=I I1=I-1 INF=M CALL Invert3dwpMatrix(M,A,I,INF,R1,R2,R) ! dcinv0 C WRITE(*,*)I1+1,I IF(INF.EQ.-1) RETURN IZ(K+1)=I I2=I+1 K=K+1 IF(I1.EQ.0) GO TO 51 D0=A(I1+1,1) A(I1,I1)=A(I1,I1)-A(I1,I1+1)*D0*A(I1+1,I1+1) D=-A(I1,I1+1)/R2(I1) DO 52 L=I1+1,I A(I1,L)=A(I1+1,L)*D 52 R1(L)=-A(L,I1+1)*D0 A(I1-1,I1)=A(I1-1,I1)/R1(I1) A(I1,I1-1)=A(I1,I1-1)/R2(I1) A(I1,I1)=A(I1,I1)/R1(I1)/R2(I1) A(I,I1)=R2(I1) 51 IF(I.EQ.M) GOTO 18 D=-A(I2,1)/R2(I2) D1=A(I2,3) DO 53 L=I1+1,I A(I2,L)=A(I,L)*D 53 R2(L)=-A(L,I)*A(I,I2) A(I2,I2)=R(I2) IF(I1.EQ.0) GO TO 13 A(I1,I2)=A(I1,I)*A(I,I2)/R1(I2) A(I2,I1)=A(I2,I1+1)*D0/R1(I1) D=A(I1,I2)/(R(I1-1)*A(I1,I1)-A(I1-1,I1)*A(I1,I1-1)) R2(I1-1)=A(I1-1,I1)*D R2(I1)=-R(I1-1)*D R(I2)=R(I2)+A(I2,I1)*R2(I1) 13 A(I1+1,I2)=R2(I2) IF(I2.NE.M) GO TO 17 D=A(M,M) IF((DABS(D).EQ.ZERO).OR.(DABS(R(M)).EQ.ZERO)) GOTO 37 INF=0 A(M,M)=1D0/R(M) IF(I1.NE.0) GO TO 15 A(M,M)=A(M,M)/(R1(M)*R2(M)) DO 14 L=I1+1,I A(M,L)=A(M,L)/(R(M)*R1(M)) A(L,M)=R2(L)*A(M,M) DO 14 J=1,M1 14 A(J,L)=A(J,L)+R2(J)*A(M,L) GO TO 38 15 A(M,I1)=-A(M,1)/D A(I1,I1)=A(I1,I1)+A(I1,M)*A(M,I1) A(I1-1,M)=R2(I1-1)/R(M) A(I1,M)=R2(I1)/R(M) M0=1 GO TO 18 17 I=I+2 A(I2,I2+1)=D1/R2(I2) D=A(I2+1,1)/R1(I2) IF(I.EQ.M) GO TO 60 A(I2+1,I2+2)=A(I2+1,3) A(I2+1,I2+1)=A(I2+1,2) A(I2+1,I2)=D GOTO 2 60 K2=K I1=I IZ(K+1)=M+1 A(M,M)=A(M,2) GO TO 19 18 K2=K-2 19 INF=0 IF((I1.EQ.0).AND.(I.EQ.M)) GO TO 38 D1=A(I1-1,I1-1) D2=A(I1-1,I1) D3=A(I1,I1-1) D4=A(I1,I1) D=R(I1-1)*D4-D2*D3 IF(DABS(D).EQ.ZERO) GO TO 37 A(I1-1,I1-1)=D4/D A(I1-1,I1)=-D2/D A(I1,I1-1)=-D3/D A(I1,I1)=R(I1-1)/D DO 20 J=K2,K1,-2 J1=IZ(J-1) J2=IZ(J) D=A(J2+1,J1-1)/(D1*D4-D2*D3) A(J2+1,J1-1)=-D4*D A(J2+2,J1-1)=D3*D D1=A(J1-2,J1-2) D2=A(J1-2,J1-1) D3=A(J1-1,J1-2) D4=A(J1-1,J1-1)+A(J1-1,J2+1)*A(J2+1,J1-1) D=R(J1-2)*D4-D2*D3 A(J1-2,J1-2)=D4/D A(J1-2,J1-1)=-D2/D A(J1-1,J1-2)=-D3/D 20 A(J1-1,J1-1)=R(J1-2)/D * DO 21 I=K1,K2,2 DO 21 J=I,K2,2 J2=IZ(J-1) DO 21 J1=1,2 D1=A(IZ(J)+J1,J2-1)*A(J2-1,IZ(I-1)-2) A(IZ(J)+J1,IZ(I-1)-1)=A(IZ(J)+J1,J2-1)*A(J2-1,IZ(I-1)-1) 21 A(IZ(J)+J1,IZ(I-1)-2)=D1 DO 22 I=K2,K1,-2 DO 22 L=I,K1,-2 J2=IZ(L-1) DO 22 L1=1,2 A(J2-1,IZ(I)+L1)=R2(J2-1)*A(IZ(L)+1,IZ(I)+L1) 22 A(J2-2,IZ(I)+L1)=R2(J2-2)*A(IZ(L)+1,IZ(I)+L1) J1=IZ(K1-1) J2=IZ(K2) IF(M0.NE.1) GO TO 24 D=A(M,J2+2) DO 23 J=2,1,-1 A(M,J1-J)=D*A(J2+2,J1-J) DO 23 L=K2,K1,-2 A(IZ(L-1)-J,M)=R2(IZ(L-1)-J)*A(IZ(L)+1,M) A(M,IZ(L)+J)=D*A(J2+2,IZ(L)+J) 23 CONTINUE 24 DO 41 L=1,K-1,2 L1=IZ(L) L2=IZ(L+1) IF(L.EQ.1) GO TO 39 A(L1-1,J1-2)=A(L1-1,J1-2)/R1(L1-1) A(L1-1,J1-1)=A(L1-1,J1-1)/R1(L1-1) DO 32 J=K1,K2,2 A(L1-1,IZ(J)+1)=A(L1-1,IZ(J)+1)/R1(L1-1) 32 A(L1-1,IZ(J)+2)=A(L1-1,IZ(J)+2)/R1(L1-1) IF(M0.EQ.1) A(L1-1,M)=A(L1-1,M)/R1(L1-1) IF(L2.EQ.M) GO TO 41 39 A(L2+1,J1-2)=A(L2+1,J1-2)/R1(L2+1) A(L2+1,J1-1)=A(L2+1,J1-1)/R1(L2+1) DO 40 J=K1,K2,2 A(L2+1,IZ(J)+1)=A(L2+1,IZ(J)+1)/R1(L2+1) 40 A(L2+1,IZ(J)+2)=A(L2+1,IZ(J)+2)/R1(L2+1) IF(M0.EQ.1) A(L2+1,M)=A(L2+1,M)/R1(L2+1) 41 CONTINUE DO 29 I=1,K-1,2 DO 29 J=IZ(I),IZ(I+1) R1(J1-2)=0.D0 R1(J1-1)=0.D0 IF(IZ(I).EQ.1) GO TO 25 R1(J1-2)=A(J1-2,IZ(I)-1)*A(IZ(I)-1,J) R1(J1-1)=A(J1-1,IZ(I)-1)*A(IZ(I)-1,J) 25 DO 27 L=K1,K2,2 L1=IZ(L-1) R1(IZ(L)+1)=0.D0 R1(IZ(L)+2)=0.D0 IF(IZ(I).EQ.1) GO TO 26 R1(IZ(L)+1)=A(IZ(L)+1,IZ(I)-1)*A(IZ(I)-1,J) R1(IZ(L)+2)=A(IZ(L)+2,IZ(I)-1)*A(IZ(I)-1,J) IF(IZ(I+1).EQ.M) GO TO 27 26 R1(L1-2)=R1(L1-2)+A(L1-2,IZ(I+1)+1)*A(IZ(I+1)+1,J) R1(L1-1)=R1(L1-1)+A(L1-1,IZ(I+1)+1)*A(IZ(I+1)+1,J) 27 CONTINUE IF(IZ(I+1).EQ.M) GO TO 31 R1(J2+1)=R1(J2+1)+A(J2+1,IZ(I+1)+1)*A(IZ(I+1)+1,J) R1(J2+2)=R1(J2+2)+A(J2+2,IZ(I+1)+1)*A(IZ(I+1)+1,J) IF(M0.NE.1) GO TO 31 R1(M)=0.D0 IF(IZ(I).NE.1) R1(M)=A(M,IZ(I)-1)*A(IZ(I)-1,J) A(M,J)=R1(M)+A(M,IZ(I+1)+1)*A(IZ(I+1)+1,J) 31 DO 29 L1=1,2 A(J1-L1,J)=R1(J1-L1) DO 29 L=K1,K2,2 A(IZ(L)+L1,J)=R1(IZ(L)+L1) 29 CONTINUE DO 46 L=1,K-1,2 L1=IZ(L) L2=IZ(L+1) IF(L.EQ.1) GO TO 44 A(J1-2,L1-1)=A(J1-2,L1-1)/A(L2,L1-1) A(J1-1,L1-1)=A(J1-1,L1-1)/A(L2,L1-1) DO 43 J=K1,K2,2 A(IZ(J)+1,L1-1)=A(IZ(J)+1,L1-1)/A(L2,L1-1) 43 A(IZ(J)+2,L1-1)=A(IZ(J)+2,L1-1)/A(L2,L1-1) IF(M0.EQ.1) A(M,L1-1)=A(M,L1-1)/A(L2,L1-1) IF(L2.EQ.M) GO TO 46 44 A(J1-2,L2+1)=A(J1-2,L2+1)/A(L1,L2+1) A(J1-1,L2+1)=A(J1-1,L2+1)/A(L1,L2+1) DO 45 J=K1,K2,2 A(IZ(J)+1,L2+1)=A(IZ(J)+1,L2+1)/A(L1,L2+1) 45 A(IZ(J)+2,L2+1)=A(IZ(J)+2,L2+1)/A(L1,L2+1) IF(M0.EQ.1) A(M,L2+1)=A(M,L2+1)/A(L1,L2+1) 46 CONTINUE * DO 49 I=IZ(1),IZ(2) DO 49 J=IZ(2)+1,M 49 A(I,J)=0.D0 DO 30 I=1,K-1,2 IF(IZ(I).EQ.1) GO TO 50 DO 33 J=IZ(I),IZ(I+1) DO 34 L=1,IZ(I)-1 34 A(J,L)=R1(J)*A(IZ(I)-1,L) DO 35 L=IZ(I),IZ(I+1) 35 A(J,L)=A(J,L)+R1(J)*A(IZ(I)-1,L) DO 33 L=IZ(I+1)+1,M A(J,L)=R1(J)*A(IZ(I)-1,L) 33 CONTINUE IF(IZ(I+1).EQ.M) GO TO 38 50 DO 47 J=IZ(I),IZ(I+1) DO 47 L=1,M 47 A(J,L)=A(J,L)+R2(J)*A(IZ(I+1)+1,L) 30 CONTINUE GO TO 38 37 INF=-1 38 RETURN END C SUBROUTINE Invert3dwpMatrix(M,A,I,INF,R1,R2,R) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(M,*),R1(*),R2(*),R(*) C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C* * C* SUBROUTINE Invert3dwpMatrix(M,A,I,INF,R1,R2,R) * C* SEPARATES IN ARRAY A AND INVERTS WELL-POSED * C* SUBMATRICES TRIDIAGONAL MATRIX C. * C* Invert3dwpMatrix WORKS IN THE TWO REGIMES: IN THE REGIME * C* TOGETHER WITH Invert3dMatrix, THAT ITSELF CALLS TO HER, AND IN AN * C* INDEPENDENT REGIME (IF IT A PRIORI KNOWN THAT THE INITIAL MATRIX C * C* IS WELL-POSED, THEN IT CAN BE INVERTED IMMEDIATELY * C* BY Invert3dwpMatrix, WITHOUT CALL TO Invert3dMatrix). * C* * C* M, A AND R1 ARE THE SAME PARAMETERS AS FOR Invert3dMatrix. * C* I, INF (INTEGER) * C* AT THE INPUT: I=1, INF=M. HERE I AND INF ARE ASSIGNED BEFORE * C* AT THE OUTPUT: I=M AND * C* INF= 0 NORMAL COMPLETION OF THE WORK; * C* INF=-1 THE INITIAL MATRIX C IS SINGLET; * C* R1,R2 (REAL*8) ONE-DIMENSIONAL WORKING ARRAYS OF DIMENSION M. * C* * C* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * COMMON /F499_RCONST/ OVFLO,UNDFLO,EPS,EPSMIN DATA ZERO,ONE/0.D0,1.D0/ C call Init_Const ! Sap: to hide this actin from User q1=eps q2=one/q1 M1=INF+1 IF((M1.LE.0).OR.(M.LE.0).OR.(I.LE.0).OR.(I.GT.M)) GOTO 23 INF=-1 R(I)=A(I,2) K0=I IF(I.EQ.M) GOTO 38 IF(I.EQ.M1) GOTO 23 D2=ONE D4=D2 IL=0 1 I=I+1 IX=0 IF(DABS(R(I-1)).NE.ZERO) GOTO 3 R(I)=-A(I,1)*A(I-1,3) IF(DABS(R(I)).EQ.ZERO) RETURN R2(I-1)=-A(I,1) R1(I-1)=A(I-1,3) IL=1 IF(I.EQ.M) GOTO 9 IF(I.EQ.M1) GOTO 11 D4=DMAX1(ONE,D4)*DABS(A(I,3)/A(I,1)) I=I+1 D2=DMAX1(ONE,D2)*DABS(A(I,1)/A(I-2,3)) R(I)=A(I,2) IX=1 3 R2(I-1)=-A(I,1)/R(I-1) R1(I-1)=-A(I-1,3)/R(I-1) IF(IX.EQ.1) GOTO 11 IF(DABS(R1(I-1)).GT.ONE) IL=1 IF(DABS(A(I-1,3)).GT.DABS(A(I,1))) GOTO 4 R(I)=A(I,2)+R1(I-1)*A(I,1) GO TO 5 4 R(I)=A(I,2)+R2(I-1)*A(I-1,3) 5 D2=DMAX1(ONE,D2)*DABS(R2(I-1)) D4=DMAX1(ONE,D4)*DABS(R1(I-1)) IF(I.EQ.M1) GOTO 6 11 IF(D2.GE.q2.OR.D2.LE.Q1.or.D4.GE.q2.OR.D4.LE.Q1) GOTO 6 IF(I.LT.M) GOTO 1 IF(DABS(R(I)).EQ.ZERO) RETURN GOTO 38 6 R(I)=R(I)/D2/D4 R1(I)=D4 R2(I)=D2 I=I-1 A(I,I+1)=A(I,3) IF(I.EQ.K0-1) GOTO 23 38 IF(K0.EQ.1) GOTO 7 K1=K0-1 R1(K1)=ONE R2(K1)=ONE D1=ONE D3=ONE DO 40 L=K0,I-1 D1=D1*R1(L) D3=R2(L)*D3 IF(DABS(R(L)).LE.Q1) GOTO 40 R1(K1)=DMAX1(R1(K1),DABS(D3)) R2(K1)=DMAX1(R2(K1),DABS(D1)) 40 CONTINUE 7 IF(K0.NE.I) GOTO 9 IF(DABS(R(K0)).EQ.ZERO) RETURN A(I,I)=ONE/R(K0) GOTO 23 9 R2(I)=A(I,2) IF(IL.NE.0) GOTO 42 A(I,I)=ONE/R(I) DO 41 J=I-1,K0,-1 A(J,J)=(ONE-R1(J)*A(J+1,J+1)*A(J+1,1))/R(J) DO 41 L=J,K0,-1 A(L,J+1)=R1(L)*A(L+1,J+1) A(J+1,L)=A(J+1,L+1)*R2(L) 41 CONTINUE GOTO 23 42 A(I,I)=R(I) J=I 10 J=J-1 IX=0 IF(DABS(R2(J+1)).NE.ZERO) GOTO 13 IL=0 D=R2(J) R2(J)=-A(J+1,1)*A(J,3) IF(DABS(R2(J)).EQ.ZERO) RETURN A(J,J+1)=-A(J,3) A(J+1,J)=D A(J,J)=R2(J) IF(J.EQ.K0) GOTO 24 J=J-1 IF(DABS(R(J)).EQ.ZERO) RETURN IX=1 D3=ZERO 13 D0=-A(J,3)/R2(J+1) IF(IX.EQ.1) GOTO 15 R1(J+1)=-A(J+1,1)/R2(J+1) IF(DABS(R1(J+1)).GT.ONE) IL=0 IF(DABS(A(J+1,1)).LT.DABS(A(J,3))) GOTO 14 D3=D0*A(J+1,1) GOTO 15 14 D3=R1(J+1)*A(J,3) 15 A(J+1,J)=R2(J) R2(J)=A(J,2)+D3 A(J,J+1)=D0 IF(J.EQ.K0) GOTO 19 IF(DABS(R(J-1)).LT.DABS(R2(J+1))) GOTO 18 A(J,J)=R(J)+D3 GOTO 10 18 IF(R(J-1).NE.ZERO) GOTO 17 A(J,J)=R(J) GOTO 10 17 A(J,J)=R2(J) IF(J.EQ.K0+1) GOTO 16 IF(DABS(R(J-2)).EQ.ZERO) GOTO 10 16 A(J,J)=R2(J)-A(J-1,3)*A(J,1)/R(J-1) GOTO 10 19 IF(DABS(R2(K0)).EQ.ZERO) RETURN A(K0,K0)=R2(K0) C WRITE(*,111)(a(l,l),R1(L),R(L),R2(L),l=k0,i) C 111 FORMAT(4d10.3) 24 IF(IL.NE.1) GOTO 26 A(K0,K0)=ONE/R2(K0) DO 28 J=K0+1,I A(J,J)=ONE/R2(J)+R1(J)*A(J-1,J-1)*A(J-1,J) DO 28 L=J,I A(J-1,L)=A(J-1,L-1)*A(L-1,L) 28 A(L,J-1)=R1(L)*A(L-1,J-1) GOTO 23 26 DO 22 J=K0+1,I IF((DABS(R(J-2)).EQ.ZERO).AND.(J.NE.K0+1)) A(J,J-1)=ZERO DO 20 L=K0,J-1 IF(J.EQ.I) GOTO 20 A(J+1,L)=A(J+1,J)*A(J,L) IF(DABS(R2(J+1)).EQ.ZERO) A(J,L)=ZERO 20 A(J,L)=A(J,L)/A(J,J) K=I-J+K0 IF((DABS(R2(K+2)).EQ.ZERO).AND.(K.NE.I-1)) A(K,K+1)=ZERO DO 21 L=K+1,I IF(K.EQ.K0) GOTO 21 A(K-1,L)=A(K-1,K)*A(K,L) IF(DABS(R(K-1)).EQ.ZERO) A(K,L)=ZERO 21 A(K,L)=A(K,L)/A(K,K) 22 CONTINUE A(K0,K0)=ONE/A(K0,K0) IF(DABS(R2(K0+1)).EQ.ZERO) A(K0,K0)=ZERO DO 25 L=K0+1,I-1 D=ONE IF(DABS(R(L-1)).EQ.ZERO) D=ZERO IF(DABS(R2(L+1)).EQ.ZERO) D=ZERO 25 A(L,L)=D/A(L,L) A(I,I)=ONE/A(I,I) IF(DABS(R(I-1)).EQ.ZERO) A(I,I)=ZERO 23 INF=0 RETURN END