From: Toby D. Young Date: Wed, 21 Jul 2010 10:21:35 +0000 (+0000) Subject: Add solver to SLEPcWrappers and trivially update documentation X-Git-Tag: v8.0.0~5783 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8c2b42e16fd9523a0f030de3bcb78791ec177166;p=dealii.git Add solver to SLEPcWrappers and trivially update documentation git-svn-id: https://svn.dealii.org/trunk@21550 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/slepc_solver.h b/deal.II/lac/include/lac/slepc_solver.h index a5ae859aeb..bdffad50d3 100644 --- a/deal.II/lac/include/lac/slepc_solver.h +++ b/deal.II/lac/include/lac/slepc_solver.h @@ -35,13 +35,13 @@ DEAL_II_NAMESPACE_OPEN * Base class for solver classes using the SLEPc solvers which are * selected based on flags passed to the eigenvalue problem solver * context. Derived classes set the right flags to set the right - * solver. On the other hand, note that: the - * AdditionalData structure is a dummy structure that - * currently exists for backward/forward compatibility only. + * solver. Note that: the AdditionalData structure is a + * dummy structure that currently exists for backward/forward + * compatibility only. * * The SLEPc solvers are intended to be used for solving the - * generalized eigenspectrum problem $(A-\lambda M)x=0$, for $x\neq0$; - * where $A$ is a system matrix, $M$ is a mass matrix, and $\lambda, + * generalized eigenspectrum problem $(A-\lambda B)x=0$, for $x\neq0$; + * where $A$ is a system matrix, $B$ is a mass matrix, and $\lambda, * x$ are a set of eigenvalues and eigenvectors respectively. The * emphasis is on methods and techniques appropriate for problems in * which the associated matrices are sparse. Most of the methods @@ -53,11 +53,10 @@ DEAL_II_NAMESPACE_OPEN * following way: @verbatim SolverControl solver_control (1000, 1e-9); - SolverArnoldi system (solver_control, - mpi_communicator); - system.solve (A, M, lambda, x, size_of_spectrum); + SolverArnoldi system (solver_control, mpi_communicator); + system.solve (A, B, lambda, x, size_of_spectrum); @endverbatim - * for the generalized eigenvalue problem $Ax=M\lambda x$, where the + * for the generalized eigenvalue problem $Ax=B\lambda x$, where the * variable const unsigned int size_of_spectrum tells * SLEPc the number of eigenvector/eigenvalue pairs to solve for: See * also step-36 for a hands-on example. @@ -89,7 +88,7 @@ DEAL_II_NAMESPACE_OPEN * @author Toby D. Young 2008, 2009, 2010; and Rickard Armiento 2008. * * @note Various tweaks and enhancments contributed by Eloy Romero and - * Jose E. Roman 2009, 2010. + * Jose E. Roman 2009, 2010. */ namespace SLEPcWrappers @@ -320,10 +319,15 @@ namespace SLEPcWrappers */ virtual void set_solver_type (EPS &eps) const = 0; - /** + /* * Attributes that store the * relevant information for the - * eigenproblem solver context. + * eigenproblem solver + * context. These are: which + * portion of the spectrum to solve + * from; and the matrices and + * vectors that discribe the + * problem. */ EPSWhich set_which; @@ -331,6 +335,12 @@ namespace SLEPcWrappers const PETScWrappers::MatrixBase *opB; const PETScWrappers::VectorBase *initial_vector; + /** + * Pointer to an an object that + * describes transformations that + * can be applied to the eigenvalue + * problem. + */ SLEPcWrappers::TransformationBase *transformation; private: @@ -535,6 +545,56 @@ namespace SLEPcWrappers }; +/** + * An implementation of the solver interface using the SLEPc Lanczos + * solver. Usage: Largest values of spectrum only, all problem types, + * complex. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2010 + */ + class SolverPower : public SolverBase + { + public: + /** + * Standardized data struct to pipe + * additional data to the solver, + * should it be needed. + */ + struct AdditionalData + {}; + + /** + * SLEPc solvers will want to have + * an MPI communicator context over + * which computations are + * parallelized. By default, this + * carries the same behaviour has + * the PETScWrappers, but you can + * change that. + */ + SolverPower (SolverControl &cn, + const MPI_Comm &mpi_communicator = PETSC_COMM_WORLD, + const AdditionalData &data = AdditionalData()); + + protected: + + /** + * Store a copy of the flags for + * this particular solver. + */ + const AdditionalData additional_data; + + /** + * Function that takes a Eigenvalue + * Problem Solver context object, + * and sets the type of solver that + * is appropriate for this class. + */ + virtual void set_solver_type (EPS &eps) const; + + }; + // --------------------------- inline and template functions ----------- /** diff --git a/deal.II/lac/source/slepc_solver.cc b/deal.II/lac/source/slepc_solver.cc index d8da4f5b9b..930e874065 100644 --- a/deal.II/lac/source/slepc_solver.cc +++ b/deal.II/lac/source/slepc_solver.cc @@ -56,6 +56,7 @@ namespace SLEPcWrappers void SolverBase::set_matrices (const PETScWrappers::MatrixBase &A) { + // standard eigenspectrum problem opA = &A; opB = NULL; } @@ -64,6 +65,7 @@ namespace SLEPcWrappers SolverBase::set_matrices (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B) { + // generalized eigenspectrum problem opA = &A; opB = &B; } @@ -95,10 +97,12 @@ namespace SLEPcWrappers solver_data.reset (new SolverData()); // create eigensolver context and - // set operators. + // set operators ierr = EPSCreate (mpi_communicator, &solver_data->eps); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // set eigenspectrum problem type + // (general/standard) AssertThrow (opA, ExcSLEPcWrappersUsageError()); if (opB) ierr = EPSSetOperators (solver_data->eps, *opA, *opB); @@ -106,6 +110,7 @@ namespace SLEPcWrappers ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // set the initial vector(s) if any if (initial_vector && initial_vector->size() != 0) { @@ -119,12 +124,15 @@ namespace SLEPcWrappers AssertThrow (ierr == 0, ExcSLEPcError(ierr)); } + // set transformation type if any if (transformation) transformation->set_context(solver_data->eps); - // set runtime options. + // set runtime options set_solver_type (solver_data->eps); + // set which portion of the + // eigenspectrum to solve for ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); @@ -316,6 +324,33 @@ namespace SLEPcWrappers AssertThrow (ierr == 0, ExcSLEPcError(ierr)); } + /* ----------------------- Power ------------------------- */ + + SolverPower::SolverPower (SolverControl &cn, + const MPI_Comm &mpi_communicator, + const AdditionalData &data) + : + SolverBase (cn, mpi_communicator), + additional_data (data) + {} + + void + SolverPower::set_solver_type (EPS &eps) const + { + int ierr; + ierr = EPSSetType (eps, const_cast(EPSPOWER)); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // hand over the absolute + // tolerance and the maximum + // number of iteration steps to + // the SLEPc convergence + // criterion. + ierr = EPSSetTolerances(eps, this->solver_control.tolerance(), + this->solver_control.max_steps()); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + } DEAL_II_NAMESPACE_CLOSE