91 lines
2.0 KiB
Fortran
91 lines
2.0 KiB
Fortran
SUBROUTINE MATINV(A,N,NR)
|
|
C =========================
|
|
C
|
|
C Matrix inversion
|
|
C by LU decomposition
|
|
C
|
|
C A - matrix of actual size (N x N) and maximum size (NR x NR)
|
|
C to be inverted;
|
|
C Inversion is accomplished in place and the original matrix is
|
|
C replaced by its inverse
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
DIMENSION A(NR,NR)
|
|
C
|
|
DO I=2,N
|
|
IM1=I-1
|
|
DO J=1,IM1
|
|
JM1=J-1
|
|
DIV=A(J,J)
|
|
SUM=0.
|
|
IF(JM1.GE.1) THEN
|
|
DO K=1,JM1
|
|
SUM=SUM+A(I,K)*A(K,J)
|
|
END DO
|
|
END IF
|
|
A(I,J)=(A(I,J)-SUM)/DIV
|
|
END DO
|
|
DO J=I,N
|
|
SUM=0.
|
|
DO K=1,IM1
|
|
SUM=SUM+A(I,K)*A(K,J)
|
|
END DO
|
|
A(I,J)=A(I,J)-SUM
|
|
END DO
|
|
END DO
|
|
C
|
|
DO II=2,N
|
|
I=N+2-II
|
|
IM1=I-1
|
|
IF(IM1.GE.1) THEN
|
|
DO JJ=1,IM1
|
|
J=I-JJ
|
|
JP1=J+1
|
|
SUM=0.
|
|
IF(JP1.LE.IM1) THEN
|
|
DO K=JP1,IM1
|
|
SUM=SUM+A(I,K)*A(K,J)
|
|
END DO
|
|
END IF
|
|
A(I,J)=-A(I,J)-SUM
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
DO II=1,N
|
|
I=N+1-II
|
|
DIV=A(I,I)
|
|
IP1=I+1
|
|
IF(IP1.LE.N) THEN
|
|
DO JJ=IP1,N
|
|
J=N+IP1-JJ
|
|
SUM=0.
|
|
DO K=IP1,J
|
|
SUM=SUM+A(I,K)*A(K,J)
|
|
END DO
|
|
A(I,J)=-SUM/DIV
|
|
END DO
|
|
END IF
|
|
A(I,I)=1.0D0/A(I,I)
|
|
END DO
|
|
C
|
|
DO I=1,N
|
|
DO 230 J=1,N
|
|
K0=I
|
|
IF(J.GE.I) GO TO 220
|
|
SUM=0.
|
|
200 DO K=K0,N
|
|
SUM=SUM+A(I,K)*A(K,J)
|
|
END DO
|
|
GO TO 230
|
|
220 K0=J
|
|
SUM=A(I,K0)
|
|
IF(K0.EQ.N) GO TO 230
|
|
K0=K0+1
|
|
GO TO 200
|
|
230 A(I,J)=SUM
|
|
END DO
|
|
RETURN
|
|
END
|