%% Project 2 part 1
clear
M = rand(5);

A = transpose(M)*M;

b = [3; 3; 3; 3; 3];
x_0 = [0; 0; 0; 0; 0];

eps = 10^-8;

x = x_0;

r_old = b - A*x;

p = r_old;

eps = eps*norm(r_old);

while norm(r_old) > eps
    alph = (transpose(r_old)*r_old)/(transpose(p)*A*p);

    x = x + alph*p;
    r_new = r_old - alph*A*p;
    Beta = (transpose(r_new)*(r_new-r_old))/(transpose(r_old)*r_old);
    p = r_new + Beta*p;
    r_old = r_new;
end



ConjugateGradient(A, x_0, b)


M25 = rand(2^5);
A25 = transpose(M25)*M25;
x_025 = []
for i = 1:2^5
   x_025 = [ x_025; 0]
end
b_25 = []
for i = 1:2^5
   b_25 = [ b_25; 3]
end


[x1, iter1] = ConjugateGradient(A25, x_025, b_25)

M26 = rand(2^6);
A26 = transpose(M26)*M26;
x_026 = []
for i = 1:2^6
   x_026 = [ x_026; 0]
end
b_26 = []
for i = 1:2^6
   b_26 = [ b_26; 3]
end
[x2, iter2] = ConjugateGradient(A26, x_026, b_26)

M27 = rand(2^7);
A27 = transpose(M27)*M27;
x_027 = []
for i = 1:2^7
   x_027 = [ x_027; 0]
end
b_27 = []
for i = 1:2^7
   b_27 = [ b_27; 3]
end
[x3, iter3] = ConjugateGradient(A27, x_027, b_27)

M28 = rand(2^8);
A28 = transpose(M28)*M28;
x_028 = []
for i = 1:2^8
   x_028 = [ x_028; 0]
end
b_28 = []
for i = 1:2^8
   b_28 = [ b_28; 3]
end
[x4, iter4] = ConjugateGradient(A28, x_028, b_28)

M29 = rand(2^9);
A29 = transpose(M29)*M29;
x_029 = []
for i = 1:2^9
   x_029 = [ x_029; 0]
end
b_29 = []
for i = 1:2^9
   b_29 = [ b_29; 3]
end
[x5, iter5] = ConjugateGradient(A29, x_029, b_29)

M210 = rand(2^10);
A210 = transpose(M210)*M210
x_0210 = []
for i = 1:2^10
   x_0210 = [ x_0210; 0]
end
b_210 = []
for i = 1:2^10
   b_210 = [ b_210; 3]
end
[x6, iter6] = ConjugateGradient(A210, x_0210, b_210)

iters = [iter1, iter2, iter3, iter4, iter5, iter6]

iter_ratio = [iters./[2^5, 2^6, 2^7, 2^8, 2^9, 2^10]]

figure(1)
plot([2^5, 2^6, 2^7, 2^8, 2^9, 2^10],iter_ratio, '-k')
title('Ratio of Iterations to Matrix Size per Matrix Size')
xlabel('Matrix Size')
ylabel('Ratio of Iterations to Matrix Size')