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.
using System;
using Imsl.Math;
public class GenMinResEx5 : GenMinRes.IFunction
{
//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)
{
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.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 = 92
The final residual norm is 2.52648529541037E-07
Link to C# source.