Example 6: The Second Householder Implementation With Preconditioning

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. Preconditioning is used with the preconditioning matrix being a diagonal matrix with 4.0 down the main diagonal and -1.0 along the upper and lower co-diagonals. This Preconditioner method solves this tridiagonal matrix. 
using System;
using Imsl.Math;
public class GenMinResEx6 : GenMinRes.IPreconditioner
{
	private double[] precondA, precondB, precondC;
	
	// Creates a new instance of GenMinResEx6 
	public GenMinResEx6(int n)
	{
		precondA = new double[n];
		precondB = new double[n];
		precondC = new double[n];
		// Define the preconditioning matrix                    
		for (int i = 0; i < n; i++)
		{
			precondA[i] = 4.0;
			precondB[i] = - 1.0;
			precondC[i] = - 1.0;
		}
	}
	
	public void  Amultp(double[] p, double[] z)
	{
		int m = z.Length;
		int n = (int) Math.Sqrt(m);
		// Multiply by diagonal blocks           
		for (int i = 0; i < m; i++)
		{
			z[i] = 4.0 * p[i];
		}
		for (int i = 0; i < m - 2; i++)
		{
			z[i] -= p[i + 1];
		}
		for (int i = 0; i < m - 2; i++)
		{
			z[i + 1] -= p[i];
		}
		// Correct for terms not properly in block diagonal
		for (int i = n - 1; i < m - n; i = i + n)
		{
			z[i] += p[i + 1];
			z[i + 1] += p[i];
		}
		// Do the super and subdiagonal blocks, the -I's     
		for (int i = 0; i < m - n; i++)
		{
			z[i] -= p[i + n];
		}
		for (int i = 0; i < m - n; i++)
		{
			z[i + n] -= p[i];
		}
	}
	
	/// Solve the tridiagonal preconditioning matrix problem for z.
	public void  Preconditioner(double[] r, double[] z)
	{
		int n = z.Length;
		double[] w = new double[n];
		double[] v = new double[n];
		double[] u = new double[n];
		w[0] = precondA[0];
		v[0] = precondC[0] / w[0];
		u[0] = r[0] / w[0];
		for (int i = 1; i < n; i++)
		{
			w[i] = precondA[i] - precondB[i] * v[i - 1];
			v[i] = precondC[i] / w[i];
			u[i] = (r[i] - precondB[i] * u[i - 1]) / w[i];
		}
		z[n - 1] = u[n - 1];
		for (int j = n - 2; j >= 0; j--)
		{
			z[j] = u[j] - v[j] * z[j + 1];
		}
	}
	
	public static void  Main(String[] args)
	{
		int n = 400;
		double[] b = new double[n];
		double[] xguess = new double[n];
		
		GenMinResEx6 atp = new GenMinResEx6(n);
		
		// 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.Method = Imsl.Math.GenMinRes.ImplementationMethod.SecondHouseholder;
		gnmnrs.ResidualUpdating = Imsl.Math.GenMinRes.ResidualMethod.
			DirectAtRestartOnly;
		// Solve Ax = b
		gnmnrs.Solve(b);
		int iterations = gnmnrs.Iterations;
		Console.Out.WriteLine("The number of iterations used = " + iterations);
		double resnorm = gnmnrs.ResidualNorm;
		Console.Out.WriteLine("The final residual norm is " + resnorm);
	}
}

Output

The number of iterations used = 60
The final residual norm is 2.84141491724154E-07
Link to C# source.