! test_inverse.f90

subroutine inverse(n,A,AA)
  implicit none
  integer, intent(in) :: n
  double precision, dimension(n,n), intent(in)  :: A
  double precision, dimension(n,n), intent(out) :: AA
  double precision, dimension(n) :: TEMP
  integer ROW(n), COL(n)
  integer HOLD, I_PIVOT, J_PIVOT
  double precision PIVOT, ABS_PIVOT, NORM1
  integer i, j, k
  
  NORM1=0.0d0

  do i=1,n
    do j=1,n
      AA(i,j)=A(i,j)
    end do
  end do

  do I=1,n
    do J=1,n
      if(abs(AA(I,J))> NORM1) NORM1=abs(AA(I,J))
    end do 
  end do 

  do K=1,n
    ROW(K)=K 
    COL(K)=K 
  end do 

  do K=1,n
    PIVOT=AA(ROW(K),COL(K))
    I_PIVOT=K 
    J_PIVOT=K 
    do I=K,N
      do J=K,N
        ABS_PIVOT=abs(PIVOT)
        if (abs(AA(ROW(I),COL(J)))> ABS_PIVOT) then
          I_PIVOT=I 
          J_PIVOT=J 
          PIVOT=AA(ROW(I),COL(J))
        end if 
      end do 
    end do 

    if (ABS_PIVOT < EPSILON(NORM1) * NORM1) then
      print *, 'matrix is singular'
      return
    end if 

    HOLD=ROW(K)
    ROW(K)=ROW(I_PIVOT)
    ROW(I_PIVOT)=HOLD 
    HOLD=COL(K)
    COL(K)=COL(J_PIVOT)
    COL(J_PIVOT)=HOLD 

    AA(ROW(K),COL(K))=1.0d0/PIVOT 
    do J=1,n
      if (J .ne. K) then
        AA(ROW(K),COL(J))=AA(ROW(K),COL(J))* AA(ROW(K),COL(K))
      end if 
    end do 

    do I=1,n
      if (K .ne. I) then
        do J=1,n
          if (K .ne. J) then
            AA(ROW(I),COL(J))=AA(ROW(I),COL(J))- AA(ROW(I),COL(K))* AA(ROW(K),COL(J))
          end if 
        end do 
        AA(ROW(I),COL(K))=-AA(ROW(I),COL(K))* AA(ROW(K),COL(K))
      end if 
    end do 
  end do 

  do J=1,n
    do I=1,n
      TEMP(COL(I))=AA(ROW(I),J)
    end do 
    do I=1,n
      AA(I,J)=TEMP(I)
    end do 
  end do 

  do I=1,n
    do J=1,n
      TEMP(ROW(J))=AA(I,COL(J))
    end do 
    do J=1,n
      AA(I,J)=TEMP(J)
    end do 
  end do 
end subroutine inverse

subroutine matmul(n, A, B, C)  ! C = A * B
  implicit none
  integer, intent(in)                            :: n
  double precision, dimension(n,n), intent (in)  :: A
  double precision, dimension(n,n), intent (in)  :: B
  double precision, dimension(n,n), intent (out) :: C
  integer                                        :: i, j, k

  ! n  A,B,C  n by n
  do i = 1,n
    do j = 1,n
      C(i,j) = 0.0d0
      do k = 1,n 
        C(i,j) = C(i,j) + A(i,k)*B(k,j)
      end do ! k 
    end do ! j 
  end do ! i
end subroutine matmul

subroutine prtmat(nn, name, n, A) ! nn = number of characters in name, A(n,n)
  implicit none
  integer, intent(in) :: nn
  character (len=nn), intent(in) :: name
  integer, intent(in) :: n
  double precision, dimension(n,n), intent (in)  :: A
  integer i, j
  do i = 1,n
    do j = 1,n
      print *, name,'(',i,',',j,') = ',A(i,j)
    end do ! j
  end do ! i
  print *, ' '
end subroutine prtmat

program main
  implicit none
  interface
     subroutine inverse(n, A, AA)
       integer, intent(in) :: n
       double precision, dimension(n,n), intent(in)  :: A
       double precision, dimension(n,n), intent(out) :: AA
     end subroutine inverse

     subroutine matmul(n, A, B, C)  ! C = A * B
       integer, intent(in) :: n
       double precision, dimension(n,n), intent(in)  :: A
       double precision, dimension(n,n), intent(in)  :: B 
       double precision, dimension(n,n), intent(out) :: C
     end subroutine matmul

     subroutine prtmat(nn, name, n, A)
       integer, intent(in) :: nn
       character (len=nn), intent(in) :: name
       integer, intent(in) :: n
       double precision, dimension(n,n), intent (in)  :: A
     end subroutine prtmat
  end interface
  
  integer, parameter   :: n = 4
  integer i, j
  double precision, dimension(4,4) :: A
  data A / 4.0, 4.0, 4.0, 4.0,  4.0, 3.0, 3.0, 3.0, &
           4.0, 3.0, 2.0, 2.0,  4.0, 3.0, 2.0, 1.0 /

  double precision, dimension(4,4) :: A2
  data A2 / 1.0, 3.0, 2.0, 5.0,  4.0, 4.0, 3.0, 3.0, &
            5.0, 1.0, 1.0, 2.0,  1.0, 4.0, 3.0, 2.0 /

  double precision, dimension(4,4) :: AH
  ! data AH / (1.0/2.0), (1.0/3.0), (1.0/4.0), (1.0/5.0), &
  !           (1.0/3.0), (1.0/4.0), (1.0/5.0), (1.0/6.0), &
  !           (1.0/4.0), (1.0/5.0), (1.0/6.0), (1.0/7.0), &
  !           (1.0/5.0), (1.0/6.0), (1.0/7.0), (1.0/8.0)  /

  double precision, dimension(n,n) :: AI   ! computed
  double precision, dimension(n,n) :: AII  ! CHECK

  print *, 'test_inverse.f90 running'
  print *, 'A=  ', A
  call prtmat(1,'A',4,A)
  call inverse(n, A, AI)
  print *, 'AI= ', AI
  print *, ' '
  call inverse(n, AI, AII)
  print *, 'A = AII= ', AII
  print *, ' '
  call matmul(n, A, AI, AII)
  print *, 'I = AII= ', AII
  print *, ' '

  print *, 'A2= ', A2
  call inverse(n, A2, AI)
  print *, 'AI= ', AI
  print *, ' '
  call inverse(n, AI, AII)
  print *, 'A2 = AII= ', AII
  print *, ' '
  call matmul(n, A2, AI, AII)
  print *, 'I = AII= ', AII
  print *, ' '

  do i = 1,n
    do j = 1,n
      AH(i,j) = 1.0d0/((i+0.0d0)+(j+0.0d0));
    end do ! j
  end do ! i
  print *, 'AH= ', AH
  call inverse(n, AH, AI)
  print *, 'AI= ', AI
  print *, ' '
  call inverse(n, AI, AII)
  print *, 'AH = AII= ', AII
  print *, ' '
  call matmul(n, AH, AI, AII)
  print *, 'I = AII= ', AII
  print *, ' '
  call prtmat(3,'AII',4,AII)
  print *, ' '
  
  print *, 'test_inverse.f90 finished'
end