16 #ifndef dealii__parpack_solver_h 17 #define dealii__parpack_solver_h 19 #include <deal.II/base/config.h> 20 #include <deal.II/base/smartpointer.h> 21 #include <deal.II/base/memory_consumption.h> 22 #include <deal.II/lac/solver_control.h> 23 #include <deal.II/lac/vector.h> 24 #include <deal.II/base/index_set.h> 29 #ifdef DEAL_II_ARPACK_WITH_PARPACK 31 DEAL_II_NAMESPACE_OPEN
36 void pdnaupd_(MPI_Fint *comm,
int *ido,
char *bmat,
int *n,
char *which,
37 int *nev,
double *tol,
double *resid,
int *ncv,
38 double *v,
int *nloc,
int *iparam,
int *ipntr,
39 double *workd,
double *workl,
int *lworkl,
43 void pdsaupd_(MPI_Fint *comm,
int *ido,
char *bmat,
int *n,
char *which,
44 int *nev,
double *tol,
double *resid,
int *ncv,
45 double *v,
int *nloc,
int *iparam,
int *ipntr,
46 double *workd,
double *workl,
int *lworkl,
50 void pdneupd_(MPI_Fint *comm,
int *rvec,
char *howmany,
int *select,
double *d,
51 double *di,
double *z,
int *ldz,
double *sigmar,
52 double *sigmai,
double *workev,
char *bmat,
int *n,
char *which,
53 int *nev,
double *tol,
double *resid,
int *ncv,
54 double *v,
int *nloc,
int *iparam,
int *ipntr,
55 double *workd,
double *workl,
int *lworkl,
int *info);
58 void pdseupd_(MPI_Fint *comm,
int *rvec,
char *howmany,
int *select,
double *d,
59 double *z,
int *ldz,
double *sigmar,
60 char *bmat,
int *n,
char *which,
61 int *nev,
double *tol,
double *resid,
int *ncv,
62 double *v,
int *nloc,
int *iparam,
int *ipntr,
63 double *workd,
double *workl,
int *lworkl,
int *info);
143 template <
typename VectorType>
207 template <
typename MatrixType>
227 void vmult (VectorType &dst,
const VectorType &
src)
const 231 A.vmult_add(dst,src);
237 void Tvmult (VectorType &dst,
const VectorType &
src)
const 241 A.Tvmult_add(dst,src);
256 const unsigned int number_of_arnoldi_vectors;
258 const bool symmetric;
260 const unsigned int number_of_arnoldi_vectors = 15,
262 const bool symmetric =
false);
289 const std::vector<IndexSet> &partitioning);
294 void reinit(
const VectorType &distributed_vector);
306 void set_shift(
const std::complex<double> sigma);
313 template <
typename MatrixType1,
314 typename MatrixType2,
typename INVERSE>
316 (
const MatrixType1 &A,
317 const MatrixType2 &B,
318 const INVERSE &inverse,
319 std::vector<std::complex<double> > &eigenvalues,
320 std::vector<VectorType> &eigenvectors,
321 const unsigned int n_eigenvalues);
323 std::size_t memory_consumption()
const;
387 std::vector<double>
v;
410 std::vector<double>
z;
461 << arg1 <<
" eigenpairs were requested, but only " 462 << arg2 <<
" converged");
465 <<
"Number of wanted eigenvalues " << arg1
466 <<
" is larger that the size of the matrix " << arg2);
469 <<
"Number of wanted eigenvalues " << arg1
470 <<
" is larger that the size of eigenvectors " << arg2);
473 <<
"To store the real and complex parts of " << arg1
474 <<
" eigenvectors in real-valued vectors, their size (currently set to " << arg2
475 <<
") should be greater than or equal to " << arg1+1);
478 <<
"Number of wanted eigenvalues " << arg1
479 <<
" is larger that the size of eigenvalues " << arg2);
482 <<
"Number of Arnoldi vectors " << arg1
483 <<
" is larger that the size of the matrix " << arg2);
486 <<
"Number of Arnoldi vectors " << arg1
487 <<
" is too small to obtain " << arg2
491 <<
" is not supported. Check documentation of ARPACK");
494 <<
" is not supported. Check documentation of ARPACK");
497 <<
"Error with Pdnaupd, info " << arg1
498 <<
". Check documentation of ARPACK");
501 <<
"Error with Pdneupd, info " << arg1
502 <<
". Check documentation of ARPACK");
505 <<
"Maximum number " << arg1
506 <<
" of iterations reached.");
509 <<
"No shifts could be applied during implicit" 510 <<
" Arnoldi update, try increasing the number of" 511 <<
" Arnoldi vectors.");
514 template <
typename VectorType>
525 src.memory_consumption() +
526 dst.memory_consumption() +
527 tmp.memory_consumption() +
531 template <
typename VectorType>
535 const bool symmetric)
537 number_of_arnoldi_vectors(number_of_arnoldi_vectors),
538 eigenvalue_of_interest(eigenvalue_of_interest),
545 ExcMessage(
"'largest real part' can only be used for non-symmetric problems!"));
547 ExcMessage(
"'smallest real part' can only be used for non-symmetric problems!"));
549 ExcMessage(
"'largest imaginary part' can only be used for non-symmetric problems!"));
551 ExcMessage(
"'smallest imaginary part' can only be used for non-symmetric problems!"));
555 template <
typename VectorType>
562 mpi_communicator( mpi_communicator ),
575 template <
typename VectorType>
582 template <
typename VectorType>
595 template <
typename VectorType>
610 v.resize (ldv*
ncv, 0.0);
624 z.resize (
ldz*ncv, 0.);
636 template <
typename VectorType>
647 template <
typename VectorType>
653 src.reinit (distributed_vector);
654 dst.reinit (distributed_vector);
655 tmp.reinit (distributed_vector);
659 template <
typename VectorType>
661 const std::vector<IndexSet> &partitioning)
671 template <
typename VectorType>
672 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
674 (
const MatrixType1 &,
675 const MatrixType2 &mass_matrix,
676 const INVERSE &inverse,
677 std::vector<std::complex<double> > &eigenvalues,
678 std::vector<VectorType> &eigenvectors,
679 const unsigned int n_eigenvalues)
684 Assert (n_eigenvalues <= eigenvectors.size(),
688 Assert (n_eigenvalues+1 <= eigenvectors.size(),
691 Assert (n_eigenvalues <= eigenvalues.size(),
697 Assert (n_eigenvalues < eigenvectors[0].size(),
735 std::strcpy (which,
"LA");
738 std::strcpy (which,
"SA");
741 std::strcpy (which,
"LM");
744 std::strcpy (which,
"SM");
747 std::strcpy (which,
"LR");
750 std::strcpy (which,
"SR");
753 std::strcpy (which,
"LI");
756 std::strcpy (which,
"SI");
759 std::strcpy (which,
"BE");
767 std::vector<int> iparam (11, 0);
786 std::vector<int> ipntr (14, 0);
797 int nev = n_eigenvalues;
798 int n_inside_arpack =
nloc;
828 const int shift_x = ipntr[0]-1;
829 const int shift_y = ipntr[1]-1;
842 mass_matrix.vmult(tmp, src);
844 inverse.vmult(dst,tmp);
861 const int shift_x = ipntr[0]-1;
862 const int shift_y = ipntr[1]-1;
863 const int shift_b_x = ipntr[2]-1;
877 &
workd[0]+shift_b_x );
888 inverse.vmult(dst,src);
904 const int shift_x = ipntr[0]-1;
905 const int shift_y = ipntr[1]-1;
918 mass_matrix.vmult(dst, src);
951 char howmany[4] =
"All";
953 std::vector<double> eigenvalues_real (n_eigenvalues+1, 0.);
954 std::vector<double> eigenvalues_im (n_eigenvalues+1, 0.);
960 bmat, &n_inside_arpack, which, &nev, &tol,
966 &
workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
983 for (
int i=0; i<nev; ++i)
985 eigenvectors[i] = 0.0;
988 eigenvectors[i].add (
nloc,
994 for (
size_type i=0; i<n_eigenvalues; ++i)
995 eigenvalues[i] = std::complex<double> (eigenvalues_real[i],
1019 template <
typename VectorType>
1025 DEAL_II_NAMESPACE_CLOSE
#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
static::ExceptionBase & PArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
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
void set_shift(const std::complex< double > sigma)
static::ExceptionBase & ExcMessage(std::string arg1)
static::ExceptionBase & PArpackExcNoShifts(int arg1)
#define DeclException1(Exception1, type1, outsequence)
bool initial_vector_provided
void Tvmult(VectorType &dst, const VectorType &src) const
unsigned int global_dof_index
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)
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)
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)
types::global_dof_index size_type
std::vector< double > resid
std_cxx11::enable_if< std_cxx11::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static::ExceptionBase & PArpackExcIdo(int arg1)
void internal_reinit(const IndexSet &locally_owned_dofs)
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
Shift(const MatrixType &A, const MatrixType &B, const double sigma)
void reinit(const IndexSet &locally_owned_dofs)
std::vector< types::global_dof_index > local_indices
static::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
std::vector< double > workl
static::ExceptionBase & PArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
std::vector< double > workd
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & PArpackExcInfoPdneupd(int arg1)