From: Toby D. Young Date: Mon, 13 Jul 2009 13:23:50 +0000 (+0000) Subject: New to the SLEPcWrappers: A shortcut for the (unlikely?) case where the mass matrix... X-Git-Tag: v8.0.0~7489 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8eb1a61cc38dbb565eef66ba869fc259d451832b;p=dealii.git New to the SLEPcWrappers: A shortcut for the (unlikely?) case where the mass matrix is the identity matrix git-svn-id: https://svn.dealii.org/trunk@19070 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 142dee9b73..3ff6913670 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -180,6 +180,18 @@ inconvenience this causes.

lac

+
    +
  1. +

    + New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now + have a handle on the generalized eigenvalue problem where B=I. +
    + (Toby D. Young 2009/06/25) +

    +
  2. +
+ +
  1. @@ -214,14 +226,16 @@ inconvenience this causes.

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

      diff --git a/deal.II/lac/include/lac/slepc_solver.h b/deal.II/lac/include/lac/slepc_solver.h index fef1bbeffd..ac5cc49ea7 100755 --- a/deal.II/lac/include/lac/slepc_solver.h +++ b/deal.II/lac/include/lac/slepc_solver.h @@ -19,8 +19,9 @@ 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. + * solver. On the other hand, note that: the AdditionalData structure + * is a dummy structure and is there for backward/forward + * compatibility. * * SLEPcWrappers can be implemented in application codes in the * following way: @@ -30,13 +31,14 @@ DEAL_II_NAMESPACE_OPEN mpi_communicator); system.solve (A, B, lambda, x, n_eigenvectors); @endverbatim - * for the generalized eigenvalue problem $Ax=B\lambda x$. + * for the generalized eigenvalue problem $Ax=B\lambda x$. See also + * @ref step 36 "step-36" for a hands-on example. * * 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 + * rather than just one. This freedom is intended for use of the + * SLEPcWrappers that require a greater handle on the eigenvalue * problem solver context. See also: @verbatim template @@ -45,8 +47,8 @@ DEAL_II_NAMESPACE_OPEN const PETScWrappers::MatrixBase &B, std::vector &kr, std::vector &vr, - const unsigned int n_eigenvectors - ) + const unsigned int n_eigenvectors) + {code...} @endverbatim * as an example on how to do this. * @@ -92,13 +94,13 @@ namespace SLEPcWrappers /** * 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 + * eigensystem $Ax=\lambda x$. 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 @@ -113,9 +115,9 @@ namespace SLEPcWrappers * reset. * * Note that the number of - * converged eigenstates can be + * converged eigenvectors can be * larger than the number of - * eigenstates requested; this is + * eigenvectors requested; this is * due to a round off error * (success) of the eigenvalue * solver context. If this is found @@ -124,6 +126,18 @@ namespace SLEPcWrappers * requested but handle that it may * be more by ignoring any extras. */ + template + void + solve (const PETScWrappers::MatrixBase &A, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors); + + /** + * Same as above, but here is a + * composite method for solving the + * system $A x=\lambda B x$. + */ template void solve (const PETScWrappers::MatrixBase &A, @@ -134,23 +148,32 @@ namespace SLEPcWrappers /** * Initialize solver for the linear - * system $Ax=\lambda - * Bx$. (required before calling - * solve) + * system $Ax=\lambda x$. (Note: + * this is required before calling + * solve ()) + */ + void + set_matrices (const PETScWrappers::MatrixBase &A); + + /** + * Same as above, but here is a + * composite method for solving the + * system $A x=\lambda B x$. */ void set_matrices (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B); /** - * Set the initial vector for the solver. + * 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 + * to be used. */ void set_transformation (SLEPcWrappers::TransformationBase &trans); @@ -170,9 +193,10 @@ namespace SLEPcWrappers * Solve the linear system for * n_eigenvectors * eigenstates. Parameter - * n_converged contains the actual - * number of eigenstates that have - * . converged; this can be both + * 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. @@ -211,7 +235,7 @@ namespace SLEPcWrappers * Access to object that controls * convergence. */ - SolverControl & control() const; + SolverControl &control() const; /** * Exceptions. @@ -458,44 +482,70 @@ namespace SLEPcWrappers // --------------------------- inline and template functions ----------- - - /** - * This is declared here to make it - * possible to take a std::vector - * of different PETScWrappers vector - * types + /** + * This is declared here to make it possible to take a std::vector + * of different PETScWrappers vector types */ + + template + void + SolverBase::solve (const PETScWrappers::MatrixBase &A, + std::vector &kr, + std::vector &vr, + const unsigned int n_eigenvectors = 1) + { + unsigned int n_converged; + + set_matrices(A); + 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]); + } + } + + template void SolverBase::solve (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B, std::vector &kr, std::vector &vr, - const unsigned int n_eigenvectors = 0) + const unsigned int n_eigenvectors = 1) { 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 diff --git a/deal.II/lac/source/slepc_solver.cc b/deal.II/lac/source/slepc_solver.cc index f9042bece1..943cdc73b9 100644 --- a/deal.II/lac/source/slepc_solver.cc +++ b/deal.II/lac/source/slepc_solver.cc @@ -44,7 +44,14 @@ namespace SLEPcWrappers if( solver_data != 0 ) solver_data.reset (); } - + + void + SolverBase::set_matrices (const PETScWrappers::MatrixBase &A) + { + opA = &A; + opB = NULL; + } + void SolverBase::set_matrices (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B) @@ -84,19 +91,22 @@ namespace SLEPcWrappers ierr = EPSCreate (mpi_communicator, &solver_data->eps); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); - AssertThrow (opA && opB, ExcSLEPcWrappersUsageError()); - ierr = EPSSetOperators (solver_data->eps, *opA, *opB); + 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)); - - if( ini_vec && ini_vec->size() != 0 ) + + if (ini_vec && ini_vec->size() != 0) { ierr = EPSSetInitialVector(solver_data->eps, *ini_vec); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); } - if( transform ) + if (transform) transform->set_context(solver_data->eps); - + // set runtime options. set_solver_type (solver_data->eps); @@ -118,7 +128,8 @@ namespace SLEPcWrappers // get number of converged // eigenstates - ierr = EPSGetConverged (solver_data->eps, reinterpret_cast(n_converged)); + ierr = EPSGetConverged (solver_data->eps, + reinterpret_cast(n_converged)); AssertThrow (ierr == 0, ExcSLEPcError(ierr)); } @@ -277,8 +288,6 @@ namespace SLEPcWrappers } - - DEAL_II_NAMESPACE_CLOSE #else