Content |
User Defined Matrix Types
We define our own symmetric matrix type and apply the conjugated gradient method on it. We will only define two things
-
The matrix class MySymatrix and
-
a function blas::mv taht computes the matrix vector product.
That's it. After that we instantly can use our matrix type with the conjugated gradient method.
Our Wired Matrix Type
This is a minimalistic user-defined matrix type:
-
It's a symmetric type as we derive it from SymmetricMatrix
-
A FLENS matrix type needs to define at least the two typedefs IndexType and ElementType
-
Otherwise this matrix type is pretty useless so far.
: public SymmetricMatrix<MySyMatrix>
{
typedef double ElementType;
typedef int IndexType;
};
Now we give it some meaning. We define a matrix-vector product for it where the vector is of type DenseVector. The Operation we define has the form \(y \leftarrow \beta y + \alpha A x\).
We define this product as if the matrix \(A\) has the special tridiagonal form:
\[ \begin{pmatrix} 2 & -1 & 0 & & 0 & 0 \\ -1 & \ddots & \ddots & \ddots & & 0 \\ 0 & \ddots & \ddots & \ddots & \ddots & \\ & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & \ddots & \ddots & \ddots & -1 \\ 0 & 0 & & 0 & -1 & 2 \end{pmatrix} \]The matrix-vector product gets computed by this function.
void
mv(const ALPHA &alpha, const MySyMatrix &A, const DenseVector<VX> &x,
const BETA &beta, DenseVector<VY> &y)
Complete Code
Here the complete code for defining a new matrix type and a corresponding matrix-vector product:
namespace flens {
struct MySyMatrix
: public SymmetricMatrix<MySyMatrix>
{
typedef double ElementType;
typedef int IndexType;
};
} // namespace flens
namespace flens { namespace blas {
template <typename ALPHA, typename VX, typename BETA, typename VY>
void
mv(const ALPHA &alpha, const MySyMatrix &A, const DenseVector<VX> &x,
const BETA &beta, DenseVector<VY> &y)
{
typedef typename DenseVector<VX>::IndexType IndexType;
const IndexType n = x.length();
// we only resize empty left hand sides
if (y.length()==0) {
y.resize(n);
}
ASSERT(y.length()==n);
for (IndexType i=1; i<=n; ++i) {
y(i) = beta*y(i) + alpha*2*x(i);
if (i>1) {
y(i) -= alpha*x(i-1);
}
if (i<n) {
y(i) -= alpha*x(i+1);
}
}
}
} } // namespace blas, flens
Conjugated Gradient (cg) Method
The implementation of the conjugated gradient method was already used in a previous example for sparse matrices.
#define FLENS_EXAMPLES_CG_H 1
#include <flens/flens.h>
#include <limits>
template <typename MA, typename VX, typename VB>
int
cg(const MA &A, const VB &b, VX &&x,
double tol = std::numeric_limits<double>::epsilon(),
int maxIterations = std::numeric_limits<int>::max())
{
typedef typename VB::ElementType ElementType;
typedef typename VB::IndexType IndexType;
typedef typename VB::NoView VectorType;
ElementType alpha, beta, rNormSquare, rNormSquarePrev;
VectorType Ap, r, p;
r = b - A*x;
p = r;
rNormSquare = r*r;
for (int k=1; k<=maxIterations; ++k) {
if (sqrt(rNormSquare)<=tol) {
return k-1;
}
Ap = A*p;
alpha = rNormSquare/(p * Ap);
x += alpha*p;
r -= alpha*Ap;
rNormSquarePrev = rNormSquare;
rNormSquare = r*r;
beta = rNormSquare/rNormSquarePrev;
p = beta*p + r;
}
return maxIterations;
}
#endif // FLENS_EXAMPLES_CG_H
Example Code
#include <flens/examples/my_symatrix.h>
#include <flens/examples/cg.h>
#include <flens/flens.cxx>
using namespace flens;
using namespace std;
int
main()
{
const int n = 5;
MySyMatrix A;
DenseVector<Array<double> > x(n), z(n), b(n);
z = 1;
b = A*z;
cg(A, b, x);
cout << "x = " << x << endl;
x -= z;
cout << "max abs error = " << abs(x(blas::iamax(x))) << std::endl;
}
Some Notes
Solve \(Ax=b\)
cout << "x = " << x << endl;
Get the infinity norm of the error
cout << "max abs error = " << abs(x(blas::iamax(x))) << std::endl;
}
Compile
$shell> cd flens/examples $shell> g++ -Wall -std=c++11 -I../.. -o my_symatrix my_symatrix.cc
Run
$shell> cd flens/examples $shell> ./my_symatrix x = 1 1 1 1 1 max abs error = 0