clear all
clc
n = 50; % change it to 5000
% generate a set of n orthogonal vectors to use as eigenvectors
B = rand(n, n);
[Q, ~] = qr(B); %QR decomposition Q(n,n) R(n,m) A(n,n)=QR the R is niglicted
% generate a set of n eigenvalues between 1 <= lambda <= 1000
lambda = linspace(1, 1000, n);
D = diag(lambda);
% generate the SPD dense matrix explicitly
A = Q * D * Q';
% compute the Cholesky decomposition
R = chol(A); % matlab returns the *upper* triangular factor A=R'*R A suppose to be symmetric matrix mxm
% YOUR CODE GOES HERE
T=R'*R; %verification T=A
b=rand(n,1);
RT=R';
o=b; %for not messing with b vector
%forward substituation
y(1)=o(1)/RT(1,1);
for i=1:n
for j=1:i-1
o(i)=o(i)-RT(i,j)*y(j);
end
y(i,1)=o(i)/RT(i,i);
end
z=RT\b; %verify the forward substituation z=y
w=y; %for not messing with y vector
%back substituation
x(n,1) = w(n) / R(n,n);
for i=n-1:-1:1
for j=i+1:n
w(i)=w(i)-R(i,j)*x(j);
end
x(i) = w(i)/R(i,i);
end
v=R\y; %verify the back substituation v=x
r=norm(A*x-b) %checking the residual for forward and back sub.
m=A\b; %verify using backslash method m=v=x
r_m=norm(A*m-b) %checking the residual for backslash method
