From: young Date: Fri, 18 Jan 2013 18:38:44 +0000 (+0000) Subject: Improve SLEPcWrappers. Fix usage of SolverControl. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3e984f83bf07784ab7ad95f112e033b878714413;p=dealii-svn.git Improve SLEPcWrappers. Fix usage of SolverControl. git-svn-id: https://svn.dealii.org/trunk@28132 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/slepc_solver.h b/deal.II/include/deal.II/lac/slepc_solver.h index 4e1b09fbc1..7e6593810c 100644 --- a/deal.II/include/deal.II/lac/slepc_solver.h +++ b/deal.II/include/deal.II/lac/slepc_solver.h @@ -191,11 +191,11 @@ namespace SLEPcWrappers set_which_eigenpairs (EPSWhich set_which); /** - * Solve the linear system for n_eigenvectors + * 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. + * 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); @@ -218,13 +218,19 @@ namespace SLEPcWrappers * Reset the solver, and return memory for eigenvectors */ void - reset(); + reset (); /** * Retrieve the SLEPc solver object that is internally used. */ - EPS * - get_EPS (); + 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. @@ -232,15 +238,21 @@ namespace SLEPcWrappers SolverControl &control () const; /** - * Exceptions. + * Exception. Standard exception. */ DeclException0 (ExcSLEPcWrappersUsageError); + /** + * Exception. SLEPc error with error number. + */ DeclException1 (ExcSLEPcError, int, << " An error with error number " << arg1 << " occurred while calling a SLEPc function"); + /** + * Exception. Convergence failure on the number of eigenvectors. + */ DeclException2 (ExcSLEPcEigenvectorConvergenceMismatchError, int, int, << " The number of converged eigenvectors is " << arg1 @@ -298,23 +310,20 @@ namespace SLEPcWrappers private: /** - * A function that is used in SLEPc as a callback to check on - * convergence. It takes the information provided from SLEPc and - * checks it against deal.II's own SolverControl objects to see if - * convergence has been reached. + * A function that can be used in SLEPc as a callback to check on + * convergence. + * + * @note This function is redundant. */ - static + static int - convergence_test (EPS eps, -#ifdef PETSC_USE_64BIT_INDICES - const PetscInt iteration, -#else - const int iteration, -#endif - const PetscReal residual_norm, - EPSConvergedReason *reason, - void *solver_control); - + convergence_test (EPS eps, + PetscScalar kr, + PetscScalar ki, + PetscReal residual_norm, + PetscReal *estimated_error, + void *solver_control); + /** * Objects of this type are explicitly created, but are destroyed * when the surrounding solver object goes out of scope, or when @@ -335,6 +344,11 @@ namespace SLEPcWrappers * Objects for Eigenvalue Problem Solver. */ EPS eps; + + /** + * Convergence. + */ + EPSConvergedReason reason; }; /** @@ -349,6 +363,7 @@ namespace SLEPcWrappers * complex. * * @ingroup SLEPcWrappers + * * @author Toby D. Young 2008 */ class SolverKrylovSchur : public SolverBase @@ -391,6 +406,7 @@ namespace SLEPcWrappers * solver. Usage: All spectrum, all problem types, complex. * * @ingroup SLEPcWrappers + * * @author Toby D. Young 2008, 2011 */ class SolverArnoldi : public SolverBase @@ -436,7 +452,6 @@ namespace SLEPcWrappers * and sets the type of solver that is appropriate for this class. */ virtual void set_solver_type (EPS &eps) const; - }; /** @@ -444,6 +459,7 @@ namespace SLEPcWrappers * solver. Usage: All spectrum, all problem types, complex. * * @ingroup SLEPcWrappers + * * @author Toby D. Young 2009 */ class SolverLanczos : public SolverBase @@ -478,7 +494,6 @@ namespace SLEPcWrappers * and sets the type of solver that is appropriate for this class. */ virtual void set_solver_type (EPS &eps) const; - }; /** @@ -486,7 +501,8 @@ namespace SLEPcWrappers * solver. Usage: Largest values of spectrum only, all problem * types, complex. * - * @ingroup SLEPcWrappers + * @ingroup SLEPcWrappers + * * @author Toby D. Young 2010 */ class SolverPower : public SolverBase @@ -521,14 +537,15 @@ namespace SLEPcWrappers * and sets the type of solver that is appropriate for this class. */ virtual void set_solver_type (EPS &eps) const; - }; /** * An implementation of the solver interface using the SLEPc - * Generalized Davidson solver. Usage: All problem types. + * Davidson solver. Usage (incomplete/untested): All problem types. * - * @ingroup SLEPcWrappers @author Toby D. Young 2013 + * @ingroup SLEPcWrappers + * + * @author Toby D. Young 2010 */ class SolverGeneralizedDavidson : public SolverBase { @@ -562,14 +579,15 @@ namespace SLEPcWrappers * and sets the type of solver that is appropriate for this class. */ virtual void set_solver_type (EPS &eps) const; - }; /** * An implementation of the solver interface using the SLEPc * Jacobi-Davidson solver. Usage: All problem types. * - * @ingroup SLEPcWrappers @author Toby D. Young 2013 + * @ingroup SLEPcWrappers + * + * @author Toby D. Young 2013 */ class SolverJacobiDavidson : public SolverBase { @@ -603,7 +621,6 @@ namespace SLEPcWrappers * and sets the type of solver that is appropriate for this class. */ virtual void set_solver_type (EPS &eps) const; - }; // --------------------------- inline and template functions ----------- @@ -619,21 +636,24 @@ namespace SLEPcWrappers std::vector &vr, const unsigned int n_eigenvectors = 1) { - unsigned int n_converged; + unsigned int n_converged = 0; + // Set the matrices of the problem set_matrices (A); - solve (n_eigenvectors,&n_converged); + // and solve + solve (n_eigenvectors, &n_converged); + if (n_converged > n_eigenvectors) n_converged = n_eigenvectors; AssertThrow (n_converged == n_eigenvectors, ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors)); - AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError()); + AssertThrow (vr.size() != 0, ExcSLEPcWrappersUsageError()); vr.resize (n_converged, vr.front()); kr.resize (n_converged); - for (unsigned int index=0; index < n_converged; ++index) + for (unsigned int index=0; index &vr, const unsigned int n_eigenvectors = 1) { - unsigned int n_converged; + int ierr; + unsigned int n_converged = 0; + // Set the matrices of the problem set_matrices (A, B); + + // and solve solve (n_eigenvectors, &n_converged); if (n_converged >= n_eigenvectors) @@ -655,12 +679,12 @@ namespace SLEPcWrappers AssertThrow (n_converged == n_eigenvectors, ExcSLEPcEigenvectorConvergenceMismatchError(n_converged, n_eigenvectors)); - AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError()); + AssertThrow (vr.size() != 0, ExcSLEPcWrappersUsageError()); vr.resize (n_converged, vr.front()); kr.resize (n_converged); - for (unsigned int index=0; index < n_converged; ++index) + for (unsigned int index=0; indexeps, *opA, PETSC_NULL); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // set runtime options + set_solver_type (solver_data->eps); + // set the initial vector(s) if any if (initial_vector && initial_vector->size() != 0) { @@ -129,10 +133,7 @@ namespace SLEPcWrappers // set transformation type if any if (transformation) - transformation->set_context(solver_data->eps); - - // set runtime options - set_solver_type (solver_data->eps); + transformation->set_context (solver_data->eps); // set which portion of the eigenspectrum to solve for ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which); @@ -147,6 +148,15 @@ namespace SLEPcWrappers ierr = EPSSetFromOptions (solver_data->eps); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // Set convergence test to be absolute + ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // Set the convergence test function + // ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test, + // reinterpret_cast(&solver_control)); + // AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + // solve the eigensystem ierr = EPSSolve (solver_data->eps); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); @@ -154,11 +164,39 @@ namespace SLEPcWrappers // get number of converged eigenstates ierr = EPSGetConverged (solver_data->eps, #ifdef PETSC_USE_64BIT_INDICES - reinterpret_cast(n_converged)); + reinterpret_cast(n_converged) #else - reinterpret_cast(n_converged)); + reinterpret_cast(n_converged) #endif + ); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + int n_iterations = 0; + double residual_norm = 1e300; + + // @todo Investigate elaborating on some of this to act on the + // complete eigenspectrum + { + // get the number of solver iterations + ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // get the residual norm of the most extreme eigenvalue + ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // check the solver state + const SolverControl::State state + = solver_control.check (n_iterations, residual_norm); + + // get the solver state according to SLEPc + get_solver_state (state); + + // and in case of failure: throw exception + if (solver_control.last_check () != SolverControl::success) + throw SolverControl::NoConvergence (solver_control.last_step (), + solver_control.last_value ()); + } } void @@ -178,7 +216,6 @@ namespace SLEPcWrappers AssertThrow (ierr == 0, ExcSLEPcError(ierr)); } - void SolverBase::reset () { @@ -189,14 +226,39 @@ namespace SLEPcWrappers } EPS * - SolverBase::get_EPS () + SolverBase::get_eps () { - if (solver_data.get() == 0) + if (solver_data.get () == 0) return NULL; - + return &solver_data->eps; } + void + SolverBase::get_solver_state (const SolverControl::State state) + { + switch (state) + { + case ::dealii::SolverControl::iterate: + solver_data->reason = EPS_CONVERGED_ITERATING; + break; + + case ::dealii::SolverControl::success: + solver_data->reason = static_cast(1); + break; + + case ::dealii::SolverControl::failure: + if (solver_control.last_step() > solver_control.max_steps()) + solver_data->reason = EPS_DIVERGED_ITS; + else + solver_data->reason = EPS_DIVERGED_BREAKDOWN; + break; + + default: + Assert (false, ExcNotImplemented()); + } + } + /* ---------------------- SolverControls ----------------------- */ SolverControl & SolverBase::control () const @@ -205,42 +267,14 @@ namespace SLEPcWrappers } int - SolverBase::convergence_test (EPS /*eps*/, -#ifdef PETSC_USE_64BIT_INDICES - const PetscInt iteration, -#else - const int iteration, -#endif - const PetscReal residual_norm, - EPSConvergedReason *reason, - void *solver_control_x) + SolverBase::convergence_test (EPS /*eps */, + PetscScalar /*kr */, + PetscScalar /*ki */, + PetscReal /*residual_norm */, + PetscReal */*estimated_error */, + void */*solver_control_x*/) { - SolverControl &solver_control - = *reinterpret_cast(solver_control_x); - - const SolverControl::State state - = solver_control.check (iteration, residual_norm); - - switch (state) - { - case ::dealii::SolverControl::iterate: - *reason = EPS_CONVERGED_ITERATING; - break; - - case ::dealii::SolverControl::success: - *reason = static_cast(1); - break; - - case ::dealii::SolverControl::failure: - if (solver_control.last_step() > solver_control.max_steps()) - *reason = EPS_DIVERGED_ITS; - else - *reason = EPS_DIVERGED_BREAKDOWN; - break; - - default: - Assert (false, ExcNotImplemented()); - } + // This function is undefined (future reference only). // return without failure. return 0;