program lapack_test_matinv

    implicit none

    integer, parameter :: n = 3
    complex(8), dimension(n, n) :: A, LU, Ainv, E
    integer :: i, j

    integer :: clock = 0, rsize
    integer, allocatable :: rseed(:)
    real :: rharvest(2)

    complex(8) :: iwork(n)
    integer:: info_lu, info_inv, ipiv(n)

    ! Add some random elements to the matrix.
    call random_seed(size=rsize)
    allocate(rseed(rsize))
    ! call system_clock(count=clock)
    rseed = clock + 71 * (/(i, i = 0, rsize - 1)/)
    call random_seed(put=rseed)
    do i = 1, n
        do j = 1, n
            call random_number(rharvest)
            A(i, j) = rharvest(1) + rharvest(2) * (0.0, 1.0)
        end do
    end do
    write(*, fmt="(a)") "A = "
    call write_mat(A, n, n)

    ! Factorize into L\U matrix.
    LU = A
    call zgetrf(n, n, LU, n, ipiv, info_lu)

    write(*, fmt="(a)") "pivoting: "
    do i = 1, n
        write(*, fmt="(a6,i4,i4)") "swap", i, ipiv(i)
    enddo
    write(*, fmt="(a)") "L\U = "
    call write_mat(LU, n, n)

    ! Invert using LU factorization.
    Ainv = LU
    call zgetri(n, Ainv, n, ipiv, iwork, n, info_inv)

    write(*, fmt="(a)") "A^-1 = "
    call write_mat(Ainv, n, n)

    ! Check if the product is unit matrix.
    E = matmul(A, Ainv)

    write(*, fmt="(a)") "A * A^-1 = "
    call write_mat(E, n, n)

end program lapack_test_matinv


subroutine write_mat (A, n, m)

    implicit none

    integer, intent(in) :: n, m
    complex(8), intent(in) :: A(n, m)

    integer :: i, j

    do i = 1, n
        do j = 1, m
            write(*, fmt="(a,(f8.4,f8.4),a)", advance="no") " (", A(i, j), "*i)"
        end do
        write(*,*)
    end do

end subroutine write_mat