Reference documentation for deal.II version GIT relicensing-216-gcda843ca70 2024-03-28 12:20:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_parpack_solver_h
16#define dealii_parpack_solver_h
17
18#include <deal.II/base/config.h>
19
23
26
27#include <cstring>
28
29
30#ifdef DEAL_II_ARPACK_WITH_PARPACK
31
33
34extern "C"
35{
36 // http://www.mathkeisan.com/usersguide/man/pdnaupd.html
37 void
38 pdnaupd_(MPI_Fint *comm,
39 int *ido,
40 char *bmat,
41 int *n,
42 char *which,
43 int *nev,
44 double *tol,
45 double *resid,
46 int *ncv,
47 double *v,
48 int *nloc,
49 int *iparam,
50 int *ipntr,
51 double *workd,
52 double *workl,
53 int *lworkl,
54 int *info);
55
56 // http://www.mathkeisan.com/usersguide/man/pdsaupd.html
57 void
58 pdsaupd_(MPI_Fint *comm,
59 int *ido,
60 char *bmat,
61 int *n,
62 char *which,
63 int *nev,
64 double *tol,
65 double *resid,
66 int *ncv,
67 double *v,
68 int *nloc,
69 int *iparam,
70 int *ipntr,
71 double *workd,
72 double *workl,
73 int *lworkl,
74 int *info);
75
76 // http://www.mathkeisan.com/usersguide/man/pdneupd.html
77 void
78 pdneupd_(MPI_Fint *comm,
79 int *rvec,
80 char *howmany,
81 int *select,
82 double *d,
83 double *di,
84 double *z,
85 int *ldz,
86 double *sigmar,
87 double *sigmai,
88 double *workev,
89 char *bmat,
90 int *n,
91 char *which,
92 int *nev,
93 double *tol,
94 double *resid,
95 int *ncv,
96 double *v,
97 int *nloc,
98 int *iparam,
99 int *ipntr,
100 double *workd,
101 double *workl,
102 int *lworkl,
103 int *info);
104
105 // http://www.mathkeisan.com/usersguide/man/pdseupd.html
106 void
107 pdseupd_(MPI_Fint *comm,
108 int *rvec,
109 char *howmany,
110 int *select,
111 double *d,
112 double *z,
113 int *ldz,
114 double *sigmar,
115 char *bmat,
116 int *n,
117 char *which,
118 int *nev,
119 double *tol,
120 double *resid,
121 int *ncv,
122 double *v,
123 int *nloc,
124 int *iparam,
125 int *ipntr,
126 double *workd,
127 double *workl,
128 int *lworkl,
129 int *info);
130
131 // other resources:
132 // http://acts.nersc.gov/superlu/example5/pnslac.c.html
133 // https://github.com/phpisciuneri/tijo/blob/master/dvr_parpack.cpp
134}
135
209template <typename VectorType>
211{
212public:
217
269
275 {
276 const unsigned int number_of_arnoldi_vectors;
278 const bool symmetric;
279 const int mode;
281 const unsigned int number_of_arnoldi_vectors = 15,
283 const bool symmetric = false,
284 const int mode = 3);
285 };
286
291 control() const;
292
298 const AdditionalData &data = AdditionalData());
299
303 void
304 reinit(const IndexSet &locally_owned_dofs);
305
312 void
313 reinit(const IndexSet &locally_owned_dofs,
314 const std::vector<IndexSet> &partitioning);
315
319 void
320 reinit(const VectorType &distributed_vector);
321
325 void
326 set_initial_vector(const VectorType &vec);
327
336 void
337 set_shift(const std::complex<double> sigma);
338
348 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
349 void
350 solve(const MatrixType1 &A,
351 const MatrixType2 &B,
352 const INVERSE &inverse,
353 std::vector<std::complex<double>> &eigenvalues,
354 std::vector<VectorType> &eigenvectors,
355 const unsigned int n_eigenvalues);
356
360 template <typename MatrixType1, typename MatrixType2, typename INVERSE>
361 void
362 solve(const MatrixType1 &A,
363 const MatrixType2 &B,
364 const INVERSE &inverse,
365 std::vector<std::complex<double>> &eigenvalues,
366 std::vector<VectorType *> &eigenvectors,
367 const unsigned int n_eigenvalues);
368
372 std::size_t
373 memory_consumption() const;
374
375protected:
381
386
387 // keep MPI communicator non-const as Arpack functions are not const either:
388
393
398
399 // C++98 guarantees that the elements of a vector are stored contiguously
400
405
409 std::vector<double> workl;
410
414 std::vector<double> workd;
415
419 int nloc;
420
424 int ncv;
425
426
430 int ldv;
431
436 std::vector<double> v;
437
442
447 std::vector<double> resid;
448
452 int ldz;
453
459 std::vector<double> z;
460
465
469 std::vector<double> workev;
470
474 std::vector<int> select;
475
479 VectorType src, dst, tmp;
480
484 std::vector<types::global_dof_index> local_indices;
485
489 double sigmar;
490
494 double sigmai;
495
496private:
503 void
504 internal_reinit(const IndexSet &locally_owned_dofs);
505
510 int,
511 int,
512 << arg1 << " eigenpairs were requested, but only " << arg2
513 << " converged");
514
516 int,
517 int,
518 << "Number of wanted eigenvalues " << arg1
519 << " is larger that the size of the matrix " << arg2);
520
522 int,
523 int,
524 << "Number of wanted eigenvalues " << arg1
525 << " is larger that the size of eigenvectors " << arg2);
526
529 int,
530 int,
531 << "To store the real and complex parts of " << arg1
532 << " eigenvectors in real-valued vectors, their size (currently set to "
533 << arg2 << ") should be greater than or equal to " << arg1 + 1);
534
536 int,
537 int,
538 << "Number of wanted eigenvalues " << arg1
539 << " is larger that the size of eigenvalues " << arg2);
540
542 int,
543 int,
544 << "Number of Arnoldi vectors " << arg1
545 << " is larger that the size of the matrix " << arg2);
546
548 int,
549 int,
550 << "Number of Arnoldi vectors " << arg1
551 << " is too small to obtain " << arg2 << " eigenvalues");
552
554 int,
555 << "This ido " << arg1
556 << " is not supported. Check documentation of ARPACK");
557
559 int,
560 << "This mode " << arg1
561 << " is not supported. Check documentation of ARPACK");
562
564 int,
565 << "Error with Pdnaupd, info " << arg1
566 << ". Check documentation of ARPACK");
567
569 int,
570 << "Error with Pdneupd, info " << arg1
571 << ". Check documentation of ARPACK");
572
574 int,
575 << "Maximum number " << arg1 << " of iterations reached.");
576
578 int,
579 << "No shifts could be applied during implicit"
580 << " Arnoldi update, try increasing the number of"
581 << " Arnoldi vectors.");
582};
583
584
585
586template <typename VectorType>
587std::size_t
589{
591 (workl.size() + workd.size() + v.size() + resid.size() + z.size() +
592 workev.size()) +
593 src.memory_consumption() + dst.memory_consumption() +
594 tmp.memory_consumption() +
596 local_indices.size();
597}
598
599
600
601template <typename VectorType>
603 const unsigned int number_of_arnoldi_vectors,
604 const WhichEigenvalues eigenvalue_of_interest,
605 const bool symmetric,
606 const int mode)
607 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
608 , eigenvalue_of_interest(eigenvalue_of_interest)
609 , symmetric(symmetric)
610 , mode(mode)
611{
612 // Check for possible options for symmetric problems
613 if (symmetric)
614 {
615 Assert(
618 "'largest real part' can only be used for non-symmetric problems!"));
619 Assert(
622 "'smallest real part' can only be used for non-symmetric problems!"));
623 Assert(
626 "'largest imaginary part' can only be used for non-symmetric problems!"));
627 Assert(
630 "'smallest imaginary part' can only be used for non-symmetric problems!"));
631 }
632 Assert(mode >= 1 && mode <= 3,
633 ExcMessage("Currently, only modes 1, 2 and 3 are supported."));
634}
635
636
637
638template <typename VectorType>
641 const AdditionalData &data)
643 , additional_data(data)
646 , lworkl(0)
647 , nloc(0)
648 , ncv(0)
649 , ldv(0)
651 , ldz(0)
652 , lworkev(0)
653 , sigmar(0.0)
654 , sigmai(0.0)
655{}
656
657
658
659template <typename VectorType>
660void
661PArpackSolver<VectorType>::set_shift(const std::complex<double> sigma)
662{
663 sigmar = sigma.real();
664 sigmai = sigma.imag();
665}
666
667
668
669template <typename VectorType>
670void
672{
673 initial_vector_provided = true;
674 Assert(resid.size() == local_indices.size(),
675 ExcDimensionMismatch(resid.size(), local_indices.size()));
676 vec.extract_subvector_to(local_indices.begin(),
677 local_indices.end(),
678 resid.data());
679}
680
681
682
683template <typename VectorType>
684void
686{
687 // store local indices to write to vectors
688 local_indices = locally_owned_dofs.get_index_vector();
689
690 // scalars
691 nloc = locally_owned_dofs.n_elements();
692 ncv = additional_data.number_of_arnoldi_vectors;
693
694 AssertDimension(local_indices.size(), nloc);
695
696 // vectors
697 ldv = nloc;
698 v.resize(ldv * ncv, 0.0);
699
700 resid.resize(nloc, 1.0);
701
702 // work arrays for ARPACK
703 workd.resize(3 * nloc, 0.0);
704
705 lworkl =
706 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
707 workl.resize(lworkl, 0.);
708
709 ldz = nloc;
710 z.resize(ldz * ncv, 0.); // TODO we actually need only ldz*nev
711
712 // WORKEV Double precision work array of dimension 3*NCV.
713 lworkev = additional_data.symmetric ? 0 /*not used in symmetric case*/
714 :
715 3 * ncv;
716 workev.resize(lworkev, 0.);
717
718 select.resize(ncv, 0);
719}
720
721
722
723template <typename VectorType>
724void
726{
727 internal_reinit(locally_owned_dofs);
728
729 // deal.II vectors:
730 src.reinit(locally_owned_dofs, mpi_communicator);
731 dst.reinit(locally_owned_dofs, mpi_communicator);
732 tmp.reinit(locally_owned_dofs, mpi_communicator);
733}
734
735
736
737template <typename VectorType>
738void
739PArpackSolver<VectorType>::reinit(const VectorType &distributed_vector)
740{
741 internal_reinit(distributed_vector.locally_owned_elements());
742
743 // deal.II vectors:
744 src.reinit(distributed_vector);
745 dst.reinit(distributed_vector);
746 tmp.reinit(distributed_vector);
747}
748
749
750
751template <typename VectorType>
752void
754 const std::vector<IndexSet> &partitioning)
755{
756 internal_reinit(locally_owned_dofs);
757
758 // deal.II vectors:
759 src.reinit(partitioning, mpi_communicator);
760 dst.reinit(partitioning, mpi_communicator);
761 tmp.reinit(partitioning, mpi_communicator);
762}
763
764
765
766template <typename VectorType>
767template <typename MatrixType1, typename MatrixType2, typename INVERSE>
768void
770 const MatrixType2 &B,
771 const INVERSE &inverse,
772 std::vector<std::complex<double>> &eigenvalues,
773 std::vector<VectorType> &eigenvectors,
774 const unsigned int n_eigenvalues)
775{
776 std::vector<VectorType *> eigenvectors_ptr(eigenvectors.size());
777 for (unsigned int i = 0; i < eigenvectors.size(); ++i)
778 eigenvectors_ptr[i] = &eigenvectors[i];
779 solve(A, B, inverse, eigenvalues, eigenvectors_ptr, n_eigenvalues);
780}
781
782
783
784template <typename VectorType>
785template <typename MatrixType1, typename MatrixType2, typename INVERSE>
786void
787PArpackSolver<VectorType>::solve(const MatrixType1 &system_matrix,
788 const MatrixType2 &mass_matrix,
789 const INVERSE &inverse,
790 std::vector<std::complex<double>> &eigenvalues,
791 std::vector<VectorType *> &eigenvectors,
792 const unsigned int n_eigenvalues)
793{
794 if (additional_data.symmetric)
795 {
796 Assert(n_eigenvalues <= eigenvectors.size(),
797 PArpackExcInvalidEigenvectorSize(n_eigenvalues,
798 eigenvectors.size()));
799 }
800 else
801 Assert(n_eigenvalues + 1 <= eigenvectors.size(),
802 PArpackExcInvalidEigenvectorSizeNonsymmetric(n_eigenvalues,
803 eigenvectors.size()));
804
805 Assert(n_eigenvalues <= eigenvalues.size(),
806 PArpackExcInvalidEigenvalueSize(n_eigenvalues, eigenvalues.size()));
807
808
809 // use eigenvectors to get the problem size so that it is possible to
810 // employ LinearOperator for mass_matrix.
811 Assert(n_eigenvalues < eigenvectors[0]->size(),
812 PArpackExcInvalidNumberofEigenvalues(n_eigenvalues,
813 eigenvectors[0]->size()));
814
815 Assert(additional_data.number_of_arnoldi_vectors < eigenvectors[0]->size(),
816 PArpackExcInvalidNumberofArnoldiVectors(
817 additional_data.number_of_arnoldi_vectors, eigenvectors[0]->size()));
818
819 Assert(additional_data.number_of_arnoldi_vectors > 2 * n_eigenvalues + 1,
820 PArpackExcSmallNumberofArnoldiVectors(
821 additional_data.number_of_arnoldi_vectors, n_eigenvalues));
822
823 int mode = additional_data.mode;
824
825 // reverse communication parameter
826 // must be zero on the first call to pdnaupd
827 int ido = 0;
828
829 // 'G' generalized eigenvalue problem
830 // 'I' standard eigenvalue problem
831 char bmat[2];
832 bmat[0] = (mode == 1) ? 'I' : 'G';
833 bmat[1] = '\0';
834
835 // Specify the eigenvalues of interest, possible parameters:
836 // "LA" algebraically largest
837 // "SA" algebraically smallest
838 // "LM" largest magnitude
839 // "SM" smallest magnitude
840 // "LR" largest real part
841 // "SR" smallest real part
842 // "LI" largest imaginary part
843 // "SI" smallest imaginary part
844 // "BE" both ends of spectrum simultaneous
845 char which[3];
846 switch (additional_data.eigenvalue_of_interest)
847 {
848 case algebraically_largest:
849 std::strcpy(which, "LA");
850 break;
851 case algebraically_smallest:
852 std::strcpy(which, "SA");
853 break;
854 case largest_magnitude:
855 std::strcpy(which, "LM");
856 break;
857 case smallest_magnitude:
858 std::strcpy(which, "SM");
859 break;
860 case largest_real_part:
861 std::strcpy(which, "LR");
862 break;
863 case smallest_real_part:
864 std::strcpy(which, "SR");
865 break;
866 case largest_imaginary_part:
867 std::strcpy(which, "LI");
868 break;
869 case smallest_imaginary_part:
870 std::strcpy(which, "SI");
871 break;
872 case both_ends:
873 std::strcpy(which, "BE");
874 break;
875 }
876
877 // tolerance for ARPACK
878 double tol = control().tolerance();
879
880 // information to the routines
881 std::vector<int> iparam(11, 0);
882
883 iparam[0] = 1;
884 // shift strategy: exact shifts with respect to the current Hessenberg matrix
885 // H.
886
887 // maximum number of iterations
888 iparam[2] = control().max_steps();
889
890 // Parpack currently works only for NB = 1
891 iparam[3] = 1;
892
893 // Sets the mode of dsaupd:
894 // 1 is A*x=lambda*x, OP = A, B = I
895 // 2 is A*x = lambda*M*x, OP = inv[M]*A, B = M
896 // 3 is shift-invert mode, OP = inv[A-sigma*M]*M, B = M
897 // 4 is buckling mode,
898 // 5 is Cayley mode.
899
900 iparam[6] = mode;
901 std::vector<int> ipntr(14, 0);
902
903 // information out of the iteration
904 // If INFO .EQ. 0, a random initial residual vector is used.
905 // If INFO .NE. 0, RESID contains the initial residual vector,
906 // possibly from a previous run.
907 // Typical choices in this situation might be to use the final value
908 // of the starting vector from the previous eigenvalue calculation
909 int info = initial_vector_provided ? 1 : 0;
910
911 // Number of eigenvalues of OP to be computed. 0 < NEV < N.
912 int nev = n_eigenvalues;
913 int n_inside_arpack = nloc;
914
915 // IDO = 99: done
916 while (ido != 99)
917 {
918 // call of ARPACK pdnaupd routine
919 if (additional_data.symmetric)
920 pdsaupd_(&mpi_communicator_fortran,
921 &ido,
922 bmat,
923 &n_inside_arpack,
924 which,
925 &nev,
926 &tol,
927 resid.data(),
928 &ncv,
929 v.data(),
930 &ldv,
931 iparam.data(),
932 ipntr.data(),
933 workd.data(),
934 workl.data(),
935 &lworkl,
936 &info);
937 else
938 pdnaupd_(&mpi_communicator_fortran,
939 &ido,
940 bmat,
941 &n_inside_arpack,
942 which,
943 &nev,
944 &tol,
945 resid.data(),
946 &ncv,
947 v.data(),
948 &ldv,
949 iparam.data(),
950 ipntr.data(),
951 workd.data(),
952 workl.data(),
953 &lworkl,
954 &info);
955
956 AssertThrow(info == 0, PArpackExcInfoPdnaupd(info));
957
958 // if we converge, we shall not modify anything in work arrays!
959 if (ido == 99)
960 break;
961
962 // IPNTR(1) is the pointer into WORKD for X,
963 // IPNTR(2) is the pointer into WORKD for Y.
964 const int shift_x = ipntr[0] - 1;
965 const int shift_y = ipntr[1] - 1;
966 Assert(shift_x >= 0, ::ExcInternalError());
967 Assert(shift_x + nloc <= static_cast<int>(workd.size()),
969 Assert(shift_y >= 0, ::ExcInternalError());
970 Assert(shift_y + nloc <= static_cast<int>(workd.size()),
972
973 src = 0.;
974
975 // switch based on both ido and mode
976 if ((ido == -1) || (ido == 1 && mode < 3))
977 // compute Y = OP * X
978 {
979 src.add(nloc, local_indices.data(), workd.data() + shift_x);
980 src.compress(VectorOperation::add);
981
982 if (mode == 3)
983 // OP = inv[K - sigma*M]*M
984 {
985 mass_matrix.vmult(tmp, src);
986 inverse.vmult(dst, tmp);
987 }
988 else if (mode == 2)
989 // OP = inv[M]*K
990 {
991 system_matrix.vmult(tmp, src);
992 // store M*X in X
993 tmp.extract_subvector_to(local_indices.begin(),
994 local_indices.end(),
995 workd.data() + shift_x);
996 inverse.vmult(dst, tmp);
997 }
998 else if (mode == 1)
999 {
1000 system_matrix.vmult(dst, src);
1001 }
1002 else
1003 AssertThrow(false, PArpackExcMode(mode));
1004 }
1005 else if (ido == 1 && mode >= 3)
1006 // compute Y = OP * X for mode 3, 4 and 5, where
1007 // the vector B * X is already available in WORKD(ipntr(3)).
1008 {
1009 const int shift_b_x = ipntr[2] - 1;
1010 Assert(shift_b_x >= 0, ::ExcInternalError());
1011 Assert(shift_b_x + nloc <= static_cast<int>(workd.size()),
1013
1014 // B*X
1015 src.add(nloc, local_indices.data(), workd.data() + shift_b_x);
1016 src.compress(VectorOperation::add);
1017
1018 // solving linear system
1019 Assert(mode == 3, ExcNotImplemented());
1020 inverse.vmult(dst, src);
1021 }
1022 else if (ido == 2)
1023 // compute Y = B * X
1024 {
1025 src.add(nloc, local_indices.data(), workd.data() + shift_x);
1026 src.compress(VectorOperation::add);
1027
1028 // Multiplication with mass matrix M
1029 if (mode == 1)
1030 {
1031 dst = src;
1032 }
1033 else
1034 // mode 2,3 and 5 have B=M
1035 {
1036 mass_matrix.vmult(dst, src);
1037 }
1038 }
1039 else
1040 AssertThrow(false, PArpackExcIdo(ido));
1041 // Note: IDO = 3 does not appear to be required for currently
1042 // implemented modes
1043
1044 // store the result
1045 dst.extract_subvector_to(local_indices.begin(),
1046 local_indices.end(),
1047 workd.data() + shift_y);
1048 } // end of pd*aupd_ loop
1049
1050 // 1 - compute eigenvectors,
1051 // 0 - only eigenvalues
1052 int rvec = 1;
1053
1054 // which eigenvectors
1055 char howmany[4] = "All";
1056
1057 std::vector<double> eigenvalues_real(n_eigenvalues + 1, 0.);
1058 std::vector<double> eigenvalues_im(n_eigenvalues + 1, 0.);
1059
1060 // call of ARPACK pdneupd routine
1061 if (additional_data.symmetric)
1062 pdseupd_(&mpi_communicator_fortran,
1063 &rvec,
1064 howmany,
1065 select.data(),
1066 eigenvalues_real.data(),
1067 z.data(),
1068 &ldz,
1069 &sigmar,
1070 bmat,
1071 &n_inside_arpack,
1072 which,
1073 &nev,
1074 &tol,
1075 resid.data(),
1076 &ncv,
1077 v.data(),
1078 &ldv,
1079 iparam.data(),
1080 ipntr.data(),
1081 workd.data(),
1082 workl.data(),
1083 &lworkl,
1084 &info);
1085 else
1086 pdneupd_(&mpi_communicator_fortran,
1087 &rvec,
1088 howmany,
1089 select.data(),
1090 eigenvalues_real.data(),
1091 eigenvalues_im.data(),
1092 v.data(),
1093 &ldz,
1094 &sigmar,
1095 &sigmai,
1096 workev.data(),
1097 bmat,
1098 &n_inside_arpack,
1099 which,
1100 &nev,
1101 &tol,
1102 resid.data(),
1103 &ncv,
1104 v.data(),
1105 &ldv,
1106 iparam.data(),
1107 ipntr.data(),
1108 workd.data(),
1109 workl.data(),
1110 &lworkl,
1111 &info);
1112
1113 if (info == 1)
1114 {
1115 AssertThrow(false, PArpackExcInfoMaxIt(control().max_steps()));
1116 }
1117 else if (info == 3)
1118 {
1119 AssertThrow(false, PArpackExcNoShifts(1));
1120 }
1121 else if (info != 0)
1122 {
1123 AssertThrow(false, PArpackExcInfoPdneupd(info));
1124 }
1125
1126 for (int i = 0; i < nev; ++i)
1127 {
1128 (*eigenvectors[i]) = 0.0;
1129 AssertIndexRange(i * nloc + nloc, v.size() + 1);
1130
1131 eigenvectors[i]->add(nloc, local_indices.data(), &v[i * nloc]);
1132 eigenvectors[i]->compress(VectorOperation::add);
1133 }
1134
1135 for (size_type i = 0; i < n_eigenvalues; ++i)
1136 eigenvalues[i] =
1137 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
1138
1139 // Throw an error if the solver did not converge.
1140 AssertThrow(iparam[4] >= static_cast<int>(n_eigenvalues),
1141 PArpackExcConvergedEigenvectors(n_eigenvalues, iparam[4]));
1142
1143 // both PDNAUPD and PDSAUPD compute eigenpairs of inv[A - sigma*M]*M
1144 // with respect to a semi-inner product defined by M.
1145
1146 // resid likely contains residual with respect to M-norm.
1147 {
1148 tmp = 0.0;
1149 tmp.add(nloc, local_indices.data(), resid.data());
1150 tmp.compress(VectorOperation::add);
1151 solver_control.check(iparam[2], tmp.l2_norm());
1152 }
1153}
1154
1155
1156
1157template <typename VectorType>
1160{
1161 return solver_control;
1162}
1163
1165
1166
1167#endif
1168#endif
size_type n_elements() const
Definition index_set.h:1923
std::vector< size_type > get_index_vector() const
Definition index_set.cc:923
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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:539
#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:516
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_v< T >, std::size_t > memory_consumption(const T &t)
unsigned int global_dof_index
Definition types.h:81
*braid_SplitCommworld & comm
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)