function [b,A,F]=LU_piv(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.
% 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. 
    

% OBSERVACIÓN: Utilizando las posibilidades de cálculo matricial de MATLAB, esta función se puede programar
% de forma más concisa y simple (a más "alto nivel"). Al ser un curso de Álgebra Lineal y no de MATLAB,
% lo que interesa es que se entienda la factorización LU y, por ejemplo, se sepa contar el número
% de operaciones en coma flotante que consume el método. Considero, por ello que es más ilustrativa la programación
% "elemento a elemento" (además de ser más eficiente computacionalmente). Para los que estén interesados,
% la función "gauss_MATLAB.m" es una versión de "alto nivel" de "gauss.m".


% 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);
        for j=k+1:n
            A(i,j)=A(i,j)-A(i,k)*A(k,j);
        end
    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
        suma=0;
        for k=1:i-1
            suma=suma+A(i,k)*b(k,j);
        end
        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
        suma=0;
        for k=(i+1):n
            suma=suma+A(i,k)*b(k,j);
        end
        b(i,j)=(b(i,j)-suma)/A(i,i);
    end
end

%Se para el cronómetro
toc