lin_sol_gen_min_residual
Solves a linear system Ax = b using the restarted generalized minimum residual (GMRES) method.
Synopsis
#include <imsl.h>
float *imsl_f_lin_sol_gen_min_residual (int n, void amultp (float *p, float *z), float *b, ..., 0)
The type double function is imsl_d_lin_sol_gen_min_residual.
Required Arguments
int n (Input)
Number of rows in the matrix.
void amultp (float *p, float *z) (Input)
User-supplied function which computes z = Ap.
float *b (Input)
Vector of length n containing the right-hand side.
Return Value
A pointer to the solution x of the linear system Ax = b. To release this space, use imsl_free. If no solution was computed, then NULL is returned.
Synopsis with Optional Arguments
#include <imsl.h>
float *imsl_f_lin_sol_gen_min_residual (int n, void amultp(), float *b,
IMSL_RETURN_USER, float x[],
IMSL_MAX_ITER, int *maxit,
IMSL_REL_ERR, float tolerance,
IMSL_PRECOND, void precond(),
IMSL_MAX_KRYLOV_SUBSPACE_DIM, int kdmax,
IMSL_HOUSEHOLDER_REORTHOG,
IMSL_FCN_W_DATA, void amultp(), void *data,
IMSL_PRECOND_W_DATA, void precond(), void *data,
0)
Optional Arguments
IMSL_RETURN_USER, float x[] (Output)
A user-allocated array of length n containing the solution x.
IMSL_MAX_ITER, int *maxit (Input/Output)
A pointer to an integer, initially set to the maximum number of GMRES iterations allowed. On exit, the number of iterations used is returned.
Default: maxit = 1000
IMSL_REL_ERR, float tolerance (Input)
The algorithm attempts to generate x such that ∥b - Ax∥2 ≤ τ∥b∥2, where τ = tolerance.
Default: tolerance = sqrt(imsl_f_machine(4))
IMSL_PRECOND, void precond (float *r, float *z) (Input)
User supplied function which sets z = M -1r, where M is the preconditioning matrix.
IMSL_MAX_KRYLOV_SUBSPACE_DIM, int kdmax, (Input)
The maximum Krylov subspace dimension, i.e., the maximum allowable number of GMRES iterations allowed before restarting.
Default: kdmax = imsl_i_min(n, 20)
IMSL_HOUSEHOLDER_REORTHOG,
Perform orthogonalization by Householder transformations, replacing the Gram-Schmidt process.
IMSL_FCN_W_DATA, void amultp (float *p, float *z, void *data), void *data, (Input)
User supplied function which computes z = Ap, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Passing Data to User-Supplied Functions in the introduction to this manual for more details.
IMSL_PRECOND_W_DATA, void precond (float *r, float *z, void *data), void *data (Input)
User supplied function which sets z = M -1r, where M is the preconditioning matrix, which also accepts a pointer to data that is supplied by the user. data is a pointer to the data to be passed to the user-supplied function. See Passing Data to User-Supplied Functions section in the introduction to this manual for more details.
Description
The function imsl_f_lin_sol_gen_min_residual, based on the FORTRAN subroutine GMRES by H.F. Walker, solves the linear system Ax = b using the GMRES method. This method is described in detail by Saad and Schultz (1986) and Walker (1988).
The GMRES method begins with an approximate solution x0 and an initial residual r0 = b - Ax0. At iteration m, a correction zm is determined in the Krylov subspace
κm (v) = span (v, Av, …, Am-1v)
v = r0 which solves the least-squares problem
Then at iteration m, xm = x0 + zm.
Orthogonalization by Householder transformations requires less storage but more arithmetic than Gram-Schmidt. However, Walker (1988) reports numerical experiments which suggest the Householder approach is more stable, especially as the limits of residual reduction are reached.
Examples
Example 1
As an example, consider the following matrix:
Let xT = (1, 2, 3, 4, 5, 6) so that Ax = (10, 7, 45, 33, -34, 31)T. The function imsl_f_mat_mul_rect_coordinate is used to form the product Ax.
#include <imsl.h>
void amultp (float*, float*);
int main()
{
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
float *x;
x = imsl_f_lin_sol_gen_min_residual (n, amultp, b,
0);
imsl_f_write_matrix ("Solution, x, to Ax = b", 1, n, x, 0);
}
void amultp (float *p, float *z)
{
Imsl_f_sparse_elem a[] = {0, 0, 10.0,
1, 1, 10.0,
1, 2, -3.0,
1, 3, -1.0,
2, 2, 15.0,
3, 0, -2.0,
3, 3, 10.0,
3, 4, -1.0,
4, 0, -1.0,
4, 3, -5.0,
4, 4, 1.0,
4, 5, -3.0,
5, 0, -1.0,
5, 1, -2.0,
5, 5, 6.0};
int n = 6;
int nz = 15;
imsl_f_mat_mul_rect_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, p,
IMSL_RETURN_USER_VECTOR, z,
0);
}
Output
Solution, x, to Ax = b
1 2 3 4 5 6
1 2 3 4 5 6
Example 2
In this example, the same system given in the first example is solved. This time a preconditioner is provided. The preconditioned matrix is chosen as the diagonal of A.
#include <imsl.h>
#include <stdio.h>
void amultp (float*, float*);
void precond (float*, float*);
int main()
{
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
float *x;
int maxit = 1000;
x = imsl_f_lin_sol_gen_min_residual (n, amultp, b,
IMSL_MAX_ITER, &maxit,
IMSL_PRECOND, precond,
0);
imsl_f_write_matrix ("Solution, x, to Ax = b", 1, n, x, 0);
printf ("\nNumber of iterations taken = %d\n", maxit);
}
/* Set z = Ap */
void amultp (float *p, float *z)
{
static Imsl_f_sparse_elem a[] =
{0, 0, 10.0,
1, 1, 10.0,
1, 2, -3.0,
1, 3, -1.0,
2, 2, 15.0,
3, 0, -2.0,
3, 3, 10.0,
3, 4, -1.0,
4, 0, -1.0,
4, 3, -5.0,
4, 4, 1.0,
4, 5, -3.0,
5, 0, -1.0,
5, 1, -2.0,
5, 5, 6.0};
int n = 6;
int nz = 15;
imsl_f_mat_mul_rect_coordinate ("A*x",
IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, p,
IMSL_RETURN_USER_VECTOR, z,
0);
}
/* Solve Mz = r */
void precond (float *r, float *z)
{
static float diagonal_inverse[] =
{0.1, 0.1, 1.0/15.0, 0.1, 1.0, 1.0/6.0};
int n = 6;
int i;
for (i=0; i<n; i++)
z[i] = diagonal_inverse[i]*r[i];
}
Output
Solution, x, to Ax = b
1 2 3 4 5 6
1 2 3 4 5 6
Number of iterations taken = 5
Fatal Errors
|
IMSL_STOP_USER_FCN |
Request from user supplied function to stop algorithm. |
