# [Dune] Implementing tri-diagonal matrices in DUNE

Eike Mueller em459 at bath.ac.uk
Sun Nov 27 19:12:16 CET 2011

```Dear dune-list,

I have a question about the efficient implementation of tridiagonal
matrices in Dune and using them in block-iterative smoothers/
preconditioners.

The reason I'm interested in this is that I'm solving a PDE on a
regular three dimensional grid (size nx*ny*nz), with a very strong
coupling in the vertical direction. I want to exploit this by using a
vertical line smoother as a preconditioner as follows: arrange the
fields in a BlockVector structure of dimension nx*ny, where each entry
is a FieldVector of length nz=BLOCKSIZE (=70 in my case), i.e. store
my fields in a structure of this type:

typedef Dune::BlockVector<Dune::FieldVector<double,BLOCKSIZE>> Vector;

I then construct a corresponding (block-)matrix like this:

typedef
Dune::BCRSMatrix<Dune::FieldMatrix<double,BLOCKSIZE,BLOCKSIZE>> Matrix;

(I should say that I do not use PDELab, but my matrix is stored
explicitly in a matrix-market file, which I read in.) This all works
fine (including the read from the matrixmarket file), If I use, for
example, BiCGStab with a SOR preconditioner, the method converges in
very few iterations, so is very beneficial for my problem (using a
point-smoother as a preconditioner is absolutely hopeless). Using the
line-smoother with an AMG preconditioner gives even better results
(and AMG with a point-smoother gives very large iteration counts as
well).

Now, I know that my BLOCKSIZE x BLOCKSIZE matrices are tridiagonal,
but as I understand the block-recursive ISTL implementation, at the
lowest level a solve() for the dense FieldMatrix will be called, which
is of course very slow as it inverts the entire BLOCKSIZE x BLOCKSIZE
matrix, so I'm sure the solution time can be reduced significantly.

I partially managed to get around this by defining a derived type

template<class K, int ROWS>
class TriDiagFieldMatrix : public Dune::FieldMatrix<K, ROWS, ROWS>
[...]

and implemented a tridiagonal solve in the solve() methods, which just
ignores all the zero entries. Although this works, I strongly suspect
it is still not the best way of doing it because

(1) I waste a lof memory by storing all the zero elements in the
BLOCKSIZE x BLOCKSIZE matrix
(2) The sparsity of the block matrices is not exploited at all in the
matrix-vector operations like y = A.x, etc.

Could anyone point me to the best way forward in this case? Would I
have to write an entirely new class, which implements the
functionality of densematrix (and is there any description of the
matrix interface that the class has to implement?, to make sure it can
be used withing the ISTL framework). I thought about deriving my class
from DenseMatrix in the same way as FieldMatrix is derived from
DenseMatrix, to avoid implementing the y = A.x etc. operations, but
looking at the DenseMatrix source code this might not work. Or is the
best way forward to derive a new class TriDiagBCRSMatrix from
BCRSMatrix, add a solve() method for this class and use a three level
matrix structure like this:

BCRSMatrix<TriDiagBCRSMatrix<Dune::FieldMatrix<double,1,1>>>

?

Thank you very much for any ideas,

Eike

```