Example 5: The Second Householder Implementation

The coefficient matrix in this example corresponds to the five-point discretization of the 2-d Poisson equation with the Dirichlet boundary condition. Assuming the natural ordering of the unknowns, and moving all boundary terms to the right hand side, we obtain a block tridiagonal matrix. (Consider the tridiagonal matrix T which has the value 4.0 down the main diagonal and -1.0 along the upper and lower co-diagonals. Then the coefficient matrix is the block tridiagonal matrix consisting of T's down the main diagonal and -I along the upper and lower codiagonals where I is the identity matrix.) Discretizing on a 20 x 20 grid implies that the coefficient matrix is 400 x 400. In the solution, the second Householder implementation is selected and we choose to update the residual vector by direct evaluation. 
import com.imsl.math.*;
public class GenMinResEx5 implements GenMinRes.Function {
    //Creates a new instance of GenMinResEx5
    public GenMinResEx5() {
    }
    public void amultp(double p[], double z[]) {
        int n = z.length;
        int k = (int)Math.sqrt(n);
        // Multiply by diagonal blocks
        for (int i = 0;  i < n;  i++) {
            z[i] = 4.0 * p[i];
        }
        for (int i = 0;  i < n-2;  i++) {
            z[i] = -1.0 * p[i+1] + z[i];
        }
        for (int i = 0;  i < n-2;  i++) {
            z[i+1] = -1.0 * p[i] + z[i+1];
        }
        // Correct for terms not properly in block diagonal
        for (int i = k-1;  i < n-k;  i=i+k) {
            z[i] += p[i+1];
            z[i+1] += p[i];
        }
        // Do the super and subdiagonal blocks, the -I's
        for (int i = 0;  i < n-k;  i++) {
            z[i] = -1.0 * p[i+k] + z[i];
        }
        for (int i = 0;  i < n-k;  i++) {
            z[i+k] = -1.0 * p[i] + z[i+k];
        }
    }
    
    
    public static void main(String args[]) throws Exception {
        int       n = 400;
        double  b[] = new double[n];
        double  xguess[] = new double[n];
        
        GenMinResEx5 atp = new GenMinResEx5();
        
        // Construct a GenMinRes object
        GenMinRes gnmnrs = new GenMinRes(n, atp);
        // Set right hand side and initial guess to ones    */
        for(int i=0; i<n; i++){
            b[i] = 1.0;
            xguess[i] = 1.0;
        }
        gnmnrs.setGuess(xguess);
        gnmnrs.setMethod(gnmnrs.SECOND_HOUSEHOLDER);
        gnmnrs.setResidualUpdating(gnmnrs.DIRECT_AT_RESTART_ONLY);
        // Solve Ax = b
        gnmnrs.solve(b);
        int iterations = gnmnrs.getIterations();
        System.out.println("The number of iterations used = "+iterations);
        double resnorm = gnmnrs.getResidualNorm();
        System.out.println("The final residual norm is "+resnorm);
    }
}

Output

The number of iterations used = 92
The final residual norm is 2.5264852954103667E-7
Link to Java source.