!     Last change:  J     8 Oct 2000   11:48 am
!                                       *** copyright 2000 ***
! *** filename mirse.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PMIRSE
!  SOLVE SYSTEM OF LINEAR EQUATIONS  A * X = B  AND COMPUTE INVERSE
!   OF MATRIX A BY FULL ELIMINATION WITH COMPLETE PIVOTING
!  MIRSE DOES BOTH
!
!                     DECLARATIONS
      IMPLICIT NONE
      INTEGER, PARAMETER :: LENGTH = 10
      REAL, DIMENSION(LENGTH,LENGTH) :: A
      REAL, DIMENSION(LENGTH) :: B
      REAL D
      INTEGER I,J,N                                                    
!                     FORMAT STATEMENTS
  20  FORMAT(2X,4F6.2)                                                  
  21  FORMAT(1X,10F13.8)                                                
  22  FORMAT(2X,'DETERMINANT ',F15.8)                                    
  26  FORMAT('  INPUT MATRIX ')
  27  FORMAT('  INVERSE MATRIX  ')
  28  FORMAT('  RIGHT HAND SIDE ')
  29  FORMAT('  SOLUTION VECTOR ')
!                   INTERFACE BLOCK
      INTERFACE
      SUBROUTINE MIRSE(A,N,DET,V)
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, DIMENSION(:), INTENT(IN OUT) :: V
      REAL, INTENT(OUT) :: DET
      INTEGER, INTENT(IN) :: N
      END SUBROUTINE MIRSE
      END INTERFACE
!
      OPEN( UNIT=5, FILE='gauspp.dat')
      OPEN (UNIT=6, FILE='mirse.out' )
!                     READ THE MATRIX A
      N = 4                                                             
      DO I = 1,N
      READ(5,20) (A(I,J),J=1,N)                                         
      ENDDO
      WRITE(6,26)
      DO I = 1,N
      WRITE(6,21) (A(I,J),J=1,N)                                        
      ENDDO
!                     READ IN THE RHS B
      READ(5,20) (B(J),J=1,N)                                           
      WRITE(6,28) 
      WRITE(6,21) (B(J),J=1,N)                                          
!
!                     CALL MIRSE TO SOLVE AND COMPUTE INVERSE
      CALL MIRSE(A,N,D,B)
!
      WRITE(6,22) D
      WRITE(6,27) 
      DO I = 1,N
      WRITE(6,21) (A(I,J),J=1,N)                                        
      ENDDO
!
      WRITE(6,29) 
      WRITE(6,21) (B(J),J=1,N)                                          
!
      STOP                                                              
      END PROGRAM PMIRSE
      SUBROUTINE MIRSE(A,N,DET,V)                                       
!   OVERWRITES THE N BY N MATRIX A WITH ITS INVERSE AND
!   OVERWRITES THE RHS VECTOR B WITH THE SOLUTION OF A * X = B
!   THE DETERMINANT OF THE MATRIX A IS RETURNED IN DET
!
!   THE METHOD IS FULL ROW ELIMINATION WITH FULL PIVOTING
!
!   THE MATRIX AND VECTOR ARE DIMENSIONED BELOW FOR CONVENIENCE ONLY
!
!   J F MONAHAN   DEPT OF STATISTICS, N C S U, RALEIGH, N C 27695 - 8203
!            CODE    CIRCA 1973
!            COMMENTS  SEPTEMBER  1984
!            RECODED FOR FORTRAN 95  SEPTEMBER 1999
!
!                           DECLARATIONS
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, DIMENSION(:), INTENT(IN OUT) :: V
      REAL, INTENT(OUT) :: DET
      INTEGER, INTENT(IN) :: N
      INTEGER I,J,L,Q,R
      INTEGER SR(N),SC(N)
      REAL T,BIG
!                          INITIALIZE PERMUTATION VECTORS AND DET
      DO I = 1,N
      SR(I) = 0
      SC(I) = 0
      ENDDO
      DET = 1.
!                          DO N PIVOTS COUNTED BY THE INDEX L
      DO L=1,N
!
!               FIND THE LARGEST PIVOTABLE ELEMENT
!
      BIG = 0.
      T = 0.
      DO I = 1,N
!                          HAVE WE PIVOTED IN THIS ROW YET?
      IF( SR(I) .EQ. 0 ) THEN
        DO J = 1,N
!                          HAVE WE PIVOTED IN THIS COLUMN YET?
        IF( (SC(J).EQ.0) .AND. (ABS(A(I,J)).GT.BIG) ) THEN
!                          PIVOT
          T = A(I,J)
          BIG = ABS(T)
          Q = I
          R = J
          ENDIF
        ENDDO       ! LOOP ON J
      ENDIF
      ENDDO         ! LOOP ON I
!
!                           PIVOT IN ROW Q AND COLUMN R
!
!                           UPDATE DETERMINANT  ** OVERFLOW DANGER HERE
      DET = DET*T
!                           IF PIVOT IS ZERO THEN MTHE MATRIX IS
!                           COMPUTATIONALLY SINGULAR SO RETURN
!                           WITH A ZERO DETERMINANT
  21  format(/' step  ',i4,'  q,r',2i3,'   pivot',f12.6/'matrix ')
      write(6,21) l,q,r,a(q,r)
      IF( BIG .LE. 0. ) RETURN
      SR(Q) = R
      SC(R) = Q
      T = 1./T
!                           ADJUST ROW Q
      DO J = 1,N
      A(Q,J) = A(Q,J)*T
      ENDDO
      A(Q,R) = T
      V(Q) = V(Q)*T
!                           UPDATE MATRIX, INVERSE, AND SOLUTION
      DO J = 1,N
      IF( Q .NE. J ) THEN
        T = A(J,R)
        A(J,R) = 0.
        DO I = 1,N
        A(J,I) = A(J,I) - T*A(Q,I)
        ENDDO         ! LOOP ON I
        V(J) = V(J) - T*V(Q)
        ENDIF
      ENDDO         ! LOOP ON J
!
!                           HAVE REDUCED ORIGINAL MATRIX TO PERMUTATION
!                           MATRIX SO PERMUTE BACK BY COLUMNS FIRST
!
  22  format(6f12.6)
      do i = 1,n
      write(6,22) (a(i,j),j=1,n)
      enddo
!
      ENDDO         ! LOOP ON L
      DO I = 1,N
      DO WHILE ( SC(I) .NE. I)
      Q = SC(I)
      DO J = 1,N
        T = A(J,I)
        A(J,I) = A(J,Q)
        A(J,Q) = T
      ENDDO         ! LOOP ON J
      SC(I) = SC(Q)
      SC(Q) = Q
      ENDDO         ! WHILE
      ENDDO         ! LOOP ON I

!                           NOW PERMUTE MATRIX AND SOLUTION BY ROWS
      DO I = 1,N
      DO WHILE ( SR(I) .NE. I )
      Q = SR(I)
      DO J = 1,N
        T = A(I,J)
        A(I,J) = A(Q,J)
        A(Q,J) = T
      ENDDO         ! LOOP ON J
      T = V(I)
      V(I) = V(Q)
      V(Q) = T
!                           FIX THE SIGN OF THE DETERMINANT
      DET = -DET
      SR(I) = SR(Q)
      SR(Q) = Q
      ENDDO         ! WHILE
      ENDDO         ! LOOP ON I
      RETURN
      END SUBROUTINE MIRSE
! *** end of filename mirse.f95         *********************