      SUBROUTINE MINV (A,N,NP)
C***********************************************************************
C***
C***  Inversion of matrix A
C***
C***  N = rank of matrix to be inverted
C***  NP = physical size of the matrix 
C***********************************************************************

      IMPLICIT NONE
      INTEGER I,J,K,L,N,NP,IMAX,IK,JK,ISAMAX
      DOUBLE PRECISION A,AMAX,ABSAMAX,SAVE,AIK
      DIMENSION A(NP,NP),IK(NP),JK(NP)

      DO 100  K=1,N
C
C     search largest matrix element
         AMAX=0.
         ABSAMAX=0.
         DO 30 J=K,N
            L=N+1-K
            IMAX=ISAMAX(L,A(K,J),1)+K-1
            IF(ABSAMAX.GT.ABS(A(IMAX,J))) GOTO 30
                AMAX=A(IMAX,J)
                ABSAMAX=ABS(AMAX)
                IK(K)=IMAX
                JK(K)=J
   30     CONTINUE
C
C***  stop in case of singularity (pivot element = 0 )
          IF (AMAX .EQ. .0) THEN
               STOP 'SINGULAR'
          ENDIF
C
C     exchange rows I and K
          I=IK(K)
          IF(I.EQ.K) GOTO 51
          DO 50 J=1,N
             SAVE =A(K,J)
             A(K,J)=A(I,J)
   50     A(I,J)=-SAVE
C
C     exchange columns J and K
   51     J=JK(K)
          IF(J.EQ.K) GOTO 61
          DO 60 I=1,N
             SAVE = A(I,K)
             A(I,K)=A(I,J)
   60     A(I,J)=-SAVE
C
C     division of column K by AMAX
   61     DO 70  I=1,N
              A(I,K)=-A(I,K)/AMAX
   70     CONTINUE
          A(K,K)=0.
C
C     transform rows and columns
C***  row I = row  I + A(I,K) * row K
          DO 81 I=1,N
             IF (I .EQ. K) GOTO 81
             AIK=A(I,K)
             DO 80 J=1,N
                 A(I,J)=A(I,J)+AIK*A(K,J)
   80        CONTINUE
   81      CONTINUE
C
C***  row K: division by AMAX
           DO 90 J=1,N
              A(K,J)=A(K,J)/AMAX
  90       CONTINUE
C
C***  diagonal element
           A(K,K)=1./AMAX
  100 CONTINUE
C
C
C     row and columns back-substitution
      DO 130 L=1,N
         K=N-L+1
         J=IK(K)
         IF(J.LE.K) GOTO 111
C***     exchange columns J and K
         DO 110 I=1,N
            SAVE=A(I,K)
            A(I,K)=-A(I,J)
  110       A(I,J)=SAVE
  111     I=JK(K)
          IF(I .LE. K) GOTO 130
C***      exchange rows I and K
          DO 120 J=1,N
            SAVE=A(K,J)
            A(K,J)=-A(I,J)
  120     A(I,J)=SAVE
  130 CONTINUE
C
      RETURN
      END
      INTEGER FUNCTION ISAMAX(N,A,INC)
C *** RETURNS THE ARRAY INDEX OF THE ELEMENT WITH MAXIMUM ABSOLUTE VALUE
      IMPLICIT NONE
      DOUBLE PRECISION A
      INTEGER N,INC,I
      DIMENSION A(N)
      ISAMAX=1
      DO 1 I=2,N,INC
        IF(ABS(A(ISAMAX)).GE.ABS(A(I)))                       GOTO 1
        ISAMAX=I
1       CONTINUE
      RETURN
      END
