function [b,A,F]=LU_piv_MATLAB(A,b)
% Función que resuelve ns sistemas lineales de ecuaciones compatibles determinados,
% todos con la misma matriz de coeficientes, mediante la factorización LU con pivote parcial:
% FA=LU y utilizando las posibilidades que ofrece el calculo matricial de MATLAB.
% Ultano Kindelán, Álgebra ETSIME, UPM. 09/2022

% Variables de entrada:
    % A = Matriz de coeficientes. 
    % b=  Matriz de dimensión n X ns con los ns vectores de términos independientes.

% Variables de salida
    % b = Matriz de dimensión n X ns con las ns soluciones de los ns sistemas.
    % A = Matriz formada por:
        % parte triangular superior, incluyendo la diagonal: coincide con la parte triangular superior de U. 
        % parte triangular inferior, sin incluir la diagonal: coincide con la parte triangular inferior de L 
        % por debajo de la diagonal.
    % F = Matriz de cambio de filas de la factorizacion FA=LU. 

% Se pone a cero el cronómetro.
tic

% Obtención de la dimensión del sistema.
n=size(A,1);
F=eye(n);
ifil=zeros(n,1);
ns=size(b,2);
%Paso 1 (triangularización)

for k=1:n-1
    %Búsqueda del pivote
    pvt=abs(A(k,k));
    fpvt=k;
    for j=k+1:n
        if (abs(A(j,k))>pvt)
            pvt = A(j,k);
            fpvt = j;
        end
    end
    if (abs(pvt) < eps)
        disp("El sistema no es compatible determinado");
        return;
    elseif (fpvt ~= k)
        A = intercambio_de_filas(A,k,fpvt);
        b = intercambio_de_filas(b,k,fpvt);
        F = intercambio_de_filas(F,k,fpvt);
    end
    % Vector que almacena los cambios de fila
    ifil(k)=fpvt;
    % Pivoteo hacia abajo
    for i=k+1:n
        A(i,k)=A(i,k)/A(k,k);
        % Se suma la fila k multiplicada por A(i,k) a la fila i, sin necesidad 
        % de utilizar un bucle
        A(i,k+1:n)=A(i,k+1:n)-A(i,k)*A(k,k+1:n);
    end
end
if(abs(A(n,n))<eps) 
   disp("El sistema no es compatible determinado");
   return;
end

%Paso 2 (resolución de los ns sistemas)

for j=1:ns
    % Paso 2.1 (para el term. indpte. j, resolución del sistema triangular 
    % inferior por descenso)
    for i=2:n
        % Se realiza la suma sin necesidad de bucle:
        suma=sum(A(i,1:i-1)*b(1:i-1,j));
        b(i,j)=b(i,j)-suma;
    end

    %Paso 2.2 (para el term. indpte. j, resolución del sistema triangular
    %inferior por remonte)

    b(n,j)=b(n,j)/A(n,n);
    for i=(n-1):-1:1
        % Se realiza la suma sin necesidad de bucle:
        suma=sum(A(i,(i+1):n)*b((i+1):n,j));
        b(i,j)=(b(i,j)-suma)/A(i,i);
    end
end

%Se para el cronómetro
toc