16#ifndef dealii_arpack_solver_h
17#define dealii_arpack_solver_h
28#ifdef DEAL_II_WITH_ARPACK
271 template <
typename VectorType>
281 set_shift(
const std::complex<double> sigma);
332 template <
typename VectorType,
333 typename MatrixType1,
334 typename MatrixType2,
337 solve(
const MatrixType1 & A,
338 const MatrixType2 & B,
339 const INVERSE & inverse,
342 const unsigned int n_eigenvalues = 0);
380 <<
"Number of wanted eigenvalues " << arg1
381 <<
" is larger that the size of the matrix " << arg2);
386 <<
"Number of wanted eigenvalues " << arg1
387 <<
" is larger that the size of eigenvectors " << arg2);
393 <<
"To store the real and complex parts of " << arg1
394 <<
" eigenvectors in real-valued vectors, their size (currently set to "
395 << arg2 <<
") should be greater than or equal to " << arg1 + 1);
400 <<
"Number of wanted eigenvalues " << arg1
401 <<
" is larger that the size of eigenvalues " << arg2);
406 <<
"Number of Arnoldi vectors " << arg1
407 <<
" is larger that the size of the matrix " << arg2);
412 <<
"Number of Arnoldi vectors " << arg1
413 <<
" is too small to obtain " << arg2 <<
" eigenvalues");
417 <<
"This ido " << arg1
418 <<
" is not supported. Check documentation of ARPACK");
422 <<
"This mode " << arg1
423 <<
" is not supported. Check documentation of ARPACK");
427 <<
"Error with dsaupd, info " << arg1
428 <<
". Check documentation of ARPACK");
432 <<
"Error with dnaupd, info " << arg1
433 <<
". Check documentation of ARPACK");
437 <<
"Error with dseupd, info " << arg1
438 <<
". Check documentation of ARPACK");
442 <<
"Error with dneupd, info " << arg1
443 <<
". Check documentation of ARPACK");
447 <<
"Maximum number " << arg1 <<
" of iterations reached.");
450 "No shifts could be applied during implicit"
451 " Arnoldi update, try increasing the number of"
452 " Arnoldi vectors.");
457 const unsigned int number_of_arnoldi_vectors,
459 const bool symmetric)
460 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
461 , eigenvalue_of_interest(eigenvalue_of_interest)
462 , symmetric(symmetric)
470 "'largest real part' can only be used for non-symmetric problems!"));
474 "'smallest real part' can only be used for non-symmetric problems!"));
478 "'largest imaginary part' can only be used for non-symmetric problems!"));
482 "'smallest imaginary part' can only be used for non-symmetric problems!"));
490 "'largest algebraic part' can only be used for symmetric problems!"));
494 "'smallest algebraic part' can only be used for symmetric problems!"));
497 "'both ends' can only be used for symmetric problems!"));
522template <
typename VectorType>
527 resid.resize(vec.size());
528 for (
size_type i = 0; i < vec.size(); ++i)
533template <
typename VectorType,
534 typename MatrixType1,
535 typename MatrixType2,
539 const MatrixType2 & mass_matrix,
540 const INVERSE & inverse,
543 const unsigned int n_eigenvalues)
549 const unsigned int nev_const =
550 (n_eigenvalues == 0) ?
eigenvalues.size() : n_eigenvalues;
552 unsigned int nev = nev_const;
598 std::strcpy(which,
"LA");
601 std::strcpy(which,
"SA");
604 std::strcpy(which,
"LM");
607 std::strcpy(which,
"SM");
610 std::strcpy(which,
"LR");
613 std::strcpy(which,
"SR");
616 std::strcpy(which,
"LI");
619 std::strcpy(which,
"SI");
622 std::strcpy(which,
"BE");
638 std::vector<double> v(ldv * ncv, 0.0);
641 std::vector<int> iparam(11, 0);
652 std::vector<int> ipntr(14, 0);
655 std::vector<double> workd(3 * n, 0.);
658 std::vector<double> workl(lworkl, 0.);
712 VectorType src, dst, tmp;
718 for (
size_type i = 0; i < src.size(); ++i)
719 src(i) = workd[ipntr[0] - 1 + i];
722 mass_matrix.vmult(tmp, src);
724 inverse.vmult(dst, tmp);
726 for (
size_type i = 0; i < dst.size(); ++i)
727 workd[ipntr[1] - 1 + i] = dst(i);
733 VectorType src, dst, tmp, tmp2;
739 for (
size_type i = 0; i < src.size(); ++i)
741 src(i) = workd[ipntr[2] - 1 + i];
742 tmp(i) = workd[ipntr[0] - 1 + i];
745 inverse.vmult(dst, src);
747 for (
size_type i = 0; i < dst.size(); ++i)
748 workd[ipntr[1] - 1 + i] = dst(i);
758 for (
size_type i = 0; i < src.size(); ++i)
759 src(i) = workd[ipntr[0] - 1 + i];
762 mass_matrix.vmult(dst, src);
764 for (
size_type i = 0; i < dst.size(); ++i)
765 workd[ipntr[1] - 1 + i] = dst(i);
801 std::vector<int> select(ncv, 1);
805 std::vector<double> eigenvalues_real(nev + 1, 0.);
806 std::vector<double> eigenvalues_im(nev + 1, 0.);
811 std::vector<double> z(ldz * nev, 0.);
815 eigenvalues_real.data(),
837 std::vector<double> workev(3 * ncv, 0.);
841 eigenvalues_real.data(),
842 eigenvalues_im.data(),
883 for (
unsigned int i = 0; i < nev; ++i)
884 for (
unsigned int j = 0; j < n; ++j)
887 for (
unsigned int i = 0; i < nev_const; ++i)
889 std::complex<double>(eigenvalues_real[i], eigenvalues_im[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)
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 > >, std::integral_constant< int, dim >::value > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)