16 #ifndef dealii_parpack_solver_h
17 #define dealii_parpack_solver_h
31 #ifdef DEAL_II_ARPACK_WITH_PARPACK
212 template <
typename VectorType>
317 const std::vector<IndexSet> &partitioning);
340 set_shift(
const std::complex<double> sigma);
351 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
353 solve(
const MatrixType1 &
A,
354 const MatrixType2 & B,
355 const INVERSE & inverse,
358 const unsigned int n_eigenvalues);
363 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
365 solve(
const MatrixType1 &
A,
366 const MatrixType2 & B,
367 const INVERSE & inverse,
370 const unsigned int n_eigenvalues);
439 std::vector<double>
v;
462 std::vector<double>
z;
515 << arg1 <<
" eigenpairs were requested, but only " << arg2
521 <<
"Number of wanted eigenvalues " << arg1
522 <<
" is larger that the size of the matrix " << arg2);
527 <<
"Number of wanted eigenvalues " << arg1
528 <<
" is larger that the size of eigenvectors " << arg2);
534 <<
"To store the real and complex parts of " << arg1
535 <<
" eigenvectors in real-valued vectors, their size (currently set to "
536 << arg2 <<
") should be greater than or equal to " << arg1 + 1);
541 <<
"Number of wanted eigenvalues " << arg1
542 <<
" is larger that the size of eigenvalues " << arg2);
547 <<
"Number of Arnoldi vectors " << arg1
548 <<
" is larger that the size of the matrix " << arg2);
553 <<
"Number of Arnoldi vectors " << arg1
554 <<
" is too small to obtain " << arg2 <<
" eigenvalues");
558 <<
"This ido " << arg1
559 <<
" is not supported. Check documentation of ARPACK");
563 <<
"This mode " << arg1
564 <<
" is not supported. Check documentation of ARPACK");
568 <<
"Error with Pdnaupd, info " << arg1
569 <<
". Check documentation of ARPACK");
573 <<
"Error with Pdneupd, info " << arg1
574 <<
". Check documentation of ARPACK");
578 <<
"Maximum number " << arg1 <<
" of iterations reached.");
582 <<
"No shifts could be applied during implicit"
583 <<
" Arnoldi update, try increasing the number of"
584 <<
" Arnoldi vectors.");
589 template <
typename VectorType>
594 (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
596 src.memory_consumption() + dst.memory_consumption() +
597 tmp.memory_consumption() +
599 local_indices.size();
604 template <
typename VectorType>
606 const unsigned int number_of_arnoldi_vectors,
610 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
611 , eigenvalue_of_interest(eigenvalue_of_interest)
621 "'largest real part' can only be used for non-symmetric problems!"));
625 "'smallest real part' can only be used for non-symmetric problems!"));
629 "'largest imaginary part' can only be used for non-symmetric problems!"));
633 "'smallest imaginary part' can only be used for non-symmetric problems!"));
636 ExcMessage(
"Currently, only modes 1, 2 and 3 are supported."));
641 template <
typename VectorType>
662 template <
typename VectorType>
666 sigmar = sigma.real();
667 sigmai = sigma.imag();
672 template <
typename VectorType>
676 initial_vector_provided =
true;
677 Assert(resid.size() == local_indices.size(),
679 vec.extract_subvector_to(local_indices.begin(),
686 template <
typename VectorType>
695 ncv = additional_data.number_of_arnoldi_vectors;
701 v.resize(ldv * ncv, 0.0);
703 resid.resize(nloc, 1.0);
706 workd.resize(3 * nloc, 0.0);
709 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
710 workl.resize(lworkl, 0.);
713 z.resize(ldz * ncv, 0.);
716 lworkev = additional_data.symmetric ? 0
719 workev.resize(lworkev, 0.);
721 select.resize(ncv, 0);
726 template <
typename VectorType>
730 internal_reinit(locally_owned_dofs);
733 src.reinit(locally_owned_dofs, mpi_communicator);
734 dst.reinit(locally_owned_dofs, mpi_communicator);
735 tmp.reinit(locally_owned_dofs, mpi_communicator);
740 template <
typename VectorType>
744 internal_reinit(distributed_vector.locally_owned_elements());
747 src.reinit(distributed_vector);
748 dst.reinit(distributed_vector);
749 tmp.reinit(distributed_vector);
754 template <
typename VectorType>
757 const std::vector<IndexSet> &partitioning)
759 internal_reinit(locally_owned_dofs);
762 src.reinit(partitioning, mpi_communicator);
763 dst.reinit(partitioning, mpi_communicator);
764 tmp.reinit(partitioning, mpi_communicator);
769 template <
typename VectorType>
770 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
773 const MatrixType2 & B,
774 const INVERSE & inverse,
777 const unsigned int n_eigenvalues)
779 std::vector<VectorType *> eigenvectors_ptr(
eigenvectors.size());
782 solve(
A, B, inverse,
eigenvalues, eigenvectors_ptr, n_eigenvalues);
787 template <
typename VectorType>
788 template <
typename MatrixType1,
typename MatrixType2,
typename INVERSE>
792 const INVERSE & inverse,
795 const unsigned int n_eigenvalues)
797 if (additional_data.symmetric)
800 PArpackExcInvalidEigenvectorSize(n_eigenvalues,
805 PArpackExcInvalidEigenvectorSizeNonsymmetric(n_eigenvalues,
809 PArpackExcInvalidEigenvalueSize(n_eigenvalues,
eigenvalues.size()));
815 PArpackExcInvalidNumberofEigenvalues(n_eigenvalues,
819 PArpackExcInvalidNumberofArnoldiVectors(
820 additional_data.number_of_arnoldi_vectors,
eigenvectors[0]->size()));
822 Assert(additional_data.number_of_arnoldi_vectors > 2 * n_eigenvalues + 1,
823 PArpackExcSmallNumberofArnoldiVectors(
824 additional_data.number_of_arnoldi_vectors, n_eigenvalues));
826 int mode = additional_data.mode;
835 bmat[0] = (mode == 1) ?
'I' :
'G';
849 switch (additional_data.eigenvalue_of_interest)
851 case algebraically_largest:
852 std::strcpy(which,
"LA");
854 case algebraically_smallest:
855 std::strcpy(which,
"SA");
857 case largest_magnitude:
858 std::strcpy(which,
"LM");
860 case smallest_magnitude:
861 std::strcpy(which,
"SM");
863 case largest_real_part:
864 std::strcpy(which,
"LR");
866 case smallest_real_part:
867 std::strcpy(which,
"SR");
869 case largest_imaginary_part:
870 std::strcpy(which,
"LI");
872 case smallest_imaginary_part:
873 std::strcpy(which,
"SI");
876 std::strcpy(which,
"BE");
881 double tol = control().tolerance();
884 std::vector<int> iparam(11, 0);
891 iparam[2] = control().max_steps();
904 std::vector<int> ipntr(14, 0);
912 int info = initial_vector_provided ? 1 : 0;
915 int nev = n_eigenvalues;
916 int n_inside_arpack = nloc;
922 if (additional_data.symmetric)
959 AssertThrow(info == 0, PArpackExcInfoPdnaupd(info));
967 const int shift_x = ipntr[0] - 1;
968 const int shift_y = ipntr[1] - 1;
970 Assert(shift_x + nloc <=
static_cast<int>(workd.size()),
973 Assert(shift_y + nloc <=
static_cast<int>(workd.size()),
979 if ((ido == -1) || (ido == 1 && mode < 3))
982 src.add(nloc, local_indices.data(), workd.data() + shift_x);
989 inverse.vmult(dst, tmp);
994 system_matrix.vmult(tmp, src);
996 tmp.extract_subvector_to(local_indices.begin(),
998 workd.data() + shift_x);
999 inverse.vmult(dst, tmp);
1003 system_matrix.vmult(dst, src);
1008 else if (ido == 1 && mode >= 3)
1012 const int shift_b_x = ipntr[2] - 1;
1014 Assert(shift_b_x + nloc <=
static_cast<int>(workd.size()),
1018 src.add(nloc, local_indices.data(), workd.data() + shift_b_x);
1023 inverse.vmult(dst, src);
1028 src.add(nloc, local_indices.data(), workd.data() + shift_x);
1048 dst.extract_subvector_to(local_indices.begin(),
1049 local_indices.end(),
1050 workd.data() + shift_y);
1058 char howmany[4] =
"All";
1060 std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1061 std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1064 if (additional_data.symmetric)
1065 pdseupd_(&mpi_communicator_fortran,
1069 eigenvalues_real.data(),
1089 pdneupd_(&mpi_communicator_fortran,
1093 eigenvalues_real.data(),
1094 eigenvalues_im.data(),
1118 AssertThrow(
false, PArpackExcInfoMaxIt(control().max_steps()));
1129 for (
int i = 0; i < nev; ++i)
1134 eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1138 for (
size_type i = 0; i < n_eigenvalues; ++i)
1140 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1143 AssertThrow(iparam[4] >=
static_cast<int>(n_eigenvalues),
1144 PArpackExcConvergedEigenvectors(n_eigenvalues, iparam[4]));
1152 tmp.add(nloc, local_indices.data(), resid.data());
1153 solver_control.check(iparam[2], tmp.l2_norm());
1159 template <
typename VectorType>
1163 return solver_control;