C
C SUBROUTINE 'INVMAT'
C ===================
C
C INVERT A SQUARE MATRIX
C
C (INVERSION BY ELIMINATION WITH PARTIAL PIVOTING)
C
      SUBROUTINE INVMAT(A,AI,N,IDIM,MSG,DET,IERR)
C
C PARAMETERS
C
C       A (I/O) MATRIX TO BE INVERTED A(I,J), I= 1 TO N, AND J = 1 TO N
C               (DOUBLE PRECISION)   NOTE: THE ARRAY IS MODIFIED DURING
C               THE COMPUTATION.
C      AI (O)   INVERTED MATRIX (DOUBLE PRECISION)
C       N (I)   SIZE OF SQUARE MATRIX
C    IDIM (I)   FIRST DIMENSION OF ARRAYS A AND AI (MUST BE THE SAME FOR
C               EACH ARRAY)
C     MSG (I)   LOGICAL UNIT NUMBER FOR ERROR MESSAGES ,0 IF NONE DESIRED
C     DET (O)   RETURNS DETERMINANT VALUE (DOUBLE PRECISION)
C    IERR (O)   ERROR FLAG =0 OK
C                          =1 SINGULAR MATRIX
C
C
C SPECIFICATION STATEMENTS
C ------------------------
C
      DOUBLE PRECISION A(IDIM,N),AI(IDIM,N),DET
      DOUBLE PRECISION EPS,AMAX,ATMP,AITMP,DIV,AMULT
C
C PERFORM MATRIX INVERSION
C ------------------------
C
C INITIALISATIONS, CONSTRUCT IDENTITY MATRIX
C
      IERR=0
      EPS=1.0E-12
      DO 20 I=1,N
      DO 10 J=1,N
      AI(I,J)=0.0
10    CONTINUE
      AI(I,I)=1.0
20    CONTINUE
      DET=1.0
      DO 200 K=1,N
C
C LOCATE MAX MAGNITUDE ON OR BELOW MAIN DIAGONAL
C
      IF(K.LT.N)THEN
        IMAX=K
        AMAX=ABS(A(K,K))
        DO 40 I=K+1,N
        IF(AMAX.LT.ABS(A(I,K)))THEN
          IMAX=I
          AMAX=ABS(A(I,K))
        ENDIF
40      CONTINUE
C
C INTERCHANGE ROWS IMAX AND K IF IMAX.NE.K
C
        IF(IMAX.NE.K)THEN
          DO 50 J=1,N
          ATMP=A(IMAX,J)
          A(IMAX,J)=A(K,J)
          A(K,J)=ATMP
          AITMP=AI(IMAX,J)
          AI(IMAX,J)=AI(K,J)
          AI(K,J)=AITMP
50        CONTINUE
          DET=-DET
        ENDIF
      ENDIF
C
C TEST FOR SINGULAR MATRIX
C
      IF(ABS(A(K,K)).LT.EPS)GO TO 800
      DET=A(K,K)*DET
C
C DIVIDE PIVOT ROW BY ITS MAIN DIAGONAL ELEMENT
C
      DIV=A(K,K)
      DO 60 J=1,N
      A(K,J)=A(K,J)/DIV
      AI(K,J)=AI(K,J)/DIV
60    CONTINUE
C
C REPLACE EACH ROW BY LINEAR COMBINATION WITH PIVOTAL ROW
C
      DO 100 I=1,N
      AMULT=A(I,K)
      IF(I.EQ.K)GO TO 100
      DO 80 J=1,N
      A(I,J)=A(I,J)-AMULT*A(K,J)
      AI(I,J)=AI(I,J)-AMULT*AI(K,J)
80    CONTINUE
100   CONTINUE
C
C END OF MAIN LOOP
C
200   CONTINUE
      RETURN
C
C ERROR CONDITIONS
C ----------------
C
800   IF(MSG.GT.0)WRITE(MSG,2001)K
      IERR=1
      RETURN
C
C FORMAT STATEMENTS
C -----------------
C
2001  FORMAT(/,'**SINGULAR MATRIX FOR ELEMENT',I5,' **')
      END