*  MATINV    MATrix INVersion subroutine
*   Larry G. Hegi          September 23, 1989

*  Note:  This subroutine is called by the main program in the files:
*            MCP309.FOR
*            MATINVC.FOR

*  This is a subprogram for inverting a matrix or solving a system of
*    simultaneous linear equations with REAL coefficients by the
*    Gauss-Jordan method.  At each stage the largest element of the
*    submatrix under consideration is found, rows and columns are
*    interchanged to make this the pivot element, and all elements
*    are divided by it.

*  This subroutine is given on page 306 of
*    John M. McCormick  &  Mario G. Salvadori
*    Numerical Methods in FORTRAN
*    Prentice Hall,   1964

*    Input:  A(N1,N1) = square matrix of coefficients  ( N1 x N1 )
*            N1 = rank of matrix A
*            B(N1,M1) = a constant-matrix of M1 column-vectors.
*                        ( b in the equation  A x = b )
*                        For simplicity, M1 can be set equal to 1.
*            M1 = the number of columns in the B-matrix of constants.
*            ND = Dimensioned number of rows in Matrix A(ND,ND)
*               = Dimensioned number of columns in Matrix A(ND,ND)
*               = Dimensioned number of rows in Matrix B(ND,MD)
*            MD = Dimensioned number of columns in Matrix B(ND,MD)

*    Output: A(N1,N1) = inverse of the input matrix A
*            B(N1,M1) = solution vector                -1
*                             ( x  in the equation x = A  b )
*            DETERM = Determinant of A
*            ID = Determinant Flag   ID = 1 if Determinant is non-zero
*                  ID = 2 if DETERM = 0  ( singular matrix A )

*  Note:  The array, INDEX, is included in the list of arguments for
*    MATINV for only 1 reason:  the size of the array is adjusted
*    according to the value of ND.  (The size of INDEX cannot be
*    adjusted unless it is included as an argument in the subroutine.)

*  General form of the Dimension Statement:
*  Dimension   A( , ),  B( , ),  INDEX( ,3)

      SUBROUTINE MATINV( A,N1,B,M1,DETERM,ID,ND,MD,INDEX )

      REAL  A(ND,ND), B(ND,MD)
      INTEGER  INDEX(ND,3)
      EQUIVALENCE  (IROW,JROW), (ICOLUM,JCOLUM), (AMAX, T, SWAP)


*****  Initialization  *****

      M = M1
      N = N1
  10  DETERM = 1.0

  15  DO 20  J = 1, N
  20     INDEX(J,3) = 0

  30  DO 550  I = 1, N

*****        Search for Pivot Element      *****

  40     AMAX = 0.0

  45     DO 105  J = 1, N
            IF ( INDEX(J,3) - 1 )  60, 105, 60

  60        DO 100  K = 1, N

* Note:  There may be a bug in the following line.  Check it carefully.

               IF ( INDEX(K,3) -1 )  80, 100, 80
  80           IF ( AMAX - ABS( A(J,K) ) )  85, 100, 100
  85           IROW = J
  90           ICOLUM = K
               AMAX = ABS( A(J,K) )
 100        CONTINUE
 105     CONTINUE

         INDEX( ICOLUM,3 ) = INDEX( ICOLUM,3 ) + 1
 260     INDEX( I,1 ) = IROW
 270     INDEX( I,2 ) = ICOLUM

*****       Interchange rows to put pivot element on diagonal    *****

 130     IF ( IROW-ICOLUM ) 140, 310, 140
 140        DETERM = -DETERM

 150     DO 200  L = 1, N
 160        SWAP = A( IROW,L )
 170        A( IROW,L ) = A( ICOLUM,L )
 200     A( ICOLUM,L ) = SWAP

         IF (M)  310, 310, 210

 210     DO 250  L = 1, M
 220        SWAP = B( IROW,L )
 230        B( IROW,L ) = B( ICOLUM,L )
 250     B( ICOLUM,L ) = SWAP

*****       Divide pivot row by pivot element     *****

 310     PIVOT = A( ICOLUM,ICOLUM )
         DETERM = DETERM * PIVOT
 330     A( ICOLUM,ICOLUM ) = 1.0

 340     DO 350  L = 1, N
 350        A( ICOLUM,L ) = A( ICOLUM,L ) / PIVOT

 355     IF (M) 380, 380, 360

 360     DO 370  L = 1, M
 370        B( ICOLUM,L ) = B( ICOLUM,L ) / PIVOT

*****      Reduce non-pivot rows     *****

 380     DO 550  L1 = 1, N
 390        IF ( L1 - ICOLUM ) 400, 550, 400
 400        T = A( L1,ICOLUM )
 420        A( L1,ICOLUM ) = 0.0

 430        DO 450  L = 1, N
 450           A( L1,L ) = A( L1,L ) - A( ICOLUM,L ) * T

 455        IF (M)  550, 550, 460

 460        DO 500  L = 1, M
 500           B( L1,L ) = B( L1,L ) - B( ICOLUM,L ) * T
 550  CONTINUE

*****     Interchange columns    *****

 600  DO 710  I = 1, N
 610     L = N + 1 - I
 620     IF ( INDEX( L,1 ) - INDEX( L,2 ) )  630, 710, 630
 630     JROW = INDEX( L,1 )
 640     JCOLUM = INDEX( L,2 )

 650     DO 705  K = 1, N
 660        SWAP = A( K,JROW )
 670        A( K,JROW ) = A( K,JCOLUM )
 700        A( K,JCOLUM ) = SWAP
 705     CONTINUE
 710  CONTINUE

      DO 730  K = 1, N
         IF ( INDEX(K,3) - 1 )  715, 720, 715
 715     ID = 2
         GO TO 740

 720     CONTINUE
 730  CONTINUE

      ID = 1

 740  RETURN

      END
