[Dune] Implementing tri-diagonal matrices in DUNE

[Oliver Sander]

Hi Eike,
I think what you are asking is a pdelab question.  Since I am not a 
pdelab user
I suggest you ask on the pdelab mailing list.
cheers,
Oliver


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