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]