From 3b7941814b102bf9a72c05fb2ea3a1bd080684aa Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 9 Jan 2017 11:14:39 +0100 Subject: [PATCH] Make Trilinos Solver classes compatible with LinearOperators. This patch adds some core functionality to the Trilinos Solver that is later necessary for them to be compatible with LinearOperators. Specifically, this is two additional solve() functions that accept Trilinos EpetraOperators as the input Matrix (with one accomodating wrapped Trilinos preconditioners, and the other another EpetraOperator as a preconditioner.). --- .../changes/minor/20170104Jean-PaulPelteret_2 | 5 + include/deal.II/lac/trilinos_solver.h | 63 +++++++- source/lac/trilinos_solver.cc | 134 ++++++++++++++++-- 3 files changed, 185 insertions(+), 17 deletions(-) create mode 100644 doc/news/changes/minor/20170104Jean-PaulPelteret_2 diff --git a/doc/news/changes/minor/20170104Jean-PaulPelteret_2 b/doc/news/changes/minor/20170104Jean-PaulPelteret_2 new file mode 100644 index 0000000000..2b2e6b066e --- /dev/null +++ b/doc/news/changes/minor/20170104Jean-PaulPelteret_2 @@ -0,0 +1,5 @@ +New: Additional solve() functions have been added to the Trilinos solver +classes, which provide the necessary extension for them to be compatible +with LinearOperators. +
+(Jean-Paul Pelteret, 2017/01/04) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index 44798f1a98..2ec95a4a85 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -61,7 +61,8 @@ namespace TrilinosWrappers * ParameterList) and is similar to the deal.II class SolverSelector. * * @ingroup TrilinosWrappers - * @author Martin Kronbichler, 2008, 2009 + * @author Martin Kronbichler, 2008, 2009; extension for full compatibility + * with LinearOperator class: Jean-Paul Pelteret, 2015 */ class SolverBase { @@ -164,11 +165,59 @@ namespace TrilinosWrappers * Trilinos is chosen. */ void - solve (Epetra_Operator &A, + solve (const Epetra_Operator &A, VectorBase &x, const VectorBase &b, const PreconditionBase &preconditioner); + /** + * Solve the linear system Ax=b where both A and its + * @p precondtioner are an operator. + * This function can be used when both A and the @p preconditioner + * are LinearOperators derived from a TrilinosPayload. + * Depending on the information provided by derived classes and the object + * passed as a preconditioner, one of the linear solvers and preconditioners + * of Trilinos is chosen. + */ + void + solve (const Epetra_Operator &A, + VectorBase &x, + const VectorBase &b, + const Epetra_Operator &preconditioner); + + /** + * Solve the linear system Ax=b where A is an operator, + * and the vectors @p x and @p b are native Trilinos vector types. + * This function can be used when A is a LinearOperators derived + * from a TrilinosPayload. + * Depending on the information provided by derived classes and the object + * passed as a preconditioner, one of the linear solvers and preconditioners + * of Trilinos is chosen. + */ + void + solve (const Epetra_Operator &A, + Epetra_MultiVector &x, + const Epetra_MultiVector &b, + const PreconditionBase &preconditioner); + + /** + * Solve the linear system Ax=b where both A and its + * @p precondtioner are an operator, and the vectors @p x and @p b are + * native Trilinos vector types. + * This function can be used when both A and the @p preconditioner + * are LinearOperators derived from a TrilinosPayload. + * Depending on the information provided by derived classes and the object + * passed as a preconditioner, one of the linear solvers and preconditioners + * of Trilinos is chosen. + */ + void + solve (const Epetra_Operator &A, + Epetra_MultiVector &x, + const Epetra_MultiVector &b, + const Epetra_Operator &preconditioner); + + + /** * Solve the linear system Ax=b. Depending on the information * provided by derived classes and the object passed as a preconditioner, @@ -257,7 +306,15 @@ namespace TrilinosWrappers * The solve function is used to set properly the Epetra_LinearProblem, * once it is done this function solves the linear problem. */ - void do_solve(const PreconditionBase &preconditioner); + template + void do_solve(const Preconditioner &preconditioner); + + /** + * A function that sets the preconditioner that the solver will apply + */ + template + void set_preconditioner (AztecOO &solver, + const Preconditioner &preconditioner); /** * A structure that collects the Trilinos sparse matrix, the right hand diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc index 3c343844f7..400b7df1ee 100644 --- a/source/lac/trilinos_solver.cc +++ b/source/lac/trilinos_solver.cc @@ -88,8 +88,10 @@ namespace TrilinosWrappers + // Note: "A" is set as a constant reference so that all patterns for ::solve + // can be used by the inverse_operator of LinearOperator void - SolverBase::solve (Epetra_Operator &A, + SolverBase::solve (const Epetra_Operator &A, VectorBase &x, const VectorBase &b, const PreconditionBase &preconditioner) @@ -99,7 +101,7 @@ namespace TrilinosWrappers // We need an Epetra_LinearProblem object to let the AztecOO solver know // about the matrix and vectors. linear_problem.reset - (new Epetra_LinearProblem(&A, + (new Epetra_LinearProblem(const_cast(&A), &x.trilinos_vector(), const_cast(&b.trilinos_vector()))); @@ -108,6 +110,72 @@ namespace TrilinosWrappers + // Note: "A" is set as a constant reference so that all patterns for ::solve + // can be used by the inverse_operator of LinearOperator + void + SolverBase::solve (const Epetra_Operator &A, + VectorBase &x, + const VectorBase &b, + const Epetra_Operator &preconditioner) + { + linear_problem.reset(); + + // We need an Epetra_LinearProblem object to let the AztecOO solver know + // about the matrix and vectors. + linear_problem.reset + (new Epetra_LinearProblem(const_cast(&A), + &x.trilinos_vector(), + const_cast(&b.trilinos_vector()))); + + do_solve(preconditioner); + } + + + + // Note: "A" is set as a constant reference so that all patterns for ::solve + // can be used by the inverse_operator of LinearOperator + void + SolverBase::solve (const Epetra_Operator &A, + Epetra_MultiVector &x, + const Epetra_MultiVector &b, + const PreconditionBase &preconditioner) + { + linear_problem.reset(); + + // We need an Epetra_LinearProblem object to let the AztecOO solver know + // about the matrix and vectors. + linear_problem.reset + (new Epetra_LinearProblem(const_cast(&A), + &x, + const_cast(&b))); + + do_solve(preconditioner); + } + + + + // Note: "A" is set as a constant reference so that all patterns for ::solve + // can be used by the inverse_operator of LinearOperator + void + SolverBase::solve (const Epetra_Operator &A, + Epetra_MultiVector &x, + const Epetra_MultiVector &b, + const Epetra_Operator &preconditioner) + { + linear_problem.reset(); + + // We need an Epetra_LinearProblem object to let the AztecOO solver know + // about the matrix and vectors. + linear_problem.reset + (new Epetra_LinearProblem(const_cast(&A), + &x, + const_cast(&b))); + + do_solve(preconditioner); + } + + + void SolverBase::solve (const SparseMatrix &A, dealii::Vector &x, @@ -214,9 +282,9 @@ namespace TrilinosWrappers } - + template void - SolverBase::do_solve(const PreconditionBase &preconditioner) + SolverBase::do_solve(const Preconditioner &preconditioner) { int ierr; @@ -246,16 +314,8 @@ namespace TrilinosWrappers Assert (false, ExcNotImplemented()); } - // Introduce the preconditioner, if the identity preconditioner is used, - // the precondioner is set to none, ... - if (preconditioner.preconditioner.use_count()!=0) - { - ierr = solver.SetPrecOperator (const_cast - (preconditioner.preconditioner.get())); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - } - else - solver.SetAztecOption(AZ_precond,AZ_none); + // Set the preconditioner + set_preconditioner(solver, preconditioner); // ... set some options, ... solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? @@ -299,6 +359,32 @@ namespace TrilinosWrappers + template<> + void + SolverBase::set_preconditioner(AztecOO &solver, + const PreconditionBase &preconditioner) + { + // Introduce the preconditioner, if the identity preconditioner is used, + // the precondioner is set to none, ... + if (preconditioner.preconditioner.use_count()!=0) + { + const int ierr = solver.SetPrecOperator (const_cast + (preconditioner.preconditioner.get())); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + else + solver.SetAztecOption(AZ_precond,AZ_none); + } + + + template<> + void + SolverBase::set_preconditioner(AztecOO &solver, + const Epetra_Operator &preconditioner) + { + const int ierr = solver.SetPrecOperator (const_cast(&preconditioner)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } /* ---------------------- SolverCG ------------------------ */ @@ -637,6 +723,26 @@ namespace TrilinosWrappers } + +// explicit instantiations +// TODO: put these instantiations into generic file +namespace TrilinosWrappers +{ + template void + SolverBase::do_solve(const PreconditionBase &preconditioner); + + template void + SolverBase::do_solve(const Epetra_Operator &preconditioner); + + template void + SolverBase::set_preconditioner(AztecOO &solver, + const PreconditionBase &preconditioner); + + template void + SolverBase::set_preconditioner(AztecOO &solver, + const Epetra_Operator &preconditioner); +} + DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_WITH_PETSC -- 2.39.5