function [x, r_norms] = sgcr(M,b,tol,maxiters)
% Generalized conjugate residual method for solving Mx = b
% INPUTS
% M - matrix
% b - right hand side
% tol - convergence tolerance, terminate on norm(b - Mx) < tol * norm(b)
% maxiters - maximum number of iterations before giving up
% OUTPUTS
% x - computed solution, returns null if no convergence
% r_norms - the scaled norm of the residual at each iteration (r_norms(1) = 1)
% Generate the initial guess for x (zero)
x = zeros(size(b));
% Set the initial residual to b - Mx^0 = b
r = b;
% Determine the norm of the initial residual
r_norms(1) = norm(r,2);
for iter = 1:maxiters
% Use the residual as the first guess for the new
% search direction and multiply by M
p(:,iter) = r;
Mp(:, iter) = M * p(:,iter);
% Make the new Mp vector orthogonal to the previous Mp vectors,
% and the p vectors M^TM orthogonal to the previous p vectors
Mp_initial = Mp(:,iter);
for j=1:iter-1,
beta = Mp_initial' * Mp(:,j);
p(:,iter) = p(:,iter) - beta * p(:,j);
Mp(:,iter) = Mp(:,iter) - beta * Mp(:,j);
end;
% Make the orthogonal Mp vector of unit length, and scale the
% p vector so that M * p is of unit length
norm_Mp = norm(Mp(:,iter),2);
Mp(:,iter) = Mp(:,iter)/norm_Mp;
p(:,iter) = p(:,iter)/norm_Mp;
% Determine the optimal amount to change x in the p direction
% by projecting r onto Mp
alpha = r' * Mp(:,iter);
% Update x and r
x = x + alpha * p(:,iter);
r = r - alpha * Mp(:,iter);
% Save the norm of r
r_norms(iter+1) = norm(r,2);
% Print the norm during the iteration
% fprintf('||r||=%g i=%d\n', norms(iter+1), iter+1);
% Check convergence.
if r_norms(iter+1) < (tol * r_norms(1))
break;
end;
end;
% Notify user of convergence
if r_norms(iter+1) > (tol * r_norms(1))
fprintf(1, 'GCR NONCONVERGENCE!!!\n');
x = [];
else
fprintf(1, 'GCR converged in %d iterations\n', iter);
end;
% Scale the r_norms with respect to the initial residual norm
r_norms = r_norms / r_norms(1);