From 784caee7cb8d25ed44967140ff6cafec9bca5199 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Nov 2017 14:50:09 +0100 Subject: [PATCH] Enable Trilinos direct solvers with LA::distributed::Vector. --- include/deal.II/lac/trilinos_solver.h | 9 +++++++ source/lac/trilinos_solver.cc | 35 ++++++++++++++++++++++++--- 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index 17e672d342..38d63ac5de 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -645,6 +645,15 @@ namespace TrilinosWrappers */ void solve (MPI::Vector &x, const MPI::Vector &b); + /** + * Solve the linear system Ax=b based on the package set in + * intialize() for deal.II's own parallel vectors. Note the matrix is not + * refactorized during this call. + */ + void + solve (dealii::LinearAlgebra::distributed::Vector &x, + const dealii::LinearAlgebra::distributed::Vector &b); + /** * Solve the linear system Ax=b. Creates a factorization of the * matrix with the package chosen from the additional data structure and diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc index 2839e10dc2..79c19a658e 100644 --- a/source/lac/trilinos_solver.cc +++ b/source/lac/trilinos_solver.cc @@ -780,17 +780,45 @@ namespace TrilinosWrappers // Assign the RHS vector to the Epetra_LinearProblem object linear_problem->SetRHS(const_cast(&b.trilinos_vector())); + // First set whether we want to print the solver information to screen or + // not. + ConditionalOStream verbose_cout (std::cout, + additional_data.output_solver_details); + + + verbose_cout << "Starting solve" << std::endl; + // Fetch return value of Amesos Solver functions - int ierr; + int ierr = solver->Solve (); + AssertThrow (ierr == 0, ExcTrilinosError (ierr)); + + // Finally, force the SolverControl object to report convergence + solver_control.check (0, 0); + } + + + + void SolverDirect::solve (dealii::LinearAlgebra::distributed::Vector &x, + const dealii::LinearAlgebra::distributed::Vector &b) + { + Epetra_Vector ep_x (View, linear_problem->GetOperator()->OperatorDomainMap(), x.begin()); + Epetra_Vector ep_b (View, linear_problem->GetOperator()->OperatorRangeMap(), const_cast(b.begin())); + + // Assign the empty LHS vector to the Epetra_LinearProblem object + linear_problem->SetLHS(&ep_x); + + // Assign the RHS vector to the Epetra_LinearProblem object + linear_problem->SetRHS(&ep_b); // First set whether we want to print the solver information to screen or // not. ConditionalOStream verbose_cout (std::cout, additional_data.output_solver_details); - verbose_cout << "Starting solve" << std::endl; - ierr = solver->Solve (); + + // Fetch return value of Amesos Solver functions + int ierr = solver->Solve (); AssertThrow (ierr == 0, ExcTrilinosError (ierr)); // Finally, force the SolverControl object to report convergence @@ -917,7 +945,6 @@ namespace TrilinosWrappers do_solve(); } - } -- 2.39.5