Solver

The conjugate-gradient solver that estimates the unblurred image vector ${\bf {\hat x}}$ implements the Jest solver interface. The solver interface defines a single method,
solve(operator B, vector y, vector x). The arguments stem from the equation ${\bf B x = y}$. The solve() invocation requires a starting solution for ${\bf x}$ which it then will update to its best estimate. The solver's termination criteria are usually set by the solver's constructor.

The solver program  literally translates the conjugate-gradient algorithm into Java code. Gill et al. 1981 list pseudo-code for the conjugate-gradient algorithm as:

$$\begin{eqnarray} p_k &=& - r_k + \beta_{k-1} p_{k-1} \\ \alpha_k &=& \frac{\vert... ... \\ \beta_k &=& \frac{\vert r_{k+1}\vert^2_2}{\vert r_k\vert^2_2}.\end{eqnarray}$$ (1) (2) (3) (4) (5)

Each pseudo-code step corresponds to one or two lines of Jest instructions, as the program's comments indicate. More importantly, the program's instructions are as general as the algorithm, since the program's objects - vectors and operators - are of the same abstraction as the algorithm's mathematical symbols.

  

package juice.solver;
import  jam.solver.*;
import  jam.vector.*;
import  jam.operator.*;
/** Conjugate Gradient Solver. */
public class CGSimple implements LinearSolver {
  private int maxIter; 
  /** @param maxIter maximum number of iterations. */ 
  public CGSimple(int maxIter) { this.maxIter = maxIter; }
  /** 
   * @param A operator to be inverted 
   * @param b known vector 
   * @param x starting estimate of unknown vector
   */
  public void solve(Operator A, Vector b, Vector x ) { 
    if (! (A instanceof LinearOperator)) 
      System.err.println("Error: A is not a LinearOperator");
    Vector r = A.getRange().newMember(); // init r = b - Ax 
    A.residual(x,b,r); 
    r.neg();
    Vector p = (Vector) r.copy();        // init p = r 
    Vector v = A.getRange ().newMember();// init v = 0 
    float rtrNew = r.norm2();

for (int iter=0; iter< maxIter; iter++) {
      A.image(p, v);                     //  v = Ap 
      float ptv   = p.dot( v);           //ptv = p^t Ap
      float alpha = rtrNew/ptv;          //  a = |r|^2 / p^t Ap
      if (ptv == 0f) System.err.println("Error: ptv=0.");
      x.addScale(  alpha, p);            //  x = x + a p 
      r.addScale( -alpha, v);            //  r = r - a Ap 
      float rtrOld = rtrNew;               
      rtrNew       = r.norm2();
      float beta = rtrNew/rtrOld;        //  b = |r|^2 / |r|^2  
      if (rtrOld == 0f) System.err.println("Error: rtrOld=0.");
      p.addScale( beta, r, p);           //  p = -r + b p 
  }}}

CGSimple is a primitive member of Jest's solver family. Its entire interaction with a program is restricted to setting the number of iterations and invoking its solve routine. Standard iterative Jest solvers accept an iterator object that contains application specific initializations, termination criteria, report instructions, and final clean up code.