function R=Cholesky(A)
%Based on Cholesky factorization algorithm in Textbook page 126
n = size(A);
R=zeros(n);
for k=1:n
if(A(k,k)<0)
error('Matrix A is not SPD!')
end
R(k,k) =sqrt(A(k,k));
u=A(k,k+1:n)/R(k,k);
R(k,k+1:n)=u;
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-u'*u;
end