%
% Conjugate Gradient Method
%

clear all;
close all;

tol = 1e-4;
itmax = 30;

A = [3 2; 2 6];
b = [ 2; -8];
x = [-2; -2];  % initial guess
X = x;  % 션쩌딱덜법넋

r = b - A*x; 
res0 = norm(r);
p = r;
rho = r'*r;
for k = 1 : itmax
    q = A*p;
    alpha = rho / (p'*q);
    x = x + alpha * p;
    X = [X,x];
    r = r - alpha*q;
    relres = norm(r) / res0;
    fprintf('k=%d, relres=%.2e\n', k, relres);
    if relres < tol
        break
    end
    rho_old = rho;
    rho = r'*r;
    beta = rho / rho_old;
    p = r + beta * p;
end


% contour plot
f = @(x1,x2) [x1,x2]*A*[x1;x2]/2 - [x1,x2]*b;
x1 = linspace(-4,6,200);
x2 = linspace(-6,4,200);
[X1,X2] = meshgrid(x1,x2);
[m,n] = size(X1);
F = zeros(m,n);
for i = 1 : m
    for j = 1 : n
        F(i,j) = f(X1(i,j),X2(i,j));
    end
end
% figure(1); mesh(X1,X2,F)
figure(2); contour(X1,X2,F,30); axis square;
hold on;
plot(X(1,:),X(2,:),'or-','LineWidth',1)
hold off
% axis([-2.2,2.2,-3.2,1.2]); % 렴댕
% axis([1.8,2.2,-2.2,-1.8]); % 렴댕

    
