15#ifndef dealii_arpack_solver_h
16#define dealii_arpack_solver_h
27#ifdef DEAL_II_WITH_ARPACK
270 template <
typename VectorType>
280 set_shift(
const std::complex<double> sigma);
331 template <
typename VectorType,
332 typename MatrixType1,
333 typename MatrixType2,
336 solve(
const MatrixType1 &A,
337 const MatrixType2 &B,
338 const INVERSE &inverse,
341 const unsigned int n_eigenvalues = 0);
379 <<
"Number of wanted eigenvalues " << arg1
380 <<
" is larger that the size of the matrix " << arg2);
385 <<
"Number of wanted eigenvalues " << arg1
386 <<
" is larger that the size of eigenvectors " << arg2);
392 <<
"To store the real and complex parts of " << arg1
393 <<
" eigenvectors in real-valued vectors, their size (currently set to "
394 << arg2 <<
") should be greater than or equal to " << arg1 + 1);
399 <<
"Number of wanted eigenvalues " << arg1
400 <<
" is larger that the size of eigenvalues " << arg2);
405 <<
"Number of Arnoldi vectors " << arg1
406 <<
" is larger that the size of the matrix " << arg2);
411 <<
"Number of Arnoldi vectors " << arg1
412 <<
" is too small to obtain " << arg2 <<
" eigenvalues");
416 <<
"This ido " << arg1
417 <<
" is not supported. Check documentation of ARPACK");
421 <<
"This mode " << arg1
422 <<
" is not supported. Check documentation of ARPACK");
426 <<
"Error with dsaupd, info " << arg1
427 <<
". Check documentation of ARPACK");
431 <<
"Error with dnaupd, info " << arg1
432 <<
". Check documentation of ARPACK");
436 <<
"Error with dseupd, info " << arg1
437 <<
". Check documentation of ARPACK");
441 <<
"Error with dneupd, info " << arg1
442 <<
". Check documentation of ARPACK");
446 <<
"Maximum number " << arg1 <<
" of iterations reached.");
449 "No shifts could be applied during implicit"
450 " Arnoldi update, try increasing the number of"
451 " Arnoldi vectors.");
456 const unsigned int number_of_arnoldi_vectors,
458 const bool symmetric)
459 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
460 , eigenvalue_of_interest(eigenvalue_of_interest)
461 , symmetric(symmetric)
469 "'largest real part' can only be used for non-symmetric problems!"));
473 "'smallest real part' can only be used for non-symmetric problems!"));
477 "'largest imaginary part' can only be used for non-symmetric problems!"));
481 "'smallest imaginary part' can only be used for non-symmetric problems!"));
489 "'largest algebraic part' can only be used for symmetric problems!"));
493 "'smallest algebraic part' can only be used for symmetric problems!"));
496 "'both ends' can only be used for symmetric problems!"));
521template <
typename VectorType>
526 resid.resize(vec.size());
527 for (
size_type i = 0; i < vec.size(); ++i)
532template <
typename VectorType,
533 typename MatrixType1,
534 typename MatrixType2,
538 const MatrixType2 &mass_matrix,
539 const INVERSE &inverse,
542 const unsigned int n_eigenvalues)
548 const unsigned int nev_const =
549 (n_eigenvalues == 0) ?
eigenvalues.size() : n_eigenvalues;
551 unsigned int nev = nev_const;
597 std::strcpy(which,
"LA");
600 std::strcpy(which,
"SA");
603 std::strcpy(which,
"LM");
606 std::strcpy(which,
"SM");
609 std::strcpy(which,
"LR");
612 std::strcpy(which,
"SR");
615 std::strcpy(which,
"LI");
618 std::strcpy(which,
"SI");
621 std::strcpy(which,
"BE");
637 std::vector<double> v(ldv * ncv, 0.0);
640 std::vector<int> iparam(11, 0);
651 std::vector<int> ipntr(14, 0);
654 std::vector<double> workd(3 * n, 0.);
657 std::vector<double> workl(lworkl, 0.);
711 VectorType src, dst, tmp;
717 for (
size_type i = 0; i < src.size(); ++i)
718 src(i) = workd[ipntr[0] - 1 + i];
721 mass_matrix.vmult(tmp, src);
723 inverse.vmult(dst, tmp);
725 for (
size_type i = 0; i < dst.size(); ++i)
726 workd[ipntr[1] - 1 + i] = dst(i);
732 VectorType src, dst, tmp, tmp2;
738 for (
size_type i = 0; i < src.size(); ++i)
740 src(i) = workd[ipntr[2] - 1 + i];
741 tmp(i) = workd[ipntr[0] - 1 + i];
744 inverse.vmult(dst, src);
746 for (
size_type i = 0; i < dst.size(); ++i)
747 workd[ipntr[1] - 1 + i] = dst(i);
757 for (
size_type i = 0; i < src.size(); ++i)
758 src(i) = workd[ipntr[0] - 1 + i];
761 mass_matrix.vmult(dst, src);
763 for (
size_type i = 0; i < dst.size(); ++i)
764 workd[ipntr[1] - 1 + i] = dst(i);
800 std::vector<int> select(ncv, 1);
804 std::vector<double> eigenvalues_real(nev + 1, 0.);
805 std::vector<double> eigenvalues_im(nev + 1, 0.);
810 std::vector<double> z(ldz * nev, 0.);
814 eigenvalues_real.data(),
836 std::vector<double> workev(3 * ncv, 0.);
840 eigenvalues_real.data(),
841 eigenvalues_im.data(),
882 for (
unsigned int i = 0; i < nev; ++i)
883 for (
unsigned int j = 0; j < n; ++j)
886 for (
unsigned int i = 0; i < nev_const; ++i)
888 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)
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)