Programa 3.3 (PA=LU: factorizacion con pivoteo - pagina 166).

%Programa 3.3 (PA=LU: factorizacion con pivoteo - pagina 166).
% Calculo de la solucion del sistema lineal AX=B cuando la matriz es invertible.
%
%Sintaxis
% X = lufact(A,B)
%
%Entrada
% - A es una matriz de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
%Resultados
% - X es la matriz de orden N x 1 que contiene la solucin de AX=B.

%Inicializamos X, Y, la matriz de almacenamiento temporal C y
%la matriz fila R donde se registran los intercambios de filas
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
C=zeros(1,N);
R=1:N;

for p=1:N-1

%Determinacion de la fila pivote para la columna q-esima
[max1,j]=max(abs(A(p:N,p)));

%Intercambio de las filas q-esima y j-esima
C=A(p,:);
A(p,:)=A(j+p-1,:);
A(j+p-1,:)=C;
d=R(p);
R(p)=R(j+p-1);
R(j+p-1)=d;

if A(p,p)==0
'A es singular, no hay solucion o no es unica'
break
end

%Calculo del multiplicador,
%que se guarda en la parte subdiagonal de A
for k=p+1:N
mult=A(k,p)/A(p,p);
A(k,p) = mult;
A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);
end
end

%Resolucion para hallar Y
Y(1) = B(R(1));
for k=2:N
Y(k)= B(R(k))-A(k,1:k-1)*Y(1:k-1);
end

%Resolucion para hallar X
X(N)=Y(N)/A(N,N);
for k=N-1:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end