Returns the inverse of a matrix calculated by finding the LU decomposition. Depends on LAPACK.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=rk), | intent(in), | dimension(:,:) | :: | A |
Matrix to invert |
|
integer, | intent(out), | optional | :: | errCode |
If error code is present return error code and do not abort |
inverse of A
function invert_matrix(A, errCode) result(Ainv)
! ---------------------------------------------------------------------------
!> Matrix to invert
real(kind=rk), dimension(:,:), intent(in) :: A
!> If error code is present return error code and do not abort
integer, optional, intent(out) :: errCode
!> inverse of A
real(kind=rk), dimension(size(A,1),size(A,2)) :: Ainv
! ---------------------------------------------------------------------------
real(kind=rk), dimension(size(A,1)) :: work ! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info
! External procedures defined in LAPACK
external DGETRF
external DGETRI
! ---------------------------------------------------------------------------
! Store A in Ainv to prevent it from being overwritten by LAPACK
Ainv = A
n = size(A,1)
! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
call DGETRF(n, n, Ainv, n, ipiv, info)
if (info /= 0) then
if (present(errCode)) then
errCode = info
write(logUnit(5),*) 'WARNING: Matrix is numerically singular!'
return
else
call tem_abort('Matrix is numerically singular!')
end if
end if
! DGETRI computes the inverse of a matrix using the LU factorization
! computed by DGETRF.
call DGETRI(n, Ainv, n, ipiv, work, n, info)
if (info /= 0) then
if (present(errCode)) then
errCode = info
write(logUnit(5),*) 'WARNING: Matrix inversion failed!'
return
else
call tem_abort('Matrix inversion failed!')
end if
end if
! successfull
if (present(errCode)) errCode = 0
end function invert_matrix