A = [1 -3 -1 1; -3 13 7 5; -1 7 14 -2; 1 5 -2 30]
c = [1;-3;2;1]
L = choleskydecomp(A)
norm(A-L*L’,2)
R = LowerTriangularSolver(L,c)
X = UpperTriangularSolver(L’,R)
A\c
norm(A*X-c,2)
function L=choleskydecomp(A)
[n,n]=size(A);
L = zeros(n,n);
for i=1:n
for j=1:n
if i==j
L(i,i) = sqrt(A(i,i)-sum(L(j,1:j-1).^2));
elseif i>j
L(i,j) = (1/L(j,j))*(A(i,j)-sum(L(i,1:j-1).*L(j,1:j-1)));
end
end
end
end
function x = UpperTriangularSolver(A,c)
[n,n] = size(A);
x = zeros(n,1);
x(n) = c(n)/A(n,n);
for i=n-1:-1:1
sum=0;
for j=i+1:n
sum=sum+A(i,j)*x(j);
end
x(i)=(c(i)-sum)/A(i,i);
end
end
function x = LowerTriangularSolver(A,c)
[n,n]=size(A);
x = zeros(n,1);
x(1)=c(1)/A(1,1);
for i =2:n
sum=0;
for j=1:i-1
sum=sum+A(i,j)*x(j);
end
x(i)=(c(i)-sum)/A(i,i);
end
end