I'm having trouble implementing this algorithm in octave even though the psuedocode for this algorithm looks really simple. The book covering this algorithm is only 1 page long, so there are not a lot of information regarding this algorithm, so I'm just going to post the psedocode:
Compute r = b - Ax and p = A r
Until convergence, Do:
a <- (r,r) / (p,r)
x <- x + a r
r <- r - r p
compute p := A r
End do
Here is my attempt at implementing this in Octave. I use an example in the book to test the program:
A = [5,2,-1;3,7,3;1,-4,6];
b = [2;-1;1];
x0 = [0;0;0];
Tol = 0.00001;
x=x0;
r = b-A*x;
p = A*r;
while true,
a = (r')*(r)/((p)*(r'));
disp(a);
x = x + a * r;
r = r - a * p;
p = A*r;
if norm(r) < Tol,
break
end
end
When I run this, i get an error saying that the first matrix im dividing with is 1x1 and the second matrix is 3v3, so I'm not able to do that and I understand that. I thought about using the ./ operator instead, but to my understanding this does not yield the result that I'm looking for and this example should be eligble to division. Did i screw up my implementation or is my understanding of this algorithm wrong? Not sure whether to post this here or math.stackexchange, but I tried here.
My first thought is the error message: You have (r') * (r) / ((p) * (r')); should the denominator be (p) * (r') or (p') * (r) (notice where the ' are)?