15#ifndef dealii_parpack_solver_h
16#define dealii_parpack_solver_h
30#ifdef DEAL_II_ARPACK_WITH_PARPACK
209template <
typename VectorType>
314 const std::vector<IndexSet> &partitioning);
320 reinit(
const VectorType &distributed_vector);
337 set_shift(
const std::complex<double> sigma);
348 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
350 solve(
const MatrixType1 &A,
351 const MatrixType2 &B,
352 const INVERSE &inverse,
355 const unsigned int n_eigenvalues);
360 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
362 solve(
const MatrixType1 &A,
363 const MatrixType2 &B,
364 const INVERSE &inverse,
367 const unsigned int n_eigenvalues);
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>
593 src.memory_consumption() +
dst.memory_consumption() +
594 tmp.memory_consumption() +
601template <
typename VectorType>
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>
669template <
typename VectorType>
683template <
typename VectorType>
723template <
typename VectorType>
737template <
typename VectorType>
744 src.reinit(distributed_vector);
745 dst.reinit(distributed_vector);
746 tmp.reinit(distributed_vector);
751template <
typename VectorType>
754 const std::vector<IndexSet> &partitioning)
766template <
typename VectorType>
767template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
770 const MatrixType2 &B,
771 const INVERSE &inverse,
774 const unsigned int n_eigenvalues)
776 std::vector<VectorType *> eigenvectors_ptr(
eigenvectors.size());
784template <
typename VectorType>
785template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
788 const MatrixType2 &mass_matrix,
789 const INVERSE &inverse,
792 const unsigned int n_eigenvalues)
832 bmat[0] = (mode == 1) ?
'I' :
'G';
849 std::strcpy(which,
"LA");
852 std::strcpy(which,
"SA");
855 std::strcpy(which,
"LM");
858 std::strcpy(which,
"SM");
861 std::strcpy(which,
"LR");
864 std::strcpy(which,
"SR");
867 std::strcpy(which,
"LI");
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);
912 int nev = n_eigenvalues;
913 int n_inside_arpack =
nloc;
964 const int shift_x = ipntr[0] - 1;
965 const int shift_y = ipntr[1] - 1;
976 if ((ido == -1) || (ido == 1 && mode < 3))
985 mass_matrix.vmult(
tmp,
src);
991 system_matrix.vmult(
tmp,
src);
995 workd.data() + shift_x);
1000 system_matrix.vmult(
dst,
src);
1005 else if (ido == 1 && mode >= 3)
1009 const int shift_b_x = ipntr[2] - 1;
1036 mass_matrix.vmult(
dst,
src);
1047 workd.data() + shift_y);
1055 char howmany[4] =
"All";
1057 std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1058 std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1066 eigenvalues_real.data(),
1090 eigenvalues_real.data(),
1091 eigenvalues_im.data(),
1126 for (
int i = 0; i < nev; ++i)
1135 for (
size_type i = 0; i < n_eigenvalues; ++i)
1137 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1140 AssertThrow(iparam[4] >=
static_cast<int>(n_eigenvalues),
1157template <
typename VectorType>
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
types::global_dof_index size_type
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::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 > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)