#include <slepc_eigen_solver.h>

Inheritance diagram for libMesh::SlepcEigenSolver< T >:

Public Member Functions
	SlepcEigenSolver (const Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
	~SlepcEigenSolver ()
void	clear ()
void	init ()
std::pair< unsigned int, unsigned int >	solve_standard (SparseMatrix< T > &matrix_A, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	solve_standard (ShellMatrix< T > &shell_matrix, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	solve_generalized (SparseMatrix< T > &matrix_A, SparseMatrix< T > &matrix_B, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	solve_generalized (ShellMatrix< T > &matrix_A, SparseMatrix< T > &matrix_B, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	solve_generalized (SparseMatrix< T > &matrix_A, ShellMatrix< T > &matrix_B, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	solve_generalized (ShellMatrix< T > &matrix_A, ShellMatrix< T > &matrix_B, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< Real, Real >	get_eigenpair (unsigned int i, NumericVector< T > &solution_in)
std::pair< Real, Real >	get_eigenvalue (unsigned int i)
Real	get_relative_error (unsigned int i)
void	attach_deflation_space (NumericVector< T > &deflation_vector)
bool	initialized () const
EigenSolverType	eigen_solver_type () const
EigenProblemType	eigen_problem_type () const
PositionOfSpectrum	position_of_spectrum () const
void	set_eigensolver_type (const EigenSolverType est)
void	set_eigenproblem_type (EigenProblemType ept)
void	set_position_of_spectrum (PositionOfSpectrum pos)
const Parallel::Communicator &	comm () const
processor_id_type	n_processors () const
processor_id_type	processor_id () const
Static Public Member Functions
static UniquePtr< EigenSolver < T > >	build (const Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD, const SolverPackage solver_package=SLEPC_SOLVERS)
static std::string	get_info ()
static void	print_info (std::ostream &out=libMesh::out)
static unsigned int	n_objects ()
static void	enable_print_counter_info ()
static void	disable_print_counter_info ()
Protected Types
typedef std::map< std::string, std::pair< unsigned int, unsigned int > >	Counts
Protected Member Functions
void	increment_constructor_count (const std::string &name)
void	increment_destructor_count (const std::string &name)
Protected Attributes
EigenSolverType	_eigen_solver_type
EigenProblemType	_eigen_problem_type
PositionOfSpectrum	_position_of_spectrum
bool	_is_initialized
const Parallel::Communicator &	_communicator
Static Protected Attributes
static Counts	_counts
static Threads::atomic < unsigned int >	_n_objects
static Threads::spin_mutex	_mutex
static bool	_enable_print_counter = true
Private Member Functions
std::pair< unsigned int, unsigned int >	_solve_standard_helper (Mat mat, int nev, int ncv, const double tol, const unsigned int m_its)
std::pair< unsigned int, unsigned int >	_solve_generalized_helper (Mat mat_A, Mat mat_B, int nev, int ncv, const double tol, const unsigned int m_its)
void	set_slepc_solver_type ()
void	set_slepc_problem_type ()
void	set_slepc_position_of_spectrum ()
Static Private Member Functions
static PetscErrorCode	_petsc_shell_matrix_mult (Mat mat, Vec arg, Vec dest)
static PetscErrorCode	_petsc_shell_matrix_get_diagonal (Mat mat, Vec dest)
Private Attributes
EPS	_eps

Detailed Description

template<typename T>
class libMesh::SlepcEigenSolver< T >

This class provides an interface to the SLEPc eigenvalue solver library www.grycap.upv.es/slepc/.

Definition at line 51 of file slepc_eigen_solver.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts [protected, inherited]

Data structure to log the information. The log is identified by the class name.

Definition at line 113 of file reference_counter.h.

Constructor & Destructor Documentation

template<typename T >

libMesh::SlepcEigenSolver< T >::SlepcEigenSolver ( const Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD ) [inline]

Constructor. Initializes Petsc data structures

Definition at line 263 of file slepc_eigen_solver.h.

References libMesh::EigenSolver< T >::_eigen_problem_type, libMesh::EigenSolver< T >::_eigen_solver_type, libMesh::ARNOLDI, and libMesh::NHEP.

                                                                          :
  EigenSolver<T>(comm_in)
{
  this->_eigen_solver_type  = ARNOLDI;
  this->_eigen_problem_type = NHEP;
}

template<typename T >

libMesh::SlepcEigenSolver< T >::~SlepcEigenSolver ( ) [inline]

Destructor.

Definition at line 274 of file slepc_eigen_solver.h.

{
  this->clear ();
}

Member Function Documentation

template<typename T >

PetscErrorCode libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal	(	Mat	mat,
		Vec	dest
	)		`[static, private]`

Internal function if shell matrix mode is used, this just calls the shell matrix's get_diagonal function. Required in order to use Jacobi preconditioning.

Definition at line 781 of file slepc_eigen_solver.C.

References libMesh::ParallelObject::comm(), libMesh::ShellMatrix< T >::get_diagonal(), and libMesh::ierr.

{
  /* Get the matrix context.  */
  PetscErrorCode ierr=0;
  void* ctx;
  ierr = MatShellGetContext(mat,&ctx);

  Parallel::communicator comm;
  PetscObjectGetComm((PetscObject)mat,&comm);
  CHKERRABORT(comm,ierr);

  /* Get user shell matrix object.  */
  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);

  /* Make \p NumericVector instances around the vector.  */
  PetscVector<T> dest_global(dest, shell_matrix.comm());

  /* Call the user function.  */
  shell_matrix.get_diagonal(dest_global);

  return ierr;
}

template<typename T >

PetscErrorCode libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult	(	Mat	mat,
		Vec	arg,
		Vec	dest
	)		`[static, private]`

Internal function if shell matrix mode is used, this just calls the shell matrix's matrix multiplication function. See PetscLinearSolver for a similar implementation.

Definition at line 756 of file slepc_eigen_solver.C.

References libMesh::ParallelObject::comm(), libMesh::ierr, and libMesh::ShellMatrix< T >::vector_mult().

{
  /* Get the matrix context.  */
  PetscErrorCode ierr=0;
  void* ctx;
  ierr = MatShellGetContext(mat,&ctx);

  Parallel::communicator comm;
  PetscObjectGetComm((PetscObject)mat,&comm);
  CHKERRABORT(comm,ierr);

  /* Get user shell matrix object.  */
  const ShellMatrix<T>& shell_matrix = *static_cast<const ShellMatrix<T>*>(ctx);

  /* Make \p NumericVector instances around the vectors.  */
  PetscVector<T> arg_global(arg,   shell_matrix.comm());
  PetscVector<T> dest_global(dest, shell_matrix.comm());

  /* Call the user function.  */
  shell_matrix.vector_mult(dest_global,arg_global);

  return ierr;
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::_solve_generalized_helper	(	Mat	mat_A,
		Mat	mat_B,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[private]`

Helper function that actually performs the generalized eigensolve.

Definition at line 448 of file slepc_eigen_solver.C.

References libMesh::ierr, and libMesh::START_LOG().

{
  START_LOG("solve_generalized()", "SlepcEigenSolver");

  PetscErrorCode ierr=0;

  // converged eigen pairs and number of iterations
  PetscInt nconv=0;
  PetscInt its=0;

#ifdef  DEBUG
  // The relative error.
  PetscReal error, re, im;

  // Pointer to vectors of the real parts, imaginary parts.
  PetscScalar kr, ki;
#endif

  // Set operators.
  ierr = EPSSetOperators (_eps, mat_A, mat_B);
  LIBMESH_CHKERRABORT(ierr);

  //set the problem type and the position of the spectrum
  set_slepc_problem_type();
  set_slepc_position_of_spectrum();

  // Set eigenvalues to be computed.
#if SLEPC_VERSION_LESS_THAN(3,0,0)
  ierr = EPSSetDimensions (_eps, nev, ncv);
#else
  ierr = EPSSetDimensions (_eps, nev, ncv, PETSC_DECIDE);
#endif
  LIBMESH_CHKERRABORT(ierr);


  // Set the tolerance and maximum iterations.
  ierr = EPSSetTolerances (_eps, tol, m_its);
  LIBMESH_CHKERRABORT(ierr);

  // Set runtime options, e.g.,
  //      -eps_type <type>, -eps_nev <nev>, -eps_ncv <ncv>
  // Similar to PETSc, these options will override those specified
  // above as long as EPSSetFromOptions() is called _after_ any
  // other customization routines.
  ierr = EPSSetFromOptions (_eps);
  LIBMESH_CHKERRABORT(ierr);

  // Solve the eigenproblem.
  ierr = EPSSolve (_eps);
  LIBMESH_CHKERRABORT(ierr);

  // Get the number of iterations.
  ierr = EPSGetIterationNumber (_eps, &its);
  LIBMESH_CHKERRABORT(ierr);

  // Get number of converged eigenpairs.
  ierr = EPSGetConverged(_eps,&nconv);
  LIBMESH_CHKERRABORT(ierr);


#ifdef DEBUG
  // ierr = PetscPrintf(this->comm().get(),
  //         "\n Number of iterations: %d\n"
  //         " Number of converged eigenpairs: %d\n\n", its, nconv);

  // Display eigenvalues and relative errors.
  ierr = PetscPrintf(this->comm().get(),
                     "           k           ||Ax-kx||/|kx|\n"
                     "   ----------------- -----------------\n" );
  LIBMESH_CHKERRABORT(ierr);

  for(PetscInt i=0; i<nconv; i++ )
    {
      ierr = EPSGetEigenpair(_eps, i, &kr, &ki, PETSC_NULL, PETSC_NULL);
      LIBMESH_CHKERRABORT(ierr);

      ierr = EPSComputeRelativeError(_eps, i, &error);
      LIBMESH_CHKERRABORT(ierr);

#ifdef LIBMESH_USE_COMPLEX_NUMBERS
      re = PetscRealPart(kr);
      im = PetscImaginaryPart(kr);
#else
      re = kr;
      im = ki;
#endif

      if (im != .0)
        {
          ierr = PetscPrintf(this->comm().get()," %9f%+9f i %12f\n", re, im, error);
          LIBMESH_CHKERRABORT(ierr);
        }
      else
        {
          ierr = PetscPrintf(this->comm().get(),"   %12f       %12f\n", re, error);
          LIBMESH_CHKERRABORT(ierr);
        }
    }

  ierr = PetscPrintf(this->comm().get(),"\n" );
  LIBMESH_CHKERRABORT(ierr);
#endif // DEBUG

  STOP_LOG("solve_generalized()", "SlepcEigenSolver");

  // return the number of converged eigenpairs
  // and the number of iterations
  return std::make_pair(nconv, its);

}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::_solve_standard_helper	(	Mat	mat,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[private]`

Helper function that actually performs the standard eigensolve.

Definition at line 167 of file slepc_eigen_solver.C.

References libMesh::ierr, and libMesh::START_LOG().

{
  START_LOG("solve_standard()", "SlepcEigenSolver");

  PetscErrorCode ierr=0;

  // converged eigen pairs and number of iterations
  PetscInt nconv=0;
  PetscInt its=0;

#ifdef  DEBUG
  // The relative error.
  PetscReal error, re, im;

  // Pointer to vectors of the real parts, imaginary parts.
  PetscScalar kr, ki;
#endif

  // Set operators.
  ierr = EPSSetOperators (_eps, mat, PETSC_NULL);
  LIBMESH_CHKERRABORT(ierr);

  //set the problem type and the position of the spectrum
  set_slepc_problem_type();
  set_slepc_position_of_spectrum();

  // Set eigenvalues to be computed.
#if SLEPC_VERSION_LESS_THAN(3,0,0)
  ierr = EPSSetDimensions (_eps, nev, ncv);
#else
  ierr = EPSSetDimensions (_eps, nev, ncv, PETSC_DECIDE);
#endif
  LIBMESH_CHKERRABORT(ierr);
  // Set the tolerance and maximum iterations.
  ierr = EPSSetTolerances (_eps, tol, m_its);
  LIBMESH_CHKERRABORT(ierr);

  // Set runtime options, e.g.,
  //      -eps_type <type>, -eps_nev <nev>, -eps_ncv <ncv>
  // Similar to PETSc, these options will override those specified
  // above as long as EPSSetFromOptions() is called _after_ any
  // other customization routines.
  ierr = EPSSetFromOptions (_eps);
  LIBMESH_CHKERRABORT(ierr);

  // Solve the eigenproblem.
  ierr = EPSSolve (_eps);
  LIBMESH_CHKERRABORT(ierr);

  // Get the number of iterations.
  ierr = EPSGetIterationNumber (_eps, &its);
  LIBMESH_CHKERRABORT(ierr);

  // Get number of converged eigenpairs.
  ierr = EPSGetConverged(_eps,&nconv);
  LIBMESH_CHKERRABORT(ierr);


#ifdef DEBUG
  // ierr = PetscPrintf(this->comm().get(),
  //         "\n Number of iterations: %d\n"
  //         " Number of converged eigenpairs: %d\n\n", its, nconv);

  // Display eigenvalues and relative errors.
  ierr = PetscPrintf(this->comm().get(),
                     "           k           ||Ax-kx||/|kx|\n"
                     "   ----------------- -----------------\n" );
  LIBMESH_CHKERRABORT(ierr);

  for(PetscInt i=0; i<nconv; i++ )
    {
      ierr = EPSGetEigenpair(_eps, i, &kr, &ki, PETSC_NULL, PETSC_NULL);
      LIBMESH_CHKERRABORT(ierr);

      ierr = EPSComputeRelativeError(_eps, i, &error);
      LIBMESH_CHKERRABORT(ierr);

#ifdef LIBMESH_USE_COMPLEX_NUMBERS
      re = PetscRealPart(kr);
      im = PetscImaginaryPart(kr);
#else
      re = kr;
      im = ki;
#endif

      if (im != .0)
        {
          ierr = PetscPrintf(this->comm().get()," %9f%+9f i %12f\n", re, im, error);
          LIBMESH_CHKERRABORT(ierr);
        }
      else
        {
          ierr = PetscPrintf(this->comm().get(),"   %12f       %12f\n", re, error);
          LIBMESH_CHKERRABORT(ierr);
        }
    }

  ierr = PetscPrintf(this->comm().get(),"\n" );
  LIBMESH_CHKERRABORT(ierr);
#endif // DEBUG


  STOP_LOG("solve_standard()", "SlepcEigenSolver");

  // return the number of converged eigenpairs
  // and the number of iterations
  return std::make_pair(nconv, its);

}

template<typename T >

void libMesh::SlepcEigenSolver< T >::attach_deflation_space ( NumericVector< T > & deflation_vector ) [virtual]

Attach a deflation space defined by a single vector.

Implements libMesh::EigenSolver< T >.

Definition at line 740 of file slepc_eigen_solver.C.

References libMesh::ierr, and libMesh::TriangleWrapper::init().

{
  this->init();

  PetscErrorCode ierr = 0;
  Vec deflation_vector = (cast_ptr<PetscVector<T>*>(&deflation_vector_in))->vec();
  Vec* deflation_space = &deflation_vector;
#if SLEPC_VERSION_LESS_THAN(3,1,0)
  ierr = EPSAttachDeflationSpace(_eps, 1, deflation_space, PETSC_FALSE);
#else
  ierr = EPSSetDeflationSpace(_eps, 1, deflation_space);
#endif
  LIBMESH_CHKERRABORT(ierr);
}

template<typename T >

UniquePtr< EigenSolver< T > > libMesh::EigenSolver< T >::build	(	const Parallel::Communicator &comm_in	LIBMESH_CAN_DEFAULT_TO_COMMWORLD,
		const SolverPackage	solver_package = `SLEPC_SOLVERS`
	)		`[static, inherited]`

Builds an EigenSolver using the linear solver package specified by solver_package

Definition at line 36 of file eigen_solver.C.

References libMesh::SLEPC_SOLVERS.

{
  // Build the appropriate solver
  switch (solver_package)
    {

#ifdef LIBMESH_HAVE_SLEPC
    case SLEPC_SOLVERS:
      return UniquePtr<EigenSolver<T> >(new SlepcEigenSolver<T>(comm));
#endif

    default:
      libmesh_error_msg("ERROR:  Unrecognized eigen solver package: " << solver_package);
    }

  return UniquePtr<EigenSolver<T> >();
}

template<typename T >

void libMesh::SlepcEigenSolver< T >::clear ( ) [virtual]

Release all memory and clear data structures.

Reimplemented from libMesh::EigenSolver< T >.

Definition at line 52 of file slepc_eigen_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, libMesh::ARNOLDI, libMesh::ierr, libMesh::initialized(), and libMesh::KRYLOVSCHUR.

{
  if (this->initialized())
    {
      this->_is_initialized = false;

      PetscErrorCode ierr=0;

      ierr = LibMeshEPSDestroy(&_eps);
      LIBMESH_CHKERRABORT(ierr);

      // SLEPc default eigenproblem solver
#if SLEPC_VERSION_LESS_THAN(2,3,2)
      this->_eigen_solver_type = ARNOLDI;
#else
      // Krylov-Schur showed up as of Slepc 2.3.2
      this->_eigen_solver_type = KRYLOVSCHUR;
#endif
    }
}

const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const [inline, inherited]

Returns:: a reference to the Parallel::Communicator object used by this mesh.

Definition at line 86 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

  { return _communicator; }

void libMesh::ReferenceCounter::disable_print_counter_info ( ) [static, inherited]

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::LibMeshInit().

{
  _enable_print_counter = false;
  return;
}

template<typename T >

EigenProblemType libMesh::EigenSolver< T >::eigen_problem_type ( ) const [inline, inherited]

Returns the type of the eigen problem.

Definition at line 102 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_eigen_problem_type.

{ return _eigen_problem_type;}

template<typename T >

EigenSolverType libMesh::EigenSolver< T >::eigen_solver_type ( ) const [inline, inherited]

Returns the type of eigensolver to use.

Definition at line 97 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_eigen_solver_type.

{ return _eigen_solver_type; }

void libMesh::ReferenceCounter::enable_print_counter_info ( ) [static, inherited]

Methods to enable/disable the reference counter output from print_info()

Definition at line 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

{
  _enable_print_counter = true;
  return;
}

template<typename T >

std::pair< Real, Real > libMesh::SlepcEigenSolver< T >::get_eigenpair	(	unsigned int	i,
		NumericVector< T > &	solution_in
	)		`[virtual]`

This function returns the real and imaginary part of the ith eigenvalue and copies the respective eigenvector to the solution vector. Note that also in case of purely real matrix entries the eigenpair may be complex values.

Implements libMesh::EigenSolver< T >.

Definition at line 671 of file slepc_eigen_solver.C.

References libMesh::PetscVector< T >::close(), libMesh::ierr, and libMesh::PetscVector< T >::vec().

{
  PetscErrorCode ierr=0;

  PetscReal re, im;

  // Make sure the NumericVector passed in is really a PetscVector
  PetscVector<T>* solution = cast_ptr<PetscVector<T>*>(&solution_in);

  // real and imaginary part of the ith eigenvalue.
  PetscScalar kr, ki;

  solution->close();

  ierr = EPSGetEigenpair(_eps, i, &kr, &ki, solution->vec(), PETSC_NULL);
  LIBMESH_CHKERRABORT(ierr);

#ifdef LIBMESH_USE_COMPLEX_NUMBERS
  re = PetscRealPart(kr);
  im = PetscImaginaryPart(kr);
#else
  re = kr;
  im = ki;
#endif

  return std::make_pair(re, im);
}

template<typename T >

std::pair< Real, Real > libMesh::SlepcEigenSolver< T >::get_eigenvalue ( unsigned int i ) [virtual]

Same as above, but does not copy the eigenvector.

Implements libMesh::EigenSolver< T >.

Definition at line 702 of file slepc_eigen_solver.C.

References libMesh::ierr.

{
  PetscErrorCode ierr=0;

  PetscReal re, im;

  // real and imaginary part of the ith eigenvalue.
  PetscScalar kr, ki;

  ierr = EPSGetEigenvalue(_eps, i, &kr, &ki);
  LIBMESH_CHKERRABORT(ierr);

#ifdef LIBMESH_USE_COMPLEX_NUMBERS
  re = PetscRealPart(kr);
  im = PetscImaginaryPart(kr);
#else
  re = kr;
  im = ki;
#endif

  return std::make_pair(re, im);
}

std::string libMesh::ReferenceCounter::get_info ( ) [static, inherited]

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

{
#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)

  std::ostringstream oss;

  oss << '\n'
      << " ---------------------------------------------------------------------------- \n"
      << "| Reference count information                                                |\n"
      << " ---------------------------------------------------------------------------- \n";

  for (Counts::iterator it = _counts.begin();
       it != _counts.end(); ++it)
    {
      const std::string name(it->first);
      const unsigned int creations    = it->second.first;
      const unsigned int destructions = it->second.second;

      oss << "| " << name << " reference count information:\n"
          << "|  Creations:    " << creations    << '\n'
          << "|  Destructions: " << destructions << '\n';
    }

  oss << " ---------------------------------------------------------------------------- \n";

  return oss.str();

#else

  return "";

#endif
}

template<typename T >

Real libMesh::SlepcEigenSolver< T >::get_relative_error ( unsigned int i )

Returns:: the relative error ||A*x-lambda*x||/|lambda*x| of the ith eigenpair. (or the equivalent for a general eigenvalue problem)

Definition at line 727 of file slepc_eigen_solver.C.

References libMesh::ierr.

{
  PetscErrorCode ierr=0;
  PetscReal error;

  ierr = EPSComputeRelativeError(_eps, i, &error);
  LIBMESH_CHKERRABORT(ierr);

  return error;
}

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string & name ) [inline, protected, inherited]

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 163 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.first++;
}

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string & name ) [inline, protected, inherited]

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 176 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

{
  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
  std::pair<unsigned int, unsigned int>& p = _counts[name];

  p.second++;
}

template<typename T >

void libMesh::SlepcEigenSolver< T >::init ( ) [virtual]

Initialize data structures if not done so already.

Implements libMesh::EigenSolver< T >.

Definition at line 76 of file slepc_eigen_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, libMesh::ierr, and libMesh::initialized().

{

  PetscErrorCode ierr=0;

  // Initialize the data structures if not done so already.
  if (!this->initialized())
    {
      this->_is_initialized = true;

      // Create the eigenproblem solver context
      ierr = EPSCreate (this->comm().get(), &_eps);
      LIBMESH_CHKERRABORT(ierr);

      // Set user-specified  solver
      set_slepc_solver_type();
    }
}

template<typename T >

bool libMesh::EigenSolver< T >::initialized ( ) const [inline, inherited]

Returns:: true if the data structures are initialized, false otherwise.

Definition at line 81 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_is_initialized.

{ return _is_initialized; }

static unsigned int libMesh::ReferenceCounter::n_objects ( ) [inline, static, inherited]

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 79 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

Referenced by libMesh::LibMeshInit::~LibMeshInit().

  { return _n_objects; }

processor_id_type libMesh::ParallelObject::n_processors ( ) const [inline, inherited]

Returns:: the number of processors in the group.

Definition at line 92 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

  { return cast_int<processor_id_type>(_communicator.size()); }

template<typename T >

PositionOfSpectrum libMesh::EigenSolver< T >::position_of_spectrum ( ) const [inline, inherited]

Returns the position of the spectrum to compute.

Definition at line 107 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_position_of_spectrum.

  { return _position_of_spectrum;}

void libMesh::ReferenceCounter::print_info ( std::ostream & out = libMesh::out ) [static, inherited]

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::~LibMeshInit().

{
  if( _enable_print_counter ) out_stream << ReferenceCounter::get_info();
}

processor_id_type libMesh::ParallelObject::processor_id ( ) const [inline, inherited]

Returns:: the rank of this processor in the group.

Definition at line 98 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

  { return cast_int<processor_id_type>(_communicator.rank()); }

template<typename T >

void libMesh::EigenSolver< T >::set_eigenproblem_type ( EigenProblemType ept ) [inline, inherited]

Sets the type of the eigenproblem.

Definition at line 119 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_eigen_problem_type.

  {_eigen_problem_type = ept;}

template<typename T >

void libMesh::EigenSolver< T >::set_eigensolver_type ( const EigenSolverType est ) [inline, inherited]

Sets the type of eigensolver to use.

Definition at line 113 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_eigen_solver_type.

  { _eigen_solver_type = est; }

template<typename T >

void libMesh::EigenSolver< T >::set_position_of_spectrum ( PositionOfSpectrum pos ) [inline, inherited]

Sets the position of the spectrum.

Definition at line 125 of file eigen_solver.h.

References libMesh::EigenSolver< T >::_position_of_spectrum.

  {_position_of_spectrum= pos;}

template<typename T >

void libMesh::SlepcEigenSolver< T >::set_slepc_position_of_spectrum ( ) [private]

Tells Slepc to compute the spectrum at the position stored in _position_of_spectrum

Definition at line 640 of file slepc_eigen_solver.C.

References libMesh::ierr, libMesh::LARGEST_IMAGINARY, libMesh::LARGEST_MAGNITUDE, libMesh::LARGEST_REAL, libMesh::SMALLEST_IMAGINARY, libMesh::SMALLEST_MAGNITUDE, and libMesh::SMALLEST_REAL.

{
  PetscErrorCode ierr = 0;

  switch (this->_position_of_spectrum)
    {
    case LARGEST_MAGNITUDE:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_MAGNITUDE);  LIBMESH_CHKERRABORT(ierr); return;
    case SMALLEST_MAGNITUDE:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_MAGNITUDE); LIBMESH_CHKERRABORT(ierr); return;
    case LARGEST_REAL:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_REAL);       LIBMESH_CHKERRABORT(ierr); return;
    case SMALLEST_REAL:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_REAL);      LIBMESH_CHKERRABORT(ierr); return;
    case LARGEST_IMAGINARY:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_LARGEST_IMAGINARY);  LIBMESH_CHKERRABORT(ierr); return;
    case SMALLEST_IMAGINARY:
      ierr = EPSSetWhichEigenpairs (_eps, EPS_SMALLEST_IMAGINARY); LIBMESH_CHKERRABORT(ierr); return;


    default:
      libmesh_error_msg("ERROR:  Unsupported SLEPc position of spectrum: " << this->_position_of_spectrum);
    }
}

template<typename T >

void libMesh::SlepcEigenSolver< T >::set_slepc_problem_type ( ) [private]

Tells Slepc to deal with the type of problem stored in _eigen_problem_type

Definition at line 610 of file slepc_eigen_solver.C.

References libMesh::err, libMesh::GHEP, libMesh::GHIEP, libMesh::GNHEP, libMesh::HEP, libMesh::ierr, and libMesh::NHEP.

{
  PetscErrorCode ierr = 0;

  switch (this->_eigen_problem_type)
    {
    case NHEP:
      ierr = EPSSetProblemType (_eps, EPS_NHEP);  LIBMESH_CHKERRABORT(ierr); return;
    case GNHEP:
      ierr = EPSSetProblemType (_eps, EPS_GNHEP); LIBMESH_CHKERRABORT(ierr); return;
    case HEP:
      ierr = EPSSetProblemType (_eps, EPS_HEP);   LIBMESH_CHKERRABORT(ierr); return;
    case GHEP:
      ierr = EPSSetProblemType (_eps, EPS_GHEP);  LIBMESH_CHKERRABORT(ierr); return;
#if !SLEPC_VERSION_LESS_THAN(3,3,0)
      // EPS_GHIEP added in 3.3.0
    case GHIEP:
      ierr = EPSSetProblemType (_eps, EPS_GHIEP);  LIBMESH_CHKERRABORT(ierr); return;
#endif

    default:
      libMesh::err << "ERROR:  Unsupported SLEPc Eigen Problem: "
                   << this->_eigen_problem_type        << std::endl
                   << "Continuing with SLEPc defaults" << std::endl;
    }
}

template<typename T >

void libMesh::SlepcEigenSolver< T >::set_slepc_solver_type ( ) [private]

Tells Slepc to use the user-specified solver stored in _eigen_solver_type

Definition at line 575 of file slepc_eigen_solver.C.

References libMesh::ARNOLDI, libMesh::Utility::enum_to_string(), libMesh::err, libMesh::ierr, libMesh::KRYLOVSCHUR, libMesh::LANCZOS, libMesh::LAPACK, libMesh::POWER, and libMesh::SUBSPACE.

{
  PetscErrorCode ierr = 0;

  switch (this->_eigen_solver_type)
    {
    case POWER:
      ierr = EPSSetType (_eps, (char*) EPSPOWER);    LIBMESH_CHKERRABORT(ierr); return;
    case SUBSPACE:
      ierr = EPSSetType (_eps, (char*) EPSSUBSPACE); LIBMESH_CHKERRABORT(ierr); return;
    case LAPACK:
      ierr = EPSSetType (_eps, (char*) EPSLAPACK);   LIBMESH_CHKERRABORT(ierr); return;
    case ARNOLDI:
      ierr = EPSSetType (_eps, (char*) EPSARNOLDI);  LIBMESH_CHKERRABORT(ierr); return;
    case LANCZOS:
      ierr = EPSSetType (_eps, (char*) EPSLANCZOS);  LIBMESH_CHKERRABORT(ierr); return;
#if !SLEPC_VERSION_LESS_THAN(2,3,2)
      // EPSKRYLOVSCHUR added in 2.3.2
    case KRYLOVSCHUR:
      ierr = EPSSetType (_eps, (char*) EPSKRYLOVSCHUR);  LIBMESH_CHKERRABORT(ierr); return;
#endif
      // case ARPACK:
      // ierr = EPSSetType (_eps, (char*) EPSARPACK);   LIBMESH_CHKERRABORT(ierr); return;

    default:
      libMesh::err << "ERROR:  Unsupported SLEPc Eigen Solver: "
                   << Utility::enum_to_string(this->_eigen_solver_type) << std::endl
                   << "Continuing with SLEPc defaults" << std::endl;
    }
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_generalized	(	SparseMatrix< T > &	matrix_A,
		SparseMatrix< T > &	matrix_B,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

This function calls the SLEPc solver to compute the eigenpairs for the generalized eigenproblem defined by the matrix_A and matrix_B, which are of type SparseMatrix. The argument nev is the number of eigenpairs to be computed and ncv is the number of basis vectors to be used in the solution procedure. Return values are the number of converged eigen values and the number of the iterations carried out by the eigen solver.

Implements libMesh::EigenSolver< T >.

Definition at line 287 of file slepc_eigen_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::TriangleWrapper::init(), and libMesh::PetscMatrix< T >::mat().

{
  this->init ();

  // Make sure the data passed in are really of Petsc types
  PetscMatrix<T>* matrix_A   = cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
  PetscMatrix<T>* matrix_B   = cast_ptr<PetscMatrix<T>*>(&matrix_B_in);

  // Close the matrix and vectors in case this wasn't already done.
  matrix_A->close ();
  matrix_B->close ();


  return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), nev, ncv, tol, m_its);
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_generalized	(	ShellMatrix< T > &	matrix_A,
		SparseMatrix< T > &	matrix_B,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

Solve generalized eigenproblem when matrix_A is of type ShellMatrix, matrix_B is of type SparseMatrix.

Implements libMesh::EigenSolver< T >.

Definition at line 310 of file slepc_eigen_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::ShellMatrix< T >::m(), libMesh::PetscMatrix< T >::mat(), and libMesh::ShellMatrix< T >::n().

{
  this->init ();

  PetscErrorCode ierr=0;

  // Prepare the matrix.
  Mat mat_A;
  ierr = MatCreateShell(this->comm().get(),
                        shell_matrix_A.m(), // Specify the number of local rows
                        shell_matrix_A.n(), // Specify the number of local columns
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        const_cast<void*>(static_cast<const void*>(&shell_matrix_A)),
                        &mat_A);
  LIBMESH_CHKERRABORT(ierr);

  PetscMatrix<T>* matrix_B   = cast_ptr<PetscMatrix<T>*>(&matrix_B_in);

  // Close the matrix and vectors in case this wasn't already done.
  matrix_B->close ();

  /* Note that the const_cast above is only necessary because PETSc
     does not accept a const void*.  Inside the member function
     _petsc_shell_matrix() below, the pointer is casted back to a
     const ShellMatrix<T>*.  */

  ierr = MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
  LIBMESH_CHKERRABORT(ierr);
  ierr = MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
  LIBMESH_CHKERRABORT(ierr);

  return _solve_generalized_helper (mat_A, matrix_B->mat(), nev, ncv, tol, m_its);
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_generalized	(	SparseMatrix< T > &	matrix_A,
		ShellMatrix< T > &	matrix_B,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

Solve generalized eigenproblem when matrix_A is of type SparseMatrix, matrix_B is of type ShellMatrix. When using this function, one should use the command line options: -st_ksp_type gmres -st_pc_type none or -st_ksp_type gmres -st_pc_type jacobi or similar.

Implements libMesh::EigenSolver< T >.

Definition at line 352 of file slepc_eigen_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::ShellMatrix< T >::m(), libMesh::PetscMatrix< T >::mat(), and libMesh::ShellMatrix< T >::n().

{
  this->init ();

  PetscErrorCode ierr=0;

  PetscMatrix<T>* matrix_A = cast_ptr<PetscMatrix<T>*>(&matrix_A_in);

  // Close the matrix and vectors in case this wasn't already done.
  matrix_A->close ();

  // Prepare the matrix.
  Mat mat_B;
  ierr = MatCreateShell(this->comm().get(),
                        shell_matrix_B.m(), // Specify the number of local rows
                        shell_matrix_B.n(), // Specify the number of local columns
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        const_cast<void*>(static_cast<const void*>(&shell_matrix_B)),
                        &mat_B);
  LIBMESH_CHKERRABORT(ierr);

  /* Note that the const_cast above is only necessary because PETSc
     does not accept a const void*.  Inside the member function
     _petsc_shell_matrix() below, the pointer is casted back to a
     const ShellMatrix<T>*.  */

  ierr = MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
  LIBMESH_CHKERRABORT(ierr);
  ierr = MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
  LIBMESH_CHKERRABORT(ierr);

  return _solve_generalized_helper (matrix_A->mat(), mat_B, nev, ncv, tol, m_its);
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_generalized	(	ShellMatrix< T > &	matrix_A,
		ShellMatrix< T > &	matrix_B,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

Solve generalized eigenproblem when both matrix_A and matrix_B are of type ShellMatrix. When using this function, one should use the command line options: -st_ksp_type gmres -st_pc_type none or -st_ksp_type gmres -st_pc_type jacobi or similar.

Implements libMesh::EigenSolver< T >.

Definition at line 394 of file slepc_eigen_solver.C.

References libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::ShellMatrix< T >::m(), and libMesh::ShellMatrix< T >::n().

{
  this->init ();

  PetscErrorCode ierr=0;

  // Prepare the matrix.
  Mat mat_A;
  ierr = MatCreateShell(this->comm().get(),
                        shell_matrix_A.m(), // Specify the number of local rows
                        shell_matrix_A.n(), // Specify the number of local columns
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        const_cast<void*>(static_cast<const void*>(&shell_matrix_A)),
                        &mat_A);
  LIBMESH_CHKERRABORT(ierr);

  Mat mat_B;
  ierr = MatCreateShell(this->comm().get(),
                        shell_matrix_B.m(), // Specify the number of local rows
                        shell_matrix_B.n(), // Specify the number of local columns
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        const_cast<void*>(static_cast<const void*>(&shell_matrix_B)),
                        &mat_B);
  LIBMESH_CHKERRABORT(ierr);

  /* Note that the const_cast above is only necessary because PETSc
     does not accept a const void*.  Inside the member function
     _petsc_shell_matrix() below, the pointer is casted back to a
     const ShellMatrix<T>*.  */

  ierr = MatShellSetOperation(mat_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
  LIBMESH_CHKERRABORT(ierr);
  ierr = MatShellSetOperation(mat_A,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
  LIBMESH_CHKERRABORT(ierr);

  ierr = MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
  LIBMESH_CHKERRABORT(ierr);
  ierr = MatShellSetOperation(mat_B,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
  LIBMESH_CHKERRABORT(ierr);

  return _solve_generalized_helper (mat_A, mat_B, nev, ncv, tol, m_its);
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_standard	(	SparseMatrix< T > &	matrix_A,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

This function calls the SLEPc solver to compute the eigenpairs of the SparseMatrix matrix_A. nev is the number of eigenpairs to be computed and ncv is the number of basis vectors to be used in the solution procedure. Return values are the number of converged eigen values and the number of the iterations carried out by the eigen solver.

Implements libMesh::EigenSolver< T >.

Definition at line 99 of file slepc_eigen_solver.C.

References libMesh::PetscMatrix< T >::close(), libMesh::TriangleWrapper::init(), and libMesh::PetscMatrix< T >::mat().

{
  //   START_LOG("solve_standard()", "SlepcEigenSolver");

  this->init ();

  // Make sure the SparseMatrix passed in is really a PetscMatrix
  PetscMatrix<T>* matrix_A   = cast_ptr<PetscMatrix<T>*>(&matrix_A_in);

  // Close the matrix and vectors in case this wasn't already done.
  matrix_A->close ();

  // just for debugging, remove this
  //   char mat_file[] = "matA.petsc";
  //   PetscViewer petsc_viewer;
  //   ierr = PetscViewerBinaryOpen(this->comm().get(), mat_file, PETSC_FILE_CREATE, &petsc_viewer);
  //          LIBMESH_CHKERRABORT(ierr);
  //   ierr = MatView(matrix_A->mat(),petsc_viewer);
  //          LIBMESH_CHKERRABORT(ierr);

  return _solve_standard_helper(matrix_A->mat(), nev, ncv, tol, m_its);
}

template<typename T >

std::pair< unsigned int, unsigned int > libMesh::SlepcEigenSolver< T >::solve_standard	(	ShellMatrix< T > &	shell_matrix,
		int	nev,
		int	ncv,
		const double	tol,
		const unsigned int	m_its
	)		`[virtual]`

Same as above except that matrix_A is a ShellMatrix in this case.

Implements libMesh::EigenSolver< T >.

Definition at line 129 of file slepc_eigen_solver.C.

References libMesh::ierr, libMesh::TriangleWrapper::init(), libMesh::ShellMatrix< T >::m(), and libMesh::ShellMatrix< T >::n().

{
  this->init ();

  PetscErrorCode ierr=0;

  // Prepare the matrix.
  Mat mat;
  ierr = MatCreateShell(this->comm().get(),
                        shell_matrix.m(), // Specify the number of local rows
                        shell_matrix.n(), // Specify the number of local columns
                        PETSC_DETERMINE,
                        PETSC_DETERMINE,
                        const_cast<void*>(static_cast<const void*>(&shell_matrix)),
                        &mat);
  LIBMESH_CHKERRABORT(ierr);

  /* Note that the const_cast above is only necessary because PETSc
     does not accept a const void*.  Inside the member function
     _petsc_shell_matrix() below, the pointer is casted back to a
     const ShellMatrix<T>*.  */

  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
  LIBMESH_CHKERRABORT(ierr);
  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
  LIBMESH_CHKERRABORT(ierr);


  return _solve_standard_helper(mat, nev, ncv, tol, m_its);
}