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 'PARAMS.FOR'
      DIMENSION A(NR,NR)
      IF(N.EQ.1) GO TO 250
      DO 50 I=2,N
         IM1=I-1
         DO 20 J=1,IM1
            JM1=J-1
            DIV=A(J,J)
            SUM=0.
            IF(JM1.LT.1) GO TO 20
            DO 10 K=1,JM1
   10          SUM=SUM+A(I,K)*A(K,J)
   20       A(I,J)=(A(I,J)-SUM)/DIV
         DO 40 J=I,N
            SUM=0.
            DO 30 K=1,IM1
   30          SUM=SUM+A(I,K)*A(K,J)
   40       A(I,J)=A(I,J)-SUM
   50 CONTINUE
      DO 80 II=2,N
         I=N+2-II
         IM1=I-1
         IF(IM1.LT.1) GO TO 80
         DO 70 JJ=1,IM1
            J=I-JJ
            JP1=J+1
            SUM=0.
            IF(JP1.GT.IM1) GO TO 70
            DO 60 K=JP1,IM1
   60          SUM=SUM+A(I,K)*A(K,J)
   70       A(I,J)=-A(I,J)-SUM
   80 CONTINUE
      DO 110 II=1,N
         I=N+1-II
         DIV=A(I,I)
         IP1=I+1
         IF(IP1.GT.N) GO TO 110
         DO 100 JJ=IP1,N
            J=N+IP1-JJ
            SUM=0.
            DO 90 K=IP1,J
   90          SUM=SUM+A(I,K)*A(K,J)
            A(I,J)=-SUM/DIV
  100    CONTINUE
  110 A(I,I)=1.0D0/A(I,I)
C
      DO 240 I=1,N
         DO 230 J=1,N
            K0=I
            IF(J.GE.I) GO TO 220
            SUM=0.
  200       DO 210 K=K0,N
  210          SUM=SUM+A(I,K)*A(K,J)
            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
  240 CONTINUE
      RETURN
  250 A(1,1)=1.0D0/A(1,1)
      RETURN
      END