subroutine pivot_total(A, P, Q)
real(dp), dimension(:, :), intent(in) :: A
real(dp), dimension(size(A, 1), size(A, 1)), intent(out) :: P, Q
real(dp), dimension(size(A, 1), size(A, 1)) :: P_tmp, Q_tmp
integer, dimension(2) :: vlmax
integer :: N, lmax, cmax, k
N = size(A, 1)
P = Identity_n(N)
Q = Identity_n(N)
do k = 1, N - 1
! Find max abs element in submatrix
vlmax = maxloc(abs(A(k:N, k:N)))
lmax = vlmax(1) + k - 1
cmax = vlmax(2) + k - 1
! permute line if necessary
P_tmp = Identity_n(N)
if (k /= lmax) then
P_tmp([k, lmax], :) = P_tmp([lmax, k], :)
end if
P = matmul(P_tmp, P)
! permute column if necessary
Q_tmp = Identity_n(N)
if (cmax /= k) then
Q_tmp(:, [k, cmax]) = Q_tmp(:, [cmax, k])
end if
Q = matmul(Q_tmp, Q)
end do
end subroutine pivot_total