
function  X = gausspiv(A,B)
%-----------------------------------------------------------------------
%             -------------                             -------
%             function file                             gausspiv.m
%             -------------                             -------
%-----------------------------------------------------------------------
%                        ----------------------------
%                        Zeros Reais de Sistemas Lineares
%                        ----------------------------
%                     ---------------------------------------
%                     Metodo de Gauss com Pivotamento Parcial
%                     ---------------------------------------
%
%            Chamada do Procedimento  gausspiv
%------------------------------------------
%X = gausspiv(A,B)
%
% Resolve AX=B por eliminação Gaussiana 
%              com pivotação parcial;
%         
% Entrada: A (N x N) matriz dos coeficientes (não singular);
%          B (N x 1) matriz dos termos independentes;
% Saída: X (N x 1) solução de AX=B.
%Inicializa X e utiliza C como matriz temporária
%
%Veja Tambem: gauss jacobi cramer sor
%--------------------------------------------------------
[N N] = size(A);
X = zeros(N,1);
C = zeros(1,N+1);

%Matriz ampliada Amp=[A|B]
Amp = [A B];

for p = 1:N-1
   
   %Pivotação parcial na coluna p
   [Y,j] = max(abs(Amp(p:N,p))); 
   %Troca linha p com j
   C = Amp(p,:);
   Amp(p,:) = Amp(j+p-1,:);
   Amp(j+p-1,:) = C;
   
   if Amp(p,p) == 0
      disp('A é singular.')
      break
   end
   
   %Eliminação de Gauss na coluna p
   for k = p+1:N
      m = Amp(k,p)/Amp(p,p);
      Amp(k,p:N+1) = Amp(k,p:N+1) - m*Amp(p,p:N+1);
   end
end

%Obtém solução por substituição ascendente 
X=subasc(Amp(1:N,1:N),Amp(1:N,N+1));

function  X = subasc(U,Y)

% Resolve UX=Y por substituição ascendente
%
% Entradas: U (n x n) matriz triangular superior (não singular);
%           Y (n x 1) matriz dos termos independentes;
% Saída: X (n x 1) solução do sistema linear UX=Y.

n = length(Y);
X = zeros(n,1);
X(n) = Y(n)/U(n,n);

for k = n-1:-1:1
   X(k) = (Y(k) - U(k,k+1:n)*X(k+1:n))/U(k,k);
end
