15#ifndef dealii_parpack_solver_h
16#define dealii_parpack_solver_h
30#ifdef DEAL_II_ARPACK_WITH_PARPACK
209template <
typename VectorType>
348 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
360 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
436 std::vector<double>
v;
459 std::vector<double>
z;
512 <<
arg1 <<
" eigenpairs were requested, but only " <<
arg2
518 <<
"Number of wanted eigenvalues " <<
arg1
519 <<
" is larger that the size of the matrix " <<
arg2);
524 <<
"Number of wanted eigenvalues " <<
arg1
525 <<
" is larger that the size of eigenvectors " <<
arg2);
531 <<
"To store the real and complex parts of " <<
arg1
532 <<
" eigenvectors in real-valued vectors, their size (currently set to "
533 <<
arg2 <<
") should be greater than or equal to " <<
arg1 + 1);
538 <<
"Number of wanted eigenvalues " <<
arg1
539 <<
" is larger that the size of eigenvalues " <<
arg2);
544 <<
"Number of Arnoldi vectors " <<
arg1
545 <<
" is larger that the size of the matrix " <<
arg2);
550 <<
"Number of Arnoldi vectors " <<
arg1
551 <<
" is too small to obtain " <<
arg2 <<
" eigenvalues");
555 <<
"This ido " <<
arg1
556 <<
" is not supported. Check documentation of ARPACK");
560 <<
"This mode " <<
arg1
561 <<
" is not supported. Check documentation of ARPACK");
565 <<
"Error with Pdnaupd, info " <<
arg1
566 <<
". Check documentation of ARPACK");
570 <<
"Error with Pdneupd, info " <<
arg1
571 <<
". Check documentation of ARPACK");
575 <<
"Maximum number " <<
arg1 <<
" of iterations reached.");
579 <<
"No shifts could be applied during implicit"
580 <<
" Arnoldi update, try increasing the number of"
581 <<
" Arnoldi vectors.");
586template <
typename VectorType>
591 (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
593 src.memory_consumption() + dst.memory_consumption() +
594 tmp.memory_consumption() +
596 local_indices.size();
601template <
typename VectorType>
603 const unsigned int number_of_arnoldi_vectors,
605 const bool symmetric,
607 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
608 , eigenvalue_of_interest(eigenvalue_of_interest)
609 , symmetric(symmetric)
618 "'largest real part' can only be used for non-symmetric problems!"));
622 "'smallest real part' can only be used for non-symmetric problems!"));
626 "'largest imaginary part' can only be used for non-symmetric problems!"));
630 "'smallest imaginary part' can only be used for non-symmetric problems!"));
633 ExcMessage(
"Currently, only modes 1, 2 and 3 are supported."));
638template <
typename VectorType>
659template <
typename VectorType>
663 sigmar =
sigma.real();
664 sigmai =
sigma.imag();
669template <
typename VectorType>
673 initial_vector_provided =
true;
674 Assert(resid.size() == local_indices.size(),
676 vec.extract_subvector_to(local_indices.begin(),
683template <
typename VectorType>
692 ncv = additional_data.number_of_arnoldi_vectors;
698 v.resize(ldv * ncv, 0.0);
700 resid.resize(nloc, 1.0);
703 workd.resize(3 * nloc, 0.0);
706 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
707 workl.resize(lworkl, 0.);
710 z.resize(ldz * ncv, 0.);
713 lworkev = additional_data.symmetric ? 0
716 workev.resize(lworkev, 0.);
718 select.resize(ncv, 0);
723template <
typename VectorType>
727 internal_reinit(locally_owned_dofs);
730 src.reinit(locally_owned_dofs, mpi_communicator);
731 dst.reinit(locally_owned_dofs, mpi_communicator);
732 tmp.reinit(locally_owned_dofs, mpi_communicator);
737template <
typename VectorType>
751template <
typename VectorType>
756 internal_reinit(locally_owned_dofs);
766template <
typename VectorType>
767template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
784template <
typename VectorType>
785template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
794 if (additional_data.symmetric)
816 PArpackExcInvalidNumberofArnoldiVectors(
817 additional_data.number_of_arnoldi_vectors,
eigenvectors[0]->size()));
820 PArpackExcSmallNumberofArnoldiVectors(
823 int mode = additional_data.mode;
832 bmat[0] = (mode == 1) ?
'I' :
'G';
846 switch (additional_data.eigenvalue_of_interest)
848 case algebraically_largest:
849 std::strcpy(
which,
"LA");
851 case algebraically_smallest:
852 std::strcpy(
which,
"SA");
854 case largest_magnitude:
855 std::strcpy(
which,
"LM");
857 case smallest_magnitude:
858 std::strcpy(
which,
"SM");
860 case largest_real_part:
861 std::strcpy(
which,
"LR");
863 case smallest_real_part:
864 std::strcpy(
which,
"SR");
866 case largest_imaginary_part:
867 std::strcpy(
which,
"LI");
869 case smallest_imaginary_part:
870 std::strcpy(
which,
"SI");
873 std::strcpy(
which,
"BE");
878 double tol = control().tolerance();
881 std::vector<int>
iparam(11, 0);
888 iparam[2] = control().max_steps();
901 std::vector<int>
ipntr(14, 0);
909 int info = initial_vector_provided ? 1 : 0;
919 if (additional_data.symmetric)
976 if ((
ido == -1) || (
ido == 1 && mode < 3))
979 src.add(nloc, local_indices.data(), workd.data() +
shift_x);
985 mass_matrix.vmult(tmp, src);
986 inverse.vmult(dst, tmp);
991 system_matrix.vmult(tmp, src);
993 tmp.extract_subvector_to(local_indices.begin(),
996 inverse.vmult(dst, tmp);
1000 system_matrix.vmult(dst, src);
1005 else if (
ido == 1 && mode >= 3)
1015 src.add(nloc, local_indices.data(), workd.data() +
shift_b_x);
1020 inverse.vmult(dst, src);
1025 src.add(nloc, local_indices.data(), workd.data() +
shift_x);
1036 mass_matrix.vmult(dst, src);
1045 dst.extract_subvector_to(local_indices.begin(),
1046 local_indices.end(),
1061 if (additional_data.symmetric)
1062 pdseupd_(&mpi_communicator_fortran,
1086 pdneupd_(&mpi_communicator_fortran,
1115 AssertThrow(
false, PArpackExcInfoMaxIt(control().max_steps()));
1126 for (
int i = 0; i <
nev; ++i)
1131 eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1149 tmp.add(nloc, local_indices.data(), resid.data());
1151 solver_control.check(
iparam[2], tmp.l2_norm());
1157template <
typename VectorType>
1161 return solver_control;
size_type n_elements() const
std::vector< size_type > get_index_vector() const
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
MPI_Comm mpi_communicator
SolverControl & control() const
@ smallest_imaginary_part
void solve(const MatrixType1 &A, const MatrixType2 &B, const INVERSE &inverse, std::vector< std::complex< double > > &eigenvalues, std::vector< VectorType > &eigenvectors, const unsigned int n_eigenvalues)
PArpackSolver(SolverControl &control, const MPI_Comm mpi_communicator, const AdditionalData &data=AdditionalData())
std::vector< types::global_dof_index > local_indices
std::vector< double > workl
std::vector< int > select
std::size_t memory_consumption() const
MPI_Fint mpi_communicator_fortran
void set_shift(const std::complex< double > sigma)
void reinit(const IndexSet &locally_owned_dofs)
std::vector< double > workd
const AdditionalData additional_data
bool initial_vector_provided
void internal_reinit(const IndexSet &locally_owned_dofs)
std::vector< double > resid
std::vector< double > workev
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & PArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
static ::ExceptionBase & PArpackExcInvalidEigenvalueSize(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & PArpackExcNoShifts(int arg1)
static ::ExceptionBase & PArpackExcInvalidEigenvectorSize(int arg1, int arg2)
static ::ExceptionBase & PArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & PArpackExcInfoPdneupd(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & PArpackExcIdo(int arg1)
static ::ExceptionBase & PArpackExcConvergedEigenvectors(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & PArpackExcInfoMaxIt(int arg1)
static ::ExceptionBase & PArpackExcMode(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & PArpackExcInfoPdnaupd(int arg1)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & PArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::vector< index_type > data
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
unsigned int global_dof_index
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)
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)
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)
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)
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false, const int mode=3)
const unsigned int number_of_arnoldi_vectors
const WhichEigenvalues eigenvalue_of_interest
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)