Calculates LU decomposition
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
real(kind=double), | intent(inout), | dimension(:,:) | :: | a | ||
integer, | intent(out), | dimension(:) | :: | indx | ||
real(kind=double), | intent(out) | :: | d |
subroutine ludcmp(a,indx,d)
!+Calculates LU decomposition
real(double), dimension(:,:),intent(INOUT):: a
integer,dimension(:), intent(OUT) :: indx
real(double), intent(OUT) :: d
real(double), dimension(size(a,1)) :: vv
integer,dimension(1) :: imaxloc
integer :: j,n,imax
n=size(indx)
d=1.0
vv=maxval(abs(a),dim=2)
if(any(vv.eq.0.))stop 'singular matrix in ludcmp'
vv=1.d0/vv
do j=1,n
imaxloc=maxloc(vv(j:n)*abs(a(j:n,j)))
imax=(j-1)+imaxloc(1)
if (j /= imax) then
call swap(a(imax,:),a(j,:))
d=-d
vv(imax)=vv(j)
endif
indx(j)=imax
if (a(j,j) == 0.0) a(j,j)=1.0d-20
a(j+1:n,j)=a(j+1:n,j)/a(j,j)
a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
enddo
end subroutine ludcmp