lin_sol_def_cg
Solves a real symmetric definite linear system using a conjugate gradient method. Using optional arguments, a preconditioner can be supplied.
Synopsis
#include <imsl.h>
float *imsl_f_lin_sol_def_cg (int n, void amultp(), float *b, ..., 0)
The type double function is imsl_d_lin_sol_def_cg.
Required Arguments
int n (Input)
Number of rows in the matrix.
void amultp (float *p, float *z)
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_def_cg (int n, void amultp(), float *b,
IMSL_RETURN_USER, float x[],
IMSL_MAX_ITER, int *maxit,
IMSL_REL_ERR, float relative_error,
IMSL_PRECOND, void precond(),
IMSL_JACOBI, float *diagonal,
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 iterations allowed. On exit, the number of iterations used is returned.
IMSL_REL_ERR, float relative_error (Input)
The relative error desired.
Default: relative_error = 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_JACOBI, float diagonal[] (Input)
Use the Jacobi preconditioner, i.e. M = diag(A). The user-supplied vector diagonal should be set so that diagonal[i] = Aii.
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 in the introduction to this manual for more details.
Description
The function imsl_f_lin_sol_def_cg solves the symmetric definite linear system Ax = b using the conjugate gradient method with optional preconditioning. This method is described in detail by Golub and Van Loan (1983, Chapter 10), and in Hageman and Young (1981, Chapter 7).
The preconditioning matrix M is a matrix that approximates A, and for which the linear system Mz = r is easy to solve. These two properties are in conflict; balancing them is a topic of much current research. In the default use of imsl_f_lin_sol_def_cg, M = I. If the option IMSL_JACOBI is selected, M is set to the diagonal of A.
The number of iterations needed depends on the matrix and the error tolerance. As a rough guide,
See the references mentioned above for details.
Let M be the preconditioning matrix, let b, p, r, x, and z be vectors and let τ be the desired relative error. Then the algorithm used is as follows:
Here λ is an estimate of λmax(G), the largest eigenvalue of the iteration matrix G = I - M -1 A. The stopping criterion is based on the result (Hageman and Young 1981, pp. 148-151)
where
It is also known that
where the Tn are the symmetric, tridiagonal matrices
with μk= 1 - βk/αk-1 - 1/αk, μ1 = 1 - 1/α1 and
Usually the eigenvalue computation is needed for only a few of the iterations.
Examples
Example 1
In this example, the solution to a linear system is found. The coefficient matrix is stored as a full matrix.
#include <imsl.h>
static void amultp (float*, float*);
int main()
{
int n = 3;
float b[] = {27.0, -78.0, 64.0};
float *x;
x = imsl_f_lin_sol_def_cg (n, amultp, b, 0);
imsl_f_write_matrix ("x", 1, n, x, 0);
}
static void amultp (float *p, float *z)
{
static float a[] = {1.0, -3.0, 2.0,
-3.0, 10.0, -5.0,
2.0, -5.0, 6.0};
int n = 3;
imsl_f_mat_mul_rect ("A*x",
IMSL_A_MATRIX, n, n, a,
IMSL_X_VECTOR, n, p,
IMSL_RETURN_USER, z,
0);
}
Output
x
1 2 3
1 -4 7
Example 2
In this example, two different preconditioners are used to find the solution of a linear system which occurs in a finite difference solution of Laplace’s equation on a regular c × c grid, c = 100. The matrix is A = E (c2, c). For the first solution, select Jacobi preconditioning and supply the diagonal, so M = diag (A). The number of iterations performed and the maximum absolute error are printed. Next, use a more complicated preconditioning matrix, M, consisting of the symmetric tridiagonal part of A.
Notice that the symmetric positive definite band solver is used to factor M once, and subsequently just perform forward and back solves. Again, the number of iterations performed and the maximum absolute error are printed. Note the substantial reduction in iterations.
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
static void amultp (float*, float*);
static void precond (float*, float*);
static Imsl_f_sparse_elem *a;
static int n = 2500;
static int c = 50;
static int nz;
int main()
{
int maxit = 1000;
int i;
int index;
float *b;
float *x;
float *mod_five;
float *diagonal;
float norm;
n = c*c;
mod_five = (float*) malloc (n*sizeof(*mod_five));
diagonal = (float*) malloc (n*sizeof(*diagonal));
b = (float*) malloc (n*sizeof(*b));
/* Generate coefficient matrix */
a = imsl_f_generate_test_coordinate (n, c, &nz,
0);
/* Set a predetermined answer and diagonal */
for (i=0; i<n; i++) {
mod_five[i] = (float) (i % 5);
diagonal[i] = 4.0;
}
/* Get right hand side */
amultp (mod_five, b);
/* Solve with jacobi preconditioning */
x = imsl_f_lin_sol_def_cg (n, amultp, b,
IMSL_MAX_ITER, &maxit,
IMSL_JACOBI, diagonal,
0);
/* Find max absolute error, print results */
norm = imsl_f_vector_norm (n, x,
IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index,
0);
printf ("iterations = %d, norm = %e\n", maxit, norm);
imsl_free (x);
/* Solve same system, with different preconditioner */
x = imsl_f_lin_sol_def_cg (n, amultp, b,
IMSL_MAX_ITER, &maxit,
IMSL_PRECOND, precond,
0);
norm = imsl_f_vector_norm (n, x,
IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index,
0);
printf ("iterations = %d, norm = %e\n", maxit, norm);
}
/* Set z = Ap */
static void amultp (float *p, float *z)
{
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 */
static void precond (float *r, float *z)
{
static float *m;
static float *factor;
static int first = 1;
float *null = (float*) 0;
if (first) {
/* Factor the first time through */
m = imsl_f_generate_test_band (n, 1,
IMSL_SYMMETRIC_STORAGE,
0);
imsl_f_lin_sol_posdef_band (n, m, 1, null,
IMSL_FACTOR, &factor,
IMSL_FACTOR_ONLY,
0);
first = 1;
}
/* Perform the forward and back solves */
imsl_f_lin_sol_posdef_band (n, m, 1, r,
IMSL_FACTOR_USER, factor,
IMSL_SOLVE_ONLY,
IMSL_RETURN_USER, z,
0);
}
Output
iterations = 115, norm = 1.382828e-05
iterations = 75, norm = 7.319450e-05
Fatal Errors
|
IMSL_STOP_USER_FCN |
Request from user supplied function to stop algorithm. |
