dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

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);


>  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.

-- 
Martin