function [g, error, niter, flag] = CGmethod(sigma, g, f, maxit, tol)
% Conjugate Gradient method for solving linear systems.
%
% sigma is a parameter for the Gaussian psf
%
% We are solving "Ag = f":
% g as input is the initial guess at a solution and then
% g as output is the solution to the system.
% f is the blurred image.
%
% maxit is the maximum number of iterations to attempt
% niter is the number of iterations actually preformed
% error is the difference between the exact answer and our answer
% tol is the allowable error
% flag is set to 1 if the function returns before error <= tol

flag = 0;
niter = 0;

% r is the residual (error).
r = f - applyGaussianPSF(sigma,g,1);

% Initial search direction.
d = r;

% norm of f (as a vector: turn rows into cols and concatenate)
norm_f = sqrt(inner_product(f,f));

if(norm_f == 0.0), norm_f = 1.0; end;

error = sqrt(inner_product(r,r)) / norm_f;

for niter = 1:maxit
   if(error < tol), return, end; % bail out if the error is small enough.
   
   % multiply our search direction by the blur matrix.
   Ad = applyGaussianPSF(sigma,d,1);
   
   % the square of the blur matrix norm of d. 
   inp_Ad_d = inner_product(Ad,d);
   
   % Step size (how far to move along the gradient direction).
   alpha = inner_product(r,d) / inp_Ad_d;
   
   % Improve the 
   g = g + alpha * d;
   
   % update the residual (error).
   r = f - applyGaussianPSF(sigma,g,1);
   
   % the conjugate factor to force our search directions
   % to be orthogonal (in some sense).
   beta = -inner_product(r,Ad) / inp_Ad_d;
   
   % New search direction.
   d = r + beta * d;
   
   % measure the relative error.
   error = sqrt(inner_product(r,r)) / norm_f;
   
end;

if(error > tol) flag = 1; end;


