Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
parpack_solver.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2010 - 2023 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
24
27
28#include <cstring>
29
30
31#ifdef DEAL_II_ARPACK_WITH_PARPACK
32
34
35extern "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
210template <typename VectorType>
212{
213public:
218
229 {
269 };
270
276 {
277 const unsigned int number_of_arnoldi_vectors;
279 const bool symmetric;
280 const int mode;
282 const unsigned int number_of_arnoldi_vectors = 15,
284 const bool symmetric = false,
285 const int mode = 3);
286 };
287
292 control() const;
293
299 const AdditionalData &data = AdditionalData());
300
304 void
305 reinit(const IndexSet &locally_owned_dofs);
306
313 void
314 reinit(const IndexSet & locally_owned_dofs,
315 const std::vector<IndexSet> &partitioning);
316
320 void
321 reinit(const VectorType &distributed_vector);
322
326 void
327 set_initial_vector(const VectorType &vec);
328
337 void
338 set_shift(const std::complex<double> sigma);
339
349 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
350 void
351 solve(const MatrixType1 & A,
352 const MatrixType2 & B,
353 const INVERSE & inverse,
354 std::vector<std::complex<double>> &eigenvalues,
355 std::vector<VectorType> & eigenvectors,
356 const unsigned int n_eigenvalues);
357
361 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
362 void
363 solve(const MatrixType1 & A,
364 const MatrixType2 & B,
365 const INVERSE & inverse,
366 std::vector<std::complex<double>> &eigenvalues,
367 std::vector<VectorType *> & eigenvectors,
368 const unsigned int n_eigenvalues);
369
373 std::size_t
374 memory_consumption() const;
375
376protected:
382
387
388 // keep MPI communicator non-const as Arpack functions are not const either:
389
394
399
400 // C++98 guarantees that the elements of a vector are stored contiguously
401
406
410 std::vector<double> workl;
411
415 std::vector<double> workd;
416
420 int nloc;
421
425 int ncv;
426
427
431 int ldv;
432
437 std::vector<double> v;
438
443
448 std::vector<double> resid;
449
453 int ldz;
454
460 std::vector<double> z;
461
466
470 std::vector<double> workev;
471
475 std::vector<int> select;
476
480 VectorType src, dst, tmp;
481
485 std::vector<types::global_dof_index> local_indices;
486
490 double sigmar;
491
495 double sigmai;
496
497private:
504 void
505 internal_reinit(const IndexSet &locally_owned_dofs);
506
511 int,
512 int,
513 << arg1 << " eigenpairs were requested, but only " << arg2
514 << " converged");
515
517 int,
518 int,
519 << "Number of wanted eigenvalues " << arg1
520 << " is larger that the size of the matrix " << arg2);
521
523 int,
524 int,
525 << "Number of wanted eigenvalues " << arg1
526 << " is larger that the size of eigenvectors " << arg2);
527
530 int,
531 int,
532 << "To store the real and complex parts of " << arg1
533 << " eigenvectors in real-valued vectors, their size (currently set to "
534 << arg2 << ") should be greater than or equal to " << arg1 + 1);
535
537 int,
538 int,
539 << "Number of wanted eigenvalues " << arg1
540 << " is larger that the size of eigenvalues " << arg2);
541
543 int,
544 int,
545 << "Number of Arnoldi vectors " << arg1
546 << " is larger that the size of the matrix " << arg2);
547
549 int,
550 int,
551 << "Number of Arnoldi vectors " << arg1
552 << " is too small to obtain " << arg2 << " eigenvalues");
553
555 int,
556 << "This ido " << arg1
557 << " is not supported. Check documentation of ARPACK");
558
560 int,
561 << "This mode " << arg1
562 << " is not supported. Check documentation of ARPACK");
563
565 int,
566 << "Error with Pdnaupd, info " << arg1
567 << ". Check documentation of ARPACK");
568
570 int,
571 << "Error with Pdneupd, info " << arg1
572 << ". Check documentation of ARPACK");
573
575 int,
576 << "Maximum number " << arg1 << " of iterations reached.");
577
579 int,
580 << "No shifts could be applied during implicit"
581 << " Arnoldi update, try increasing the number of"
582 << " Arnoldi vectors.");
583};
584
585
586
587template <typename VectorType>
588std::size_t
590{
592 (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
593 workev.size()) +
594 src.memory_consumption() + dst.memory_consumption() +
595 tmp.memory_consumption() +
597 local_indices.size();
598}
599
600
601
602template <typename VectorType>
604 const unsigned int number_of_arnoldi_vectors,
605 const WhichEigenvalues eigenvalue_of_interest,
606 const bool symmetric,
607 const int mode)
608 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
609 , eigenvalue_of_interest(eigenvalue_of_interest)
610 , symmetric(symmetric)
611 , mode(mode)
612{
613 // Check for possible options for symmetric problems
614 if (symmetric)
615 {
616 Assert(
619 "'largest real part' can only be used for non-symmetric problems!"));
620 Assert(
623 "'smallest real part' can only be used for non-symmetric problems!"));
624 Assert(
627 "'largest imaginary part' can only be used for non-symmetric problems!"));
628 Assert(
631 "'smallest imaginary part' can only be used for non-symmetric problems!"));
632 }
633 Assert(mode >= 1 && mode <= 3,
634 ExcMessage("Currently, only modes 1, 2 and 3 are supported."));
635}
636
637
638
639template <typename VectorType>
642 const AdditionalData &data)
644 , additional_data(data)
647 , lworkl(0)
648 , nloc(0)
649 , ncv(0)
650 , ldv(0)
652 , ldz(0)
653 , lworkev(0)
654 , sigmar(0.0)
655 , sigmai(0.0)
656{}
657
658
659
660template <typename VectorType>
661void
662PArpackSolver<VectorType>::set_shift(const std::complex<double> sigma)
663{
664 sigmar = sigma.real();
665 sigmai = sigma.imag();
666}
667
668
669
670template <typename VectorType>
671void
673{
674 initial_vector_provided = true;
675 Assert(resid.size() == local_indices.size(),
676 ExcDimensionMismatch(resid.size(), local_indices.size()));
677 vec.extract_subvector_to(local_indices.begin(),
678 local_indices.end(),
679 resid.data());
680}
681
682
683
684template <typename VectorType>
685void
687{
688 // store local indices to write to vectors
689 local_indices = locally_owned_dofs.get_index_vector();
690
691 // scalars
692 nloc = locally_owned_dofs.n_elements();
693 ncv = additional_data.number_of_arnoldi_vectors;
694
695 AssertDimension(local_indices.size(), nloc);
696
697 // vectors
698 ldv = nloc;
699 v.resize(ldv * ncv, 0.0);
700
701 resid.resize(nloc, 1.0);
702
703 // work arrays for ARPACK
704 workd.resize(3 * nloc, 0.0);
705
706 lworkl =
707 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
708 workl.resize(lworkl, 0.);
709
710 ldz = nloc;
711 z.resize(ldz * ncv, 0.); // TODO we actually need only ldz*nev
712
713 // WORKEV Double precision work array of dimension 3*NCV.
714 lworkev = additional_data.symmetric ? 0 /*not used in symmetric case*/
715 :
716 3 * ncv;
717 workev.resize(lworkev, 0.);
718
719 select.resize(ncv, 0);
720}
721
722
723
724template <typename VectorType>
725void
727{
728 internal_reinit(locally_owned_dofs);
729
730 // deal.II vectors:
731 src.reinit(locally_owned_dofs, mpi_communicator);
732 dst.reinit(locally_owned_dofs, mpi_communicator);
733 tmp.reinit(locally_owned_dofs, mpi_communicator);
734}
735
736
737
738template <typename VectorType>
739void
740PArpackSolver<VectorType>::reinit(const VectorType &distributed_vector)
741{
742 internal_reinit(distributed_vector.locally_owned_elements());
743
744 // deal.II vectors:
745 src.reinit(distributed_vector);
746 dst.reinit(distributed_vector);
747 tmp.reinit(distributed_vector);
748}
749
750
751
752template <typename VectorType>
753void
755 const std::vector<IndexSet> &partitioning)
756{
757 internal_reinit(locally_owned_dofs);
758
759 // deal.II vectors:
760 src.reinit(partitioning, mpi_communicator);
761 dst.reinit(partitioning, mpi_communicator);
762 tmp.reinit(partitioning, mpi_communicator);
763}
764
765
766
767template <typename VectorType>
768template <typename MatrixType1, typename MatrixType2, typename INVERSE>
769void
771 const MatrixType2 & B,
772 const INVERSE & inverse,
773 std::vector<std::complex<double>> &eigenvalues,
774 std::vector<VectorType> &eigenvectors,
775 const unsigned int n_eigenvalues)
776{
777 std::vector<VectorType *> eigenvectors_ptr(eigenvectors.size());
778 for (unsigned int i = 0; i < eigenvectors.size(); ++i)
779 eigenvectors_ptr[i] = &eigenvectors[i];
780 solve(A, B, inverse, eigenvalues, eigenvectors_ptr, n_eigenvalues);
781}
782
783
784
785template <typename VectorType>
786template <typename MatrixType1, typename MatrixType2, typename INVERSE>
787void
788PArpackSolver<VectorType>::solve(const MatrixType1 &system_matrix,
789 const MatrixType2 &mass_matrix,
790 const INVERSE & inverse,
791 std::vector<std::complex<double>> &eigenvalues,
792 std::vector<VectorType *> &eigenvectors,
793 const unsigned int n_eigenvalues)
794{
795 if (additional_data.symmetric)
796 {
797 Assert(n_eigenvalues <= eigenvectors.size(),
798 PArpackExcInvalidEigenvectorSize(n_eigenvalues,
799 eigenvectors.size()));
800 }
801 else
802 Assert(n_eigenvalues + 1 <= eigenvectors.size(),
803 PArpackExcInvalidEigenvectorSizeNonsymmetric(n_eigenvalues,
804 eigenvectors.size()));
805
806 Assert(n_eigenvalues <= eigenvalues.size(),
807 PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
808
809
810 // use eigenvectors to get the problem size so that it is possible to
811 // employ LinearOperator for mass_matrix.
812 Assert(n_eigenvalues < eigenvectors[0]->size(),
813 PArpackExcInvalidNumberofEigenvalues(n_eigenvalues,
814 eigenvectors[0]->size()));
815
816 Assert(additional_data.number_of_arnoldi_vectors < eigenvectors[0]->size(),
817 PArpackExcInvalidNumberofArnoldiVectors(
818 additional_data.number_of_arnoldi_vectors, eigenvectors[0]->size()));
819
820 Assert(additional_data.number_of_arnoldi_vectors > 2 * n_eigenvalues + 1,
821 PArpackExcSmallNumberofArnoldiVectors(
822 additional_data.number_of_arnoldi_vectors, n_eigenvalues));
823
824 int mode = additional_data.mode;
825
826 // reverse communication parameter
827 // must be zero on the first call to pdnaupd
828 int ido = 0;
829
830 // 'G' generalized eigenvalue problem
831 // 'I' standard eigenvalue problem
832 char bmat[2];
833 bmat[0] = (mode == 1) ? 'I' : 'G';
834 bmat[1] = '\0';
835
836 // Specify the eigenvalues of interest, possible parameters:
837 // "LA" algebraically largest
838 // "SA" algebraically smallest
839 // "LM" largest magnitude
840 // "SM" smallest magnitude
841 // "LR" largest real part
842 // "SR" smallest real part
843 // "LI" largest imaginary part
844 // "SI" smallest imaginary part
845 // "BE" both ends of spectrum simultaneous
846 char which[3];
847 switch (additional_data.eigenvalue_of_interest)
848 {
849 case algebraically_largest:
850 std::strcpy(which, "LA");
851 break;
852 case algebraically_smallest:
853 std::strcpy(which, "SA");
854 break;
855 case largest_magnitude:
856 std::strcpy(which, "LM");
857 break;
858 case smallest_magnitude:
859 std::strcpy(which, "SM");
860 break;
861 case largest_real_part:
862 std::strcpy(which, "LR");
863 break;
864 case smallest_real_part:
865 std::strcpy(which, "SR");
866 break;
867 case largest_imaginary_part:
868 std::strcpy(which, "LI");
869 break;
870 case smallest_imaginary_part:
871 std::strcpy(which, "SI");
872 break;
873 case both_ends:
874 std::strcpy(which, "BE");
875 break;
876 }
877
878 // tolerance for ARPACK
879 double tol = control().tolerance();
880
881 // information to the routines
882 std::vector<int> iparam(11, 0);
883
884 iparam[0] = 1;
885 // shift strategy: exact shifts with respect to the current Hessenberg matrix
886 // H.
887
888 // maximum number of iterations
889 iparam[2] = control().max_steps();
890
891 // Parpack currently works only for NB = 1
892 iparam[3] = 1;
893
894 // Sets the mode of dsaupd:
895 // 1 is A*x=lambda*x, OP = A, B = I
896 // 2 is A*x = lambda*M*x, OP = inv[M]*A, B = M
897 // 3 is shift-invert mode, OP = inv[A-sigma*M]*M, B = M
898 // 4 is buckling mode,
899 // 5 is Cayley mode.
900
901 iparam[6] = mode;
902 std::vector<int> ipntr(14, 0);
903
904 // information out of the iteration
905 // If INFO .EQ. 0, a random initial residual vector is used.
906 // If INFO .NE. 0, RESID contains the initial residual vector,
907 // possibly from a previous run.
908 // Typical choices in this situation might be to use the final value
909 // of the starting vector from the previous eigenvalue calculation
910 int info = initial_vector_provided ? 1 : 0;
911
912 // Number of eigenvalues of OP to be computed. 0 < NEV < N.
913 int nev = n_eigenvalues;
914 int n_inside_arpack = nloc;
915
916 // IDO = 99: done
917 while (ido != 99)
918 {
919 // call of ARPACK pdnaupd routine
920 if (additional_data.symmetric)
921 pdsaupd_(&mpi_communicator_fortran,
922 &ido,
923 bmat,
924 &n_inside_arpack,
925 which,
926 &nev,
927 &tol,
928 resid.data(),
929 &ncv,
930 v.data(),
931 &ldv,
932 iparam.data(),
933 ipntr.data(),
934 workd.data(),
935 workl.data(),
936 &lworkl,
937 &info);
938 else
939 pdnaupd_(&mpi_communicator_fortran,
940 &ido,
941 bmat,
942 &n_inside_arpack,
943 which,
944 &nev,
945 &tol,
946 resid.data(),
947 &ncv,
948 v.data(),
949 &ldv,
950 iparam.data(),
951 ipntr.data(),
952 workd.data(),
953 workl.data(),
954 &lworkl,
955 &info);
956
957 AssertThrow(info == 0, PArpackExcInfoPdnaupd(info));
958
959 // if we converge, we shall not modify anything in work arrays!
960 if (ido == 99)
961 break;
962
963 // IPNTR(1) is the pointer into WORKD for X,
964 // IPNTR(2) is the pointer into WORKD for Y.
965 const int shift_x = ipntr[0] - 1;
966 const int shift_y = ipntr[1] - 1;
967 Assert(shift_x >= 0, ::ExcInternalError());
968 Assert(shift_x + nloc <= static_cast<int>(workd.size()),
970 Assert(shift_y >= 0, ::ExcInternalError());
971 Assert(shift_y + nloc <= static_cast<int>(workd.size()),
973
974 src = 0.;
975
976 // switch based on both ido and mode
977 if ((ido == -1) || (ido == 1 && mode < 3))
978 // compute Y = OP * X
979 {
980 src.add(nloc, local_indices.data(), workd.data() + shift_x);
981 src.compress(VectorOperation::add);
982
983 if (mode == 3)
984 // OP = inv[K - sigma*M]*M
985 {
986 mass_matrix.vmult(tmp, src);
987 inverse.vmult(dst, tmp);
988 }
989 else if (mode == 2)
990 // OP = inv[M]*K
991 {
992 system_matrix.vmult(tmp, src);
993 // store M*X in X
994 tmp.extract_subvector_to(local_indices.begin(),
995 local_indices.end(),
996 workd.data() + shift_x);
997 inverse.vmult(dst, tmp);
998 }
999 else if (mode == 1)
1000 {
1001 system_matrix.vmult(dst, src);
1002 }
1003 else
1004 AssertThrow(false, PArpackExcMode(mode));
1005 }
1006 else if (ido == 1 && mode >= 3)
1007 // compute Y = OP * X for mode 3, 4 and 5, where
1008 // the vector B * X is already available in WORKD(ipntr(3)).
1009 {
1010 const int shift_b_x = ipntr[2] - 1;
1011 Assert(shift_b_x >= 0, ::ExcInternalError());
1012 Assert(shift_b_x + nloc <= static_cast<int>(workd.size()),
1014
1015 // B*X
1016 src.add(nloc, local_indices.data(), workd.data() + shift_b_x);
1017 src.compress(VectorOperation::add);
1018
1019 // solving linear system
1020 Assert(mode == 3, ExcNotImplemented());
1021 inverse.vmult(dst, src);
1022 }
1023 else if (ido == 2)
1024 // compute Y = B * X
1025 {
1026 src.add(nloc, local_indices.data(), workd.data() + shift_x);
1027 src.compress(VectorOperation::add);
1028
1029 // Multiplication with mass matrix M
1030 if (mode == 1)
1031 {
1032 dst = src;
1033 }
1034 else
1035 // mode 2,3 and 5 have B=M
1036 {
1037 mass_matrix.vmult(dst, src);
1038 }
1039 }
1040 else
1041 AssertThrow(false, PArpackExcIdo(ido));
1042 // Note: IDO = 3 does not appear to be required for currently
1043 // implemented modes
1044
1045 // store the result
1046 dst.extract_subvector_to(local_indices.begin(),
1047 local_indices.end(),
1048 workd.data() + shift_y);
1049 } // end of pd*aupd_ loop
1050
1051 // 1 - compute eigenvectors,
1052 // 0 - only eigenvalues
1053 int rvec = 1;
1054
1055 // which eigenvectors
1056 char howmany[4] = "All";
1057
1058 std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1059 std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1060
1061 // call of ARPACK pdneupd routine
1062 if (additional_data.symmetric)
1063 pdseupd_(&mpi_communicator_fortran,
1064 &rvec,
1065 howmany,
1066 select.data(),
1067 eigenvalues_real.data(),
1068 z.data(),
1069 &ldz,
1070 &sigmar,
1071 bmat,
1072 &n_inside_arpack,
1073 which,
1074 &nev,
1075 &tol,
1076 resid.data(),
1077 &ncv,
1078 v.data(),
1079 &ldv,
1080 iparam.data(),
1081 ipntr.data(),
1082 workd.data(),
1083 workl.data(),
1084 &lworkl,
1085 &info);
1086 else
1087 pdneupd_(&mpi_communicator_fortran,
1088 &rvec,
1089 howmany,
1090 select.data(),
1091 eigenvalues_real.data(),
1092 eigenvalues_im.data(),
1093 v.data(),
1094 &ldz,
1095 &sigmar,
1096 &sigmai,
1097 workev.data(),
1098 bmat,
1099 &n_inside_arpack,
1100 which,
1101 &nev,
1102 &tol,
1103 resid.data(),
1104 &ncv,
1105 v.data(),
1106 &ldv,
1107 iparam.data(),
1108 ipntr.data(),
1109 workd.data(),
1110 workl.data(),
1111 &lworkl,
1112 &info);
1113
1114 if (info == 1)
1115 {
1116 AssertThrow(false, PArpackExcInfoMaxIt(control().max_steps()));
1117 }
1118 else if (info == 3)
1119 {
1120 AssertThrow(false, PArpackExcNoShifts(1));
1121 }
1122 else if (info != 0)
1123 {
1124 AssertThrow(false, PArpackExcInfoPdneupd(info));
1125 }
1126
1127 for (int i = 0; i < nev; ++i)
1128 {
1129 (*eigenvectors[i]) = 0.0;
1130 AssertIndexRange(i * nloc + nloc, v.size() + 1);
1131
1132 eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1133 eigenvectors[i]->compress(VectorOperation::add);
1134 }
1135
1136 for (size_type i = 0; i < n_eigenvalues; ++i)
1137 eigenvalues[i] =
1138 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1139
1140 // Throw an error if the solver did not converge.
1141 AssertThrow(iparam[4] >= static_cast<int>(n_eigenvalues),
1142 PArpackExcConvergedEigenvectors(n_eigenvalues, iparam[4]));
1143
1144 // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
1145 // with respect to a semi-inner product defined by M.
1146
1147 // resid likely contains residual with respect to M-norm.
1148 {
1149 tmp = 0.0;
1150 tmp.add(nloc, local_indices.data(), resid.data());
1151 tmp.compress(VectorOperation::add);
1152 solver_control.check(iparam[2], tmp.l2_norm());
1153 }
1154}
1155
1156
1157
1158template <typename VectorType>
1161{
1162 return solver_control;
1163}
1164
1166
1167
1168#endif
1169#endif
size_type n_elements() const
Definition index_set.h:1816
std::vector< size_type > get_index_vector() const
Definition index_set.cc:690
void set_initial_vector(const VectorType &vec)
SolverControl & solver_control
MPI_Comm mpi_communicator
std::vector< double > z
std::vector< double > v
SolverControl & control() const
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)
PArpackSolver(SolverControl &control, const MPI_Comm mpi_communicator, const AdditionalData &data=AdditionalData())
std::vector< types::global_dof_index > local_indices
std::vector< double > workl
std::vector< int > select
std::size_t memory_consumption() const
MPI_Fint mpi_communicator_fortran
void set_shift(const std::complex< double > sigma)
void reinit(const IndexSet &locally_owned_dofs)
std::vector< double > workd
const AdditionalData additional_data
bool initial_vector_provided
void internal_reinit(const IndexSet &locally_owned_dofs)
std::vector< double > resid
std::vector< double > workev
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & PArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
static ::ExceptionBase & PArpackExcInvalidEigenvalueSize(int arg1, int arg2)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & PArpackExcNoShifts(int arg1)
static ::ExceptionBase & PArpackExcInvalidEigenvectorSize(int arg1, int arg2)
static ::ExceptionBase & PArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & PArpackExcInfoPdneupd(int arg1)
#define AssertIndexRange(index, range)
static ::ExceptionBase & PArpackExcIdo(int arg1)
static ::ExceptionBase & PArpackExcConvergedEigenvectors(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & PArpackExcInfoMaxIt(int arg1)
static ::ExceptionBase & PArpackExcMode(int arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & PArpackExcInfoPdnaupd(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & PArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & PArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
unsigned int global_dof_index
Definition types.h:82
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)
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)
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)
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)
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false, const int mode=3)
const unsigned int number_of_arnoldi_vectors
const WhichEigenvalues eigenvalue_of_interest
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
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)
const MPI_Comm comm