Matrix inversion with LU-decomposition
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=double), | intent(in), | dimension(:,:) | :: | a | ||
real(kind=double), | intent(out), | dimension(:,:) | :: | b |
subroutine matinv(a, b)
!+ Matrix inversion with LU-decomposition
!====================================================
real(double), dimension(:,:), intent(IN) :: a
real(double), dimension(:,:), intent(OUT) :: b
real(double), dimension(size(a,dim=1),size(a,dim=2)) :: ah, y
integer :: i, N
integer, dimension(size(a,dim=1)) :: indx
real(double) :: d
N = size(a,dim=1)
if (N /= size(a,dim=2)) stop 'SUB matinv: ludcmp matrix must be square!'
ah = a
y = 0.
do i = 1, N
y(i,i) = 1.d0
enddo
call ludcmp(ah,indx,d)
do i = 1, N
call lubksb(ah,indx,y(:,i))
enddo
b = y
end subroutine matinv