C
C Subroutine for inverting a matrix
C
       SUBROUTINE INVERSE(N,A)
C
C CALCULATION  OF THE INVERSE MATRIX IN PLACE
C N is the dimension of the matrix A, to be inverted. When the subroutine
C returns, A will hold the inverse matrix.
C Adapted from an old scientific subroutine package.
C mvr-2011
C     
       DIMENSION X(N),A(N,N)
       PARAMETER(EPS = 1.0e-20)
C
C      ...... CALL SIMUL........
       DETER=SIMUL(N,A,X,EPS,-1,n)
C
       WRITE (*,202) DETER
       
 202   FORMAT (/' DETER=',F12.6/)
C
       RETURN
       END
C
       FUNCTION SIMUL(N,A,X,EPS,INDIC,NRC)
       DIMENSION IROW(N),JCOL(N),JORD(N),Y(N),A(NRC,NRC),X(N)
C
       MAX=N
       IF(INDIC.GE.0) MAX=N+1
C
C      .......IS N LARGER THAN 50 ...............
       IF(N.LE.50) GO TO 5
       WRITE (6,200)
       SIMUL=0.
       RETURN
C
C      ......... BEGIN ELIMINATION PROCEDURE .........
   5   DETER=1.
       DO 18 K= 1,N
       KM1=K-1
C
C      .......... SEARCH FOR THE PIVOT ELEMENT........
       PIVOT=0.
       DO 11 I=1,N
       DO 11 J=1,N
C      .......... SCAN IROW AND JCOL ARRAYS FOR INVALID PIVOT SUBS....
       IF (K.EQ.1) GO TO 9
       DO 8  ISCAN=1,KM1
       DO 8  JSCAN=1,KM1
       IF (I.EQ.IROW(ISCAN)) GO TO 11
       IF (J.EQ.JCOL(JSCAN)) GO TO 11
   8   CONTINUE
   9   IF(ABS(A(I,J)).LE.ABS(PIVOT)) GO TO 11
       PIVOT=A(I,J)
       IROW(K)=I
       JCOL(K)=J
  11   CONTINUE
C
C      .....INSURE THAT SELETED PIVOT IS LARGER THAN EPS.....
       IF(ABS(PIVOT).GT.EPS) GO TO 13
       SIMUL=0.
       RETURN
C
C      ......UPDATE THE DETERMINANT VALUE......
  13   IROWK=IROW(K)
       JCOLK=JCOL(K)
       DETER=DETER*PIVOT
C
C      .......NORMALISE PIVOT ROW ELEMENTS.......
       DO 14 J= 1,MAX
C
  14       A(IROWK,J)=A(IROWK,J)/PIVOT
C      ....... CARRY OUT ELIMINATION AND DEVELOP INVERSE......
       A(IROWK,JCOLK) =1./PIVOT
       DO 18 I= 1 , N
       AIJCK=A(I,JCOLK)
       IF(I.EQ.IROWK) GO TO 18
       A(I,JCOLK)=-AIJCK/PIVOT
       DO 17 J= 1,MAX
 17    IF(J.NE.JCOLK) A(I,J)=A(I,J)-AIJCK*A(IROWK,J)
 18    CONTINUE
C
C      ..... ORDER SOLUTION VALUES(IF ANY)AND CREATE JORD ARRAY....
       DO 20 I=1,N
       IROWI=IROW(I)
       JCOLI=JCOL(I)
       JORD(IROWI)=JCOLI
 20    IF(INDIC.GE.0)   X(JCOLI)=A(IROWI,MAX)
C
C      ..... ADJUST SIGN OF DETERMINANT........
       INTCH=0
       NM1=N-1
       DO 22 I=1,NM1
       IP1=I+1
       DO 22 J=IP1,N
       IF(JORD(J).GE.JORD(I)) GO TO 22
       JTEMP=JORD(J)
       JORD(J)=JORD(I)
       JORD(I)=JTEMP
       INTCH=INTCH+1
 22    CONTINUE
       IF (INTCH/2*2.NE.INTCH)   DETER=-DETER
C
C      .........IF INDIC IS POSITIVE RETURN WITH RESULTS.....
       IF (INDIC.LE.0) GO TO 26
       SIMUL=DETER
       RETURN
C
C      ..IF INDIC IS NEGATIVE OR 0 UNSCRAMBLE THE INVERSE 1ST.BY ROWS...
 26    DO 28 J=1,N
       DO 27 I=1,N
       IROWI=IROW(I)
       JCOLI=JCOL(I)
 27    Y(JCOLI)=A(IROWI,J)
       DO 28 I= 1,N
 28    A(I,J)=Y(I)
C      ...... THEN BY COLUMNS.....
       DO 30 I= 1,N
       DO 29 J=1,N
       IROWJ=IROW(J)
       JCOLJ=JCOL(J)
 29    Y(IROWJ)=A(I,JCOLJ)
       DO 30 J= 1,N
 30    A(I,J)=Y(J)
C
C      ....RETURN FOR INDIG NEGATIVE OR ZERO......
       SIMUL=DETER
       RETURN
C
C      ......FORMAT FOR OUTPUT STATEMENT......
 200   FORMAT('  N TOO BIG')
C
       END