function [X,M,AB,L]= T31gaussps(A,b,nfo)
% Sintaxis:
% =========
%       [X,M,AB,L]= slgauss(A,b,opt,nfo)
%       Resolución de un sistema de ecuaciones por el método de Gauss con 
%       pivote doble.
% EJEMPLOS:
%       [a,b,c,d]=slgauss([2,1;1,2],[3;3],1,1)
%       [a,b,c,d]=slgauss([2,1;1,2],[3;3],'ninguno')
%
% ENTRADA
%   A   Matriz del sistema
%   b   Segundo miembro
%   nfo Nivel de información 
%       >0 saca resultados en cada iteración
%
% SALIDA
%   X   Solución del sistema
%   M   Vector de intercambio de filas (:,1)) y columnas (:,2)
%   AB  Matriz triangular ampliada obtenida por el proceso
%   L   Matriz con los coeficientes usados para anular columnas
%
%   Nota: L*AB(:,1:5)=slpermuta(M(:,1))*A*slpermuta(M(:,2))
%

% Control sobre los argumentos de entrada
if nargin<2;error('Numero de parametros insuficientes');end
if nargin<3;nfo=0;end
% Coherencia del sistema
[n,nc]=size(A);if n~=nc;error('la matriz no es cuadrada');end
[p,nc]=size(b);if nc~=1 | p~=n;error('el sistema no es coherente');end
M=[1:n;1:n]';
A=[A,b];      % Entrada
%
%Triangularizacion del sistema
%
for i=1:n
    % Pivoteo
        [m,ifil]=max(abs( A(i:n,i)) );     %Maximo de la columna y posición
        % Control del pivote e intercambio de filas
        if m==0; error('sistema indeterminado o incompatible');end
        if ifil>1;        
            j=ifil+i-1;
            p=A(i,:);A(i,:)=A(j,:);A(j,:)=p;%Intercambio de filas 
            p=M(i,1);M(i,1)=M(j,1);M(j,1)=p; %Historico del intercambio de filas 
        end
   % Triangularización
    for j=i+1:n;    % Anulamos la columna bajo la diagonal
        m=A(j,i)/A(i,i);  % coeficiente de anulación
        A(j,i+1:end)=A(j,i+1:end)-m*A(i,i+1:end); %Combinamos desde la diagonal
        A(j,i)=m;  % Almacenamos el coeficiente en la posición "nula"
    end
    if nfo;disp([A,M]);end 
end
%
% Resolución del sistema triangular
%
X=0*b; X(n)=A(n,n+1)/A(n,n);
for i=n-1:-1:1
   X(i)=(A(i,n+1)-A(i,i+1:n)*X(i+1:n))/A(i,i);
end
X(M(:,2))=X; % Recolocamos las soluciones (solo con pivote doble)
if nargout>3;L=tril(A(1:n,1:n),-1)+eye(n);end
if nargout>2;AB=triu(A);end