dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #07078
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