From 0e9866f59fe84127a8ca5216a0e3255d747e41ed Mon Sep 17 00:00:00 2001 From: Mike Harmon Date: Sun, 10 Apr 2016 16:54:19 +0200 Subject: [PATCH] Refactored trilinos direct solver --- include/deal.II/lac/trilinos_solver.h | 15 ++++++ source/lac/trilinos_solver.cc | 77 ++++++++++++++++++++++++++- 2 files changed, 90 insertions(+), 2 deletions(-) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index 9ab13fc2aa..084bf77a41 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -567,6 +567,21 @@ namespace TrilinosWrappers */ virtual ~SolverDirect (); + /** + * Initializes the direct solver for the matrix A and creates a + * factorization for it with the package chosen from the additional + * data structure. Note that there is no need for a preconditioner + * here and solve() is not called. + */ + void initialize (const SparseMatrix &A); + + /* + * Solve the linear system Ax=b based on the + * package set in intialize(). Note the matrix is not refactorized during + * this call. + */ + void solve (VectorBase &x, const VectorBase &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 6b539705ce..b28637bc55 100644 --- a/source/lac/trilinos_solver.cc +++ b/source/lac/trilinos_solver.cc @@ -442,9 +442,15 @@ namespace TrilinosWrappers - void - SolverDirect::do_solve() + void SolverDirect::initialize (const SparseMatrix &A) { + // We need an Epetra_LinearProblem object to let the Amesos solver know + // about the matrix and vectors. + linear_problem.reset (new Epetra_LinearProblem ()); + + // Assign the matrix operator to the Epetra_LinearProblem object + linear_problem->SetOperator(const_cast(&A.trilinos_matrix())); + // Fetch return value of Amesos Solver functions int ierr; @@ -480,6 +486,73 @@ namespace TrilinosWrappers verbose_cout << "Starting numeric factorization" << std::endl; ierr = solver->NumericFactorization(); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + void SolverDirect::solve (VectorBase &x, const VectorBase &b) + { + // Assign the empty LHS vector to the Epetra_LinearProblem object + linear_problem->SetLHS(&x.trilinos_vector()); + + // Assign the RHS vector to the Epetra_LinearProblem object + linear_problem->SetRHS(const_cast(&b.trilinos_vector())); + + // Fetch return value of Amesos Solver functions + int ierr; + + // 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 (); + AssertThrow (ierr == 0, ExcTrilinosError (ierr)); + + // Finally, force the SolverControl object to report convergence + solver_control.check (0, 0); + } + + + + void + SolverDirect::do_solve() + { + // Fetch return value of Amesos Solver functions + int ierr; + + // First set whether we want to print the solver information to screen or + // not. + ConditionalOStream verbose_cout (std::cout, + additional_data.output_solver_details); + + solver.reset (); + + // Next allocate the Amesos solver, this is done in two steps, first we + // create a solver Factory and and generate with that the concrete Amesos + // solver, if possible. + Amesos Factory; + + AssertThrow (Factory.Query (additional_data.solver_type.c_str ()), + ExcMessage (std:: + string ("You tried to select the solver type <") + + additional_data.solver_type + + "> but this solver is not supported by Trilinos either " + "because it does not exist, or because Trilinos was not " + "configured for its use.")); + + solver.reset (Factory. + Create (additional_data.solver_type.c_str (), + *linear_problem)); + + verbose_cout << "Starting symbolic factorization" << std::endl; + ierr = solver->SymbolicFactorization (); + AssertThrow (ierr == 0, ExcTrilinosError (ierr)); + + verbose_cout << "Starting numeric factorization" << std::endl; + ierr = solver->NumericFactorization (); + AssertThrow (ierr == 0, ExcTrilinosError (ierr)); verbose_cout << "Starting solve" << std::endl; ierr = solver->Solve(); -- 2.39.5