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
signature.asc
Description: Digital signature
_______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
