← Back to team overview

dolfin team mailing list archive

Re: linear solvers

 

On Tue, Jun 01, 2010 at 02:09:47PM +0100, Garth N. Wells wrote:
>
>
> On 01/06/10 10:29, Anders Logg wrote:
> >On Mon, May 31, 2010 at 01:58:08PM +0100, Garth N. Wells wrote:
> >>
> >>
> >>On 31/05/10 13:50, Anders Logg wrote:
> >>>On Mon, May 31, 2010 at 01:17:15PM +0100, Garth N. Wells wrote:
> >>>>At the moment, a linear solver is not associated with a particular
> >>>>matrix. This makes it awkward and dangerous to re-use
> >>>>preconditioners, factorisations, etc. I suggest that we change the
> >>>>interface of the linear solver classes to accept a matrix at
> >>>>construction, and that the solvers maintain a (smart) pointer to the
> >>>>matrix. Opinions?
> >>>
> >>>It feels awkward to associate a linear solver with a specific system.
> >>>
> >>
> >>Why is that? All our la backends do it. It seems natural to me since
> >>solving a linear system involves a number of data structures which
> >>are specific to the matrix.
> >
> >Would would the interface look like?
> >
> >   LinearSolver solver(A);
>
> This is more or less how how PETSc goes about creating solvers
> (associating the matrix with the solver).
>
> >   solver.parameters["preconditioner"] = "amg";
> >
> >   solver.solve(x, b); ?
> >   x = solver.solve(b); ?
> >
>
> We would need to make sure that we avoid a copy constructor call for
> the latter.
>
> >The other option would be to introduce a class LinearSystem:
> >
> >   LinearSystem system(A, b);
>
> This is more or less how how AztecOO (Trilinos) goes about creating
> solvers, except LinearSystem is just a holder, which used to create
> a solver
>
>     LinearSystem system(A, x, b);
>     KrylovSolver solver(system);
>
> >   x = system.solve();
> >
> >That would be analogous to VariationalProblem:
> >
> >   VariationalProblem problem(a, L);
> >   u = problem.solve()
> >
> >Looks like there are three options:
> >
> >   1. LinearSolver / VariationalSolver
> >   2. LinearSystem / VariationalProblem
> >   3. Both 1 and 2
> >
>
> Perhaps both. LinearSystem/VariationalProblem hold data, and
> LinearSolver/VariationalSolver do something.

That's what I'm thinking too now. We should update the blueprint:

https://blueprints.launchpad.net/dolfin/+spec/solver-interfaces

We also need shortcuts such as

1. Calling LinearSolver::solve from LinearSystem::solve

  LinearSystem::solve()
  {
    LinearSolver solver(*this);
    return solver.solve();
  }

2. Simple solve() function that wraps everything

  x = solve(A, b, parameters)
  u = solve(a, L, parameters)

Concerning the creation of temporaries, didn't you figure out how to
get around that in creating dolfin::refine()?

--
Anders

Attachment: signature.asc
Description: Digital signature


References