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