Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
parpack_solver.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2018 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>
22 #include <deal.II/base/memory_consumption.h>
23 #include <deal.II/base/smartpointer.h>
24 
25 #include <deal.II/lac/solver_control.h>
26 #include <deal.II/lac/vector_operation.h>
27 
28 #include <cstring>
29 
30 
31 #ifdef DEAL_II_ARPACK_WITH_PARPACK
32 
33 DEAL_II_NAMESPACE_OPEN
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 
281  template <typename MatrixType>
282  class DEAL_II_DEPRECATED Shift : public ::Subscriptor
283  {
284  public:
288  Shift(const MatrixType &A, const MatrixType &B, const double sigma)
289  : A(A)
290  , B(B)
291  , sigma(sigma)
292  {}
293 
297  void
298  vmult(VectorType &dst, const VectorType &src) const
299  {
300  B.vmult(dst, src);
301  dst *= (-sigma);
302  A.vmult_add(dst, src);
303  }
304 
308  void
309  Tvmult(VectorType &dst, const VectorType &src) const
310  {
311  B.Tvmult(dst, src);
312  dst *= (-sigma);
313  A.Tvmult_add(dst, src);
314  }
315 
316  private:
317  const MatrixType &A;
318  const MatrixType &B;
319  const double sigma;
320  };
321 
327  {
328  const unsigned int number_of_arnoldi_vectors;
329  const WhichEigenvalues eigenvalue_of_interest;
330  const bool symmetric;
331  const int mode;
333  const unsigned int number_of_arnoldi_vectors = 15,
334  const WhichEigenvalues eigenvalue_of_interest = largest_magnitude,
335  const bool symmetric = false,
336  const int mode = 3);
337  };
338 
342  SolverControl &
343  control() const;
344 
348  PArpackSolver(SolverControl & control,
349  const MPI_Comm & mpi_communicator,
350  const AdditionalData &data = AdditionalData());
351 
355  void
356  reinit(const IndexSet &locally_owned_dofs);
357 
364  void
365  reinit(const IndexSet & locally_owned_dofs,
366  const std::vector<IndexSet> &partitioning);
367 
371  void
372  reinit(const VectorType &distributed_vector);
373 
377  void
378  set_initial_vector(const VectorType &vec);
379 
388  void
389  set_shift(const std::complex<double> sigma);
390 
400  template <typename MatrixType1, typename MatrixType2, typename INVERSE>
401  void
402  solve(const MatrixType1 & A,
403  const MatrixType2 & B,
404  const INVERSE & inverse,
405  std::vector<std::complex<double>> &eigenvalues,
406  std::vector<VectorType> & eigenvectors,
407  const unsigned int n_eigenvalues);
408 
412  template <typename MatrixType1, typename MatrixType2, typename INVERSE>
413  void
414  solve(const MatrixType1 & A,
415  const MatrixType2 & B,
416  const INVERSE & inverse,
417  std::vector<std::complex<double>> &eigenvalues,
418  std::vector<VectorType *> & eigenvectors,
419  const unsigned int n_eigenvalues);
420 
424  std::size_t
425  memory_consumption() const;
426 
427 protected:
433 
438 
439  // keep MPI communicator non-const as Arpack functions are not const either:
440 
445 
450 
451  // C++98 guarantees that the elements of a vector are stored contiguously
452 
456  int lworkl;
457 
461  std::vector<double> workl;
462 
466  std::vector<double> workd;
467 
471  int nloc;
472 
476  int ncv;
477 
478 
482  int ldv;
483 
488  std::vector<double> v;
489 
494 
499  std::vector<double> resid;
500 
504  int ldz;
505 
511  std::vector<double> z;
512 
516  int lworkev;
517 
521  std::vector<double> workev;
522 
526  std::vector<int> select;
527 
531  VectorType src, dst, tmp;
532 
536  std::vector<types::global_dof_index> local_indices;
537 
541  double sigmar;
542 
546  double sigmai;
547 
548 private:
555  void
556  internal_reinit(const IndexSet &locally_owned_dofs);
557 
562  int,
563  int,
564  << arg1 << " eigenpairs were requested, but only " << arg2
565  << " converged");
566 
568  int,
569  int,
570  << "Number of wanted eigenvalues " << arg1
571  << " is larger that the size of the matrix " << arg2);
572 
574  int,
575  int,
576  << "Number of wanted eigenvalues " << arg1
577  << " is larger that the size of eigenvectors " << arg2);
578 
581  int,
582  int,
583  << "To store the real and complex parts of " << arg1
584  << " eigenvectors in real-valued vectors, their size (currently set to "
585  << arg2 << ") should be greater than or equal to " << arg1 + 1);
586 
588  int,
589  int,
590  << "Number of wanted eigenvalues " << arg1
591  << " is larger that the size of eigenvalues " << arg2);
592 
594  int,
595  int,
596  << "Number of Arnoldi vectors " << arg1
597  << " is larger that the size of the matrix " << arg2);
598 
600  int,
601  int,
602  << "Number of Arnoldi vectors " << arg1
603  << " is too small to obtain " << arg2 << " eigenvalues");
604 
606  int,
607  << "This ido " << arg1
608  << " is not supported. Check documentation of ARPACK");
609 
611  int,
612  << "This mode " << arg1
613  << " is not supported. Check documentation of ARPACK");
614 
616  int,
617  << "Error with Pdnaupd, info " << arg1
618  << ". Check documentation of ARPACK");
619 
621  int,
622  << "Error with Pdneupd, info " << arg1
623  << ". Check documentation of ARPACK");
624 
626  int,
627  << "Maximum number " << arg1 << " of iterations reached.");
628 
630  int,
631  << "No shifts could be applied during implicit"
632  << " Arnoldi update, try increasing the number of"
633  << " Arnoldi vectors.");
634 };
635 
636 
637 
638 template <typename VectorType>
639 std::size_t
641 {
642  return MemoryConsumption::memory_consumption(double()) *
643  (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
644  workev.size()) +
645  src.memory_consumption() + dst.memory_consumption() +
646  tmp.memory_consumption() +
648  local_indices.size();
649 }
650 
651 
652 
653 template <typename VectorType>
655  const unsigned int number_of_arnoldi_vectors,
656  const WhichEigenvalues eigenvalue_of_interest,
657  const bool symmetric,
658  const int mode)
659  : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
660  , eigenvalue_of_interest(eigenvalue_of_interest)
661  , symmetric(symmetric)
662  , mode(mode)
663 {
664  // Check for possible options for symmetric problems
665  if (symmetric)
666  {
667  Assert(
668  eigenvalue_of_interest != largest_real_part,
669  ExcMessage(
670  "'largest real part' can only be used for non-symmetric problems!"));
671  Assert(
672  eigenvalue_of_interest != smallest_real_part,
673  ExcMessage(
674  "'smallest real part' can only be used for non-symmetric problems!"));
675  Assert(
676  eigenvalue_of_interest != largest_imaginary_part,
677  ExcMessage(
678  "'largest imaginary part' can only be used for non-symmetric problems!"));
679  Assert(
680  eigenvalue_of_interest != smallest_imaginary_part,
681  ExcMessage(
682  "'smallest imaginary part' can only be used for non-symmetric problems!"));
683  }
684  Assert(mode >= 1 && mode <= 3,
685  ExcMessage("Currently, only modes 1, 2 and 3 are supported."));
686 }
687 
688 
689 
690 template <typename VectorType>
692  const MPI_Comm & mpi_communicator,
693  const AdditionalData &data)
694  : solver_control(control)
695  , additional_data(data)
696  , mpi_communicator(mpi_communicator)
697  , mpi_communicator_fortran(MPI_Comm_c2f(mpi_communicator))
698  , lworkl(0)
699  , nloc(0)
700  , ncv(0)
701  , ldv(0)
702  , initial_vector_provided(false)
703  , ldz(0)
704  , lworkev(0)
705  , sigmar(0.0)
706  , sigmai(0.0)
707 {}
708 
709 
710 
711 template <typename VectorType>
712 void
713 PArpackSolver<VectorType>::set_shift(const std::complex<double> sigma)
714 {
715  sigmar = sigma.real();
716  sigmai = sigma.imag();
717 }
718 
719 
720 
721 template <typename VectorType>
722 void
724 {
726  Assert(resid.size() == local_indices.size(),
727  ExcDimensionMismatch(resid.size(), local_indices.size()));
728  vec.extract_subvector_to(local_indices.begin(),
729  local_indices.end(),
730  resid.data());
731 }
732 
733 
734 
735 template <typename VectorType>
736 void
738 {
739  // store local indices to write to vectors
740  locally_owned_dofs.fill_index_vector(local_indices);
741 
742  // scalars
743  nloc = locally_owned_dofs.n_elements();
744  ncv = additional_data.number_of_arnoldi_vectors;
745 
747 
748  // vectors
749  ldv = nloc;
750  v.resize(ldv * ncv, 0.0);
751 
752  resid.resize(nloc, 1.0);
753 
754  // work arrays for ARPACK
755  workd.resize(3 * nloc, 0.0);
756 
757  lworkl =
758  additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
759  workl.resize(lworkl, 0.);
760 
761  ldz = nloc;
762  z.resize(ldz * ncv, 0.); // TODO we actually need only ldz*nev
763 
764  // WORKEV Double precision work array of dimension 3*NCV.
765  lworkev = additional_data.symmetric ? 0 /*not used in symmetric case*/
766  :
767  3 * ncv;
768  workev.resize(lworkev, 0.);
769 
770  select.resize(ncv, 0);
771 }
772 
773 
774 
775 template <typename VectorType>
776 void
777 PArpackSolver<VectorType>::reinit(const IndexSet &locally_owned_dofs)
778 {
779  internal_reinit(locally_owned_dofs);
780 
781  // deal.II vectors:
782  src.reinit(locally_owned_dofs, mpi_communicator);
783  dst.reinit(locally_owned_dofs, mpi_communicator);
784  tmp.reinit(locally_owned_dofs, mpi_communicator);
785 }
786 
787 
788 
789 template <typename VectorType>
790 void
791 PArpackSolver<VectorType>::reinit(const VectorType &distributed_vector)
792 {
793  internal_reinit(distributed_vector.locally_owned_elements());
794 
795  // deal.II vectors:
796  src.reinit(distributed_vector);
797  dst.reinit(distributed_vector);
798  tmp.reinit(distributed_vector);
799 }
800 
801 
802 
803 template <typename VectorType>
804 void
805 PArpackSolver<VectorType>::reinit(const IndexSet &locally_owned_dofs,
806  const std::vector<IndexSet> &partitioning)
807 {
808  internal_reinit(locally_owned_dofs);
809 
810  // deal.II vectors:
811  src.reinit(partitioning, mpi_communicator);
812  dst.reinit(partitioning, mpi_communicator);
813  tmp.reinit(partitioning, mpi_communicator);
814 }
815 
816 
817 
818 template <typename VectorType>
819 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
820 void
821 PArpackSolver<VectorType>::solve(const MatrixType1 & A,
822  const MatrixType2 & B,
823  const INVERSE & inverse,
824  std::vector<std::complex<double>> &eigenvalues,
825  std::vector<VectorType> &eigenvectors,
826  const unsigned int n_eigenvalues)
827 {
828  std::vector<VectorType *> eigenvectors_ptr(eigenvectors.size());
829  for (unsigned int i = 0; i < eigenvectors.size(); ++i)
830  eigenvectors_ptr[i] = &eigenvectors[i];
831  solve(A, B, inverse, eigenvalues, eigenvectors_ptr, n_eigenvalues);
832 }
833 
834 
835 
836 template <typename VectorType>
837 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
838 void
839 PArpackSolver<VectorType>::solve(const MatrixType1 &system_matrix,
840  const MatrixType2 &mass_matrix,
841  const INVERSE & inverse,
842  std::vector<std::complex<double>> &eigenvalues,
843  std::vector<VectorType *> &eigenvectors,
844  const unsigned int n_eigenvalues)
845 {
846  if (additional_data.symmetric)
847  {
848  Assert(n_eigenvalues <= eigenvectors.size(),
849  PArpackExcInvalidEigenvectorSize(n_eigenvalues,
850  eigenvectors.size()));
851  }
852  else
853  Assert(n_eigenvalues + 1 <= eigenvectors.size(),
855  eigenvectors.size()));
856 
857  Assert(n_eigenvalues <= eigenvalues.size(),
858  PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
859 
860 
861  // use eigenvectors to get the problem size so that it is possible to
862  // employ LinearOperator for mass_matrix.
863  Assert(n_eigenvalues < eigenvectors[0]->size(),
865  eigenvectors[0]->size()));
866 
867  Assert(additional_data.number_of_arnoldi_vectors < eigenvectors[0]->size(),
869  additional_data.number_of_arnoldi_vectors, eigenvectors[0]->size()));
870 
871  Assert(additional_data.number_of_arnoldi_vectors > 2 * n_eigenvalues + 1,
873  additional_data.number_of_arnoldi_vectors, n_eigenvalues));
874 
875  int mode = additional_data.mode;
876 
877  // reverse communication parameter
878  // must be zero on the first call to pdnaupd
879  int ido = 0;
880 
881  // 'G' generalized eigenvalue problem
882  // 'I' standard eigenvalue problem
883  char bmat[2];
884  bmat[0] = (mode == 1) ? 'I' : 'G';
885  bmat[1] = '\0';
886 
887  // Specify the eigenvalues of interest, possible parameters:
888  // "LA" algebraically largest
889  // "SA" algebraically smallest
890  // "LM" largest magnitude
891  // "SM" smallest magnitude
892  // "LR" largest real part
893  // "SR" smallest real part
894  // "LI" largest imaginary part
895  // "SI" smallest imaginary part
896  // "BE" both ends of spectrum simultaneous
897  char which[3];
898  switch (additional_data.eigenvalue_of_interest)
899  {
901  std::strcpy(which, "LA");
902  break;
904  std::strcpy(which, "SA");
905  break;
906  case largest_magnitude:
907  std::strcpy(which, "LM");
908  break;
909  case smallest_magnitude:
910  std::strcpy(which, "SM");
911  break;
912  case largest_real_part:
913  std::strcpy(which, "LR");
914  break;
915  case smallest_real_part:
916  std::strcpy(which, "SR");
917  break;
919  std::strcpy(which, "LI");
920  break;
922  std::strcpy(which, "SI");
923  break;
924  case both_ends:
925  std::strcpy(which, "BE");
926  break;
927  }
928 
929  // tolerance for ARPACK
930  double tol = control().tolerance();
931 
932  // information to the routines
933  std::vector<int> iparam(11, 0);
934 
935  iparam[0] = 1;
936  // shift strategy: exact shifts with respect to the current Hessenberg matrix
937  // H.
938 
939  // maximum number of iterations
940  iparam[2] = control().max_steps();
941 
942  // Parpack currently works only for NB = 1
943  iparam[3] = 1;
944 
945  // Sets the mode of dsaupd:
946  // 1 is A*x=lambda*x, OP = A, B = I
947  // 2 is A*x = lambda*M*x, OP = inv[M]*A, B = M
948  // 3 is shift-invert mode, OP = inv[A-sigma*M]*M, B = M
949  // 4 is buckling mode,
950  // 5 is Cayley mode.
951 
952  iparam[6] = mode;
953  std::vector<int> ipntr(14, 0);
954 
955  // information out of the iteration
956  // If INFO .EQ. 0, a random initial residual vector is used.
957  // If INFO .NE. 0, RESID contains the initial residual vector,
958  // possibly from a previous run.
959  // Typical choices in this situation might be to use the final value
960  // of the starting vector from the previous eigenvalue calculation
961  int info = initial_vector_provided ? 1 : 0;
962 
963  // Number of eigenvalues of OP to be computed. 0 < NEV < N.
964  int nev = n_eigenvalues;
965  int n_inside_arpack = nloc;
966 
967  // IDO = 99: done
968  while (ido != 99)
969  {
970  // call of ARPACK pdnaupd routine
971  if (additional_data.symmetric)
972  pdsaupd_(&mpi_communicator_fortran,
973  &ido,
974  bmat,
975  &n_inside_arpack,
976  which,
977  &nev,
978  &tol,
979  resid.data(),
980  &ncv,
981  v.data(),
982  &ldv,
983  iparam.data(),
984  ipntr.data(),
985  workd.data(),
986  workl.data(),
987  &lworkl,
988  &info);
989  else
990  pdnaupd_(&mpi_communicator_fortran,
991  &ido,
992  bmat,
993  &n_inside_arpack,
994  which,
995  &nev,
996  &tol,
997  resid.data(),
998  &ncv,
999  v.data(),
1000  &ldv,
1001  iparam.data(),
1002  ipntr.data(),
1003  workd.data(),
1004  workl.data(),
1005  &lworkl,
1006  &info);
1007 
1008  AssertThrow(info == 0, PArpackExcInfoPdnaupd(info));
1009 
1010  // if we converge, we shall not modify anything in work arrays!
1011  if (ido == 99)
1012  break;
1013 
1014  // IPNTR(1) is the pointer into WORKD for X,
1015  // IPNTR(2) is the pointer into WORKD for Y.
1016  const int shift_x = ipntr[0] - 1;
1017  const int shift_y = ipntr[1] - 1;
1018  Assert(shift_x >= 0, ::ExcInternalError());
1019  Assert(shift_x + nloc <= static_cast<int>(workd.size()),
1020  ::ExcInternalError());
1021  Assert(shift_y >= 0, ::ExcInternalError());
1022  Assert(shift_y + nloc <= static_cast<int>(workd.size()),
1023  ::ExcInternalError());
1024 
1025  src = 0.;
1026 
1027  // switch based on both ido and mode
1028  if ((ido == -1) || (ido == 1 && mode < 3))
1029  // compute Y = OP * X
1030  {
1031  src.add(nloc, local_indices.data(), workd.data() + shift_x);
1032  src.compress(VectorOperation::add);
1033 
1034  if (mode == 3)
1035  // OP = inv[K - sigma*M]*M
1036  {
1037  mass_matrix.vmult(tmp, src);
1038  inverse.vmult(dst, tmp);
1039  }
1040  else if (mode == 2)
1041  // OP = inv[M]*K
1042  {
1043  system_matrix.vmult(tmp, src);
1044  // store M*X in X
1045  tmp.extract_subvector_to(local_indices.begin(),
1046  local_indices.end(),
1047  workd.data() + shift_x);
1048  inverse.vmult(dst, tmp);
1049  }
1050  else if (mode == 1)
1051  {
1052  system_matrix.vmult(dst, src);
1053  }
1054  else
1055  AssertThrow(false, PArpackExcMode(mode));
1056  }
1057  else if (ido == 1 && mode >= 3)
1058  // compute Y = OP * X for mode 3, 4 and 5, where
1059  // the vector B * X is already available in WORKD(ipntr(3)).
1060  {
1061  const int shift_b_x = ipntr[2] - 1;
1062  Assert(shift_b_x >= 0, ::ExcInternalError());
1063  Assert(shift_b_x + nloc <= static_cast<int>(workd.size()),
1064  ::ExcInternalError());
1065 
1066  // B*X
1067  src.add(nloc, local_indices.data(), workd.data() + shift_b_x);
1068  src.compress(VectorOperation::add);
1069 
1070  // solving linear system
1071  Assert(mode == 3, ExcNotImplemented());
1072  inverse.vmult(dst, src);
1073  }
1074  else if (ido == 2)
1075  // compute Y = B * X
1076  {
1077  src.add(nloc, local_indices.data(), workd.data() + shift_x);
1078  src.compress(VectorOperation::add);
1079 
1080  // Multiplication with mass matrix M
1081  if (mode == 1)
1082  {
1083  dst = src;
1084  }
1085  else
1086  // mode 2,3 and 5 have B=M
1087  {
1088  mass_matrix.vmult(dst, src);
1089  }
1090  }
1091  else
1092  AssertThrow(false, PArpackExcIdo(ido));
1093  // Note: IDO = 3 does not appear to be required for currently
1094  // implemented modes
1095 
1096  // store the result
1097  dst.extract_subvector_to(local_indices.begin(),
1098  local_indices.end(),
1099  workd.data() + shift_y);
1100  } // end of pd*aupd_ loop
1101 
1102  // 1 - compute eigenvectors,
1103  // 0 - only eigenvalues
1104  int rvec = 1;
1105 
1106  // which eigenvectors
1107  char howmany[4] = "All";
1108 
1109  std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1110  std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1111 
1112  // call of ARPACK pdneupd routine
1113  if (additional_data.symmetric)
1114  pdseupd_(&mpi_communicator_fortran,
1115  &rvec,
1116  howmany,
1117  select.data(),
1118  eigenvalues_real.data(),
1119  z.data(),
1120  &ldz,
1121  &sigmar,
1122  bmat,
1123  &n_inside_arpack,
1124  which,
1125  &nev,
1126  &tol,
1127  resid.data(),
1128  &ncv,
1129  v.data(),
1130  &ldv,
1131  iparam.data(),
1132  ipntr.data(),
1133  workd.data(),
1134  workl.data(),
1135  &lworkl,
1136  &info);
1137  else
1138  pdneupd_(&mpi_communicator_fortran,
1139  &rvec,
1140  howmany,
1141  select.data(),
1142  eigenvalues_real.data(),
1143  eigenvalues_im.data(),
1144  v.data(),
1145  &ldz,
1146  &sigmar,
1147  &sigmai,
1148  workev.data(),
1149  bmat,
1150  &n_inside_arpack,
1151  which,
1152  &nev,
1153  &tol,
1154  resid.data(),
1155  &ncv,
1156  v.data(),
1157  &ldv,
1158  iparam.data(),
1159  ipntr.data(),
1160  workd.data(),
1161  workl.data(),
1162  &lworkl,
1163  &info);
1164 
1165  if (info == 1)
1166  {
1167  AssertThrow(false, PArpackExcInfoMaxIt(control().max_steps()));
1168  }
1169  else if (info == 3)
1170  {
1171  AssertThrow(false, PArpackExcNoShifts(1));
1172  }
1173  else if (info != 0)
1174  {
1175  AssertThrow(false, PArpackExcInfoPdneupd(info));
1176  }
1177 
1178  for (int i = 0; i < nev; ++i)
1179  {
1180  (*eigenvectors[i]) = 0.0;
1181  AssertIndexRange(i * nloc + nloc, v.size() + 1);
1182 
1183  eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1184  eigenvectors[i]->compress(VectorOperation::add);
1185  }
1186 
1187  for (size_type i = 0; i < n_eigenvalues; ++i)
1188  eigenvalues[i] =
1189  std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1190 
1191  // Throw an error if the solver did not converge.
1192  AssertThrow(iparam[4] >= static_cast<int>(n_eigenvalues),
1193  PArpackExcConvergedEigenvectors(n_eigenvalues, iparam[4]));
1194 
1195  // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
1196  // with respect to a semi-inner product defined by M.
1197 
1198  // resid likely contains residual with respect to M-norm.
1199  {
1200  tmp = 0.0;
1201  tmp.add(nloc, local_indices.data(), resid.data());
1202  solver_control.check(iparam[2], tmp.l2_norm());
1203  }
1204 }
1205 
1206 
1207 
1208 template <typename VectorType>
1209 SolverControl &
1211 {
1212  return solver_control;
1213 }
1214 
1215 DEAL_II_NAMESPACE_CLOSE
1216 
1217 
1218 #endif
1219 #endif
types::global_dof_index size_type
void Tvmult(VectorType &dst, const VectorType &src) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
static ::ExceptionBase & PArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
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)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1637
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:1519
const AdditionalData additional_data
void set_shift(const std::complex< double > sigma)
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
bool initial_vector_provided
SolverControl & control() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
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)
VectorType src
static ::ExceptionBase & PArpackExcInfoMaxIt(int arg1)
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:476
std::vector< double > v
std::vector< double > resid
unsigned int global_dof_index
Definition: types.h:89
void vmult(VectorType &dst, const VectorType &src) const
static ::ExceptionBase & PArpackExcIdo(int arg1)
void internal_reinit(const IndexSet &locally_owned_dofs)
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
unsigned int max_steps() const
static ::ExceptionBase & ExcNotImplemented()
std::size_t memory_consumption() const
Shift(const MatrixType &A, const MatrixType &B, const double sigma)
void reinit(const IndexSet &locally_owned_dofs)
static ::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
std::vector< double > workl
size_type n_elements() const
Definition: index_set.h:1799
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)