C***********************************************************************
C   THIS SUBROUTINE HAS BEEN WRITTEN SPECIFICALLY FOR SYMMETRIC POSITIVE
C   DEFINITE MATRICES OF SMALL TO MODERATE SIZE.
C   IT WILL WORK ON CERTAIN OTHER SYMMETRIC MATRICES,BUT WAS NOT WRITTEN
C   WITH THIS IN MIND
C
C   CALLING ARGUMENTS:
C   A  THE M BY M ARRAY TRANSFERED TO THE SUBROUTINE.  THE ACTUAL MATRIX
C      TO BE INVERTED MAY BE OF SMALLER SIZE, AS IT MAY BE ONLY AN UPPER
C      SQUARE OF THE TRANSFERED ARRAY.
C   N  THE SIZE OF THE MATRIX TO BE INVERTED. THUS THE MATRIX TO BE
C     INVERTED IS ONLY THE UPPER N BY N SUB-MATRIX OF THE M BY M ARRAY A
C        THIS FEATURE ALLOWS MORE FLEXIBILITY IN THE USE OF ARRAYS IN
C        THE MAIN PROGRAM.
C   M    THE ACTUAL DIMENSION OF THE SQUARE ARRAY 'A'. IT IS REQUIRED
C        THAT M IS GREATER THAN OR EQUAL TO N.
C   ACC  A NUMBER RETURNED BY THE ROUTINE TO INDICATE VALIDITY OF RESULT
C          -1. MATRIX DOES NOT APPEAR SYMMETRIC
C           0. ACCURACY IS NOT ASSURED. IN THIS CASE INVERSION WAS AT-
C           TEMPTED BUT PROBLEMS WERE ENCOUNTERED. THE MATRIX MAY BE
C           SINGULAR
C          +1.  MATRIX INVERTED WITH GOOD ACCURACY.
C
C   THE BASIC ALGORITHM IS SIMPLE AND BY ITSELF IS NOT TOO ACCURATE.
C   HOWEVER, RESULTS ARE CHECKED FOR ACCURACY AND ITERATIVE IMPROVEMENTS
C   ARE DONE AS NEEDED.  THIS IMPROVES THE ACCURACY OF RESULTS CON-
C   SIDERABLY.   THE BASIC ALGORITHM IS AN ESCALATOR METHOD BASED ON
C   RESULTS IN PARTITIONED MATRIX THEORY. IT IS ALSO KNOWN AS BORDERING
C   AND INVERTING REFERENCE: HOUSEHOLDER,1957. "A SURVEY OF SOME CLOSED
C   METHODS FOR INVERTING MATRICES" J. SOC. INDUST. APPL. MATH.  VOL. 5,
C   PAGES 159-160.
C   BECAUSE A IS SYMMETRIC IT HAS ONLY M*(M+1)/2 DISTINCT ELEMENTS OUT
C   OF M*M PLACES. CONSEQUENTLY THERE IS WASTED SPACE IN A SENSE IN THE
C   ARRAYS A, B, C, D. HOWEVER CORE SPACE IS USUALLY NOT A PROBLEM,
C   ESPECIALLY FOR VALUES OF N, M LESS THAN 20. THUS THE ONLY THING THAT
C   MIGHT BE GAINED BY TAKING SYMMETRY MORE FULLY INTO CONSIDERATION IS
C   SPEED. BUT IT IS NOT CLEAR THAT SINGLY DIMENSIONING THE ARRAYS WOULD
C   IMPROVE THE SPEED OF THE SUBROUTINE. IT WOULD MAKE IT VERY CLUMSY
C   THOUGH.
C
C   SUBROUTINES CALLED: MINV
C***********************************************************************
      SUBROUTINE INVERT (A,N,M,ACC)
      INCLUDE 'PARMTR.INC'
C***********************************************************************
C     DECLARATIONS
C***********************************************************************
      DOUBLE PRECISION A, B, C, D, SCALE, Q, TEST, Z
      DIMENSION A(M,M), B(MAXPAR,MAXPAR), C(MAXPAR,MAXPAR),
     1   D(MAXPAR,MAXPAR), SCALE(MAXPAR)
      ACC=-1.
      IF (N.NE.1) GO TO 10
      ACC=1.
      A(1,1)=1./A(1,1)
      RETURN
C
C   CHECK FOR SYMMETRY
C
   10 DO 20 I=2,N
      I1=I-1
      DO 20 J=1,I1
      Q=A(I,J)
      Z=Q
      IF (Q.EQ.0.) Z=1.
      TEST=(Q-A(J,I))/Z
      IF (TEST.GT..1D-7) RETURN
   20 CONTINUE
      ACC=1.
C
C   FIRST A IS SCALED THEN SET EQUAL TO B. IT IS THE SCALED
C        MATRIX THAT IS INVERTED THEN RESCALED.
C
      DO 30 I=1,N
   30 SCALE(I)=1./DSQRT(DABS(A(I,I)))
      DO 40 I=1,N
      DO 40 J=1,N
      A(I,J)=SCALE(I)*A(I,J)*SCALE(J)
   40 B(I,J)=A(I,J)
      CALL MINV (B,N)
C
C   THIS RETURNS THE INVERSE OF B
C
      IC=0
   50 Z=0.
      IC=IC+1
      DO 70 I=1,N
      DO 70 J=1,I
      C(I,J)=0.
      DO 60 L=1,N
   60 C(I,J)=C(I,J)+A(I,L)*B(L,J)
      C(J,I)=C(I,J)
      Q=C(I,J)
      IF (I.EQ.J) Q=Q-1.
   70 Z=Z+DABS(Q)
      IF (Z.LT..000000001) GO TO 110
C
C   C IS A*B HENCE IT SHOULD BE CLOSE TO THE IDENTITY MATRIX.
C   IF NOT(AS MEASURED BY Z) FIND C INVERSE.
C
      CALL MINV (C,N)
      DO 90 I=1,N
      DO 90 J=1,I
      D(I,J)=0.
      DO 80 L=1,N
   80 D(I,J)=D(I,J)+B(I,L)*C(L,J)
   90 D(J,I)=D(I,J)
      DO 100 I=1,N
      DO 100 J=1,N
  100 B(I,J)=D(I,J)
      IF (IC.LT.2) GO TO 50
      ACC=0.
  110 DO 120 I=1,N
      DO 120 J=1,N
  120 A(I,J)=SCALE(I)*B(I,J)*SCALE(J)
      RETURN
      END