set_which_eigenpairs (EPSWhich set_which);
/**
- * Solve the linear system for <code>n_eigenvectors</code>
+ * Solve the linear system for n_eigenvectors
* eigenstates. Parameter <code>n_converged</code> contains the
- * actual number of eigenstates that have converged; this can be
- * both fewer or more than <code>n_eigenvectors</code>, 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);
* 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.
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
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
* Objects for Eigenvalue Problem Solver.
*/
EPS eps;
+
+ /**
+ * Convergence.
+ */
+ EPSConvergedReason reason;
};
/**
* complex.
*
* @ingroup SLEPcWrappers
+ *
* @author Toby D. Young 2008
*/
class SolverKrylovSchur : public SolverBase
* solver. Usage: All spectrum, all problem types, complex.
*
* @ingroup SLEPcWrappers
+ *
* @author Toby D. Young 2008, 2011
*/
class SolverArnoldi : public SolverBase
* and sets the type of solver that is appropriate for this class.
*/
virtual void set_solver_type (EPS &eps) const;
-
};
/**
* solver. Usage: All spectrum, all problem types, complex.
*
* @ingroup SLEPcWrappers
+ *
* @author Toby D. Young 2009
*/
class SolverLanczos : public SolverBase
* and sets the type of solver that is appropriate for this class.
*/
virtual void set_solver_type (EPS &eps) const;
-
};
/**
* solver. Usage: Largest values of spectrum only, all problem
* types, complex.
*
- * @ingroup SLEPcWrappers
+ * @ingroup SLEPcWrappers
+ *
* @author Toby D. Young 2010
*/
class SolverPower : public SolverBase
* 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
{
* 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
{
* and sets the type of solver that is appropriate for this class.
*/
virtual void set_solver_type (EPS &eps) const;
-
};
// --------------------------- inline and template functions -----------
std::vector<OutputVector> &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<n_converged; ++index)
get_eigenpair (index, kr[index], vr[index]);
}
std::vector<OutputVector> &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)
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<n_converged; ++index)
get_eigenpair (index, kr[index], vr[index]);
}
}
solver_control (cn),
mpi_communicator (mpi_communicator),
set_which (EPS_LARGEST_MAGNITUDE),
- opA (NULL), opB (NULL),
+ opA (NULL),
+ opB (NULL),
initial_vector (NULL),
transformation (NULL)
{}
ierr = EPSSetOperators (solver_data->eps, *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)
{
// 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);
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<void *>(&solver_control));
+ // 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,
#ifdef PETSC_USE_64BIT_INDICES
- reinterpret_cast<PetscInt *>(n_converged));
+ reinterpret_cast<PetscInt *>(n_converged)
#else
- reinterpret_cast<int *>(n_converged));
+ reinterpret_cast<int *>(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
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
-
void
SolverBase::reset ()
{
}
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<EPSConvergedReason>(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
}
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<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;
- else
- *reason = EPS_DIVERGED_BREAKDOWN;
- break;
-
- default:
- Assert (false, ExcNotImplemented());
- }
+ // This function is undefined (future reference only).
// return without failure.
return 0;