gcr.m

n = 400;
A = diag( log( [2:n+1] ) ) + randn(n) / n;
b = randn(n,1);
y = A \ b;
EPS = 1e-8;


% GCR
disp(' '); disp('GRR');
tic;
R = zeros(n);
c = zeros(n,1);
v = zeros(n); % krylov basis
u = zeros(n); % residuals

x = zeros(n,1);
r = b - A * x;
for k = 1 : n
  % Extend Krylov subspace
  u(:,k) = r;
  v(:,k) = A * u(:,k);

  % Modified Gram Schmidt for the basis
  for i = 1 : k-1 
    R(i,k) = v(:,i)' * v(:,k);
    v(:,k) = v(:,k) - R(i,k) * v(:,i);
  end
  R(k,k) = norm( v(:,k) );
  v(:,k) = v(:,k) / R(k,k);

  % Modified Gram Schmidt for the residual
  c(k) = r' * v(:,k);
  r = r - c(k) * v(:,k);
  if norm(r) < EPS
    break;
  end
end
disp(k);
y = R(1:k,1:k) \ c(1:k);
x = u(:,1:k) * y;
disp(norm(A * x - b));
toc;