To begin with the interface to cg has been completely
templated. This allows for great flexibility in terms of the types of
matrices used, as well as the other components such as
Preconditioner and Iterator. Also, the parameters
are in terms of MTL Matrix and Vector objects. This is much more convenient
and flexible than passing around data pointers and dimension numbers
separately (it is also much better software engineering).
The first line in the example pulls the scaled()
function in from the mtl namespace. Since the scaled()
function is used many times in cg this is easier than typing
mtl::scaled() each time.
using mtl::scaled;
The second line declares some scalar variables of the same
type as the elements in vector x. It is important
to do this instead of using double or some other
concrete type in order to keep the algorithm as generic
as possible.
typename VectorX::value_type rho,
rho_1, alpha, beta;
After we get the number of rows in matrix A we
declare several vectors with size N.
VectorX p(N), q(N), r(N), z(N);
The next operation is to multiple matrix A times vector
x, subtract b, and put the result in vector
r. This is carried out all in one line with the following
call to the mtl::mult() function.
mtl::mult(A, scaled(x, -1), b, r);
The scaled()
function is somewhat special. It creates a vector adaptor that wraps
up vector x. When elements inside x are accessed the
adaptor causes the element to be multiplied by -1. The other
adaptors in MTL are strided()
and trans().
The next step in the algorithm is to start the iteration,
implemented here with a while loop. The iteration is controled by the
iter object. See the ITL documentation for a full
description of the Iteration responsibilities.
The preconditioner M is then used on vectors r
and z, followed by a call the the mtl::dot_conj() operation.
If this is the first iteration of the loop vector z is
then copied into vector p. In MTL, if one wishes to perform a
deep copy use the mtl::copy() function.
mtl::copy(z, p);
To do a shallow assignment of the vector handles (all MTL objects are
just handles to the real data) use the assignment operator.
p = z;
If this is not the first iteration a call to the mtl::add()
function is made. The scaled()
function is used again, but this time with a vector instead of a
matrix.
The remainder of the algorithm makes several more calls to MTL
operations that we have already seen.
|
|
Using MTL to Create an Iterative Solver
