module fatoracao_LU

	use dados

	real*8 :: soma

	real*8,dimension(:),allocatable :: d

contains

!----------------------------------------------------------

subroutine fatoracao

	integer :: i, j, k

	allocate ( L(N,N), U(N,N) )

	L = 0; U = 0

	! Passo 1
	L(1,1) = 1.0d0
	U(1,1) = A(1,1)/L(1,1)
	if ( A(1,1) == 0.0d0 ) then
		write(*,*) "Fatoracao impossivel"
		stop
	end if

	! Passo 2
	do j = 2, N
		U(1,j) = A(1,j)/L(1,1) ! Primeira linha de U
		L(j,1) = A(j,1)/U(1,1) ! Primeira coluna de L
	end do

	! Passo 3
	do i = 2, N-1
		
		! Passo 4 (Definicao da diagonal principal de L e U)
		L(i,i) = 1.0d0
		soma = 0
		do k = 1, i-1
			soma = soma + L(i,k)*U(k,i)
		end do
		U(i,i) = ( A(i,i) - soma )/L(i,i)

		! Passo 5 
		do j = i+1, N

			soma = 0
			do k = 1, i-1
				soma = soma + L(i,k)*U(k,j)
			end do
			U(i,j) = ( A(i,j) - soma ) / L(i,i)

			soma = 0
			do k = 1, i-1
				soma = soma + L(j,k)*U(k,i)
			end do
			L(j,i) = ( A(j,i) - soma ) / U(i,i)

		end do

	end do

	! Passo 6
	L(N,N) = 1.0d0
	soma = 0
	do k = 1, N-1
		soma = soma + L(N,k)*U(k,N)
	end do
	U(N,N) = ( A(N,N) - soma ) / L(N,N)

end subroutine fatoracao

!----------------------------------------------------------

subroutine subs_progressiva

	integer :: i, j

	allocate ( d(N) )
	
	d(1) = b(1)/L(1,1)
	
	do i = 2, N
		soma = 0
		do j = 1, i-1
			soma = soma + L(i,j)*d(j)
		end do
		d(i) = ( b(i) - soma ) / L(i,i)

	end do

end subroutine subs_progressiva

!----------------------------------------------------------

subroutine subs_regressiva

	integer :: i, j

	x(N) = d(N) / U(N,N)

	do i = N-1, 1, -1
		soma = 0
		do j = i+1, N
			soma = soma + U(i,j)*x(j)
		end do
		x(i) = ( d(i) - soma ) / U(i,i)

	end do

end subroutine subs_regressiva

!----------------------------------------------------------
end module fatoracao_LU