Jacobi algorithm updates b=b-Ax rather than b=b-(L+U)x
Perhaps I am missing something, but while debugging a unrelated bug, I noticed that the jacobi method is computing { \mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -(L+U+D)\mathbf {x} ^{(k)})=D^{-1}(\mathbf {b}-A\mathbf {x} ^{(k)})}
rather than { \mathbf {x} ^{(k+1)}=D^{-1}(\mathbf {b} -(L+U)\mathbf {x} ^{(k)}).}
Here is the offending line:
for (; j.index()<i.index(); ++j)
rhs -= (*j) * x[j.index()]; // b-=Lx
coliterator diag=j; // this is the diagonal entry!
for (; j!=endj; ++j)
rhs -= (*j) * x[j.index()]; // since entry was not skipped, this is b-=(D+U)x
If I am not making an stupid mistake, it should be something like this:
for (; j.index()<i.index(); ++j)
rhs -= (*j) * x[j.index()]; // b-=Lx
coliterator diag=j++; // store diagonal entry and advance iterator
for (; j!=endj; ++j)
rhs -= (*j) * x[j.index()]; // b-=Ux
Does this make sense?
Edited by Santiago Ospina De Los Ríos