subroutine Zinv(Am)
! Inverts the general complex matrix Am
complex(8), intent(inout) :: Am(:,:) ! Matrix to be inverted
complex(8), allocatable :: Amt(:,:), work(:)
integer :: n, nb
integer :: lwork, info
integer, allocatable :: ipiv(:)
n = size(Am, 1)
call assert_shape(Am, [n, n], "inv", "Am")
nb = ilaenv(1, 'ZGETRI', "UN", n, -1, -1, -1) ! TODO: check UN param
if (nb < 1) nb = max(1, n)
lwork = n*nb
allocate(Amt(n,n), ipiv(n), work(lwork))
Amt = Am
call zgetrf(n, n, Amt, n, ipiv, info)
if (info /= 0) then
print *, "zgetrf returned info =", info
if (info < 0) then
print *, "the", -info, "-th argument had an illegal value"
else
print *, "U(", info, ",", info, ") is exactly zero; The factorization"
print *, "has been completed, but the factor U is exactly"
print *, "singular, and division by zero will occur if it is used"
print *, "to solve a system of equations."
end if
stop 'Z_mat_invert: zgetrf error'
end if
call zgetri(n, Amt, n, ipiv, work, lwork, info)
if (info /= 0) then
print *, "zgetri returned info =", info
if (info < 0) then
print *, "the", -info, "-th argument had an illegal value"
else
print *, "U(", info, ",", info, ") is exactly zero; the matrix is"
print *, "singular and its inverse could not be computed."
end if
stop 'Z_mat_invert: zgetri error'
end if
Am = Amt
end subroutine Zinv