function ShowLDLPosDef()
% Using ldl to solve a symmetric positive definite system...
clc
format short
n = 6;
A = pascal(n);
b = randn(n,1);
[L,D,P] = ldl(A)

% P'AP = LDL', i.e., A = PLDL'P'
% To solve Ax = b we want to solve P(L(D(L'(P'x)))) = b

x = P*(L'\(D\(L\(P'*b))));

% Equivalently...
% Solve P*x1 = b for x1...
x1 = P'*b;
% Solve L*x2 = x1 for x2...
x2 = L\x1;
% Solve D*x3 = x2 for x3...
x3 = D\x2;
% Solve L'*x4 = x3 for x4...
x4 = L'\x3;
% Solve P'x = x4 for x...
x = P*x4;