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);