In [1]:
"""
Return the Givens rotation for a pair of numbers.
"""
function planerot(x)
  a,b = -x[2],x[1]
  d = sqrt(b^2+a^2)
  a,b = a/d,b/d
  G = [b -a;a b]
  y = G*x
  return G,y
end

"""
An "efficient" implementation of GMRES that uses the Givens rotations
to update the solution of the least squares problem.
(Note that there are more efficiencies possible here, but this
  gets at the main ideas)
"""
function gmres_efficient(A,b,tol,maxit)
  n = size(b)
  beta = norm(b)

  Q = zeros(length(b),maxit+1)
  Q[:,1] = b/norm(b)

  g = zeros(maxit+1)
  g[1] = beta;
  Js = zeros(2,maxit) # memory to store the Js
  hist = zeros(maxit)
  H = zeros(maxit+1,maxit)

  i = 1
  for i = 1:maxit
    v = A*Q[:,i]
    H[i+1,i] = 0 # expand k
    for j = 1:i
      H[j,i]= v'*Q[:,j]
      v = v - H[j,i]*Q[:,j]
    end

    # at this point, v = H[i+1,i]*Q[:,i+1]
    H[i+1,i] = norm(v)
    #@show v
    #@show H[i+1,i]
    Q[:,i+1] = v/H[i+1,i]

    # apply J1, ... Ji-1 to z
    for j = 1:i-1
      Jj = [Js[1,j] Js[2,j]; -Js[2,j] Js[1,j]];
      H[j:j+1,i] = Jj*H[j:j+1,i];
    end
    # create Ji to zero out the final row
    Ji,d = planerot(H[i:i+1,i])
    H[i:i+1,i] = d
    Js[1,i] = Ji[1,1]
    Js[2,i] = Ji[1,2]

    # apply Ji to g
    g[i:i+1] = Ji*g[i:i+1]
    hist[i] = abs(g[i+1]);

    if (abs(g[i+1])/beta) < tol
      break
    end
  end

  # produce the solution
  y = H[1:i,1:i]\g[1:i]
  x = Q[:,1:i]*y;

  if (norm(b-A*x)/beta) < tol
    flag = 0
  elseif (hist[i]/beta) < tol
    flag = -1
  else
    flag = 1
  end

  hist = hist[1:i]
  return x,hist,flag
end
Out[1]:
gmres_efficient
In [2]:
##
include("lecture-22-gmres.jl")
Out[2]:
gmres
In [3]:
##
n = 10
A = randn(n,n)
b = randn(n)

x1,hist1 = gmres(A,b,5)

x2,hist2,flag = gmres_efficient(A,b,1.0e-8,5)

@show hist1 - hist2
hist1 - hist2 = [0.0, 0.0, -4.44089e-16, 0.0, 0.0]