[Dune] Implementing tri-diagonal matrices in DUNE

[Eike Mueller]

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