## We'll look at a Jacobi type iteration on this system
A = [1 2; -2 1.0]
b = [2.0; 1.0]
A\b
##
M = zeros(Int,2)
M[1] = 2 # use row 2 to update x_1
M[2] = 1 # use row 1 to update x_2
function update(x)
xnew = similar(x)
for i=1:length(x)
row = M[i]
rhs = A[row,:]'*x - A[row,i]*x[i] # there are ways to avoid subtraction :)
rhs = b[row] - rhs
xnew[i] = (rhs/A[row,i])
end
return xnew
end
update([1,0.5])
##
x = [1.0,0.5]
for iter=1:10
global x = update(x)
@show x
end
## The choice of 1->2 and 2->1 was arbitrary, what about 1-1, 2-2?
M[1] = 1
M[2] = 2
x = [1.0,0.5]
for iter=1:20
global x = update(x)
@show x
end
##
# I'm in a hurry... and tried this for my homework.
M[1] = 2
M[2] = 1
function update2(x)
for i=1:length(x)
row = M[i]
rhs = A[row,:]'*x - A[row,i]*x[i]
rhs = b[row] - rhs
x[i] = rhs/A[row,i]
end
return x
end
x = [1.0,0.5]
for i=1:10
global x = update2(x)
@show x
end
## This is the Gauss Seidel iteration Note that
## we didn't allocate an extra vector of memory.