16 #ifndef dealii_parpack_solver_h 17 #define dealii_parpack_solver_h 19 #include <deal.II/base/config.h> 21 #include <deal.II/base/index_set.h> 22 #include <deal.II/base/memory_consumption.h> 23 #include <deal.II/base/smartpointer.h> 25 #include <deal.II/lac/solver_control.h> 26 #include <deal.II/lac/vector_operation.h> 31 #ifdef DEAL_II_ARPACK_WITH_PARPACK 33 DEAL_II_NAMESPACE_OPEN
39 pdnaupd_(MPI_Fint *comm,
59 pdsaupd_(MPI_Fint *comm,
79 pdneupd_(MPI_Fint *comm,
108 pdseupd_(MPI_Fint *comm,
212 template <
typename VectorType>
281 template <
typename MatrixType>
288 Shift(
const MatrixType &A,
const MatrixType &B,
const double sigma)
298 vmult(VectorType &dst,
const VectorType &
src)
const 302 A.vmult_add(dst, src);
313 A.Tvmult_add(dst, src);
328 const unsigned int number_of_arnoldi_vectors;
330 const bool symmetric;
333 const unsigned int number_of_arnoldi_vectors = 15,
335 const bool symmetric =
false,
366 const std::vector<IndexSet> &partitioning);
372 reinit(
const VectorType &distributed_vector);
389 set_shift(
const std::complex<double> sigma);
400 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
402 solve(
const MatrixType1 & A,
403 const MatrixType2 & B,
404 const INVERSE & inverse,
407 const unsigned int n_eigenvalues);
412 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
414 solve(
const MatrixType1 & A,
415 const MatrixType2 & B,
416 const INVERSE & inverse,
418 std::vector<VectorType *> & eigenvectors,
419 const unsigned int n_eigenvalues);
488 std::vector<double>
v;
511 std::vector<double>
z;
564 << arg1 <<
" eigenpairs were requested, but only " << arg2
570 <<
"Number of wanted eigenvalues " << arg1
571 <<
" is larger that the size of the matrix " << arg2);
576 <<
"Number of wanted eigenvalues " << arg1
577 <<
" is larger that the size of eigenvectors " << arg2);
583 <<
"To store the real and complex parts of " << arg1
584 <<
" eigenvectors in real-valued vectors, their size (currently set to " 585 << arg2 <<
") should be greater than or equal to " << arg1 + 1);
590 <<
"Number of wanted eigenvalues " << arg1
591 <<
" is larger that the size of eigenvalues " << arg2);
596 <<
"Number of Arnoldi vectors " << arg1
597 <<
" is larger that the size of the matrix " << arg2);
602 <<
"Number of Arnoldi vectors " << arg1
603 <<
" is too small to obtain " << arg2 <<
" eigenvalues");
607 <<
"This ido " << arg1
608 <<
" is not supported. Check documentation of ARPACK");
612 <<
"This mode " << arg1
613 <<
" is not supported. Check documentation of ARPACK");
617 <<
"Error with Pdnaupd, info " << arg1
618 <<
". Check documentation of ARPACK");
622 <<
"Error with Pdneupd, info " << arg1
623 <<
". Check documentation of ARPACK");
627 <<
"Maximum number " << arg1 <<
" of iterations reached.");
631 <<
"No shifts could be applied during implicit" 632 <<
" Arnoldi update, try increasing the number of" 633 <<
" Arnoldi vectors.");
638 template <
typename VectorType>
645 src.memory_consumption() + dst.memory_consumption() +
646 tmp.memory_consumption() +
653 template <
typename VectorType>
655 const unsigned int number_of_arnoldi_vectors,
657 const bool symmetric,
659 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
660 , eigenvalue_of_interest(eigenvalue_of_interest)
661 , symmetric(symmetric)
670 "'largest real part' can only be used for non-symmetric problems!"));
674 "'smallest real part' can only be used for non-symmetric problems!"));
678 "'largest imaginary part' can only be used for non-symmetric problems!"));
682 "'smallest imaginary part' can only be used for non-symmetric problems!"));
684 Assert(mode >= 1 && mode <= 3,
685 ExcMessage(
"Currently, only modes 1, 2 and 3 are supported."));
690 template <
typename VectorType>
696 , mpi_communicator(mpi_communicator)
711 template <
typename VectorType>
721 template <
typename VectorType>
735 template <
typename VectorType>
750 v.resize(ldv *
ncv, 0.0);
762 z.resize(
ldz * ncv, 0.);
775 template <
typename VectorType>
789 template <
typename VectorType>
796 src.reinit(distributed_vector);
797 dst.reinit(distributed_vector);
798 tmp.reinit(distributed_vector);
803 template <
typename VectorType>
806 const std::vector<IndexSet> &partitioning)
818 template <
typename VectorType>
819 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
822 const MatrixType2 & B,
823 const INVERSE & inverse,
826 const unsigned int n_eigenvalues)
828 std::vector<VectorType *> eigenvectors_ptr(eigenvectors.size());
829 for (
unsigned int i = 0; i < eigenvectors.size(); ++i)
830 eigenvectors_ptr[i] = &eigenvectors[i];
836 template <
typename VectorType>
837 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
840 const MatrixType2 &mass_matrix,
841 const INVERSE & inverse,
844 const unsigned int n_eigenvalues)
848 Assert(n_eigenvalues <= eigenvectors.size(),
850 eigenvectors.size()));
853 Assert(n_eigenvalues + 1 <= eigenvectors.size(),
855 eigenvectors.size()));
863 Assert(n_eigenvalues < eigenvectors[0]->size(),
865 eigenvectors[0]->size()));
884 bmat[0] = (mode == 1) ?
'I' :
'G';
901 std::strcpy(which,
"LA");
904 std::strcpy(which,
"SA");
907 std::strcpy(which,
"LM");
910 std::strcpy(which,
"SM");
913 std::strcpy(which,
"LR");
916 std::strcpy(which,
"SR");
919 std::strcpy(which,
"LI");
922 std::strcpy(which,
"SI");
925 std::strcpy(which,
"BE");
933 std::vector<int> iparam(11, 0);
953 std::vector<int> ipntr(14, 0);
964 int nev = n_eigenvalues;
965 int n_inside_arpack =
nloc;
1016 const int shift_x = ipntr[0] - 1;
1017 const int shift_y = ipntr[1] - 1;
1026 if ((ido == -1) || (ido == 1 && mode < 3))
1035 mass_matrix.vmult(tmp, src);
1036 inverse.vmult(dst, tmp);
1041 system_matrix.vmult(tmp, src);
1045 workd.data() + shift_x);
1046 inverse.vmult(dst, tmp);
1050 system_matrix.vmult(dst, src);
1055 else if (ido == 1 && mode >= 3)
1059 const int shift_b_x = ipntr[2] - 1;
1070 inverse.vmult(dst, src);
1086 mass_matrix.vmult(dst, src);
1097 workd.data() + shift_y);
1105 char howmany[4] =
"All";
1107 std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1108 std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1116 eigenvalues_real.data(),
1140 eigenvalues_real.data(),
1141 eigenvalues_im.data(),
1176 for (
int i = 0; i < nev; ++i)
1178 (*eigenvectors[i]) = 0.0;
1185 for (
size_type i = 0; i < n_eigenvalues; ++i)
1187 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1206 template <
typename VectorType>
1213 DEAL_II_NAMESPACE_CLOSE
types::global_dof_index size_type
#define DeclException2(Exception2, type1, type2, outsequence)
static::ExceptionBase & PArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
static::ExceptionBase & PArpackExcInvalidEigenvectorSize(int arg1, int arg2)
void vmult(VectorType &dst, const VectorType &src) const
virtual State check(const unsigned int step, const double check_value)
std::vector< double > workev
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)
static::ExceptionBase & PArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
SolverControl & control() const
static::ExceptionBase & PArpackExcMode(int arg1)
size_type n_elements() const
PArpackSolver(SolverControl &control, const MPI_Comm &mpi_communicator, const AdditionalData &data=AdditionalData())
MPI_Fint mpi_communicator_fortran
#define AssertThrow(cond, exc)
const AdditionalData additional_data
unsigned long long int global_dof_index
void set_shift(const std::complex< double > sigma)
static::ExceptionBase & ExcMessage(std::string arg1)
static::ExceptionBase & PArpackExcNoShifts(int arg1)
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)
#define DeclException1(Exception1, type1, outsequence)
bool initial_vector_provided
void Tvmult(VectorType &dst, const VectorType &src) const
unsigned int max_steps() const
void fill_index_vector(std::vector< size_type > &indices) const
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< int > select
static::ExceptionBase & PArpackExcConvergedEigenvectors(int arg1, int arg2)
static::ExceptionBase & PArpackExcInfoPdnaupd(int arg1)
static::ExceptionBase & PArpackExcInfoMaxIt(int arg1)
MPI_Comm mpi_communicator
static::ExceptionBase & PArpackExcInvalidEigenvalueSize(int arg1, int arg2)
std::size_t memory_consumption() const
std::vector< double > resid
static::ExceptionBase & PArpackExcIdo(int arg1)
void internal_reinit(const IndexSet &locally_owned_dofs)
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
static::ExceptionBase & ExcNotImplemented()
Shift(const MatrixType &A, const MatrixType &B, const double sigma)
void reinit(const IndexSet &locally_owned_dofs)
static::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
std::vector< double > workl
static::ExceptionBase & PArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
std::vector< types::global_dof_index > local_indices
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::vector< double > workd
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & PArpackExcInfoPdneupd(int arg1)