function Feb13()
n = 5;
A = randn(n,n); A = A+A';
b = randn(n,1);
[L,D,piv] = ldl(A,'vector');
% P'*A*P = L*D*L' where P = I(piv,)
% Ax = b implies (P'AP)P'x = P'b,i.e., L(D(L'(P'x))) =P'b
c = zeros(n,1);
c(piv) = b;    % c = P'b
w = L\c;       % Lw = P'b
z = D\w;       % LDz = P'b
y = L'\z;      % LDL'y = P'b;
x = y(piv)     % LDL'Px = P'b
norm(A*x-b)
