gausspiv.f90

subroutine gausspiv(a, b, ip, det)

!Eliminación de Gauss con estrategia de pivote parcial
!Entrada:	a (matriz de coeficientes)
!		b (términos independientes)
!Salida:	a (matriz triangular superior salvo permutación de filas)
!		b (términos independientes permutados)
!		ip (puntero de las permutaciones)
!		det (determinante del sistema)

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

n = size(b)

!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 = 0d0
!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 + 1d0
	else
		ipk = ip(k)
	end if
	!Actualización del determinante
	det = det*piv
	!Eliminación
	do i = k+1, n
		ipi = ip(i)
		factor = a(ipi, k)/piv
		do j = k+1,n
			a(ipi, j) = a(ipi, j)-factor*a(ipk, j)
		end do
		b(ipi) = b(ipi)-factor*b(ipk)
	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*(-1d0)**index
print*, "Determinante: ", det

return
end subroutine gausspiv