C-------------------------------------------------------------     ************
C                                                                      MXM
C                                                                  ************
      SUBROUTINE MXM(A,NAR,B,NAC,C,NBC)
C
C    Routine on VAX to emulate CRAY SCILIB routine, MXM p. 4-22
C         LIBRARY manual, implemented in VAX single precision.
C
      INTEGER NAR,NAC,NBC,I,J,K
      REAL A(NAR,1),B(NAC,1),C(NAR,1)

      DO 30 J=1,NBC
         DO 20 I = 1,NAR
            C(I,J) = 0.0
            DO 10 K = 1,NAC
               C(I,J) = C(I,J) + A(I,K)*B(K,J)
   10       CONTINUE
   20    CONTINUE
   30 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                      MXV
C                                                                  ************
      SUBROUTINE MXV(A,NAR,B,NBR,C)
C
C     Computes a matrix, A, times a vector, B, and assumes a skip
C          distance between elements of the matrix to be 1.
C
      INTEGER NAR,NBR,I,K
      REAL A(NAR,1),B(1),C(1)

      DO 20 I = 1,NAR
         C(I) = 0.0
         DO 10 K = 1,NBR
   10       C(I) = C(I) + A(I,K)*B(K)
   20 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                      MXMA
C                                                                  ************
      SUBROUTINE MXMA(A,NA,IAD,B,NB,IBD,C,NC,ICD,NAR,NAC,NBC)
C
C     Computes a matrix, A, times a matrix, B, and allows for
C          arbitrary spacing of matrix elements.
C     Martin J. McBride.  8/21/85.
C     General Electric CRD, Information System Operation.
C
      INTEGER NA,IAD,NB,IBD,NC,ICD,NAR,NAC,NBC
      INTEGER AC,AR,BC,BR,CC,CR,I,J,K
      REAL A(1,1),B(1,1),C(1,1)

      BC = 1
      CC = 1
      DO 30 J = 1,NBC
         CR = 1
         AR = 1
         DO 20 I = 1,NAR
            C(CR,CC) = 0.0
            AC = 1
            BR = 1
            DO 10 K = 1,NAC
               C(CR,CC) = C(CR,CC) + A(AR,AC)*B(BR,BC)
               AC = AC + IAD
               BR = BR + NB
   10       CONTINUE
            CR = CR + NC
            AR = AR + NA
   20    CONTINUE
         CC = CC + ICD
         BC = BC + IBD
   30 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                      MXVA
C                                                                  ************
      SUBROUTINE MXVA(A,NA,IAD,B,NB,C,NC,NAR,NBR)
C
C     Computes a matrix, A, times a vector, B, and allows for
C          arbitrary spacing of matrix elements.
C     Martin J. McBride.  8/21/85.
C     General Electric CRD, Information System Operation.
C
      INTEGER NA,IAD,NB,NC,NAR,NBR
      INTEGER AC,AR,IB,IC,I,K
      REAL A(1,1),B(1),C(1)

      IC = 1
      AR = 1
      DO 20 I = 1,NAR
         C(IC) = 0.0
         AC = 1
         IB = 1
         DO 10 K = 1,NBR
            C(IC) = C(IC) + A(AR,AC)*B(IB)
            AC = AC + IAD
            IB = IB + NB
   10    CONTINUE
         IC = IC + NC
         AR = AR + NA
   20 CONTINUE
      RETURN
      END

C-------------------------------------------------------------     ************
C                                                                      MINV
C                                                                  ************
      SUBROUTINE MINV(AB,N,ND,SCRATCH,DET,EPS,M,MODE)
C
C     A subroutine that calculates the determinant and inverse of
C          a matrix, as well as solving systems of linear equations.
C     Martin J. McBride.  11/25/85.
C     General Electric CRD, Information System Operation.
C
      INTEGER N,ND,M,MODE,OUTER,ROW,COL,I,SCOL,SROW,PIVCNT
      REAL AB(ND,1),SCRATCH(1),DET,EPS,MULT,COLNUM,TEMP

C  Initialize scratch space, with 1 to N holding the diagonal of the identity
C  matrix used to compute the inverse and N+1 to 2N holding the positions of
C  the first N columns of the matrix (for use when pivot occurs).
      DO 5 I = 1,N
    5    SCRATCH(I) = 1.0
      COLNUM = 1.0
      DO 6 I = N+1,2*N
         SCRATCH(I) = COLNUM
         COLNUM = COLNUM + 1.0
    6 CONTINUE

C  Make left, square matrix an upper triangular matrix.
      DET = 0.0
      PIVCNT = 0
      DO 10 OUTER = 1,N-1
         IF (ABS(AB(OUTER,OUTER)) .LE. EPS) THEN
            CALL PIVOT(AB,N,ND,OUTER,SCRATCH,EPS)
            IF (AB(OUTER,OUTER) .EQ. 0.0) THEN
               PRINT*
               PRINT*,'*************************************'
               PRINT*,'  MINV called with singular matrix.'
               PRINT*,'*************************************'
               PRINT*
               STOP
            ENDIF
            PIVCNT = PIVCNT + 1
         ENDIF
         DO 20 ROW = OUTER+1,N
            MULT = AB(ROW,OUTER)/AB(OUTER,OUTER)
            DO 30 COL = OUTER,N+M
   30          AB(ROW,COL) = AB(ROW,COL) - AB(OUTER,COL)*MULT
            DO 25 SCOL = 1,OUTER-1
   25          AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT
            AB(ROW,OUTER) = AB(ROW,OUTER) - SCRATCH(OUTER)*MULT
   20    CONTINUE
   10 CONTINUE

C  Compute determinant.
      DET = AB(1,1)
      DO 40 I = 2,N
   40    DET = DET*AB(I,I)
      DET = (-1.0)**PIVCNT * DET

C  Return if inverse is not to be found and there are no systems of equations
C  to solve.
      IF (MODE .EQ. 0 .AND. M .EQ. 0) RETURN

C  Place ones in diagonal of square matrix A.
      DO 80 ROW = 1,N
         DIV = AB(ROW,ROW)
         DO 90 COL = 1,N+M
            AB(ROW,COL) = AB(ROW,COL)/DIV
   90    CONTINUE
         SCRATCH(ROW) = SCRATCH(ROW)/DIV
   80 CONTINUE

C  Reduce upper triangle to zeros to give matrix A = I.
      DO 50 OUTER = 2,N
         DO 60 ROW = OUTER-1,1,-1
            MULT = AB(ROW,OUTER)/AB(OUTER,OUTER)
            DO 70 COL = OUTER,N+M
   70          AB(ROW,COL) = AB(ROW,COL) - AB(OUTER,COL)*MULT
            DO 65 SCOL = 1,ROW-1
   65          AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT
            SCRATCH(ROW) = SCRATCH(ROW) - AB(OUTER,ROW)*MULT
            DO 63 SCOL = ROW+1,OUTER-1
   63          AB(ROW,SCOL) = AB(ROW,SCOL) - AB(OUTER,SCOL)*MULT
            AB(ROW,OUTER) = AB(ROW,OUTER) - SCRATCH(OUTER)*MULT
   60    CONTINUE
   50 CONTINUE

C  Move diagonals of inverse to matrix AB.
      DO 85 I = 1,N
   85    AB(I,I) = SCRATCH(I)

C  If pivot was made, switch rows corresponding to the columns that were
C  pivoted.
      IF (PIVCNT .EQ. 0) RETURN
      ROW = 1
      DO 95 I = 1,N-1
         SROW = INT(SCRATCH(ROW+N))
         IF (SROW .NE. ROW) THEN
            DO 92 COL = 1,N+M
               TEMP = AB(ROW,COL)
               AB(ROW,COL) = AB(SROW,COL)
               AB(SROW,COL) = TEMP
   92       CONTINUE
            TEMP = SCRATCH(ROW+N)
            SCRATCH(ROW+N) = SCRATCH(SROW+N)
            SCRATCH(SROW+N) = TEMP
         ELSE
            ROW = ROW + 1
         ENDIF
   95 CONTINUE

      RETURN
      END

C--------------------------------------------------------------
      SUBROUTINE PIVOT(AB,N,ND,OUTER,SCRATCH,EPS)
C
C     This subroutine switches two columns of a matrix to get
C         a nonzero entry in the diagonal.
C     Martin J. McBride.  12/04/85.
C     General Electric CRD, Information System Operation.
C
      INTEGER N,ND,COL,OUTER,I
      REAL AB(ND,1),SCRATCH(1),TEMP,EPS

C  Get first column with non-zero element in row OUTER.
      COL = OUTER + 1
   10 IF (COL .GT. N) GO TO 90
      IF (ABS(AB(OUTER,COL)) .GT. EPS) GO TO 20
         COL = COL + 1
         GO TO 10

C  Switch column OUTER with column COL, which has non-zero element in
C  row OUTER.
   20 DO 30 I = 1,N
         TEMP = AB(I,OUTER)
         AB(I,OUTER) = AB(I,COL)
         AB(I,COL) = TEMP
   30 CONTINUE
      TEMP = SCRATCH(N+OUTER)
      SCRATCH(N+OUTER) = SCRATCH(N+COL)
      SCRATCH(N+COL) = TEMP

   90 CONTINUE
      RETURN
      END
