!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Matrix Inverse by Gauss-Jordan Elimination  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Take in a matrix A, replace it with C so that the    !!
!!  original matrix is not destroyed.  Upon completion,  !!
!!  rearrange the elements of C into 1:n elements of B.  !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  ABW 2013 - 2015  !
!!!!!!!!!!!!!!!!!!!!!
module inv_mod
  implicit none
  
contains
  
  subroutine matinv(A,B,n)
    integer :: n
    real(kind=16),dimension(:,:) :: A
    real(kind=16),dimension(n,n) :: B
    real(kind=16),dimension(2*n,2*n) :: C
    integer :: i,j,k
    
    C = A
    
    do i = 1,n
       do j = 1,n
          if (i == j) then
             C(i,n+j) = 1.0q0
          else
             C(i,n+j) = 0.0q0
          endif
       enddo
    enddo
    
    do k = 1,n
       do i = 1,n
          if (i /= k) then
             do j = k+1,2*n
                C(i,j) = C(i,j) - (C(i,k)/C(k,k))*C(k,j)
             enddo
          endif
       enddo
    enddo
    
    do i = 1,n
       if (C(i,i) == 0.0q0) then
          write(*,*) 'No Inverse Exists'
          exit
       endif
    enddo
    
    do i = 1,n
       do j = n+1,2*n
          C(i,j) = C(i,j)/C(i,i)
          B(i,j-n) = C(i,j)
       enddo
    enddo
    
  end subroutine matinv
  
end module inv_mod
