lupiv.f90

subroutine lupiv(a, ip, det)

!Obtención de la factorización LU aplicando el método del pivote parcial
!Entrada:	a (matriz de coeficientes)
!Salida:	a (triangular inferior: L (sin los 1 de la diagonal) y triangular superior: U)
!		ip (puntero de las permutaciones)
!		det (determinante de la matriz)

implicit none
real(8), intent (inout) :: a(:,:)
integer, intent (out) :: ip(:)
real(8), intent (out) :: det
real(8) :: piv
integer :: i, j, k, n, ipiv, ipi, ipk, index

n = size(ip)

!Inicialización del determinante
det = 1d0
!Inicialización de la permutación de filas
ip = (/(i, i=1, n)/)
!Inicialización del contador de cambios
index = 0
!Etapa k-ésima de la eliminación
do k = 1, n-1
	!Búsqueda del pivote y de la fila en la que se encuentra
	piv = a(ip(k), k)
	ipiv = k
	do i = k+1, n
		if(dabs(piv) < dabs(a(ip(i), k)))then
		piv = a(ip(i), k)
		ipiv = i
		end if
	end do
	!Comprobación de que el k-ésimo pivote no es nulo
	if(abs(piv) < 1d-12) then
		print*, "Pivote nulo en la etapa: ", k
		print*, "¡La matriz del sistema es Singular!"
		stop
	end if
	!Puesta al día de la permutación, si el pivote no está en la fila k
	if(ipiv /= k) then
		ipk = ip(ipiv)
		ip(ipiv) = ip(k)
		ip(k) = ipk
		index = index + 1
	else
		ipk = ip(k)
	end if
	!Actualización del determinante
	det = det*piv
	!Eliminación
	do i = k+1, n
		ipi = ip(i)
		a(ipi, k) = a(ipi, k)/piv
		do j = k+1,n
			a(ipi, j) = a(ipi, j)-a(ipi, k)*a(ipk, j)
		end do
	end do
end do
!Comprobación de que el último pivote no es nulo
piv = a(ip(n), n)
if(dabs(piv) < 1d-12) then
  print*, "Pivote nulo en la etapa: ", n
  print*, "¡La matriz del sistema es Singular!"
  stop
end if
!Finalización del cálculo del determinante
det = det*piv*(-1)**index
print*, "Determinante: ", det

return
end subroutine lupiv