[Dune] Implementing tri-diagonal matrices in DUNE

[Eike Mueller]

Hi Oliver,

thank you very much, I now got the solver to work with the BTDMatrix, as 
you suggested.

I also tried to build this into the PDELab framework by writing a new 
version of ISTLVectorBackend, based on

Dune::BlockVector<Dune::BlockVector<Dune::BlockFieldVector<E,1>>>

instead of

Dune::BlockVector<BlockFieldVector<E,BLOCKSIZE>>

and a modified version of ISTLBCRSMatrixBackend, based on

Dune::BCRSMatrix<Dune::BTDMatrix<Dune::FieldMatrix<E,1,1>>>

instead of

Dune::BCRSMatrix<Dune::FieldMatrix<E,BLOCKROWSIZE,BLOCKCOLSIZE>>

and this all works fine, i.e. I really just need to replace the backends 
in my code and everything else stays the same:

[...]
// (2) Construct finite element space
     typedef Dune::PDELab::P0LocalFiniteElementMap<Coord, double, dim> FEM;
     FEM fem(Dune::GeometryType(Dune::GeometryType::cube,dim));
     typedef Dune::PDELab::OverlappingConformingDirichletConstraints CON;

     // (3) Set up grid function space
     typedef ISTLBlockVectorBackend<32> VBE;
//    typedef Dune::PDELab::ISTLVectorBackend<1> VBE;
     typedef Dune::PDELab::GridFunctionSpace<GV, FEM, CON, VBE> GFS;
     GFS gfs(gv, fem);

     // (4) Set up grid operator space
     LocalOperator lop(epsilon2);
     typedef ISTLBlockTridiagMatrixBackend<32> MBE;
//    typedef Dune::PDELab::ISTLBCRSMatrixBackend<1,1> MBE;
     typedef Dune::PDELab::EmptyTransformation CC;
     typedef Dune::PDELab::GridOperatorSpace<GFS, GFS, LocalOperator, 
CC, CC, MBE> GOS;
     GOS gos(gfs,gfs,lop);
[...]

If I solve a Poisson problem on a unit cube with a strong anisotropy in 
the x-direction and use SSOR as a preconditioner (i.e. an x-line 
smoother), the tridiagonal solve is very efficient and the solver 
converges much faster than just by using the default Vector/Matrix 
backends (i.e. a point smoother). The only issues I still have are:

* I have to specify the size of the x-direction, i.e. the block size as 
a template parameter to my backends and
* I don't seem to have any control over how the degrees of freedom are 
mapped from the grid onto the vector container, for example in my case 
all points on an x-line are stored in one block of the block vector, but 
I actually want the vertical degrees at a given (x,y) to be grouped 
together. I guess I could modify the access() methods in the backends, 
but the backends do not know anything about the grid function space.

Maybe rewriting the backends is also not the right way forward and the 
grid function space class needs to be adapted as well.

Thanks,

Eike


Oliver Sander wrote:
> 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
>