function V = arnoldi ( x0, b, atv, k )

%*****************************************************************************80
%
%% arnoldi() implements the Arnoldi orthogonormalization process.
%
%  Discussion:
%
%    The columns of the output matrix V form an orthonormal basis for Kk,
%    the affine Krylov space based at x0, comprising the span of
%    (r0,Ar0,...,A^(k-1)r0), where r0 = b - A x0.
%
%    The method can break down because of an unlikely division by zero,
%    in which case the exact solution is in the affine subspace, a
%    so-called "happy breakdown".
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 April 2021
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%    
%  Input:
%
%    real x0(n): the initial solution estimate.
%
%    real b(n): the right hand side.
%
%    function ax = atv(x):, a matrix-vector product routine computing A*x.
%
%    real params(2): iteration parameters.
%    params(1): relative residual reduction factor.
%    params(2): kmax = maximum number of iterations.
%
%  Output:
%
%    real x(n): the estimated solution.
%
%    real err(*): vector of iteration residual norms.
%
%    integer total_iters: number of iterations.
%
  b = b(:);
  x0 = x0(:);

  n = size ( b, 1 );
  V = zeros ( n, k );

  r0 = b - feval ( atv, x0 );
  V(1:n,1) = r0 / norm ( r0 );

  for i = 1 : k - 1
    avi = feval ( atv, V(1:n,i) );
    w = avi;
    for j = 1 : i
      w = w - ( avi' * V(1:n,j) ) * V(1:n,j);
    end
    V(1:n,i+1) = w / norm ( w );
  end

  return
end

