</li>
</ol>
+<ol>
+ <li>
+ <p>
+ 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.
+ <br>
+ (Toby D. Young 2009/06/25)
+ </p>
+ </li>
+</ol>
<a name="deal.II"></a>
--- /dev/null
+
+#ifndef __deal2__slepc_solver_h
+#define __deal2__slepc_solver_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+#include <lac/slepc_spectral_transformation.h>
+#include <boost/shared_ptr.hpp>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#include <petscksp.h>
+#include <slepceps.h>
+
+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 <typename OutputVector>
+ void
+ SolverBase::solve (const PETScWrappers::MatrixBase &A,
+ const PETScWrappers::MatrixBase &B,
+ std::vector<double> &kr,
+ std::vector<OutputVector> &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 <tt>vr</tt> 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 <typename OutputVector>
+ void
+ solve (const PETScWrappers::MatrixBase &A,
+ const PETScWrappers::MatrixBase &B,
+ std::vector<double> &kr,
+ std::vector<OutputVector> &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<SolverData> 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 <typename OutputVector>
+ void
+ SolverBase::solve (const PETScWrappers::MatrixBase &A,
+ const PETScWrappers::MatrixBase &B,
+ std::vector<double> &kr,
+ std::vector<OutputVector> &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 ---------------------------*/
+
--- /dev/null
+
+#ifndef __deal2__slepc_spectral_transformation_h
+#define __deal2__slepc_spectral_transformation_h
+
+#include <base/config.h>
+#include <lac/exceptions.h>
+#include <lac/solver_control.h>
+#include <lac/slepc_solver.h>
+#include <boost/shared_ptr.hpp>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#include <petscksp.h>
+#include <slepceps.h>
+
+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<TransformationData> 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 ------------------*/
--- /dev/null
+
+#include <lac/petsc_matrix_base.h>
+#include <lac/petsc_vector_base.h>
+#include <lac/petsc_vector.h>
+#include <lac/slepc_solver.h>
+#include <lac/slepc_spectral_transformation.h>
+
+#include <cmath>
+#include <vector>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
+#include <petscsles.h>
+#endif
+#include <petscversion.h>
+
+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<int *>(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<SolverControl*>(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<EPSConvergedReason>(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<char *>(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<char *>(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<char *>(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
+
--- /dev/null
+
+#include <lac/petsc_matrix_base.h>
+#include <lac/petsc_vector_base.h>
+#include <lac/petsc_vector.h>
+#include <lac/slepc_solver.h>
+#include <lac/slepc_spectral_transformation.h>
+
+#include <cmath>
+#include <vector>
+
+#ifdef DEAL_II_USE_SLEPC
+
+#if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR < 2)
+# include <petscsles.h>
+#endif
+#include <petscversion.h>
+
+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<char *>(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<char *>(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<char *>(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<char *>(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
+