% Nicolae Tecu, small changes by M.O.
% Numerical Methods I
% Problem Set 5, Problem 5

% this function solves Ax=b using the conjugate gradient method.
% The parameters are A, b, and a starting vector.
% The function returns the iteration history (the 2-norm of the residual
% at each step), the number of iterations needed, and the solution x.

function [norm2res, it, sol] = conjgradmeth (A,b,st)

  error      = 10^(-9);     % Desired tolerance
  [n,n]      = size(A);
  iterations = n;           % Max number or iterations 
  x          = st; 
  res        = zeros(n, iterations);
  p          = zeros(n, iterations);

  res(:,1) = b - A*x;
  norm2res(1,1) = norm(res(:,1));
  p(:,1) = res(:,1);
  alpha  = (p(:,1))' * res(:,1)/((p(:,1))' *A* p(:,1));
  i = 1;
  while (i <= iterations & norm2res(i,1)>error)
      x = x + alpha * p(:,i); 
      i = i + 1;
      res(:,i) = b - A*x;
      norm2res(i,1) = norm(res(:,i));
      p(:,i) = res(:,i) - ...
	       p(:,i-1) * (res(:,i))'*A*p(:,i-1)/((p(:,i-1))'*A*p(:,i-1));
      alpha = (p(:,i))'*res(:,i)/((p(:,i))' * A * p(:,i));
  end

  it = i;
  sol = x;

end