[Dune] Applying an ISTL InverseOperator to multiple right-hand sides

[Bård Skaflestad]

Dear Christian,

Thank you for your comments.  I'm simply trying to solve multiple linear systems with the exact same coefficient matrix.  The right-hand sides may be calculated independently.  On the other hand, calculating the right-hand side values in our case amounts (almost) to a complete matrix assembly (one pass over the grid's codim=0 entities for each right-hand side, and with a fairly expensive calculation for each entity).  Thus, my desire to assemble all components at the same time and have the solver automagically DTRT[tm] when presented with the resulting matrix system.

Still, you are right of course, that operating on multiple vectors within a Krylov method is probably not particularly well defined.  I'll restructure our solver to follow the

    krylov-method
      preconditioner
        foreach rhs

format.  Thank you, again, for your help.


Best regards,
-- 
Bård Skaflestad <Bard.Skaflestad at sintef.no>
SINTEF ICT, Applied Mathematics

________________________________________
From: Christian Engwer [christi at uni-hd.de]
Sent: 02 May 2010 21:27
To: Bård Skaflestad
Cc: dune at dune-project.org
Subject: Re: [Dune] Applying an ISTL InverseOperator to multiple right-hand sides

Dear Bard,

> In one of our projects we need to solve sparse systems of linear equations of the form
>
>     Sx = r  (*)
>
> with 'x' and 'r' being ordinary one-dimensional vectors of scalars.  The operator 'S' is symmetric and positive definite, so we've been using the CGSolver template from ISTL (with the Dune::Amg::AMG template as a preconditioner).  Now a need has arisen where we  would like to extend the system above to the form
>
>     SX = R
>
> where the number of right-hand sides (and solutions) is variable--though fixed (typically at three) during any one run.  So far we've assembled the system matrix once for each right hand side.  I'm probably missing something simple, but is there way in ISTL to have the solvers (e.g., the CGSolver or others) accept multiple right-hand sides (and solutions)?
>
> For reference, these are the types/classes we're currently using to represent the system (*)
>
>     BCRSMatrix< FieldMatrix<double, 1, 1> > S
>     BlockVector< FieldVector<double, 1> > x, r
>
> whence we implement the linear solver as a suitable MatrixAdapter of 'S'.
>
> Any and all information and comments is greatly appreciated.

perhaps I didn't get your point, but what exactly are you trying to
achieve? As you are solving your system using an iterative solver, the
computational costs will be same, no matter if you do
foreach rhs
  krylov-method
    preconditioner
or
krylov-method
  preconditioner
    foreach rhs
further more I'm not quite sure what the Krylov method should do when
you try to operate on multiple vectors, as each rhs will lead to a
different defect and thus to different search direction.

What do want to make sure is that you assemble the system matrix
and construct the coarse space hierarchy in the AMG only once, but
there is no reason why you shouldn't be able to do so.

Christian