Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
parpack_solver.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_parpack_solver_h
17 #define dealii_parpack_solver_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/index_set.h>
24 
27 
28 #include <cstring>
29 
30 
31 #ifdef DEAL_II_ARPACK_WITH_PARPACK
32 
34 
35 extern "C"
36 {
37  // http://www.mathkeisan.com/usersguide/man/pdnaupd.html
38  void
39  pdnaupd_(MPI_Fint *comm,
40  int * ido,
41  char * bmat,
42  int * n,
43  char * which,
44  int * nev,
45  double * tol,
46  double * resid,
47  int * ncv,
48  double * v,
49  int * nloc,
50  int * iparam,
51  int * ipntr,
52  double * workd,
53  double * workl,
54  int * lworkl,
55  int * info);
56 
57  // http://www.mathkeisan.com/usersguide/man/pdsaupd.html
58  void
59  pdsaupd_(MPI_Fint *comm,
60  int * ido,
61  char * bmat,
62  int * n,
63  char * which,
64  int * nev,
65  double * tol,
66  double * resid,
67  int * ncv,
68  double * v,
69  int * nloc,
70  int * iparam,
71  int * ipntr,
72  double * workd,
73  double * workl,
74  int * lworkl,
75  int * info);
76 
77  // http://www.mathkeisan.com/usersguide/man/pdneupd.html
78  void
79  pdneupd_(MPI_Fint *comm,
80  int * rvec,
81  char * howmany,
82  int * select,
83  double * d,
84  double * di,
85  double * z,
86  int * ldz,
87  double * sigmar,
88  double * sigmai,
89  double * workev,
90  char * bmat,
91  int * n,
92  char * which,
93  int * nev,
94  double * tol,
95  double * resid,
96  int * ncv,
97  double * v,
98  int * nloc,
99  int * iparam,
100  int * ipntr,
101  double * workd,
102  double * workl,
103  int * lworkl,
104  int * info);
105 
106  // http://www.mathkeisan.com/usersguide/man/pdseupd.html
107  void
108  pdseupd_(MPI_Fint *comm,
109  int * rvec,
110  char * howmany,
111  int * select,
112  double * d,
113  double * z,
114  int * ldz,
115  double * sigmar,
116  char * bmat,
117  int * n,
118  char * which,
119  int * nev,
120  double * tol,
121  double * resid,
122  int * ncv,
123  double * v,
124  int * nloc,
125  int * iparam,
126  int * ipntr,
127  double * workd,
128  double * workl,
129  int * lworkl,
130  int * info);
131 
132  // other resources:
133  // http://acts.nersc.gov/superlu/example5/pnslac.c.html
134  // https://github.com/phpisciuneri/tijo/blob/master/dvr_parpack.cpp
135 }
136 
212 template <typename VectorType>
214 {
215 public:
220 
231  {
271  };
272 
278  {
279  const unsigned int number_of_arnoldi_vectors;
281  const bool symmetric;
282  const int mode;
284  const unsigned int number_of_arnoldi_vectors = 15,
285  const WhichEigenvalues eigenvalue_of_interest = largest_magnitude,
286  const bool symmetric = false,
287  const int mode = 3);
288  };
289 
293  SolverControl &
294  control() const;
295 
300  const MPI_Comm & mpi_communicator,
301  const AdditionalData &data = AdditionalData());
302 
306  void
307  reinit(const IndexSet &locally_owned_dofs);
308 
315  void
316  reinit(const IndexSet & locally_owned_dofs,
317  const std::vector<IndexSet> &partitioning);
318 
322  void
323  reinit(const VectorType &distributed_vector);
324 
328  void
329  set_initial_vector(const VectorType &vec);
330 
339  void
340  set_shift(const std::complex<double> sigma);
341 
351  template <typename MatrixType1, typename MatrixType2, typename INVERSE>
352  void
353  solve(const MatrixType1 & A,
354  const MatrixType2 & B,
355  const INVERSE & inverse,
356  std::vector<std::complex<double>> &eigenvalues,
357  std::vector<VectorType> & eigenvectors,
358  const unsigned int n_eigenvalues);
359 
363  template <typename MatrixType1, typename MatrixType2, typename INVERSE>
364  void
365  solve(const MatrixType1 & A,
366  const MatrixType2 & B,
367  const INVERSE & inverse,
368  std::vector<std::complex<double>> &eigenvalues,
369  std::vector<VectorType *> & eigenvectors,
370  const unsigned int n_eigenvalues);
371 
375  std::size_t
376  memory_consumption() const;
377 
378 protected:
384 
389 
390  // keep MPI communicator non-const as Arpack functions are not const either:
391 
396 
401 
402  // C++98 guarantees that the elements of a vector are stored contiguously
403 
407  int lworkl;
408 
412  std::vector<double> workl;
413 
417  std::vector<double> workd;
418 
422  int nloc;
423 
427  int ncv;
428 
429 
433  int ldv;
434 
439  std::vector<double> v;
440 
445 
450  std::vector<double> resid;
451 
455  int ldz;
456 
462  std::vector<double> z;
463 
467  int lworkev;
468 
472  std::vector<double> workev;
473 
477  std::vector<int> select;
478 
483 
487  std::vector<types::global_dof_index> local_indices;
488 
492  double sigmar;
493 
497  double sigmai;
498 
499 private:
506  void
507  internal_reinit(const IndexSet &locally_owned_dofs);
508 
513  int,
514  int,
515  << arg1 << " eigenpairs were requested, but only " << arg2
516  << " converged");
517 
519  int,
520  int,
521  << "Number of wanted eigenvalues " << arg1
522  << " is larger that the size of the matrix " << arg2);
523 
525  int,
526  int,
527  << "Number of wanted eigenvalues " << arg1
528  << " is larger that the size of eigenvectors " << arg2);
529 
532  int,
533  int,
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);
537 
539  int,
540  int,
541  << "Number of wanted eigenvalues " << arg1
542  << " is larger that the size of eigenvalues " << arg2);
543 
545  int,
546  int,
547  << "Number of Arnoldi vectors " << arg1
548  << " is larger that the size of the matrix " << arg2);
549 
551  int,
552  int,
553  << "Number of Arnoldi vectors " << arg1
554  << " is too small to obtain " << arg2 << " eigenvalues");
555 
557  int,
558  << "This ido " << arg1
559  << " is not supported. Check documentation of ARPACK");
560 
562  int,
563  << "This mode " << arg1
564  << " is not supported. Check documentation of ARPACK");
565 
567  int,
568  << "Error with Pdnaupd, info " << arg1
569  << ". Check documentation of ARPACK");
570 
572  int,
573  << "Error with Pdneupd, info " << arg1
574  << ". Check documentation of ARPACK");
575 
577  int,
578  << "Maximum number " << arg1 << " of iterations reached.");
579 
581  int,
582  << "No shifts could be applied during implicit"
583  << " Arnoldi update, try increasing the number of"
584  << " Arnoldi vectors.");
585 };
586 
587 
588 
589 template <typename VectorType>
590 std::size_t
592 {
593  return MemoryConsumption::memory_consumption(double()) *
594  (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
595  workev.size()) +
596  src.memory_consumption() + dst.memory_consumption() +
597  tmp.memory_consumption() +
599  local_indices.size();
600 }
601 
602 
603 
604 template <typename VectorType>
606  const unsigned int number_of_arnoldi_vectors,
608  const bool symmetric,
609  const int mode)
610  : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
611  , eigenvalue_of_interest(eigenvalue_of_interest)
612  , symmetric(symmetric)
613  , mode(mode)
614 {
615  // Check for possible options for symmetric problems
616  if (symmetric)
617  {
618  Assert(
619  eigenvalue_of_interest != largest_real_part,
620  ExcMessage(
621  "'largest real part' can only be used for non-symmetric problems!"));
622  Assert(
623  eigenvalue_of_interest != smallest_real_part,
624  ExcMessage(
625  "'smallest real part' can only be used for non-symmetric problems!"));
626  Assert(
627  eigenvalue_of_interest != largest_imaginary_part,
628  ExcMessage(
629  "'largest imaginary part' can only be used for non-symmetric problems!"));
630  Assert(
631  eigenvalue_of_interest != smallest_imaginary_part,
632  ExcMessage(
633  "'smallest imaginary part' can only be used for non-symmetric problems!"));
634  }
635  Assert(mode >= 1 && mode <= 3,
636  ExcMessage("Currently, only modes 1, 2 and 3 are supported."));
637 }
638 
639 
640 
641 template <typename VectorType>
643  const MPI_Comm & mpi_communicator,
644  const AdditionalData &data)
645  : solver_control(control)
646  , additional_data(data)
647  , mpi_communicator(mpi_communicator)
648  , mpi_communicator_fortran(MPI_Comm_c2f(mpi_communicator))
649  , lworkl(0)
650  , nloc(0)
651  , ncv(0)
652  , ldv(0)
653  , initial_vector_provided(false)
654  , ldz(0)
655  , lworkev(0)
656  , sigmar(0.0)
657  , sigmai(0.0)
658 {}
659 
660 
661 
662 template <typename VectorType>
663 void
664 PArpackSolver<VectorType>::set_shift(const std::complex<double> sigma)
665 {
666  sigmar = sigma.real();
667  sigmai = sigma.imag();
668 }
669 
670 
671 
672 template <typename VectorType>
673 void
675 {
677  Assert(resid.size() == local_indices.size(),
678  ExcDimensionMismatch(resid.size(), local_indices.size()));
679  vec.extract_subvector_to(local_indices.begin(),
680  local_indices.end(),
681  resid.data());
682 }
683 
684 
685 
686 template <typename VectorType>
687 void
689 {
690  // store local indices to write to vectors
691  locally_owned_dofs.fill_index_vector(local_indices);
692 
693  // scalars
694  nloc = locally_owned_dofs.n_elements();
696 
698 
699  // vectors
700  ldv = nloc;
701  v.resize(ldv * ncv, 0.0);
702 
703  resid.resize(nloc, 1.0);
704 
705  // work arrays for ARPACK
706  workd.resize(3 * nloc, 0.0);
707 
708  lworkl =
709  additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
710  workl.resize(lworkl, 0.);
711 
712  ldz = nloc;
713  z.resize(ldz * ncv, 0.); // TODO we actually need only ldz*nev
714 
715  // WORKEV Double precision work array of dimension 3*NCV.
716  lworkev = additional_data.symmetric ? 0 /*not used in symmetric case*/
717  :
718  3 * ncv;
719  workev.resize(lworkev, 0.);
720 
721  select.resize(ncv, 0);
722 }
723 
724 
725 
726 template <typename VectorType>
727 void
728 PArpackSolver<VectorType>::reinit(const IndexSet &locally_owned_dofs)
729 {
730  internal_reinit(locally_owned_dofs);
731 
732  // deal.II vectors:
733  src.reinit(locally_owned_dofs, mpi_communicator);
734  dst.reinit(locally_owned_dofs, mpi_communicator);
735  tmp.reinit(locally_owned_dofs, mpi_communicator);
736 }
737 
738 
739 
740 template <typename VectorType>
741 void
743 {
744  internal_reinit(distributed_vector.locally_owned_elements());
745 
746  // deal.II vectors:
747  src.reinit(distributed_vector);
748  dst.reinit(distributed_vector);
749  tmp.reinit(distributed_vector);
750 }
751 
752 
753 
754 template <typename VectorType>
755 void
756 PArpackSolver<VectorType>::reinit(const IndexSet &locally_owned_dofs,
757  const std::vector<IndexSet> &partitioning)
758 {
759  internal_reinit(locally_owned_dofs);
760 
761  // deal.II vectors:
762  src.reinit(partitioning, mpi_communicator);
763  dst.reinit(partitioning, mpi_communicator);
764  tmp.reinit(partitioning, mpi_communicator);
765 }
766 
767 
768 
769 template <typename VectorType>
770 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
771 void
773  const MatrixType2 & B,
774  const INVERSE & inverse,
775  std::vector<std::complex<double>> &eigenvalues,
776  std::vector<VectorType> &eigenvectors,
777  const unsigned int n_eigenvalues)
778 {
779  std::vector<VectorType *> eigenvectors_ptr(eigenvectors.size());
780  for (unsigned int i = 0; i < eigenvectors.size(); ++i)
781  eigenvectors_ptr[i] = &eigenvectors[i];
782  solve(A, B, inverse, eigenvalues, eigenvectors_ptr, n_eigenvalues);
783 }
784 
785 
786 
787 template <typename VectorType>
788 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
789 void
790 PArpackSolver<VectorType>::solve(const MatrixType1 &system_matrix,
791  const MatrixType2 &mass_matrix,
792  const INVERSE & inverse,
793  std::vector<std::complex<double>> &eigenvalues,
794  std::vector<VectorType *> &eigenvectors,
795  const unsigned int n_eigenvalues)
796 {
798  {
799  Assert(n_eigenvalues <= eigenvectors.size(),
800  PArpackExcInvalidEigenvectorSize(n_eigenvalues,
801  eigenvectors.size()));
802  }
803  else
804  Assert(n_eigenvalues + 1 <= eigenvectors.size(),
806  eigenvectors.size()));
807 
808  Assert(n_eigenvalues <= eigenvalues.size(),
809  PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
810 
811 
812  // use eigenvectors to get the problem size so that it is possible to
813  // employ LinearOperator for mass_matrix.
814  Assert(n_eigenvalues < eigenvectors[0]->size(),
816  eigenvectors[0]->size()));
817 
818  Assert(additional_data.number_of_arnoldi_vectors < eigenvectors[0]->size(),
820  additional_data.number_of_arnoldi_vectors, eigenvectors[0]->size()));
821 
822  Assert(additional_data.number_of_arnoldi_vectors > 2 * n_eigenvalues + 1,
825 
826  int mode = additional_data.mode;
827 
828  // reverse communication parameter
829  // must be zero on the first call to pdnaupd
830  int ido = 0;
831 
832  // 'G' generalized eigenvalue problem
833  // 'I' standard eigenvalue problem
834  char bmat[2];
835  bmat[0] = (mode == 1) ? 'I' : 'G';
836  bmat[1] = '\0';
837 
838  // Specify the eigenvalues of interest, possible parameters:
839  // "LA" algebraically largest
840  // "SA" algebraically smallest
841  // "LM" largest magnitude
842  // "SM" smallest magnitude
843  // "LR" largest real part
844  // "SR" smallest real part
845  // "LI" largest imaginary part
846  // "SI" smallest imaginary part
847  // "BE" both ends of spectrum simultaneous
848  char which[3];
850  {
852  std::strcpy(which, "LA");
853  break;
855  std::strcpy(which, "SA");
856  break;
857  case largest_magnitude:
858  std::strcpy(which, "LM");
859  break;
860  case smallest_magnitude:
861  std::strcpy(which, "SM");
862  break;
863  case largest_real_part:
864  std::strcpy(which, "LR");
865  break;
866  case smallest_real_part:
867  std::strcpy(which, "SR");
868  break;
870  std::strcpy(which, "LI");
871  break;
873  std::strcpy(which, "SI");
874  break;
875  case both_ends:
876  std::strcpy(which, "BE");
877  break;
878  }
879 
880  // tolerance for ARPACK
881  double tol = control().tolerance();
882 
883  // information to the routines
884  std::vector<int> iparam(11, 0);
885 
886  iparam[0] = 1;
887  // shift strategy: exact shifts with respect to the current Hessenberg matrix
888  // H.
889 
890  // maximum number of iterations
891  iparam[2] = control().max_steps();
892 
893  // Parpack currently works only for NB = 1
894  iparam[3] = 1;
895 
896  // Sets the mode of dsaupd:
897  // 1 is A*x=lambda*x, OP = A, B = I
898  // 2 is A*x = lambda*M*x, OP = inv[M]*A, B = M
899  // 3 is shift-invert mode, OP = inv[A-sigma*M]*M, B = M
900  // 4 is buckling mode,
901  // 5 is Cayley mode.
902 
903  iparam[6] = mode;
904  std::vector<int> ipntr(14, 0);
905 
906  // information out of the iteration
907  // If INFO .EQ. 0, a random initial residual vector is used.
908  // If INFO .NE. 0, RESID contains the initial residual vector,
909  // possibly from a previous run.
910  // Typical choices in this situation might be to use the final value
911  // of the starting vector from the previous eigenvalue calculation
912  int info = initial_vector_provided ? 1 : 0;
913 
914  // Number of eigenvalues of OP to be computed. 0 < NEV < N.
915  int nev = n_eigenvalues;
916  int n_inside_arpack = nloc;
917 
918  // IDO = 99: done
919  while (ido != 99)
920  {
921  // call of ARPACK pdnaupd routine
924  &ido,
925  bmat,
926  &n_inside_arpack,
927  which,
928  &nev,
929  &tol,
930  resid.data(),
931  &ncv,
932  v.data(),
933  &ldv,
934  iparam.data(),
935  ipntr.data(),
936  workd.data(),
937  workl.data(),
938  &lworkl,
939  &info);
940  else
942  &ido,
943  bmat,
944  &n_inside_arpack,
945  which,
946  &nev,
947  &tol,
948  resid.data(),
949  &ncv,
950  v.data(),
951  &ldv,
952  iparam.data(),
953  ipntr.data(),
954  workd.data(),
955  workl.data(),
956  &lworkl,
957  &info);
958 
959  AssertThrow(info == 0, PArpackExcInfoPdnaupd(info));
960 
961  // if we converge, we shall not modify anything in work arrays!
962  if (ido == 99)
963  break;
964 
965  // IPNTR(1) is the pointer into WORKD for X,
966  // IPNTR(2) is the pointer into WORKD for Y.
967  const int shift_x = ipntr[0] - 1;
968  const int shift_y = ipntr[1] - 1;
969  Assert(shift_x >= 0, ::ExcInternalError());
970  Assert(shift_x + nloc <= static_cast<int>(workd.size()),
971  ::ExcInternalError());
972  Assert(shift_y >= 0, ::ExcInternalError());
973  Assert(shift_y + nloc <= static_cast<int>(workd.size()),
974  ::ExcInternalError());
975 
976  src = 0.;
977 
978  // switch based on both ido and mode
979  if ((ido == -1) || (ido == 1 && mode < 3))
980  // compute Y = OP * X
981  {
982  src.add(nloc, local_indices.data(), workd.data() + shift_x);
983  src.compress(VectorOperation::add);
984 
985  if (mode == 3)
986  // OP = inv[K - sigma*M]*M
987  {
988  mass_matrix.vmult(tmp, src);
989  inverse.vmult(dst, tmp);
990  }
991  else if (mode == 2)
992  // OP = inv[M]*K
993  {
994  system_matrix.vmult(tmp, src);
995  // store M*X in X
996  tmp.extract_subvector_to(local_indices.begin(),
997  local_indices.end(),
998  workd.data() + shift_x);
999  inverse.vmult(dst, tmp);
1000  }
1001  else if (mode == 1)
1002  {
1003  system_matrix.vmult(dst, src);
1004  }
1005  else
1006  AssertThrow(false, PArpackExcMode(mode));
1007  }
1008  else if (ido == 1 && mode >= 3)
1009  // compute Y = OP * X for mode 3, 4 and 5, where
1010  // the vector B * X is already available in WORKD(ipntr(3)).
1011  {
1012  const int shift_b_x = ipntr[2] - 1;
1013  Assert(shift_b_x >= 0, ::ExcInternalError());
1014  Assert(shift_b_x + nloc <= static_cast<int>(workd.size()),
1015  ::ExcInternalError());
1016 
1017  // B*X
1018  src.add(nloc, local_indices.data(), workd.data() + shift_b_x);
1019  src.compress(VectorOperation::add);
1020 
1021  // solving linear system
1022  Assert(mode == 3, ExcNotImplemented());
1023  inverse.vmult(dst, src);
1024  }
1025  else if (ido == 2)
1026  // compute Y = B * X
1027  {
1028  src.add(nloc, local_indices.data(), workd.data() + shift_x);
1029  src.compress(VectorOperation::add);
1030 
1031  // Multiplication with mass matrix M
1032  if (mode == 1)
1033  {
1034  dst = src;
1035  }
1036  else
1037  // mode 2,3 and 5 have B=M
1038  {
1039  mass_matrix.vmult(dst, src);
1040  }
1041  }
1042  else
1043  AssertThrow(false, PArpackExcIdo(ido));
1044  // Note: IDO = 3 does not appear to be required for currently
1045  // implemented modes
1046 
1047  // store the result
1048  dst.extract_subvector_to(local_indices.begin(),
1049  local_indices.end(),
1050  workd.data() + shift_y);
1051  } // end of pd*aupd_ loop
1052 
1053  // 1 - compute eigenvectors,
1054  // 0 - only eigenvalues
1055  int rvec = 1;
1056 
1057  // which eigenvectors
1058  char howmany[4] = "All";
1059 
1060  std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1061  std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1062 
1063  // call of ARPACK pdneupd routine
1066  &rvec,
1067  howmany,
1068  select.data(),
1069  eigenvalues_real.data(),
1070  z.data(),
1071  &ldz,
1072  &sigmar,
1073  bmat,
1074  &n_inside_arpack,
1075  which,
1076  &nev,
1077  &tol,
1078  resid.data(),
1079  &ncv,
1080  v.data(),
1081  &ldv,
1082  iparam.data(),
1083  ipntr.data(),
1084  workd.data(),
1085  workl.data(),
1086  &lworkl,
1087  &info);
1088  else
1090  &rvec,
1091  howmany,
1092  select.data(),
1093  eigenvalues_real.data(),
1094  eigenvalues_im.data(),
1095  v.data(),
1096  &ldz,
1097  &sigmar,
1098  &sigmai,
1099  workev.data(),
1100  bmat,
1101  &n_inside_arpack,
1102  which,
1103  &nev,
1104  &tol,
1105  resid.data(),
1106  &ncv,
1107  v.data(),
1108  &ldv,
1109  iparam.data(),
1110  ipntr.data(),
1111  workd.data(),
1112  workl.data(),
1113  &lworkl,
1114  &info);
1115 
1116  if (info == 1)
1117  {
1118  AssertThrow(false, PArpackExcInfoMaxIt(control().max_steps()));
1119  }
1120  else if (info == 3)
1121  {
1122  AssertThrow(false, PArpackExcNoShifts(1));
1123  }
1124  else if (info != 0)
1125  {
1126  AssertThrow(false, PArpackExcInfoPdneupd(info));
1127  }
1128 
1129  for (int i = 0; i < nev; ++i)
1130  {
1131  (*eigenvectors[i]) = 0.0;
1132  AssertIndexRange(i * nloc + nloc, v.size() + 1);
1133 
1134  eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1135  eigenvectors[i]->compress(VectorOperation::add);
1136  }
1137 
1138  for (size_type i = 0; i < n_eigenvalues; ++i)
1139  eigenvalues[i] =
1140  std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1141 
1142  // Throw an error if the solver did not converge.
1143  AssertThrow(iparam[4] >= static_cast<int>(n_eigenvalues),
1144  PArpackExcConvergedEigenvectors(n_eigenvalues, iparam[4]));
1145 
1146  // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
1147  // with respect to a semi-inner product defined by M.
1148 
1149  // resid likely contains residual with respect to M-norm.
1150  {
1151  tmp = 0.0;
1152  tmp.add(nloc, local_indices.data(), resid.data());
1153  solver_control.check(iparam[2], tmp.l2_norm());
1154  }
1155 }
1156 
1157 
1158 
1159 template <typename VectorType>
1160 SolverControl &
1162 {
1163  return solver_control;
1164 }
1165 
1167 
1168 
1169 #endif
1170 #endif
void pdnaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *nloc, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
static ::ExceptionBase & PArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
static ::ExceptionBase & PArpackExcInvalidEigenvectorSize(int arg1, int arg2)
virtual State check(const unsigned int step, const double check_value)
std::vector< double > z
std::vector< double > workev
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)
static ::ExceptionBase & PArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
const unsigned int number_of_arnoldi_vectors
VectorType dst
const WhichEigenvalues eigenvalue_of_interest
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
static ::ExceptionBase & PArpackExcMode(int arg1)
PArpackSolver(SolverControl &control, const MPI_Comm &mpi_communicator, const AdditionalData &data=AdditionalData())
MPI_Fint mpi_communicator_fortran
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
void pdseupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d, double *z, int *ldz, double *sigmar, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *nloc, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
const AdditionalData additional_data
void set_shift(const std::complex< double > sigma)
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false, const int mode=3)
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & PArpackExcNoShifts(int arg1)
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)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:518
void mass_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const double factor=1.)
Definition: l2.h:63
bool initial_vector_provided
SolverControl & control() const
#define Assert(cond, exc)
Definition: exceptions.h:1419
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< int > select
static ::ExceptionBase & PArpackExcConvergedEigenvectors(int arg1, int arg2)
static ::ExceptionBase & PArpackExcInfoPdnaupd(int arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType src
static ::ExceptionBase & PArpackExcInfoMaxIt(int arg1)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
MPI_Comm mpi_communicator
static ::ExceptionBase & PArpackExcInvalidEigenvalueSize(int arg1, int arg2)
double tolerance() const
void fill_index_vector(std::vector< size_type > &indices) const
Definition: index_set.cc:507
std::vector< double > v
static const char A
std::vector< double > resid
unsigned int global_dof_index
Definition: types.h:76
void pdneupd_(MPI_Fint *comm, int *rvec, char *howmany, int *select, double *d, double *di, double *z, int *ldz, double *sigmar, double *sigmai, double *workev, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *nloc, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
static ::ExceptionBase & PArpackExcIdo(int arg1)
void internal_reinit(const IndexSet &locally_owned_dofs)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
unsigned int max_steps() const
VectorType tmp
void pdsaupd_(MPI_Fint *comm, int *ido, char *bmat, int *n, char *which, int *nev, double *tol, double *resid, int *ncv, double *v, int *nloc, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
static ::ExceptionBase & ExcNotImplemented()
std::size_t memory_consumption() const
void reinit(const IndexSet &locally_owned_dofs)
static ::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
std::vector< double > workl
Eigenvalue vector is filled.
size_type n_elements() const
Definition: index_set.h:1833
static ::ExceptionBase & PArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
std::vector< types::global_dof_index > local_indices
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
std::vector< double > workd
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & PArpackExcInfoPdneupd(int arg1)