* @ingroup SLEPcWrappers
*
* @author Toby D. Young 2008, 2009, 2010, 2011, 2013; and Rickard Armiento
- * 2008.
+ * 2008; and Denis Davydov 2015.
*
* @note Various tweaks and enhancements contributed by Eloy Romero and Jose
* E. Roman 2009, 2010.
*/
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;
-
/**
* 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 <code>n_eigenpairs</code> eigenstates.
* Parameter <code>n_converged</code> contains the actual number of
set_matrices (const PETScWrappers::MatrixBase &A,
const PETScWrappers::MatrixBase &B);
- /**
- * Target eigenvalue to solve for.
- */
- std_cxx11::shared_ptr<PetscScalar> target_eigenvalue;
-
- /**
- * Which portion of the spectrum to solve from.
- */
- EPSWhich set_which;
+ protected:
/**
- * Set the eigenspectrum problem type.
+ * Objects for Eigenvalue Problem Solver.
*/
- EPSProblemType set_problem;
+ EPS eps;
private:
-
- /**
- * The matrix $A$ of the generalized eigenvalue problem $Ax=B\lambda x$,
- * or the standard eigenvalue problem $Ax=\lambda x$.
- */
- const PETScWrappers::MatrixBase *opA;
-
/**
- * The matrix $B$ of the generalized eigenvalue problem $Ax=B\lambda x$.
+ * Convergence.
*/
- const PETScWrappers::MatrixBase *opB;
+ EPSConvergedReason reason;
- /**
- * An initial vector space used in SLEPc solvers.
- */
- std::vector<Vec> initial_space;
-
- /**
- * Pointer to an an object that describes transformations that can be
- * applied to the eigenvalue problem.
- */
- SLEPcWrappers::TransformationBase *transformation;
/**
* A function that can be used in SLEPc as a callback to check on
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 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;
-
- /**
- * Convergence.
- */
- EPSConvergedReason reason;
- };
-
- /**
- * Pointer to the <code>SolverData</code> object.
- */
- std_cxx11::shared_ptr<SolverData> solver_data;
};
/**
* 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;
};
/**
* 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;
};
/**
*
* @ingroup SLEPcWrappers
*
- * @author Toby D. Young 2009; Denis Davydov 2015;
+ * @author Toby D. Young 2009; and Denis Davydov 2015;
*/
class SolverLanczos : public SolverBase
{
* 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;
};
/**
* 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;
};
/**
* 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;
};
/**
* 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;
};
* 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 -----------
void
SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
{
- initial_space.resize(this_initial_space.size());
- for (unsigned int i = 0; i < initial_space.size(); i++)
+ int ierr;
+ std::vector<Vec> vecs(this_initial_space.size());
+
+ for (unsigned int i = 0; i < this_initial_space.size(); i++)
{
- initial_space[i] = this_initial_space[i];
+ vecs[i] = this_initial_space[i];
}
+
+ // if the eigensolver supports only a single initial vector, but several
+ // guesses are provided, then all except the first one will be discarded.
+ // One could still build a vector that is rich in the directions of all guesses,
+ // by taking a linear combination of them. (TODO: make function virtual?)
+
+#if DEAL_II_PETSC_VERSION_LT(3,1,0)
+ ierr = EPSSetInitialVector (eps, &vecs[0]);
+#else
+ ierr = EPSSetInitialSpace (eps, vecs.size(), &vecs[0]);
+#endif
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
}
namespace SLEPcWrappers
{
- SolverBase::SolverData::~SolverData ()
- {
- // Destroy the solver object.
-#if DEAL_II_PETSC_VERSION_LT(3,2,0)
- int ierr = EPSDestroy (eps);
-#else
- int ierr = EPSDestroy (&eps);
-#endif
- 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),
- set_problem (EPS_GNHEP),
- opA (NULL),
- opB (NULL),
- transformation (NULL)
- {}
+ mpi_communicator (mpi_communicator)
+ {
+ // create eigensolver context
+ int ierr = EPSCreate (mpi_communicator, &eps);
+ 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));
+
+ // default values:
+ set_which_eigenpairs(EPS_LARGEST_MAGNITUDE);
+ set_problem_type(EPS_GNHEP);
+
+ // TODO:
+ // By default, EPS initializes the starting vector or the initial subspace randomly.
+ }
SolverBase::~SolverBase ()
- {}
+ {
+ if (eps != NULL)
+ {
+ // Destroy the solver object.
+#if DEAL_II_PETSC_VERSION_LT(3,2,0)
+ int ierr = EPSDestroy (eps);
+#else
+ int ierr = EPSDestroy (&eps);
+#endif
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+ }
+ }
void
SolverBase::set_matrices (const PETScWrappers::MatrixBase &A)
{
// standard eigenspectrum problem
- opA = &A;
- opB = NULL;
+ int ierr = EPSSetOperators (eps, A, PETSC_NULL);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
void
const PETScWrappers::MatrixBase &B)
{
// generalized eigenspectrum problem
- opA = &A;
- opB = &B;
- }
-
- void
- SolverBase::set_initial_vector (const PETScWrappers::VectorBase &this_initial_vector)
- {
- initial_space.resize(1);
- initial_space[0] = this_initial_vector;
+ int ierr = EPSSetOperators (eps, A, B);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
void
- SolverBase::set_transformation (SLEPcWrappers::TransformationBase &this_transformation)
+ SolverBase::set_transformation (SLEPcWrappers::TransformationBase &transformation)
{
- transformation = &this_transformation;
+ // set transformation type if any
+ // STSetShift is called inside
+ int ierr = EPSSetST(eps,transformation.st);
+ AssertThrow (ierr == 0, SolverBase::ExcSLEPcError(ierr));
}
void
SolverBase::set_target_eigenvalue (const PetscScalar &this_target)
{
- target_eigenvalue.reset (new PetscScalar(this_target));
+ // set target eigenvalues to solve for
+ // in all transformation except STSHIFT there is a direct connection between
+ // the target and the shift, read more on p41 of SLEPc manual.
+ int ierr = EPSSetTarget (eps, this_target );
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
void
SolverBase::set_which_eigenpairs (const EPSWhich eps_which)
{
- set_which = eps_which;
+ // set which portion of the eigenspectrum to solve for
+ int ierr = EPSSetWhichEigenpairs (eps, eps_which);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
void
SolverBase::set_problem_type (const EPSProblemType eps_problem)
{
- set_problem = eps_problem;
+ int ierr = EPSSetProblemType (eps, eps_problem);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
void
{
int ierr;
- // create a solver object if this is necessary
- if (solver_data.get() == 0)
- {
- // reset solver dtaa
- solver_data.reset (new SolverData());
-
- // create eigensolver context and 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);
- else
- ierr = EPSSetOperators (solver_data->eps, *opA, PETSC_NULL);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
- // set runtime options
- set_solver_type (solver_data->eps);
- }
-
- // set the problem type
- ierr = EPSSetProblemType (solver_data->eps, set_problem);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
-
- // set the initial vector(s) if any
- if (initial_space.size() != 0)
- {
-
-#if DEAL_II_PETSC_VERSION_LT(3,1,0)
- ierr = EPSSetInitialVector (solver_data->eps, &initial_space[0]);
-#else
- ierr = EPSSetInitialSpace (solver_data->eps, initial_space.size(), &initial_space[0]);
-#endif
-
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
- }
-
- // if a spectral transformation is to be used, set the
- // transformation and target the wanted eigenvalues
- if (transformation)
- {
- // set transformation type if any
- // STSetShift is called inside
- transformation->set_context (solver_data->eps);
- }
-
- // set target eigenvalues to solve for
- // in all transformation except STSHIFT there is a direct connection between
- // the target and the shift, read more on p41 of SLEPc manual.
- if (target_eigenvalue.get() != 0)
- {
- int ierr = EPSSetTarget (solver_data->eps, *(target_eigenvalue.get()) );
- 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));
-
// set number of eigenvectors to compute
- ierr = EPSSetDimensions (solver_data->eps, n_eigenpairs,
+ ierr = EPSSetDimensions (eps, n_eigenpairs,
PETSC_DECIDE, PETSC_DECIDE);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// set the solve options to the eigenvalue problem solver context
- ierr = EPSSetFromOptions (solver_data->eps);
+ ierr = EPSSetFromOptions (eps);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// TODO breaks step-36
- // force solvers to use true residual
+ // force Krylov solver to use true residual instead of an estimate.
//EPSSetTrueResidual(solver_data->eps, PETSC_TRUE);
//AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// Set convergence test to be absolute
- ierr = EPSSetConvergenceTest (solver_data->eps, EPS_CONV_ABS);
+ ierr = EPSSetConvergenceTest (eps, EPS_CONV_ABS);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
- // Set the convergence test function
+ // TODO Set the convergence test function
// ierr = EPSSetConvergenceTestFunction (solver_data->eps, &convergence_test,
// reinterpret_cast<void *>(&solver_control));
// AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// solve the eigensystem
- ierr = EPSSolve (solver_data->eps);
+ ierr = EPSSolve (eps);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// get number of converged eigenstates
- ierr = EPSGetConverged (solver_data->eps,
+ ierr = EPSGetConverged (eps,
reinterpret_cast<PetscInt *>(n_converged));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// complete eigenspectrum
{
// get the number of solver iterations
- ierr = EPSGetIterationNumber (solver_data->eps, &n_iterations);
+ ierr = EPSGetIterationNumber (eps, &n_iterations);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// get the maximum of residual norm among converged eigenvectors.
// used during the solution process.
// Yet, this is the norm which gives error bounds (Saad, 1992, ch3):
// | \lambda - \widehat\lambda | <= ||r||_2
- ierr = EPSComputeResidualNorm (solver_data->eps, i, &residual_norm_i);
+ ierr = EPSComputeResidualNorm (eps, i, &residual_norm_i);
// EPSComputeRelativeError may not be consistent with the stopping criteria
// used during the solution process. Given EPS_CONV_ABS set above,
PetscScalar &eigenvalues,
PETScWrappers::VectorBase &eigenvectors)
{
- AssertThrow (solver_data.get() != 0, ExcSLEPcWrappersUsageError());
-
// get converged eigenpair
- int ierr = EPSGetEigenpair (solver_data->eps, index,
+ int ierr = EPSGetEigenpair (eps, index,
&eigenvalues, PETSC_NULL,
eigenvectors, PETSC_NULL);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
PETScWrappers::VectorBase &imag_eigenvectors)
{
#ifndef PETSC_USE_COMPLEX
- AssertThrow (solver_data.get() != 0, ExcSLEPcWrappersUsageError());
-
// get converged eigenpair
- int ierr = EPSGetEigenpair (solver_data->eps, index,
+ int ierr = EPSGetEigenpair (eps, index,
&real_eigenvalues, &imag_eigenvalues,
real_eigenvectors, imag_eigenvectors);
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
#endif
}
-
void
SolverBase::reset ()
{
- AssertThrow (solver_data.get() != 0, ExcSLEPcWrappersUsageError());
-
- // destroy solver object.
- solver_data.reset ();
+ //TODO
}
- EPS *
- SolverBase::get_eps ()
- {
- 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;
+ reason = EPS_CONVERGED_ITERATING;
break;
case ::dealii::SolverControl::success:
- solver_data->reason = static_cast<EPSConvergedReason>(1);
+ reason = static_cast<EPSConvergedReason>(1);
break;
case ::dealii::SolverControl::failure:
if (solver_control.last_step() > solver_control.max_steps())
- solver_data->reason = EPS_DIVERGED_ITS;
+ reason = EPS_DIVERGED_ITS;
else
- solver_data->reason = EPS_DIVERGED_BREAKDOWN;
+ reason = EPS_DIVERGED_BREAKDOWN;
break;
default:
SolverBase::convergence_test (EPS /*eps */,
PetscScalar /*real_eigenvalue */,
PetscScalar /*imag_eigenvalue */,
- PetscReal /*residual_norm */,
- PetscReal */*estimated_error */,
+ PetscReal /*residual norm associated to the eigenpair */,
+ PetscReal */*(output) computed error estimate */,
void */*solver_control_x*/)
{
+ // If the error estimate returned by the convergence test function is less
+ // than the tolerance, then the eigenvalue is accepted as converged.
// This function is undefined (future reference only).
// return without failure.
:
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());
+ int ierr = EPSSetType (eps, const_cast<char *>(EPSKRYLOVSCHUR));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
:
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());
+ int ierr = EPSSetType (eps, const_cast<char *>(EPSARNOLDI));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
// if requested, set delayed reorthogonalization in the Arnoldi
}
}
+
/* ---------------------- Lanczos ------------------------ */
SolverLanczos::AdditionalData::
AdditionalData(const EPSLanczosReorthogType r)
:
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());
+ int ierr = EPSSetType (eps, const_cast<char *>(EPSLANCZOS));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
ierr = EPSLanczosSetReorthog(eps,additional_data.reorthog);
:
SolverBase (cn, mpi_communicator),
additional_data (data)
- {}
-
- void
- SolverPower::set_solver_type (EPS &eps) const
{
- int ierr;
- ierr = EPSSetType (eps, const_cast<char *>(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());
+ int ierr = EPSSetType (eps, const_cast<char *>(EPSPOWER));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
:
SolverBase (cn, mpi_communicator),
additional_data (data)
- {}
-
- void
- SolverGeneralizedDavidson::set_solver_type (EPS &eps) const
{
#if DEAL_II_PETSC_VERSION_GTE(3,1,0)
- int ierr;
- ierr = EPSSetType (eps, const_cast<char *>(EPSGD));
- 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());
+ int ierr = EPSSetType (eps, const_cast<char *>(EPSGD));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
if (additional_data.double_expansion)
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
#else
- // Suppress compiler warnings about unused parameters.
- (void) eps;
-
// PETSc/SLEPc version must be > 3.1.0.
Assert ((false),
ExcMessage ("Your SLEPc installation does not include a copy of the "
:
SolverBase (cn, mpi_communicator),
additional_data (data)
- {}
-
- void
- SolverJacobiDavidson::set_solver_type (EPS &eps) const
{
#if DEAL_II_PETSC_VERSION_GTE(3,1,0)
int ierr;
ierr = EPSSetType (eps, const_cast<char *>(EPSJD));
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));
#else
- // Suppress compiler warnings about unused parameters.
- (void) eps;
-
// PETSc/SLEPc version must be > 3.1.0.
Assert ((false),
ExcMessage ("Your SLEPc installation does not include a copy of the "
#endif
}
-
/* ---------------------- LAPACK ------------------------- */
SolverLAPACK::SolverLAPACK (SolverControl &cn,
const MPI_Comm &mpi_communicator,
:
SolverBase (cn, mpi_communicator),
additional_data (data)
- {}
-
- void
- SolverLAPACK::set_solver_type (EPS &eps) const
{
// 'Tis overwhelmingly likely that PETSc/SLEPc *always* has
// BLAS/LAPACK, but let's be defensive.
int ierr;
ierr = EPSSetType (eps, const_cast<char *>(EPSLAPACK));
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));
#else
- // Suppress compiler warnings about unused parameters.
- (void) eps;
-
Assert ((false),
ExcMessage ("Your PETSc/SLEPc installation was not configured with BLAS/LAPACK "
"but this is needed to use the LAPACK solver."));
#endif
}
-
}
DEAL_II_NAMESPACE_CLOSE
projects / dealii.git / commitdiff