From: Toby D. Young Date: Wed, 24 Jun 2009 12:24:08 +0000 (+0000) Subject: Added SLEPcWrappers which give a handle on SLEPc: the resultant two base classes... X-Git-Tag: v8.0.0~7593 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=83e3f3ef12ca0bbc68ba8baecbcc4d4ed70be54f;p=dealii.git Added SLEPcWrappers which give a handle on SLEPc: the resultant two base classes (sover base and transformation base) can be used for solving the generalized eigenvalue problem git-svn-id: https://svn.dealii.org/trunk@18957 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index eed45d43f2..e7b452b24c 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -129,6 +129,15 @@ inconvenience this causes. +
    +
  1. +

    + New: Based on work with Rickard Armiento, Francisco Alvaro, and Jose E. Roman, SLEPcWrappers that give a handle on some of the features of SLEPc (Scalable Library for Eigenvalue Problem Computations): (1) The SLEPcWrappers::SolverBase class can be used for specifying an eigenvalue problem, either in standard or generalized form, on serial or parallel architectures with support for a few solver types; and (2) The SLEPcWrappers::TransformationBase class encapsulates a variety of spectral transformations providing some functionality required for acceleration techniques based on the transformation of the spectrum. +
    + (Toby D. Young 2009/06/25) +

    +
  2. +
diff --git a/deal.II/lac/include/lac/slepc_solver.h b/deal.II/lac/include/lac/slepc_solver.h new file mode 100755 index 0000000000..fef1bbeffd --- /dev/null +++ b/deal.II/lac/include/lac/slepc_solver.h @@ -0,0 +1,510 @@ + +#ifndef __deal2__slepc_solver_h +#define __deal2__slepc_solver_h + +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_SLEPC + +#include +#include + +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 -- in particular, note that the AdditionalData structure is + * a dummy structure and is there for backward compatibility. + * + * SLEPcWrappers can be implemented in application codes in the + * following way: + @verbatim + SolverControl solver_control (1000, 1e-10); + SolverArnoldi system (solver_control, + mpi_communicator); + system.solve (A, B, lambda, x, n_eigenvectors); + @endverbatim + * for the generalized eigenvalue problem $Ax=B\lambda x$. + * + * An alternative implementation to the one above is to use the API + * internals directly within the application code. In this way the + * calling sequence requires calling several of SolverBase functions + * rather than just one. This freedom is intended for use by users of + * the SLEPcWrappers that require a greater handle on the eigenvalue + * problem solver context. See also: + @verbatim + template + void + SolverBase::solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors + ) + @endverbatim + * as an example on how to do this. + * + * See also the @ref PETScWrappers "PETScWrappers", on which the + * SLEPcWrappers depend, for additional explanations. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2008, 2009 + */ +namespace SLEPcWrappers +{ + + /** + * Base class for solver classes + * using the SLEPc solvers. Since + * solvers in SLEPc are selected + * based on flags passed to a + * generic solver object, basically + * all the actual solver calls + * happen in this class, and + * derived classes simply set the + * right flags to select one solver + * or another, or to set certain + * parameters for individual + * solvers. + */ + class SolverBase + { + public: + + /** + * Constructor. Takes the MPI + * communicator over which parallel + * computations are to happen. + */ + SolverBase (SolverControl &cn, + const MPI_Comm &mpi_communicator); + + /** + * Destructor. + */ + virtual ~SolverBase (); + + /** + * Composite method that solves the + * linear system $Ax=\lambda + * Bx$. The eigenvector sent in has + * to have at least one element + * that we can use as a template + * when resizing, since we do not + * know the parameters of the + * specific vector class used + * (i.e. local_dofs for MPI + * vectors). However, while copying + * eigenvectors, at least twice the + * memory size of vr is + * being used (and can be more). To + * avoid doing this, use instead + * the calling sequence used here + * is fairly standard: Initialise; + * set up matrices for solving; + * actually solve the system; + * gather the solution(s); and + * reset. + * + * Note that the number of + * converged eigenstates can be + * larger than the number of + * eigenstates requested; this is + * due to a round off error + * (success) of the eigenvalue + * 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 by ignoring any extras. + */ + template + void + solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors); + + /** + * Initialize solver for the linear + * system $Ax=\lambda + * Bx$. (required before calling + * solve) + */ + void + set_matrices (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B); + + /** + * Set the initial vector for the solver. + */ + void + set_initial_vector (const PETScWrappers::VectorBase &initial_vec); + + /** + * Set the spectral transformation + * to be used. By default SLEPc + */ + void + set_transformation (SLEPcWrappers::TransformationBase &trans); + + /** + * Indicate which part of the + * spectrum is to be computed. By + * default largest magnitude + * eigenvalues are computed. For + * other allowed values see the + * SLEPc documentation. + */ + void + set_which_eigenpairs (EPSWhich set_which); + + /** + * 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, index = 0 + * ... n_converged-1 + */ + void + get_eigenpair (const unsigned int index, + double &kr, + PETScWrappers::VectorBase &vr); + + /** + * Reset the solver, and return + * memory for eigenvectors + */ + void + reset(); + + /** + * Retrieve the SLEPc solver object + * used internally. + */ + EPS * + get_EPS (); + + + /** + * Access to object that controls + * convergence. + */ + SolverControl & control() const; + + /** + * Exceptions. + */ + DeclException0 (ExcSLEPcWrappersUsageError); + DeclException1 (ExcSLEPcError, + int, + << "An error with error number " << arg1 + << " occured while calling a SLEPc function"); + + protected: + + /** + * Reference to the object that + * controls convergence of the + * iterative solver. + */ + SolverControl &solver_control; + + /** + * Copy of the MPI communicator + * object to be used for the + * solver. + */ + const MPI_Comm mpi_communicator; + + /** + * Function that takes an + * Eigenvalue Problem Solver + * context object, and sets the + * type of solver that is requested + * by the derived class. + */ + virtual void set_solver_type (EPS &eps) const = 0; + + /** + * Attributes that store the + * relevant information for the + * eigenproblem solver context. + */ + EPSWhich set_which; + const PETScWrappers::MatrixBase *opA; + const PETScWrappers::MatrixBase *opB; + const PETScWrappers::VectorBase *ini_vec; + SLEPcWrappers::TransformationBase *transform; + + 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. + */ + static + int + convergence_test (EPS eps, + const int iteration, + const PetscScalar residual_norm, + EPSConvergedReason *reason, + void *solver_control); + + /** + * Objects of this type are + * explicitly created, but are + * destroyed when the surrounding + * solver object goes out of scope, + * or when we assign a new value to + * the pointer to this object. The + * respective Destroy functions are + * therefore written into the + * destructor of this object, even + * though the object does not have + * a constructor. + */ + struct SolverData + { + + /** + * Destructor. + */ + ~SolverData (); + + /** + * Objects for Eigenvalue Problem + * Solver. + */ + EPS eps; + }; + + boost::shared_ptr solver_data; + }; + +/** + * An implementation of the solver interface using the SLEPc + * Krylov-Schur solver. Usage: All spectrum, all problem types, + * complex. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2008 + */ + class SolverKrylovSchur : 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. + */ + SolverKrylovSchur (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; + }; + +/** + * An implementation of the solver interface using the SLEPc Arnoldi + * solver. Usage: All spectrum, all problem types, complex. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2008 + */ + class SolverArnoldi : 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. + */ + SolverArnoldi (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; + + }; + +/** + * An implementation of the solver interface using the SLEPc Lanczos + * solver. Usage: All spectrum, all problem types, complex. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + */ + class SolverLanczos : 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. + */ + SolverLanczos (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 ----------- + + /** + * This is declared here to make it + * possible to take a std::vector + * of different PETScWrappers vector + * types + */ + template + void + SolverBase::solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors = 0) + { + unsigned int n_converged; + + set_matrices(A,B); + + solve(n_eigenvectors,&n_converged); + + if (n_converged > n_eigenvectors) + { + n_converged = n_eigenvectors; + } + + AssertThrow (vr.size() >= 1, ExcSLEPcWrappersUsageError()); + vr.resize(n_converged, vr.front()); + kr.resize(n_converged); + + for (unsigned int index=0; index < n_converged; + ++index) + { + get_eigenpair(index, kr[index], vr[index]); + } + } + + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_SLEPC + +/*---------------------------- slepc_solver.h ---------------------------*/ + +#endif + +/*---------------------------- slepc_solver.h ---------------------------*/ + diff --git a/deal.II/lac/include/lac/slepc_spectral_transformation.h b/deal.II/lac/include/lac/slepc_spectral_transformation.h new file mode 100755 index 0000000000..0b229a3386 --- /dev/null +++ b/deal.II/lac/include/lac/slepc_spectral_transformation.h @@ -0,0 +1,317 @@ + +#ifndef __deal2__slepc_spectral_transformation_h +#define __deal2__slepc_spectral_transformation_h + +#include +#include +#include +#include +#include + +#ifdef DEAL_II_USE_SLEPC + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace SLEPcWrappers +{ + +/** + * Base class for spectral transformation classes using the SLEPc + * solvers which are selected based on flags passed to the spectral + * transformation. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + **/ + class TransformationBase + { + public: + + /** + * Constructor. Takes the MPI + * communicator over which parallel + * computations are to happen. + */ + TransformationBase (); + + /** + * Destructor. + */ + virtual ~TransformationBase (); + + /** + * Record the EPS object that is associated + * to the spectral transformation + */ + void set_context (EPS &eps); + + protected: + + virtual void set_transformation_type (ST &st) const = 0; + + private: + + /** + * Objects of this type are + * explicitly created, but are + * destroyed when the surrounding + * solver object goes out of scope, + * or when we assign a new value to + * the pointer to this object. The + * respective Destroy functions are + * therefore written into the + * destructor of this object, even + * though the object does not have + * a constructor. + */ + struct TransformationData + { + + /** + * Destructor. + */ + ~TransformationData (); + + /** + * Objects for Eigenvalue Problem + * Solver. + */ + ST st; + }; + + boost::shared_ptr transformation_data; + }; + +/** + * An implementation of the transformation interface using the SLEPc + * Shift. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + */ + class TransformationShift : public TransformationBase + { + public: + + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + struct AdditionalData + { + + /** + * Constructor. By default, set the + * shift parameter to zero. + */ + AdditionalData (const double shift_parameter = 0); + + /** + * Shift parameter. + */ + const double shift_parameter; + }; + + + /** + * Constructor. + */ + TransformationShift (const AdditionalData &data = AdditionalData()); + + + protected: + + /** + * Store a copy of the flags for this + * particular solver. + */ + const AdditionalData additional_data; + + /** + * Function that takes a Spectral + * Transformation context object, + * and sets the type of spectral + * transformationthat is + * appropriate for this class. + */ + virtual void set_transformation_type (ST &st) const; + }; + +/** + * An implementation of the transformation interface using the SLEPc + * Shift and Invert. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + */ + class TransformationShiftInvert : public TransformationBase + { + public: + + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + struct AdditionalData + { + /** + * Constructor. By default, set the + * shift parameter to zero. + */ + AdditionalData (const double shift_parameter = 0); + + /** + * Shift parameter. + */ + const double shift_parameter; + }; + + + /** + * Constructor. + */ + TransformationShiftInvert (const AdditionalData &data = AdditionalData()); + + protected: + + /** + * Store a copy of the flags for this + * particular solver. + */ + const AdditionalData additional_data; + + /** + * Function that takes a Spectral + * Transformation context object, + * and sets the type of spectral + * transformationthat is + * appropriate for this class. + */ + virtual void set_transformation_type (ST &st) const; + }; + +/** + * An implementation of the transformation interface using the SLEPc + * Spectrum Folding. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + */ + class TransformationSpectrumFolding : public TransformationBase + { + public: + + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + struct AdditionalData + { + /** + * Constructor. By default, set the + * shift parameter to zero. + */ + AdditionalData (const double shift_parameter = 0); + + /** + * Shift parameter. + */ + const double shift_parameter; + }; + + + /** + * Constructor. + */ + TransformationSpectrumFolding (const AdditionalData &data = AdditionalData()); + + protected: + + /** + * Store a copy of the flags for this + * particular solver. + */ + const AdditionalData additional_data; + + /** + * Function that takes a Spectral + * Transformation context object, + * and sets the type of spectral + * transformationthat is + * appropriate for this class. + */ + virtual void set_transformation_type (ST &st) const; + }; + +/** + * An implementation of the transformation interface using the SLEPc + * Cayley. + * + * @ingroup SLEPcWrappers + * @author Toby D. Young 2009 + */ + class TransformationCayley : public TransformationBase + { + public: + + /** + * Standardized data struct to + * pipe additional data to the + * solver. + */ + struct AdditionalData + { + /** + * Constructor. Requires two shift parameters + */ + AdditionalData (const double shift_parameter = 0, const double antishift_parameter = 0); + + /** + * Shift and antishift parameter. + */ + const double shift_parameter; + const double antishift_parameter; + }; + + + /** + * Constructor. + */ + TransformationCayley (const double shift, const double antishift); + + protected: + + /** + * Store a copy of the flags for this + * particular solver. + */ + const AdditionalData additional_data; + + /** + * Function that takes a Spectral + * Transformation context object, + * and sets the type of spectral + * transformationthat is + * appropriate for this class. + */ + virtual void set_transformation_type (ST &st) const; + }; + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_SLEPC + +/*-------------------- slepc_spectral_transformation.h ------------------*/ + +#endif + +/*-------------------- slepc_spectral_transformation.h ------------------*/ diff --git a/deal.II/lac/source/slepc_solver.cc b/deal.II/lac/source/slepc_solver.cc new file mode 100644 index 0000000000..f9042bece1 --- /dev/null +++ b/deal.II/lac/source/slepc_solver.cc @@ -0,0 +1,289 @@ + +#include +#include +#include +#include +#include + +#include +#include + +#ifdef DEAL_II_USE_SLEPC + +#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2) +#include +#endif +#include + +DEAL_II_NAMESPACE_OPEN + +namespace SLEPcWrappers +{ + + SolverBase::SolverData::~SolverData () + { + // Destroy the solver object. + int ierr = EPSDestroy (eps); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + + SolverBase::SolverBase (SolverControl &cn, + const MPI_Comm &mpi_communicator) + : + solver_control (cn), + mpi_communicator (mpi_communicator), + set_which (EPS_LARGEST_MAGNITUDE), + opA (NULL), opB (NULL), + ini_vec (NULL), + transform (NULL) + { + } + + SolverBase::~SolverBase () + { + if( solver_data != 0 ) + solver_data.reset (); + } + + void + SolverBase::set_matrices (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B) + { + opA = &A; + opB = &B; + } + + void + SolverBase::set_initial_vector (const PETScWrappers::VectorBase &initial_vec) + { + ini_vec = &initial_vec; + } + + void + SolverBase::set_transformation (SLEPcWrappers::TransformationBase &trans) + { + transform = &trans; + } + + void + SolverBase::set_which_eigenpairs (const EPSWhich eps_which) + { + set_which = eps_which; + } + + void + SolverBase::solve (const unsigned int n_eigenvectors, unsigned int *n_converged) + { + int ierr; + + AssertThrow (solver_data == 0, ExcSLEPcWrappersUsageError()); + solver_data.reset (new SolverData()); + + // create eigensolver context and + // set operators. + ierr = EPSCreate (mpi_communicator, &solver_data->eps); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + AssertThrow (opA && opB, ExcSLEPcWrappersUsageError()); + ierr = EPSSetOperators (solver_data->eps, *opA, *opB); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + if( ini_vec && ini_vec->size() != 0 ) + { + ierr = EPSSetInitialVector(solver_data->eps, *ini_vec); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + + if( transform ) + transform->set_context(solver_data->eps); + + // set runtime options. + set_solver_type (solver_data->eps); + + ierr = EPSSetWhichEigenpairs (solver_data->eps, set_which); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // set number of eigenvectors to + // compute + ierr = EPSSetDimensions (solver_data->eps, n_eigenvectors, + PETSC_DECIDE, PETSC_DECIDE); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + ierr = EPSSetFromOptions (solver_data->eps); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // solve the eigensystem + ierr = EPSSolve (solver_data->eps); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + + // get number of converged + // eigenstates + ierr = EPSGetConverged (solver_data->eps, reinterpret_cast(n_converged)); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + + void + SolverBase::get_eigenpair (const unsigned int index, + double &kr, + PETScWrappers::VectorBase &vr) + { + AssertThrow (solver_data != 0, ExcSLEPcWrappersUsageError()); + + // get converged eigenpair + int ierr = EPSGetEigenpair(solver_data->eps, index, + &kr, PETSC_NULL, vr, PETSC_NULL); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + + + void + SolverBase::reset () + { + AssertThrow (solver_data != 0, ExcSLEPcWrappersUsageError()); + + // destroy solver object. + solver_data.reset (); + } + + EPS * + SolverBase::get_EPS () + { + if( solver_data == 0 ) + return NULL; + return &solver_data->eps; + } + + /* ---------------------- SolverControls ----------------------- */ + SolverControl & + SolverBase::control () const + { + return solver_control; + } + + int + SolverBase::convergence_test (EPS /*eps*/, + const int iteration, + const PetscScalar residual_norm, + EPSConvergedReason *reason, + 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; + break; + + default: + Assert (false, ExcNotImplemented()); + } + + // return without failure. + return 0; + } + + /* ---------------------- SolverKrylovSchur ------------------------ */ + SolverKrylovSchur::SolverKrylovSchur (SolverControl &cn, + const MPI_Comm &mpi_communicator, + const AdditionalData &data) + : + SolverBase (cn, mpi_communicator), + additional_data (data) + {} + + void + SolverKrylovSchur::set_solver_type (EPS &eps) const + { + int ierr; + ierr = EPSSetType (eps, const_cast(EPSKRYLOVSCHUR)); + 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)); + } + + /* ---------------------- SolverArnoldi ------------------------ */ + SolverArnoldi::SolverArnoldi (SolverControl &cn, + const MPI_Comm &mpi_communicator, + const AdditionalData &data) + : + SolverBase (cn, mpi_communicator), + additional_data (data) + {} + + void + SolverArnoldi::set_solver_type (EPS &eps) const + { + int ierr; + ierr = EPSSetType (eps, const_cast(EPSARNOLDI)); + 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)); + } + + /* ---------------------- Lanczos ------------------------ */ + + SolverLanczos::SolverLanczos (SolverControl &cn, + const MPI_Comm &mpi_communicator, + const AdditionalData &data) + : + SolverBase (cn, mpi_communicator), + additional_data (data) + {} + + void + SolverLanczos::set_solver_type (EPS &eps) const + { + int ierr; + ierr = EPSSetType (eps, const_cast(EPSLANCZOS)); + 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 + +#else +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {} } +#endif // DEAL_II_USE_SLEPC + diff --git a/deal.II/lac/source/slepc_spectral_transformation.cc b/deal.II/lac/source/slepc_spectral_transformation.cc new file mode 100644 index 0000000000..6b4839c39f --- /dev/null +++ b/deal.II/lac/source/slepc_spectral_transformation.cc @@ -0,0 +1,154 @@ + +#include +#include +#include +#include +#include + +#include +#include + +#ifdef DEAL_II_USE_SLEPC + +#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2) +# include +#endif +#include + +DEAL_II_NAMESPACE_OPEN + +namespace SLEPcWrappers +{ + TransformationBase::TransformationData::~TransformationData () + { + } + + TransformationBase::TransformationBase () + { + } + + TransformationBase::~TransformationBase () + { + } + + void TransformationBase::set_context (EPS &eps) + { + AssertThrow (transformation_data == 0, + SolverBase::ExcSLEPcWrappersUsageError()); + transformation_data.reset (new TransformationData()); + + int ierr = EPSGetST(eps, &transformation_data->st); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + set_transformation_type(transformation_data->st); + } + + /* ------------------- TransformationShift --------------------- */ + TransformationShift::AdditionalData:: + AdditionalData (const double shift_parameter) + : + shift_parameter (shift_parameter) + {} + + TransformationShift::TransformationShift (const AdditionalData &data) + : + additional_data (data) + {} + + void + TransformationShift::set_transformation_type (ST &st) const + { + int ierr; + ierr = STSetType (st, const_cast(STSHIFT)); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + ierr = STSetShift (st, additional_data.shift_parameter); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + } + + /* ---------------- TransformationShiftInvert ------------------ */ + TransformationShiftInvert::AdditionalData:: + AdditionalData (const double shift_parameter) + : + shift_parameter (shift_parameter) + {} + + TransformationShiftInvert::TransformationShiftInvert (const AdditionalData &data) + : + additional_data (data) + {} + + void + TransformationShiftInvert::set_transformation_type (ST &st) const + { + int ierr; + ierr = STSetType (st, const_cast(STSINV)); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + ierr = STSetShift (st, additional_data.shift_parameter); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + } + + /* --------------- TransformationSpectrumFolding ----------------- */ + TransformationSpectrumFolding::AdditionalData:: + AdditionalData (const double shift_parameter) + : + shift_parameter (shift_parameter) + {} + + TransformationSpectrumFolding::TransformationSpectrumFolding (const AdditionalData &data) + : + additional_data (data) + {} + + + void + TransformationSpectrumFolding::set_transformation_type (ST &st) const + { + int ierr; + ierr = STSetType (st, const_cast(STFOLD)); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + ierr = STSetShift (st, additional_data.shift_parameter); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + } + + /* ------------------- TransformationCayley --------------------- */ + TransformationCayley::AdditionalData:: + AdditionalData (const double shift_parameter, + const double antishift_parameter) + : + shift_parameter (shift_parameter), + antishift_parameter (antishift_parameter) + { + } + + TransformationCayley::TransformationCayley (const double shift, + const double antishift) + : + additional_data (shift, antishift) + {} + + void + TransformationCayley::set_transformation_type (ST &st) const + { + int ierr = STSetType (st, const_cast(STCAYLEY)); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + ierr = STSetShift (st, additional_data.shift_parameter); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + + ierr = STCayleySetAntishift (st, additional_data.antishift_parameter); + AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr)); + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#else +// On gcc2.95 on Alpha OSF1, the native assembler does not like empty +// files, so provide some dummy code +namespace { void dummy () {} } +#endif // DEAL_II_USE_SLEPC +