function [V, H] = arnoldi_gs(A, v, m)
%
% function [V, H] = arnoldi_gs(A, v, m)
% Performs m iterations of Arnoldi's procedure (Alg. 6.1 in Saad)
% Returns: V: An n by (m + 1) matrx, orthonormal basis for K_{m + 1}(A, v)
%          H: An upper Hessenberg (m + 1) by m matrix.

[nr, nc] = size(A);
H = zeros(m + 1, m);
V = zeros(nr, m + 1);

V(:, 1) = v / norm(v, 2);
for j = 1 : m
    w = A * V(:, j);
    for i = 1 : j
        H(i, j) = V(:, i)' * w;
    end
    for i = 1 : j;
        w = w - H(i, j) * V(:, i);
    end
    Ht = norm(w, 2);
    if abs(Ht) < 1.e-12;
        warning('Returns because h(j + 1, j) = 0');
        return
    end
    V(:, j + 1) = w / Ht;
    H(j + 1, j) = Ht;
end
