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(); }