sistuspf.f90

subroutine sistuspf(a, b, ip)

!Resolución del sistema triangular superior salvo permutación de filas ip
!Entrada:	a (triangular inferior: L y triangular superior (sin la diagonal de 1): U)
!		b (solución de L*u=Pb, donde ahora u=b)
!		ip (puntero de las permutaciones)
!Salida:	b (solución del sistema)


implicit none
real(8), intent (in) :: a(:,:)
real(8), intent (inout) :: b(:)
integer, intent (in) :: ip(:)
real(8) :: sum
integer :: i, n, j

n = size(b)

do i = n, 1, -1
	sum = 0d0
	do j = i+1, n
		sum = sum+a(ip(i), j)*b(j)
	end do
	b(i) = (b(i)-sum)/a(ip(i), i)
end do

return
end subroutine sistuspf