From 9a03bc439277d85502accef56d8f9048558ae925 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Sat, 24 Jan 2015 17:37:34 +0100 Subject: [PATCH] added an PArpack eigensolver class --- cmake/configure/configure_arpack.cmake | 1 + cmake/modules/FindARPACK.cmake | 31 +- include/deal.II/base/config.h.in | 2 + include/deal.II/lac/parpack_solver.h | 870 +++++++++++++++++++++++++ 4 files changed, 902 insertions(+), 2 deletions(-) create mode 100644 include/deal.II/lac/parpack_solver.h diff --git a/cmake/configure/configure_arpack.cmake b/cmake/configure/configure_arpack.cmake index 7cffb145bd..5ac73fed24 100644 --- a/cmake/configure/configure_arpack.cmake +++ b/cmake/configure/configure_arpack.cmake @@ -20,3 +20,4 @@ SET(FEATURE_ARPACK_DEPENDS LAPACK) CONFIGURE_FEATURE(ARPACK) +SET(DEAL_II_ARPACK_WITH_PARPACK ARPACK_WITH_PARPACK) diff --git a/cmake/modules/FindARPACK.cmake b/cmake/modules/FindARPACK.cmake index 0e185cbc8b..509f490ce7 100644 --- a/cmake/modules/FindARPACK.cmake +++ b/cmake/modules/FindARPACK.cmake @@ -20,6 +20,7 @@ # # ARPACK_LIBRARIES # ARPACK_LINKER_FLAGS +# ARPACK_WITH_PARPACK # # @@ -35,8 +36,34 @@ DEAL_II_FIND_LIBRARY(ARPACK_LIBRARY PATH_SUFFIXES lib${LIB_SUFFIX} lib64 lib ) +# +# Sanity check: Only search the parpack library in the same directory as +# the arpack library... +# +GET_FILENAME_COMPONENT(_path "${ARPACK_LIBRARY}" PATH) +DEAL_II_FIND_LIBRARY(PARPACK_LIBRARY + NAMES parpack + HINTS ${_path} + NO_DEFAULT_PATH + NO_CMAKE_ENVIRONMENT_PATH + NO_CMAKE_PATH + NO_SYSTEM_ENVIRONMENT_PATH + NO_CMAKE_SYSTEM_PATH + NO_CMAKE_FIND_ROOT_PATH + ) + +IF(NOT PARPACK_LIBRARY MATCHES "-NOTFOUND") + SET(ARPACK_WITH_PARPACK TRUE) +ELSE() + SET(ARPACK_WITH_PARPACK FALSE) +ENDIF() + + DEAL_II_PACKAGE_HANDLE(ARPACK - LIBRARIES REQUIRED ARPACK_LIBRARY LAPACK_LIBRARIES + LIBRARIES + OPTIONAL PARPACK_LIBRARY + REQUIRED ARPACK_LIBRARY LAPACK_LIBRARIES + OPTIONAL MPI_C_LIBRARIES LINKER_FLAGS OPTIONAL LAPACK_LINKER_FLAGS - CLEAR ARPACK_LIBRARY + CLEAR ARPACK_LIBRARY PARPACK_LIBRARY ) diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in index fda7109f58..dc5f57ccfa 100644 --- a/include/deal.II/base/config.h.in +++ b/include/deal.II/base/config.h.in @@ -62,6 +62,8 @@ #cmakedefine DEAL_II_WITH_UMFPACK #cmakedefine DEAL_II_WITH_ZLIB +#cmakedefine DEAL_II_ARPACK_WITH_PARPACK + /*********************************************************************** * Compiler bugs: diff --git a/include/deal.II/lac/parpack_solver.h b/include/deal.II/lac/parpack_solver.h new file mode 100644 index 0000000000..10861b29bd --- /dev/null +++ b/include/deal.II/lac/parpack_solver.h @@ -0,0 +1,870 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2010 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef __deal2__parpack_solver_h +#define __deal2__parpack_solver_h + +#include +#include +#include +#include +#include +#include + +#include + + +#ifdef DEAL_II_ARPACK_WITH_PARPACK + +DEAL_II_NAMESPACE_OPEN + +extern "C" { + + // http://www.mathkeisan.com/usersguide/man/pdnaupd.html + void pdnaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which, + int *nev, double *tol, double *resid, int *ncv, + double *v, int *nloc, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, + int *info); + + // http://www.mathkeisan.com/usersguide/man/pdsaupd.html + void pdsaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which, + int *nev, double *tol, double *resid, int *ncv, + double *v, int *nloc, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, + int *info); + + // http://www.mathkeisan.com/usersguide/man/pdneupd.html + void pdneupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d, + double *di, double *z, int *ldz, double *sigmar, + double *sigmai, double *workev, char *bmat, int *n, char *which, + int *nev, double *tol, double *resid, int *ncv, + double *v, int *nloc, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, int *info); + + // http://www.mathkeisan.com/usersguide/man/pdseupd.html + void pdseupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d, + double *z, int *ldz, double *sigmar, + char *bmat, int *n, char *which, + int *nev, double *tol, double *resid, int *ncv, + double *v, int *nloc, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, int *info); + + // other resources: + // http://acts.nersc.gov/superlu/example5/pnslac.c.html + // https://github.com/phpisciuneri/tijo/blob/master/dvr_parpack.cpp + +} + +/** + * Interface for using PARPACK. PARPACK is a collection of Fortran77 + * subroutines designed to solve large scale eigenvalue problems. + * Here we interface to the routines pdneupd, + * pdseupd, pdnaupd, pdsaupd of + * PARPACK. The package is designed to compute a few eigenvalues and + * corresponding eigenvectors of a general n by n matrix A. It is most + * appropriate for large sparse matrices A. + * + * In this class we make use of the method applied to the generalized + * eigenspectrum problem $(A-\lambda B)x=0$, for $x\neq0$; where $A$ + * is a system matrix, $B$ is a mass matrix, and $\lambda, x$ are a + * set of eigenvalues and eigenvectors respectively. + * + * The ArpackSolver can be used in application codes in the + * following way: + @code + SolverControl solver_control (1000, 1e-9); + const unsigned int num_arnoldi_vectors = 2*size_of_spectrum + 2; + PArpackSolver::AdditionalData + additional_data(num_arnoldi_vectors, + dealii::PArpackSolver::largest_magnitude, + true); + + PArpackSolver eigensolver (solver_control, + mpi_communicator, + additional_data); + eigensolver.set_shift(sigma); + eigensolver.reinit(locally_owned_dofs); + eigensolver.solve (A, + B, + OP, + lambda, + x, + size_of_spectrum); + @endcode + * for the generalized eigenvalue problem $Ax=B\lambda x$, where the + * variable size_of_spectrum tells PARPACK the number of + * eigenvector/eigenvalue pairs to solve for. Here, + * lambda is a vector that will contain the eigenvalues + * computed, x a vector of objects of type V + * that will contain the eigenvectors computed. OP is an + * inverse operation for the matrix A - sigma * B, where + * sigma is a shift value, set to zero by default. + * + * Through the AdditionalData the user can specify some of the + * parameters to be set. + * + * The class is intended to be used with MPI and can work on arbitrary + * vector and matrix distributed classes. Both symmetric and + * non-symmetric A are supported. + * + * For further information on how the PARPACK routines + * pdneupd, pdseupd, pdnaupd, + * pdsaupd work and also how to set the parameters + * appropriately please take a look into the PARPACK manual. + * + * @author Denis Davydov, 2014. + */ +template +class PArpackSolver : public Subscriptor +{ +public: + /** + * Declare the type for container size. + */ + typedef types::global_dof_index size_type; + + /** + * An enum that lists the possible choices for which eigenvalues to + * compute in the solve() function. + * + * A particular choice is limited based on symmetric or + * non-symmetric matrix A considered. + */ + enum WhichEigenvalues + { + algebraically_largest, + algebraically_smallest, + largest_magnitude, + smallest_magnitude, + largest_real_part, + smallest_real_part, + largest_imaginary_part, + smallest_imaginary_part, + both_ends + }; + + /** + * Auxiliary class to represent A-sigma*B operator. + */ + template + class Shift : public dealii::Subscriptor + { + public: + + /** + * Constructor. + */ + Shift (const MATRIX &A, + const MATRIX &B, + const double sigma) + : + A(A), + B(B), + sigma(sigma) + {} + + /** + * Apply A-sigma * B + */ + void vmult (VECTOR &dst, const VECTOR &src) const + { + B.vmult(dst,src); + dst *= (-sigma); + A.vmult_add(dst,src); + } + + /** + * Apply A^T-sigma * B^T + */ + void Tvmult (VECTOR &dst, const VECTOR &src) const + { + B.Tvmult(dst,src); + dst *= (-sigma); + A.Tvmult_add(dst,src); + } + + private: + const MATRIX &A; + const MATRIX &B; + const double sigma; + }; + + /** + * Standardized data struct to pipe additional data to the solver, + * should it be needed. + */ + struct AdditionalData + { + const unsigned int number_of_arnoldi_vectors; + const WhichEigenvalues eigenvalue_of_interest; + const bool symmetric; + AdditionalData( + const unsigned int number_of_arnoldi_vectors = 15, + const WhichEigenvalues eigenvalue_of_interest = largest_magnitude, + const bool symmetric = false); + }; + + /** + * Access to the object that controls convergence. + */ + SolverControl &control () const; + + /** + * Constructor. + */ + PArpackSolver(SolverControl &control, + const MPI_Comm &mpi_communicator, + const AdditionalData &data = AdditionalData()); + + /** + * Initialise internal variables. + */ + void reinit(const dealii::IndexSet &locally_owned_dofs ); + + /** + * Set desired shift value. + */ + void set_shift(const double s ); + + /** + * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by + * calling the pd(n/s)eupd and pd(n/s)aupd + * functions of PARPACK. + */ + template + void solve( + const MATRIX1 &A, + const MATRIX2 &B, + const INVERSE &inverse, + std::vector > &eigenvalues, + std::vector &eigenvectors, + const unsigned int n_eigenvalues); + + std::size_t memory_consumption() const; + +protected: + + /** + * Reference to the object that controls convergence of the + * iterative solver. + */ + SolverControl &solver_control; + + /** + * Store a copy of the flags for this particular solver. + */ + const AdditionalData additional_data; + + // keep MPI communicator non-const as Arpack functions are not const either: + + /** + * C++ MPI communicator. + */ + MPI_Comm mpi_communicator; + + /** + * Fortran MPI communicator. + */ + MPI_Fint mpi_communicator_fortran; + + // C++98 guarantees that the elements of a vector are stored contiguously + + /** + * Length of the work array workl. + */ + int lworkl; + + /** + * Double precision work array of length lworkl + */ + std::vector workl; + + /** + * Double precision work array of length 3*N + */ + std::vector workd; + + /** + * Number of local degrees of freedom. + */ + int nloc; + + /** + * Number of Arnoldi basis vectors specified in additional_data + */ + int ncv; + + + /** + * The leading dimension of the array v + */ + int ldv; + + /** + * Double precision vector of size ldv by NCV. Will contains the + * final set of Arnoldi basis vectors. + */ + std::vector v; + + /** + * The initial residual vector, possibly from a previous run. On + * output, it contains the final residual vector. + */ + std::vector resid; + + /** + * The leading dimension of the array Z equal to nloc. + */ + int ldz; + + /** + * A vector of minimum size of nloc by NEV+1. Z contains the + * B-orthonormal Ritz vectors of the eigensystem A*z = lambda*B*z + * corresponding to the Ritz value approximations. + */ + std::vector z; + + /** + * The size of the workev array. + */ + int lworkev; + + /** + * Double precision work array of dimension 3* NCV. + */ + std::vector workev; + + /** + * A vector of dimension NCV. + */ + std::vector select; + + /** + * Temporary vectors used between Arpack and deal.II + */ + VECTOR src,dst,tmp; + + /** + * Indices of local degrees of freedom. + */ + std::vector< types::global_dof_index > local_indices; + + /** + * The shift value to be applied during solution + */ + double shift_value; + +private: + + /** + * PArpackExcInfoPdnaupds. + */ + DeclException2 (PArpackExcConvergedEigenvectors, int, int, + << arg1 << "eigenpairs were requested, but only" + << arg2 << " converged"); + + DeclException2 (PArpackExcInvalidNumberofEigenvalues, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of the matrix " << arg2); + + DeclException2 (PArpackExcInvalidEigenvectorSize, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of eigenvectors " << arg2); + + DeclException2 (PArpackExcInvalidEigenvalueSize, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of eigenvalues " << arg2); + + DeclException2 (PArpackExcInvalidNumberofArnoldiVectors, int, int, + << "Number of Arnoldi vectors " << arg1 + << " is larger that the size of the matrix " << arg2); + + DeclException2 (PArpackExcSmallNumberofArnoldiVectors, int, int, + << "Number of Arnoldi vectors " << arg1 + << " is too small to obtain " << arg2 + << " eigenvalues"); + + DeclException1 (PArpackExcIdo, int, << "This ido " << arg1 + << " is not supported. Check documentation of ARPACK"); + + DeclException1 (PArpackExcMode, int, << "This mode " << arg1 + << " is not supported. Check documentation of ARPACK"); + + DeclException1 (PArpackExcInfoPdnaupd, int, + << "Error with Pdnaupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (PArpackExcInfoPdneupd, int, + << "Error with Pdneupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (PArpackExcInfoMaxIt, int, + << "Maximum number " << arg1 + << " of iterations reached."); + + DeclException1 (PArpackExcNoShifts, int, + << "No shifts could be applied during implicit" + << " Arnoldi update, try increasing the number of" + << " Arnoldi vectors."); +}; + +template +std::size_t +PArpackSolver::memory_consumption() const +{ + return MemoryConsumption::memory_consumption (double()) * + (workl.size() + + workd.size() + + v.size() + + resid.size() + + z.size() + + workev.size() ) + + src.memory_consumption() + + dst.memory_consumption() + + tmp.memory_consumption() + + MemoryConsumption::memory_consumption (types::global_dof_index()) * local_indices.size(); +} + +template +PArpackSolver::AdditionalData:: +AdditionalData (const unsigned int number_of_arnoldi_vectors, + const WhichEigenvalues eigenvalue_of_interest, + const bool symmetric) + : + number_of_arnoldi_vectors(number_of_arnoldi_vectors), + eigenvalue_of_interest(eigenvalue_of_interest), + symmetric(symmetric) +{} + +template +PArpackSolver::PArpackSolver (SolverControl &control, + const MPI_Comm &mpi_communicator, + const AdditionalData &data) + : + solver_control (control), + additional_data (data), + mpi_communicator( mpi_communicator ), + mpi_communicator_fortran ( MPI_Comm_c2f( mpi_communicator ) ), + shift_value(0.0) + +{} + +template +void PArpackSolver::set_shift(const double s ) +{ + shift_value = s; +} + +template +void PArpackSolver::reinit(const dealii::IndexSet &locally_owned_dofs) +{ + // store local indices to write to vectors + locally_owned_dofs.fill_index_vector(local_indices); + + // scalars + nloc = locally_owned_dofs.n_elements (); + ncv = additional_data.number_of_arnoldi_vectors; + + Assert (local_indices.size() == nloc, ExcInternalError() ); + + // vectors + ldv = nloc; + v.resize (ldv*ncv, 0.0); + + // TODO: add optional input for resid + resid.resize(nloc, 1.0); + + // work arrays for ARPACK + workd.resize(3*nloc,0.0); + + lworkl = additional_data.symmetric ? + ncv*ncv + 8*ncv + : + 3*ncv*ncv+6*ncv; + workl.resize (lworkl, 0.); + + ldz = nloc; + z.resize (ldz*ncv, 0.); // TODO we actually need only ldz*nev + + // WORKEV Double precision work array of dimension 3*NCV. + lworkev = additional_data.symmetric ? + 0 /*not used in symmetric case*/ + : + 3*ncv; + workev.resize (lworkev, 0.); + + select.resize (ncv, 0); + + // deal.II vectors: + src.reinit (locally_owned_dofs,mpi_communicator); + dst.reinit (locally_owned_dofs,mpi_communicator); + tmp.reinit (locally_owned_dofs,mpi_communicator); + +} + +template +template +void PArpackSolver::solve ( + const MATRIX1 &/*system_matrix*/, + const MATRIX2 &mass_matrix, + const INVERSE &inverse, + std::vector > &eigenvalues, + std::vector &eigenvectors, + const unsigned int n_eigenvalues) +{ + + Assert (n_eigenvalues <= eigenvectors.size(), + PArpackExcInvalidEigenvectorSize(n_eigenvalues, eigenvectors.size())); + + Assert (n_eigenvalues <= eigenvalues.size(), + PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size())); + + + Assert (n_eigenvalues < mass_matrix.m(), + PArpackExcInvalidNumberofEigenvalues(n_eigenvalues, mass_matrix.m())); + + Assert (additional_data.number_of_arnoldi_vectors < mass_matrix.m(), + PArpackExcInvalidNumberofArnoldiVectors( + additional_data.number_of_arnoldi_vectors, mass_matrix.m())); + + Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1, + PArpackExcSmallNumberofArnoldiVectors( + additional_data.number_of_arnoldi_vectors, n_eigenvalues)); + // ARPACK mode for dnaupd, here only mode 3 + int mode = 3; + + // reverse communication parameter + // must be zero on the first call to pdnaupd + int ido = 0; + + /** + * 'G' generalized eigenvalue problem + * 'I' standard eigenvalue problem + */ + char bmat[2] = "G"; + + /** Specify the eigenvalues of interest, + * possible parameters + * "LA" algebraically largest + * "SA" algebraically smallest + * "LM" largest magnitude + * "SM" smallest magnitude + * "LR" largest real part + * "SR" smallest real part + * "LI" largest imaginary part + * "SI" smallest imaginary part + * "BE" both ends of spectrum simultaneous + */ + char which[3]; + switch (additional_data.eigenvalue_of_interest) + { + case algebraically_largest: + std::strcpy (which, "LA"); + break; + case algebraically_smallest: + std::strcpy (which, "SA"); + break; + case largest_magnitude: + std::strcpy (which, "LM"); + break; + case smallest_magnitude: + std::strcpy (which, "SM"); + break; + case largest_real_part: + std::strcpy (which, "LR"); + break; + case smallest_real_part: + std::strcpy (which, "SR"); + break; + case largest_imaginary_part: + std::strcpy (which, "LI"); + break; + case smallest_imaginary_part: + std::strcpy (which, "SI"); + break; + case both_ends: + std::strcpy (which, "BE"); + break; + } + + // tolerance for ARPACK + double tol = control().tolerance(); + + //information to the routines + std::vector iparam (11, 0); + + iparam[0] = 1; + // shift strategy: exact shifts with respect to the current Hessenberg matrix H. + + // maximum number of iterations + iparam[2] = control().max_steps(); + + /** Sets the mode of dsaupd. + * 1 is exact shifting, + * 2 is user-supplied shifts, + * 3 is shift-invert mode, + * 4 is buckling mode, + * 5 is Cayley mode. + */ + + iparam[6] = mode; + std::vector ipntr (14, 0); + + //information out of the iteration + // If INFO .EQ. 0, a random initial residual vector is used. + // If INFO .NE. 0, RESID contains the initial residual vector, + // possibly from a previous run. + // Typical choices in this situation might be to use the final value + // of the starting vector from the previous eigenvalue calculation + int info = 1; + + // Number of eigenvalues of OP to be computed. 0 < NEV < N. + int nev = n_eigenvalues; + int n_inside_arpack = nloc; + + while (ido != 99) + { + // call of ARPACK pdnaupd routine + if (additional_data.symmetric) + pdsaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], + &workd[0], &workl[0], &lworkl, &info); + else + pdnaupd_(&mpi_communicator_fortran,&ido, bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], + &workd[0], &workl[0], &lworkl, &info); + + if (ido == 99) + break; + + switch (mode) + { +// OP = (inv[K - sigma*M])*M + case 3: + { + switch (ido) + { +// compute Y = OP * X where +// IPNTR(1) is the pointer into WORKD for X, +// IPNTR(2) is the pointer into WORKD for Y. + case -1: + { + const int shift_x = ipntr[0]-1; + const int shift_y = ipntr[1]-1; + Assert (shift_x>=0, dealii::ExcInternalError() ); + Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() ); + Assert (shift_y>=0, dealii::ExcInternalError() ); + Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() ); + + src = 0.0; + src.add (nloc, + &local_indices[0], + &workd[0]+shift_x ); + src.compress (VectorOperation::add); + + // multiplication with mass matrix M + mass_matrix.vmult(tmp, src); + // solving linear system + inverse.vmult(dst,tmp); + + // store the result + dst.extract_subvector_to (local_indices.begin(), + local_indices.end(), + &workd[0]+shift_y ); + } + break; + +// compute Y = OP * X where +// IPNTR(1) is the pointer into WORKD for X, +// IPNTR(2) is the pointer into WORKD for Y. +// In mode 3,4 and 5, the vector B * X is already +// available in WORKD(ipntr(3)). It does not +// need to be recomputed in forming OP * X. + case 1: + { + const int shift_x = ipntr[0]-1; + const int shift_y = ipntr[1]-1; + const int shift_b_x = ipntr[2]-1; + + Assert (shift_x>=0, dealii::ExcInternalError() ); + Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() ); + Assert (shift_y>=0, dealii::ExcInternalError() ); + Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() ); + Assert (shift_b_x>=0, dealii::ExcInternalError() ); + Assert (shift_b_x+nloc <= workd.size(), dealii::ExcInternalError() ); + Assert (shift_y>=0, dealii::ExcInternalError() ); + Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() ); + + src = 0.0; // B*X + src.add (nloc, + &local_indices[0], + &workd[0]+shift_b_x ); + + tmp = 0.0; // X + tmp.add (nloc, + &local_indices[0], + &workd[0]+shift_x); + + src.compress (VectorOperation::add); + tmp.compress (VectorOperation::add); + + // solving linear system + inverse.vmult(dst,src); + + // store the result + dst.extract_subvector_to (local_indices.begin(), + local_indices.end(), + &workd[0]+shift_y ); + + } + break; + +// compute Y = B * X where +// IPNTR(1) is the pointer into WORKD for X, +// IPNTR(2) is the pointer into WORKD for Y. + case 2: + { + + const int shift_x = ipntr[0]-1; + const int shift_y = ipntr[1]-1; + Assert (shift_x>=0, dealii::ExcInternalError() ); + Assert (shift_x+nloc <= workd.size(), dealii::ExcInternalError() ); + Assert (shift_y>=0, dealii::ExcInternalError() ); + Assert (shift_y+nloc <= workd.size(), dealii::ExcInternalError() ); + + src = 0.0; + src.add (nloc, + &local_indices[0], + &workd[0]+shift_x ); + src.compress (VectorOperation::add); + + // Multiplication with mass matrix M + mass_matrix.vmult(dst, src); + + // store the result + dst.extract_subvector_to (local_indices.begin(), + local_indices.end(), + &workd[0]+shift_y); + + } + break; + + default: + Assert (false, PArpackExcIdo(ido)); + break; + } + } + break; + default: + Assert (false, PArpackExcMode(mode)); + break; + } + } + + if (info<0) + { + Assert (false, PArpackExcInfoPdnaupd(info)); + } + else + { + /** 1 - compute eigenvectors, + * 0 - only eigenvalues + */ + int rvec = 1; + + // which eigenvectors + char howmany[4] = "All"; + + double sigmar = shift_value; // real part of the shift + double sigmai = 0.0; // imaginary part of the shift + + std::vector eigenvalues_real (n_eigenvalues, 0.); + std::vector eigenvalues_im (n_eigenvalues, 0.); + + // call of ARPACK pdneupd routine + if (additional_data.symmetric) + pdseupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0], + &z[0], &ldz, &sigmar, + bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, + &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); + else + pdneupd_(&mpi_communicator_fortran, &rvec, howmany, &select[0], &eigenvalues_real[0], + &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai, + &workev[0], bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, + &iparam[0], &ipntr[0], &workd[0], &workl[0], &lworkl, &info); + + if (info == 1) + { + Assert (false, PArpackExcInfoMaxIt(control().max_steps())); + } + else if (info == 3) + { + Assert (false, PArpackExcNoShifts(1)); + } + else if (info!=0) + { + Assert (false, PArpackExcInfoPdneupd(info)); + } + + for (size_type i=0; i (eigenvalues_real[i], + eigenvalues_im[i]); + } + + Assert (iparam[4] == n_eigenvalues, + PArpackExcConvergedEigenvectors(iparam[4], n_eigenvalues)); + + // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M + // with respect to a semi-inner product defined by M. + + // resid likely contains residual with respect to M-norm. + { + + tmp = 0.0; + tmp.add (nloc, + &local_indices[0], + &resid[0]); + solver_control.check ( iparam[2], tmp.l2_norm() ); + } + + +} + +template +SolverControl &PArpackSolver::control () const +{ + return solver_control; +} + +DEAL_II_NAMESPACE_CLOSE + + +#endif +#endif -- 2.39.5