From 613ca24e852d6aae59d0e1f5ff42744140d83120 Mon Sep 17 00:00:00 2001 From: young Date: Tue, 23 Apr 2013 11:08:30 +0000 Subject: [PATCH] New function set_target_eigenvalue and minor API tweaks. git-svn-id: https://svn.dealii.org/trunk@29370 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/slepc_solver.h | 150 +++++++++++---------- deal.II/source/lac/slepc_solver.cc | 11 ++ 2 files changed, 92 insertions(+), 69 deletions(-) diff --git a/deal.II/include/deal.II/lac/slepc_solver.h b/deal.II/include/deal.II/lac/slepc_solver.h index 7985623114..d28d46c028 100644 --- a/deal.II/include/deal.II/lac/slepc_solver.h +++ b/deal.II/include/deal.II/lac/slepc_solver.h @@ -63,7 +63,7 @@ DEAL_II_NAMESPACE_OPEN * SLEPc solvers before calling solve(). For example, if * the matrices of the general eigenspectrum problem are not hermitian * and the lower eigenvalues are wanted only, the following code can - * be implemented: + * be implemented before calling solve(): * @code * system.set_problem_type (EPS_NHEP); * system.set_which_eigenpairs (EPS_SMALLEST_REAL); @@ -102,7 +102,6 @@ DEAL_II_NAMESPACE_OPEN * @note Various tweaks and enhancments contributed by Eloy Romero and * Jose E. Roman 2009, 2010. */ - namespace SLEPcWrappers { @@ -142,9 +141,9 @@ namespace SLEPcWrappers * Set up matrices for solving; Actually solve the system; Gather * the solution(s); and reset. * - * Note that the number of converged eigenvectors can be larger - * than the number of eigenvectors requested; this is due to a - * round off error (success) of the eigenproblem solver + * @note Note that the number of converged eigenvectors can be + * larger than the number of eigenvectors requested; this is due + * to a round off error (success) of the eigenproblem solver * context. If this is found to be the case we simply do not * bother with more eigenpairs than requested, but handle that it * may be more than specified by ignoring any extras. By default @@ -169,21 +168,6 @@ namespace SLEPcWrappers std::vector &vr, const unsigned int n_eigenvectors = 1); - /** - * Initialize solver for the linear system $Ax=\lambda x$. (Note: - * this is required before calling solve ()) - */ - void - set_matrices (const PETScWrappers::MatrixBase &A); - - /** - * Same as above, but here initialize solver for the linear system - * $A x=\lambda B x$. - */ - void - set_matrices (const PETScWrappers::MatrixBase &A, - const PETScWrappers::MatrixBase &B); - /** * Set the initial vector for the solver. */ @@ -197,6 +181,13 @@ namespace SLEPcWrappers void set_transformation (SLEPcWrappers::TransformationBase &this_transformation); + /** + * Set target eigenvalues in the spectrum to be computed. By + * default, no target is set. + */ + void + set_target_eigenvalue (const double &this_target); + /** * Indicate which part of the spectrum is to be computed. By * default largest magnitude eigenvalues are computed. @@ -217,52 +208,13 @@ namespace SLEPcWrappers void set_problem_type (EPSProblemType set_problem); - /** - * Solve the linear system for n_eigenvectors - * eigenstates. Parameter n_converged contains the - * actual number of eigenstates that have . converged; this can - * be both fewer or more than n_eigenvectors, depending on the - * SLEPc eigensolver used. - */ - void - solve (const unsigned int n_eigenvectors, unsigned int *n_converged); - - /** - * Access the solutions for a solved eigenvector problem, pair - * index solutions, $\text{index}\,\in\,0\hdots - * \text{n\_converged}-1$. - */ - void - get_eigenpair (const unsigned int index, -#ifndef PETSC_USE_COMPLEX - double &kr, -#else - std::complex &kr, -#endif - PETScWrappers::VectorBase &vr); - - /** - * Reset the solver, and return memory for eigenvectors - */ - void - reset (); - - /** - * Retrieve the SLEPc solver object that is internally used. - */ - EPS *get_eps (); - /** * Take the information provided from SLEPc and checks it against * deal.II's own SolverControl objects to see if convergence has * been reached. */ - void get_solver_state (const SolverControl::State state); - - /** - * Access to the object that controls convergence. - */ - SolverControl &control () const; + void + get_solver_state (const SolverControl::State state); /** * Exception. Standard exception. @@ -285,6 +237,11 @@ namespace SLEPcWrappers << " The number of converged eigenvectors is " << arg1 << " but " << arg2 << " were requested. "); + /** + * Access to the object that controls convergence. + */ + SolverControl &control () const; + protected: /** @@ -305,6 +262,61 @@ namespace SLEPcWrappers */ virtual void set_solver_type (EPS &eps) const = 0; + /** + * Reset the solver, and return memory for eigenvectors + */ + void + reset (); + + /** + * Retrieve the SLEPc solver object that is internally used. + */ + EPS *get_eps (); + + /** + * Solve the linear system for n_eigenvectors + * eigenstates. Parameter n_converged contains the + * actual number of eigenstates that have . converged; this can + * be both fewer or more than n_eigenvectors, depending on the + * SLEPc eigensolver used. + */ + void + solve (const unsigned int n_eigenvectors, unsigned int *n_converged); + + /** + * Access the solutions for a solved eigenvector problem, pair + * index solutions, $\text{index}\,\in\,0\hdots + * \text{n\_converged}-1$. + */ + void + get_eigenpair (const unsigned int index, +#ifndef PETSC_USE_COMPLEX + double &kr, +#else + std::complex &kr, +#endif + PETScWrappers::VectorBase &vr); + + /** + * Initialize solver for the linear system $Ax=\lambda x$. (Note: + * this is required before calling solve ()) + */ + void + set_matrices (const PETScWrappers::MatrixBase &A); + + /** + * Same as above, but here initialize solver for the linear system + * $A x=\lambda B x$. + */ + void + set_matrices (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B); + + /** + * Target eigenvalue to solve for. + */ + PetscScalar target_eigenvalue; + /** * Which portion of the spectrum to solve from. */ @@ -315,6 +327,8 @@ namespace SLEPcWrappers */ EPSProblemType set_problem; + private: + /** * The matrix $A$ of the generalized eigenvalue problem * $Ax=B\lambda x$, or the standard eigenvalue problem $Ax=\lambda @@ -339,8 +353,6 @@ namespace SLEPcWrappers */ SLEPcWrappers::TransformationBase *transformation; - private: - /** * A function that can be used in SLEPc as a callback to check on * convergence. @@ -718,7 +730,7 @@ namespace SLEPcWrappers unsigned int n_converged = 0; // Set the matrices of the problem - set_matrices (A); + set_matrices (A); // and solve solve (n_eigenvectors, &n_converged); @@ -738,11 +750,11 @@ namespace SLEPcWrappers template void - SolverBase::solve (const PETScWrappers::MatrixBase &A, - const PETScWrappers::MatrixBase &B, - std::vector &kr, - std::vector &vr, - const unsigned int n_eigenvectors) + SolverBase::solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors) { // Panic if no eigenpairs are wanted. AssertThrow (n_eigenvectors != 0, ExcSLEPcWrappersUsageError()); diff --git a/deal.II/source/lac/slepc_solver.cc b/deal.II/source/lac/slepc_solver.cc index c76b6c4331..ded535c6bc 100644 --- a/deal.II/source/lac/slepc_solver.cc +++ b/deal.II/source/lac/slepc_solver.cc @@ -49,6 +49,7 @@ namespace SLEPcWrappers : solver_control (cn), mpi_communicator (mpi_communicator), + target_eigenvalue (PETSC_NULL), set_which (EPS_LARGEST_MAGNITUDE), set_problem (EPS_NHEP), opA (NULL), @@ -89,6 +90,12 @@ namespace SLEPcWrappers transformation = &this_transformation; } + void + SolverBase::set_target_eigenvalue (const double &this_target) + { + target_eigenvalue = this_target; + } + void SolverBase::set_which_eigenpairs (const EPSWhich eps_which) { @@ -147,6 +154,10 @@ namespace SLEPcWrappers if (transformation) transformation->set_context (solver_data->eps); + // set target eigenvalues to solve for + ierr = EPSSetTarget (solver_data->eps, target_eigenvalue); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // set which portion of the eigenspectrum to solve for ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); -- 2.39.5