function [L,D,P,rho] = Diagpiv_1(A)
%DIAGPIV   Factorización de Pivotación diagonal con pivotación parcial.
%          Matriz de entrada, A (simétrica). Matrices de salida: permutación P,
%          triangular inferior unitaria L, D matriz diagonal en bloques 1x1 o 2x2.
%          P*A*P' = L*D*L'.
%          Ref. J.R. Bunch and L. Kaufman, Some stable methods for calculating
%          inertia and solving symmetric linear systems, Math. Comp. 1977.
%
 if norm(triu(A,1)'-tril(A,-1),1), error('La matriz A debe ser simétrica.'), end
 n = max(size(A)); k = 1; pp = 1:n;
 D = eye(n); L = eye(n); 
 normA = norm(A(:),inf); rho = normA;
 alpha = (1 + sqrt(17))/8;
 while k < n
    [lambda,r] = max( abs(A(k+1:n,k)) );
    r = r(1) + k;
    if lambda > 0
        swap = 0;
        if abs(A(k,k))>=alpha*lambda
           s = 1;
        else
           temp  = A(k:n,r); temp(r-k+1) = 0;
           sigma = norm(temp, inf);
           if alpha*lambda^2 <= abs(A(k,k))*sigma
              s = 1;
           elseif abs(A(r,r)) >= alpha*sigma
              swap = 1;
              m1 = k; m2 = r;
              s = 1;
           else
              swap = 1;
              m1 = k+1; m2 = r;
              s = 2;
           end
        end
        % -----
        if swap
           A( [m1, m2],: ) = A( [m2, m1],: ); A( :,[m1, m2] ) = A( :,[m2, m1] );
           L( [m1, m2],: ) = L( [m2, m1],: ); L( :,[m1, m2] ) = L( :,[m2, m1] );
           pp( [m1, m2] ) = pp( [m2, m1] );
        end
        % -----
        if s == 1                   % Bloque diagonal 1x1
           D(k,k) = A(k,k); i = k+1:n;
           A(i,k) = A(i,k)/A(k,k);
           L(i,k) = A(i,k);
           A(i,i) = A(i,i) - A(i,k) * A(k,i);
        elseif s == 2               % Bloque diagonal 2x2
           E = A(k:k+1,k:k+1);
           D(k:k+1,k:k+1) = E;
           C = A(k+2:n,k:k+1);
           temp = C/E;
           L(k+2:n,k:k+1) = temp;
           A(k+2:n,k+2:n) = A(k+2:n,k+2:n) - temp*C';
       end
       if k+s <= n
          rho = max( rho, max(max(abs(A(k+s:n,k+s:n)))) );
       end
    else                        % Nada, elemento en diagonal 0
       s = 1; D(k,k) = A(k,k);
    end
    k = k + s;
    if k == n, D(n,n) = A(n,n); break, end
 end
 if nargout>=3, P = eye(n); P = P(pp,:); end
 rho = rho/normA; 
end
