conjugate gradient algorithm
for j = 2:ny for i = 2:nx
q[i,j] = (p[i+1,j] - 2.0*p[i,j] + p[i-1,j])/(dx^2) +
(p[i,j+1] - 2.0*p[i,j] + p[i,j-1])/(dy^2)
end~end
aa, bb = 0.0, 0.0
for j = 2:ny for i = 2:nx
aa = aa + r[i,j]*r[i,j] # <r,r>
bb = bb + q[i,j]*p[i,j] # <q,p>
end end
cc = aa/(bb + tiny) #alpha
for j = 2:ny for i = 2:nx
un[i,j] = un[i,j] + cc*p[i,j] # update the numerical solution
end~end
bb, aa = aa, 0.0
for j = 2:ny for i = 2:nx
r[i,j] = r[i,j] - cc*q[i,j] # update the residual
aa = aa + r[i,j]*r[i,j]
end end
cc = aa/(bb+tiny) # beta
for j = 1:ny for i = 1:nx
p[i,j] = r[i,j] + cc*p[i,j] # update the conjugate vector
end~end