The solve() method of Jag's Steepest Decent Solver is almost literally translated from Shewchuk excellent paper 1994.
public void solve(Operator A, Vector b, Vector x ) {
if (! (A instanceof LinearOperator)) {
System.err.println("Error in SteepDescentSolver.solve(): " +
"operator A must be a LinearOperator");
}
this.A = A;
this.b = b;
this.x = x;
// initialize the solver
Vector v = A.getRange().newMember();
Vector res = A.getRange().newMember();
A.residual(x,b,res);
res.neg();
rtr = res.norm2();
// iterator is supplied to the solver by the programmer
iterator.start();
while(iterator.hasMoreIterations()) { // iterator decides about stopping
doOneStep(v);
iterator.nextIteration(); // iterator reports about iteration
iter++;
}
iterator.finish();
}
protected void doOneStep(Vector v) {
A.image(res, v);
// compute step length alpha
float rtv = res.dot(v);
float alpha = rtr/rtv;
if (rtv == 0f) System.err.println("Error in SteepDescentSlv: " +
"divide by rtv (=0).");
// update x and r
// every reset-th iteration I recalculate the exact residual to
// remove the accumulated floating point error.
x.addScale( alpha, res);
if (((iter+1) % reset) == 0) {
A.residual(x,b,res);
res.neg();
}
else {
res.addScale( -alpha, v);
}
rtr = res.norm2();
}