From f32a5dc61fa23d30f4536e40e817e547bf5eef89 Mon Sep 17 00:00:00 2001 From: turcksin Date: Wed, 12 Jun 2013 19:34:39 +0000 Subject: [PATCH] TrilinosWrappers::SolverBase::Solve support Epetra_Operator for matrix free calculation. git-svn-id: https://svn.dealii.org/trunk@29811 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/lac/petsc_matrix_free.h | 2 +- deal.II/include/deal.II/lac/trilinos_solver.h | 52 ++++++ deal.II/source/lac/trilinos_solver.cc | 157 +++++++++--------- 3 files changed, 130 insertions(+), 81 deletions(-) diff --git a/deal.II/include/deal.II/lac/petsc_matrix_free.h b/deal.II/include/deal.II/lac/petsc_matrix_free.h index c7217178d3..f5580a0cce 100644 --- a/deal.II/include/deal.II/lac/petsc_matrix_free.h +++ b/deal.II/include/deal.II/lac/petsc_matrix_free.h @@ -34,7 +34,7 @@ namespace PETScWrappers * Implementation of a parallel matrix class based on PETSc MatShell matrix-type. * This base class implements only the interface to the PETSc matrix object, * while all the functionality is contained in the matrix-vector - * multiplication which must be reimplmented in derived classes. + * multiplication which must be reimplemented in derived classes. * * This interface is an addition to the dealii::MatrixFree class to realize * user-defined matrix-classes together with PETSc solvers and functionalities. diff --git a/deal.II/include/deal.II/lac/trilinos_solver.h b/deal.II/include/deal.II/lac/trilinos_solver.h index a0ed043623..54a72beaf1 100644 --- a/deal.II/include/deal.II/lac/trilinos_solver.h +++ b/deal.II/include/deal.II/lac/trilinos_solver.h @@ -155,6 +155,26 @@ namespace TrilinosWrappers const VectorBase &b, const PreconditionBase &preconditioner); + /** + * Solve the linear system + * Ax=b where A + * is an operator. This function + * can be used for matrix free + * computation. 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 (Epetra_Operator &A, + VectorBase &x, + const VectorBase &b, + const PreconditionBase &preconditioner); + /** * Solve the linear system * Ax=b. Depending on the @@ -179,6 +199,32 @@ namespace TrilinosWrappers const dealii::Vector &b, const PreconditionBase &preconditioner); + /** + * Solve the linear system + * Ax=b where A + * is an operator. This function can + * be used for matric free. 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. This class + * works with matrices according to + * the TrilinosWrappers format, but + * can take deal.II vectors as + * argument. Since deal.II are serial + * vectors (not distributed), this + * function does only what you expect + * in case the matrix is locally + * owned. Otherwise, an exception + * will be thrown. + */ + void + solve (Epetra_Operator &A, + dealii::Vector &x, + const dealii::Vector &b, + const PreconditionBase &preconditioner); + /** * Access to object that controls * convergence. @@ -211,6 +257,12 @@ namespace TrilinosWrappers private: + /** + * The solve function is used to set properly the Epetra_LinearProblem, + * once it is done this function solves the linear problem. + */ + void execute_solve(const PreconditionBase &preconditioner); + /** * A structure that collects * the Trilinos sparse matrix, diff --git a/deal.II/source/lac/trilinos_solver.cc b/deal.II/source/lac/trilinos_solver.cc index 60ae66c5de..c315f19139 100644 --- a/deal.II/source/lac/trilinos_solver.cc +++ b/deal.II/source/lac/trilinos_solver.cc @@ -73,8 +73,6 @@ namespace TrilinosWrappers const VectorBase &b, const PreconditionBase &preconditioner) { - int ierr; - linear_problem.reset(); // We need an @@ -87,86 +85,30 @@ namespace TrilinosWrappers &x.trilinos_vector(), const_cast(&b.trilinos_vector()))); - // Next we can allocate the - // AztecOO solver... - solver.SetProblem(*linear_problem); + execute_solve(preconditioner); + } - // ... and we can specify the - // solver to be used. - switch (solver_name) - { - case cg: - solver.SetAztecOption(AZ_solver, AZ_cg); - break; - case cgs: - solver.SetAztecOption(AZ_solver, AZ_cgs); - break; - case gmres: - solver.SetAztecOption(AZ_solver, AZ_gmres); - solver.SetAztecOption(AZ_kspace, additional_data.gmres_restart_parameter); - break; - case bicgstab: - solver.SetAztecOption(AZ_solver, AZ_bicgstab); - break; - case tfqmr: - solver.SetAztecOption(AZ_solver, AZ_tfqmr); - break; - default: - Assert (false, ExcNotImplemented()); - } - // Introduce the - // preconditioner, ... - ierr = solver.SetPrecOperator (const_cast - (preconditioner.preconditioner.get())); - AssertThrow (ierr == 0, ExcTrilinosError(ierr)); - // ... set some options, ... - solver.SetAztecOption (AZ_output, additional_data.output_solver_details ? - AZ_all : AZ_none); - solver.SetAztecOption (AZ_conv, AZ_noscaled); - - // ... and then solve! - ierr = solver.Iterate (solver_control.max_steps(), - solver_control.tolerance()); - - // report errors in more detail - // than just by checking whether - // the return status is zero or - // greater. the error strings are - // taken from the implementation - // of the AztecOO::Iterate - // function - switch (ierr) - { - case -1: - AssertThrow (false, ExcMessage("AztecOO::Iterate error code -1: " - "option not implemented")); - case -2: - AssertThrow (false, ExcMessage("AztecOO::Iterate error code -2: " - "numerical breakdown")); - case -3: - AssertThrow (false, ExcMessage("AztecOO::Iterate error code -3: " - "loss of precision")); - case -4: - AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: " - "GMRES hessenberg ill-conditioned")); - default: - AssertThrow (ierr >= 0, ExcTrilinosError(ierr)); - } + void + SolverBase::solve (Epetra_Operator &A, + VectorBase &x, + const VectorBase &b, + const PreconditionBase &preconditioner) + { + linear_problem.reset(); - // Finally, let the deal.II - // SolverControl object know - // what has happened. If the - // solve succeeded, the status - // of the solver control will - // turn into - // SolverControl::success. - solver_control.check (solver.NumIters(), solver.TrueResidual()); + // We need an + // Epetra_LinearProblem object + // to let the AztecOO solver + // know about the matrix and + // vectors. + linear_problem.reset + (new Epetra_LinearProblem(&A, + &x.trilinos_vector(), + const_cast(&b.trilinos_vector()))); - if (solver_control.last_check() != SolverControl::success) - throw SolverControl::NoConvergence (solver_control.last_step(), - solver_control.last_value()); + execute_solve(preconditioner); } @@ -177,8 +119,6 @@ namespace TrilinosWrappers const dealii::Vector &b, const PreconditionBase &preconditioner) { - int ierr; - linear_problem.reset(); // In case we call the solver with @@ -205,6 +145,39 @@ namespace TrilinosWrappers (const_cast(&A.trilinos_matrix()), &ep_x, &ep_b)); + execute_solve(preconditioner); + } + + + + void + SolverBase::solve (Epetra_Operator &A, + dealii::Vector &x, + const dealii::Vector &b, + const PreconditionBase &preconditioner) + { + linear_problem.reset(); + + Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.begin()); + Epetra_Vector ep_b (View, A.OperatorRangeMap(), const_cast(b.begin())); + + // We need an + // Epetra_LinearProblem object + // to let the AztecOO solver + // know about the matrix and + // vectors. + linear_problem.reset (new Epetra_LinearProblem(&A,&ep_x, &ep_b)); + + execute_solve(preconditioner); + } + + + + void + SolverBase::execute_solve(const PreconditionBase &preconditioner) + { + int ierr; + // Next we can allocate the // AztecOO solver... solver.SetProblem(*linear_problem); @@ -247,7 +220,31 @@ namespace TrilinosWrappers // ... and then solve! ierr = solver.Iterate (solver_control.max_steps(), solver_control.tolerance()); - AssertThrow (ierr >= 0, ExcTrilinosError(ierr)); + + // report errors in more detail + // than just by checking whether + // the return status is zero or + // greater. the error strings are + // taken from the implementation + // of the AztecOO::Iterate + // function + switch (ierr) + { + case -1: + AssertThrow (false, ExcMessage("AztecOO::Iterate error code -1: " + "option not implemented")); + case -2: + AssertThrow (false, ExcMessage("AztecOO::Iterate error code -2: " + "numerical breakdown")); + case -3: + AssertThrow (false, ExcMessage("AztecOO::Iterate error code -3: " + "loss of precision")); + case -4: + AssertThrow (false, ExcMessage("AztecOO::Iterate error code -4: " + "GMRES hessenberg ill-conditioned")); + default: + AssertThrow (ierr >= 0, ExcTrilinosError(ierr)); + } // Finally, let the deal.II // SolverControl object know -- 2.39.5