15#ifndef dealii_arpack_solver_h
16#define dealii_arpack_solver_h
30#ifdef DEAL_II_WITH_ARPACK
273 template <
typename VectorType>
334 template <
typename VectorType,
382 <<
"Number of wanted eigenvalues " <<
arg1
383 <<
" is larger that the size of the matrix " <<
arg2);
388 <<
"Number of wanted eigenvalues " <<
arg1
389 <<
" is larger that the size of eigenvectors " <<
arg2);
395 <<
"To store the real and complex parts of " <<
arg1
396 <<
" eigenvectors in real-valued vectors, their size (currently set to "
397 <<
arg2 <<
") should be greater than or equal to " <<
arg1 + 1);
402 <<
"Number of wanted eigenvalues " <<
arg1
403 <<
" is larger that the size of eigenvalues " <<
arg2);
408 <<
"Number of Arnoldi vectors " <<
arg1
409 <<
" is larger that the size of the matrix " <<
arg2);
414 <<
"Number of Arnoldi vectors " <<
arg1
415 <<
" is too small to obtain " <<
arg2 <<
" eigenvalues");
419 <<
"This ido " <<
arg1
420 <<
" is not supported. Check documentation of ARPACK");
424 <<
"This mode " <<
arg1
425 <<
" is not supported. Check documentation of ARPACK");
429 <<
"Error with dsaupd, info " <<
arg1
430 <<
". Check documentation of ARPACK");
434 <<
"Error with dnaupd, info " <<
arg1
435 <<
". Check documentation of ARPACK");
439 <<
"Error with dseupd, info " <<
arg1
440 <<
". Check documentation of ARPACK");
444 <<
"Error with dneupd, info " <<
arg1
445 <<
". Check documentation of ARPACK");
449 <<
"Maximum number " <<
arg1 <<
" of iterations reached.");
452 "No shifts could be applied during implicit"
453 " Arnoldi update, try increasing the number of"
454 " Arnoldi vectors.");
459 const unsigned int number_of_arnoldi_vectors,
461 const bool symmetric)
462 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
463 , eigenvalue_of_interest(eigenvalue_of_interest)
464 , symmetric(symmetric)
472 "'largest real part' can only be used for non-symmetric problems!"));
476 "'smallest real part' can only be used for non-symmetric problems!"));
480 "'largest imaginary part' can only be used for non-symmetric problems!"));
484 "'smallest imaginary part' can only be used for non-symmetric problems!"));
492 "'largest algebraic part' can only be used for symmetric problems!"));
496 "'smallest algebraic part' can only be used for symmetric problems!"));
499 "'both ends' can only be used for symmetric problems!"));
524template <
typename VectorType>
535template <
typename VectorType,
600 std::strcpy(
which,
"LA");
603 std::strcpy(
which,
"SA");
606 std::strcpy(
which,
"LM");
609 std::strcpy(
which,
"SM");
612 std::strcpy(
which,
"LR");
615 std::strcpy(
which,
"SR");
618 std::strcpy(
which,
"LI");
621 std::strcpy(
which,
"SI");
624 std::strcpy(
which,
"BE");
640 std::vector<double> v(ldv * ncv, 0.0);
643 std::vector<int>
iparam(11, 0);
654 std::vector<int>
ipntr(14, 0);
657 std::vector<double> workd(3 * n, 0.);
659 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
660 std::vector<double> workl(lworkl, 0.);
714 VectorType src, dst, tmp;
720 for (
size_type i = 0; i < src.size(); ++i)
721 src(i) = workd[
ipntr[0] - 1 + i];
724 mass_matrix.vmult(tmp, src);
726 inverse.vmult(dst, tmp);
728 for (
size_type i = 0; i < dst.size(); ++i)
729 workd[
ipntr[1] - 1 + i] = dst(i);
735 VectorType src, dst, tmp,
tmp2;
741 for (
size_type i = 0; i < src.size(); ++i)
743 src(i) = workd[
ipntr[2] - 1 + i];
744 tmp(i) = workd[
ipntr[0] - 1 + i];
747 inverse.vmult(dst, src);
749 for (
size_type i = 0; i < dst.size(); ++i)
750 workd[
ipntr[1] - 1 + i] = dst(i);
760 for (
size_type i = 0; i < src.size(); ++i)
761 src(i) = workd[
ipntr[0] - 1 + i];
764 mass_matrix.vmult(dst, src);
766 for (
size_type i = 0; i < dst.size(); ++i)
767 workd[
ipntr[1] - 1 + i] = dst(i);
803 std::vector<int> select(ncv, 1);
813 std::vector<double> z(ldz *
nev, 0.);
839 std::vector<double> workev(3 * ncv, 0.);
885 for (
unsigned int i = 0; i <
nev; ++i)
886 for (
unsigned int j = 0;
j < n; ++
j)
889 for (
unsigned int i = 0; i <
nev_const; ++i)
void dseupd_(int *rvec, char *howmany, int *select, double *d, double *z, int *ldz, double *sigmar, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dsaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dnaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, const double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dneupd_(int *rvec, char *howmany, int *select, double *d, double *di, double *z, int *ldz, double *sigmar, double *sigmai, double *workev, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
SolverControl & control() const
void set_initial_vector(const VectorType &vec)
bool initial_vector_provided
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
std::vector< double > resid
SolverControl & solver_control
@ 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=0)
void set_shift(const std::complex< double > sigma)
const AdditionalData additional_data
virtual State check(const unsigned int step, const double check_value)
unsigned int max_steps() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ArpackExcInvalidEigenvalueSize(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackInfoMaxIt(int arg1)
static ::ExceptionBase & ArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & ArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackNoShifts()
static ::ExceptionBase & ArpackExcInvalidEigenvectorSize(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackMode(int arg1)
static ::ExceptionBase & ArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ArpackExcArpackInfodneupd(int arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
static ::ExceptionBase & ArpackExcArpackInfodsaupd(int arg1)
static ::ExceptionBase & ArpackExcArpackInfodnaupd(int arg1)
#define DeclExceptionMsg(Exception, defaulttext)
static ::ExceptionBase & ArpackExcArpackInfodseupd(int arg1)
static ::ExceptionBase & ArpackExcArpackIdo(int arg1)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
unsigned int global_dof_index
const WhichEigenvalues eigenvalue_of_interest
const unsigned int number_of_arnoldi_vectors
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false)
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)