!***********************************************************************************
!	MODULE:	MatInv
!
!	PURPOSE: Contains subroutines that do matrix inversion
!
!	AUTHOR:	Ben Malin ... taken from Tao Pang
!
!	DATE:	8/14/07
!       REVISED: 
!***********************************************************************************

module MatInv

  implicit none
  
contains
  
  !****************************************************************************
  !
  !	SUBROUTINE:  ELGS
  !
  !	PURPOSE: 
  ! Subroutine to perform the partial-pivoting Gaussian elimination.
  ! A(N,N) is the original matrix in the input and transformed matrix
  ! plus the pivoting element ratios below the diagonal in the output.
  ! INDX(N) records the pivoting order.  Copyright (c) Tao Pang 2001.
  !
  !	AUTHOR:	 Tao Pang
  !
  !	DATE:    10/24/2001
  !     
  !****************************************************************************
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                       !
  ! Please Note:                                                          !
  !                                                                       !
  ! (1) This computer program is written by Tao Pang in conjunction with  !
  !     his book, "An Introduction to Computational Physics," published   !
  !     by Cambridge University Press in 1997.                            !
  !                                                                       !
  ! (2) No warranties, express or implied, are made for this program.     !
  !                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE ELGS (A,N,INDX)

    IMPLICIT NONE
    INTEGER, INTENT (IN) :: N
    INTEGER :: I,J,K,ITMP
    INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
    REAL :: C1,PI,PI1,PJ
    REAL, INTENT (INOUT), DIMENSION (N,N) :: A
    REAL, DIMENSION (N) :: C
    !
    ! Initialize the index
    !
    DO I = 1, N
       INDX(I) = I
    END DO
    !
    ! Find the rescaling factors, one from each row
    !
    DO I = 1, N
       C1= 0.0
       DO J = 1, N
          C1 = AMAX1(C1,ABS(A(I,J)))
       END DO
       C(I) = C1
    END DO
    !
    ! Search the pivoting (largest) element from each column
    !
    DO J = 1, N-1
       PI1 = 0.0
       DO I = J, N
          PI = ABS(A(INDX(I),J))/C(INDX(I))
          IF (PI.GT.PI1) THEN
             PI1 = PI
             K   = I
          ENDIF
       END DO
       !
       ! Interchange the rows via INDX(N) to record pivoting order
       !
       ITMP    = INDX(J)
       INDX(J) = INDX(K)
       INDX(K) = ITMP
       DO I = J+1, N
          PJ  = A(INDX(I),J)/A(INDX(J),J)
          !
          ! Record pivoting ratios below the diagonal
          !
          A(INDX(I),J) = PJ
          !
          ! Modify other elements accordingly
          !
          DO K = J+1, N
             A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
          END DO
       END DO
    END DO
    
  END SUBROUTINE ELGS

  !****************************************************************************
  !
  !	SUBROUTINE:  MIGS
  !
  !	PURPOSE: 
  ! Subroutine to invert matrix A(N,N) with the inverse stored
  ! in X(N,N) in the output.  Copyright (c) Tao Pang 2001.
  !
  !	AUTHOR:	 Tao Pang
  !
  !	DATE:    10/24/2001
  !     
  !****************************************************************************
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !                                                                       !
  ! Please Note:                                                          !
  !                                                                       !
  ! (1) This computer program is written by Tao Pang in conjunction with  !
  !     his book, "An Introduction to Computational Physics," published   !
  !     by Cambridge University Press in 1997.                            !
  !                                                                       !
  ! (2) No warranties, express or implied, are made for this program.     !
  !                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE MIGS (A,N,X,INDX)

    IMPLICIT NONE
    INTEGER, INTENT (IN) :: N
    INTEGER :: I,J,K
    INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
    REAL, INTENT (IN), DIMENSION (N,N):: A
    REAL, INTENT (OUT), DIMENSION (N,N):: X
    REAL, DIMENSION (N,N) :: B
    real, dimension(N,N) :: A1   !Added by bmalin
    !
    A1 = A
    !
    DO I = 1, N
       DO J = 1, N
          B(I,J) = 0.0
       END DO
    END DO
    DO I = 1, N
       B(I,I) = 1.0
    END DO
    !
    CALL ELGS (A1,N,INDX)
    !
    DO I = 1, N-1
       DO J = I+1, N
          DO K = 1, N
             B(INDX(J),K) = B(INDX(J),K)-A1(INDX(J),I)*B(INDX(I),K)
          END DO
       END DO
    END DO
    !
    DO I = 1, N
       X(N,I) = B(INDX(N),I)/A1(INDX(N),N)
       DO J = N-1, 1, -1
          X(J,I) = B(INDX(J),I)
          DO K = J+1, N
             X(J,I) = X(J,I)-A1(INDX(J),K)*X(K,I)
          END DO
          X(J,I) =  X(J,I)/A1(INDX(J),J)
       END DO
    END DO

  END SUBROUTINE MIGS
  
end module MatInv
