[Dune] Implementing tri-diagonal matrices in DUNE

[Oliver Sander]

Hi Eike,
check out BTDMatrix from dune-istl.  It is an implementation of a
dynamic tridiagonal matrix, including the Thomas algorithm for
direct solving in O(n).  That's what I use for line-smoothing,
maybe you can use it, too.
cheers,
Oliver

PS: The implementation of BTDMatrix can still be improved, but
it does already give you optimal space and time complexity.

Am 27.11.2011 19:12, schrieb 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
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune