dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Apr 01, 2008 at 03:48:38PM +0200, Martin Sandve Alnæs wrote:
> 2008/4/1, Ola Skavhaug <skavhaug@xxxxxxxxx>:
> > Anders Logg skrev den 01/04-2008 følgende:
> >
> > > On Tue, Apr 01, 2008 at 11:23:29AM +0200, Martin Sandve Alnæs wrote:
> >  > > 2008/4/1, Kent-Andre Mardal <kent-and@xxxxxxxxx>:
> >  > > >
> >  > > >  tir, 01.04.2008 kl. 11.08 +0200, skrev Martin Sandve Alnæs:
> >  > > >
> >  > > > > I don't see what's the point of that. The only reasons for
> >  > > >  > Matrix/Vector are related to having a single LA backend in each
> >  > > >  > application run, but a variety of solvers may be used in the same
> >  > > >  > application.
> >  > > >  >
> >  > > >
> >  > > >
> >  > > > I thought the reason for having Matrix/Vector was to ensure that only
> >  > > >  GenericMatrix/GenericVector functionality was used. I guess the PETSc
> >  > > >  and the uBlas families now live seperate lives.
> >  > >
> >  > > Then why don't you come up with an example of where this makes sense
> >  > > for solvers.
> >  >
> >  > Here's an example:
> >  >
> >  > Choosing the type of solver based on some criterion and then
> >  > reusing the solver:
> >  >
> >  >    // At startup
> >  >    LinearSolver* solver = 0;
> >  >    if (solver_type == "direct")
> >  >      solver = new LUSolver();
> >  >    else
> >  >      solver = new KrylovSolver()
> >  >
> >  >    // Then reuse it
> >  >    while (t < T)
> >  >    {
> >  >      ...
> >  >      solver->solve(A, x, b);
> >  >    }
> >  >
> >  > This is used in the ODE solvers, where one may specify the type of
> >  > solver used to solve the linearized equations in each iteration based
> >  > on a parameter.
> >  >
> >  > When calling LUSolver::solve(), it needs to choose a specific backend
> >  > for the solution, depending on the backend used for A, x, b. Same for
> >  > KrylovSolver.
> >  >
> >  > Still not sure if it should be LinearSolver or GenericLinearSolver
> >  > above.
> >
> >
> > What about a solution something like the following. It only uses the top level
> >  generic types during solve:
> >
> >  GenericSolver* solver = 0;
> >
> >  if (solver_type == "direct")
> >     solver = A.factory()->directSolver();
> >  else
> >     solver = A.factory()->iterativeSolver();
> 
> 
> Why all the if's?
> 
>    LAFactory * fac = A.factory();
> 
>    LinearSolver* solver = fac->createLinearSolver(solver_type);
>    solver->setParameters(linsolveParams);
> 
>    Preconditioner* B = fac->preconditioner(preconditioner_type);
>    B->setParameters(precParams);
>    B->update(A);

Looks good to me.

> >  while (t < T)
> >  {
> >     if (B != 0)
> >         solver->solve(B, A, x, b);
> >     else
> >         solver->solve(A, x, b);
> >  }
> 
> solve(B, A, x, b) should be able to deal with no preconditioner.
> Alternatively we collect B, A, x, b in LinearSystem and pass that to solve.

I think we should collect all optional arguments at the end:

   solve(A, x, b)
   solve(A, x, b, options)
   solve(A, x, b, pc)
   solve(A, x, b, pc, options)

where options is a dictionary (in Python) and a ParameterList in C++.

-- 
Anders