← Back to team overview

dolfin team mailing list archive

Re: Linear algebra

 

Ola Skavhaug skrev den 01/04-2008 følgende:
> Martin Sandve Alnæs skrev den 01/04-2008 følgende:
> > 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);
> 
> Aha, I see what your're getting at. Moving the if into the factory ;-)

Just kidding. It does make perfect sense to do that.

Ola
 
> >    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.
> 
> Agree.
> 
> Ola
>  
> > -- 
> > Martin


References