From 0f641f2ccf73491903a50121bb909093e3310dde Mon Sep 17 00:00:00 2001 From: young Date: Fri, 10 May 2013 11:51:21 +0000 Subject: [PATCH] Get imaginary parts of eigenpairs. git-svn-id: https://svn.dealii.org/trunk@29491 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/lac/slepc_solver.h | 227 ++++++++++++++------- deal.II/source/lac/slepc_solver.cc | 16 ++ 2 files changed, 169 insertions(+), 74 deletions(-) diff --git a/deal.II/include/deal.II/lac/slepc_solver.h b/deal.II/include/deal.II/lac/slepc_solver.h index 402f4ebf69..fb74ab55f7 100644 --- a/deal.II/include/deal.II/lac/slepc_solver.h +++ b/deal.II/include/deal.II/lac/slepc_solver.h @@ -84,7 +84,7 @@ DEAL_II_NAMESPACE_OPEN const PETScWrappers::MatrixBase &B, std::vector &kr, std::vector &vr, - const unsigned int n_eigenvectors) + const unsigned int n_eigenpairs) { ... } @endcode * as an example on how to do this. @@ -134,11 +134,11 @@ namespace SLEPcWrappers * 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 vr is - * being used (and can be more). To avoid doing this, the fairly - * standard calling sequence executed here is used: Initialise; - * Set up matrices for solving; Actually solve the system; Gather - * the solution(s); and reset. + * eigenvectors, at least twice the memory size of + * eigenvectors is being used (and can be more). To avoid + * doing this, the fairly standard calling sequence executed here + * is used: Initialise; Set up matrices for solving; Actually + * solve the system; Gather the solution(s); and reset. * * @note Note that the number of converged eigenvectors can be * larger than the number of eigenvectors requested; this is due @@ -151,8 +151,8 @@ namespace SLEPcWrappers template void solve (const PETScWrappers::MatrixBase &A, - std::vector &r_eigenvalues, - std::vector &r_eigenvectors = std::vector (), + std::vector &eigenvalues, + std::vector &eigenvectors, const unsigned int n_eigenvectors = 1); /** @@ -163,8 +163,23 @@ namespace SLEPcWrappers void solve (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B, - std::vector &r_eigenvalues, - std::vector &r_eigenvectors = std::vector (), + std::vector &eigenvalues, + std::vector &eigenvectors, + const unsigned int n_eigenvectors = 1); + + /** + * Same as above, but here a composite method for solving the + * system $A x=\lambda B x$ with real and imaginary parts to the + * eigenvalues/eigenvectors. + */ + template + void + solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &real_eigenvalues, + std::vector &imag_eigenvalues, + std::vector &real_eigenvectors, + std::vector &imag_eigenvectors, const unsigned int n_eigenvectors = 1); /** @@ -283,14 +298,26 @@ namespace SLEPcWrappers solve (const unsigned int n_eigenvectors, unsigned int *n_converged); /** - * Access the solutions for a solved eigenvector problem, pair - * index solutions, $\text{index}\,\in\,0\hdots + * Access the real parts of solutions for a solved eigenvector + * problem, pair index solutions, $\text{index}\,\in\,0\hdots * \text{n\_converged}-1$. */ void get_eigenpair (const unsigned int index, - double &kr, - PETScWrappers::VectorBase &vr); + double &eigenvalues, + PETScWrappers::VectorBase &eigenvectors); + + /** + * Access the real and imaginary parts of solutions for a solved + * eigenvector problem, pair index solutions, + * $\text{index}\,\in\,0\hdots \text{n\_converged}-1$. + */ + void + get_eigenpair (const unsigned int index, + double &real_eigenvalues, + double &imag_eigenvalues, + PETScWrappers::VectorBase &real_eigenvectors, + PETScWrappers::VectorBase &imag_eigenvectors); /** * Initialize solver for the linear system $Ax=\lambda x$. (Note: @@ -711,74 +738,126 @@ namespace SLEPcWrappers * This is declared here to make it possible to take a std::vector * of different PETScWrappers vector types */ + // todo: The logic of these functions can be simplified without breaking backward compatibility... template - void - SolverBase::solve (const PETScWrappers::MatrixBase &A, - std::vector &kr, - std::vector &vr, - const unsigned int n_eigenvectors) - { - // Panic if the number of eigenpairs wanted is out of bounds. - AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), - ExcSLEPcWrappersUsageError()); - - // Set the matrices of the problem - set_matrices (A); - - // and solve - unsigned int n_converged = 0; - 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() != 0, ExcSLEPcWrappersUsageError()); - vr.resize (n_converged, vr.front()); - kr.resize (n_converged); - - for (unsigned int index=0; index &eigenvalues, + std::vector &eigenvectors, + const unsigned int n_eigenvectors) + { + // Panic if the number of eigenpairs wanted is out of bounds. + AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), + ExcSLEPcWrappersUsageError()); + + // Set the matrices of the problem + set_matrices (A); + + // and solve + unsigned int n_converged = 0; + 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 (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError()); + eigenvectors.resize (n_converged, eigenvectors.front()); + eigenvalues.resize (n_converged); + + for (unsigned int index=0; index - void + void SolverBase::solve (const PETScWrappers::MatrixBase &A, const PETScWrappers::MatrixBase &B, - std::vector &kr, - std::vector &vr, + std::vector &eigenvalues, + std::vector &eigenvectors, const unsigned int n_eigenvectors) - { - // Guard against incompatible matrix sizes: - AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m())); - AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n())); + { + // Guard against incompatible matrix sizes: + AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m())); + AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n())); + + // Panic if the number of eigenpairs wanted is out of bounds. + AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), + ExcSLEPcWrappersUsageError()); + + // Set the matrices of the problem + set_matrices (A, B); + + // and solve + unsigned int n_converged = 0; + 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 (eigenvectors.size() != 0, ExcSLEPcWrappersUsageError()); + + eigenvectors.resize (n_converged, eigenvectors.front()); + eigenvalues.resize (n_converged); + + for (unsigned int index=0; index 0) && (n_eigenvectors <= A.m ()), - ExcSLEPcWrappersUsageError()); - - // Set the matrices of the problem - set_matrices (A, B); - - // and solve - unsigned int n_converged = 0; - 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() != 0, ExcSLEPcWrappersUsageError()); - - vr.resize (n_converged, vr.front()); - kr.resize (n_converged); + template + void + SolverBase::solve (const PETScWrappers::MatrixBase &A, + const PETScWrappers::MatrixBase &B, + std::vector &real_eigenvalues, + std::vector &imag_eigenvalues, + std::vector &real_eigenvectors, + std::vector &imag_eigenvectors, + const unsigned int n_eigenvectors) + { + // Guard against incompatible matrix sizes: + AssertThrow (A.m() == B.m (), ExcDimensionMismatch(A.m(), B.m())); + AssertThrow (A.n() == B.n (), ExcDimensionMismatch(A.n(), B.n())); + + // and incompatible eigenvalue/eigenvector sizes + AssertThrow (real_eigenvalues.size() == imag_eigenvalues.size(), + ExcDimensionMismatch(real_eigenvalues.size(), imag_eigenvalues.size())); + AssertThrow (real_eigenvectors.size() == imag_eigenvectors.n (), + ExcDimensionMismatch(real_eigenvectors.size(), imag_eigenvectors.size())); + + // Panic if the number of eigenpairs wanted is out of bounds. + AssertThrow ((n_eigenvectors > 0) && (n_eigenvectors <= A.m ()), + ExcSLEPcWrappersUsageError()); + + // Set the matrices of the problem + set_matrices (A, B); + + // and solve + unsigned int n_converged = 0; + 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 ((real_eigenvectors.size()!=0) && (imag_eigenvectors.size()!=0), + ExcSLEPcWrappersUsageError()); + + real_eigenvectors.resize (n_converged, real_eigenvectors.front()); + imag_eigenvectors.resize (n_converged, imag_eigenvectors.front()); + real_eigenvalues.resize (n_converged); + imag_eigenvalues.resize (n_converged); + + for (unsigned int index=0; indexeps, index, + &real_eigenvalues, &imag_eigenvalues, + real_eigenvectors, imag_eigenvectors); + AssertThrow (ierr == 0, ExcSLEPcError(ierr)); + } + void SolverBase::reset () { -- 2.39.5