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
arpack_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_arpack_solver_h
17#define dealii_arpack_solver_h
18
19#include <deal.II/base/config.h>
20
22
24
25#include <cstring>
26
27
28#ifdef DEAL_II_WITH_ARPACK
29
31
32
33extern "C" void
34dnaupd_(int * ido,
35 char * bmat,
36 unsigned int *n,
37 char * which,
38 unsigned int *nev,
39 const double *tol,
40 double * resid,
41 int * ncv,
42 double * v,
43 int * ldv,
44 int * iparam,
45 int * ipntr,
46 double * workd,
47 double * workl,
48 int * lworkl,
49 int * info);
50
51extern "C" void
52dsaupd_(int * ido,
53 char * bmat,
54 unsigned int *n,
55 char * which,
56 unsigned int *nev,
57 double * tol,
58 double * resid,
59 int * ncv,
60 double * v,
61 int * ldv,
62 int * iparam,
63 int * ipntr,
64 double * workd,
65 double * workl,
66 int * lworkl,
67 int * info);
68
69extern "C" void
70dneupd_(int * rvec,
71 char * howmany,
72 int * select,
73 double * d,
74 double * di,
75 double * z,
76 int * ldz,
77 double * sigmar,
78 double * sigmai,
79 double * workev,
80 char * bmat,
81 unsigned int *n,
82 char * which,
83 unsigned int *nev,
84 double * tol,
85 double * resid,
86 int * ncv,
87 double * v,
88 int * ldv,
89 int * iparam,
90 int * ipntr,
91 double * workd,
92 double * workl,
93 int * lworkl,
94 int * info);
95
96extern "C" void
97dseupd_(int * rvec,
98 char * howmany,
99 int * select,
100 double * d,
101 double * z,
102 int * ldz,
103 double * sigmar,
104 char * bmat,
105 unsigned int *n,
106 char * which,
107 unsigned int *nev,
108 double * tol,
109 double * resid,
110 int * ncv,
111 double * v,
112 int * ldv,
113 int * iparam,
114 int * ipntr,
115 double * workd,
116 double * workl,
117 int * lworkl,
118 int * info);
119
168{
169public:
174
175
181 {
221 };
222
227 {
233 explicit AdditionalData(
234 const unsigned int number_of_arnoldi_vectors = 15,
236 const bool symmetric = false);
237
243 const unsigned int number_of_arnoldi_vectors;
244
249
253 const bool symmetric;
254 };
255
260 control() const;
261
266 const AdditionalData &data = AdditionalData());
267
271 template <typename VectorType>
272 void
273 set_initial_vector(const VectorType &vec);
274
280 void
281 set_shift(const std::complex<double> sigma);
282
332 template <typename VectorType,
333 typename MatrixType1,
334 typename MatrixType2,
335 typename INVERSE>
336 void
337 solve(const MatrixType1 & A,
338 const MatrixType2 & B,
339 const INVERSE & inverse,
340 std::vector<std::complex<double>> &eigenvalues,
341 std::vector<VectorType> & eigenvectors,
342 const unsigned int n_eigenvalues = 0);
343
344protected:
350
355
360 std::vector<double> resid;
361
365 double sigmar;
366
370 double sigmai;
371
372
373private:
378 int,
379 int,
380 << "Number of wanted eigenvalues " << arg1
381 << " is larger that the size of the matrix " << arg2);
382
384 int,
385 int,
386 << "Number of wanted eigenvalues " << arg1
387 << " is larger that the size of eigenvectors " << arg2);
388
391 int,
392 int,
393 << "To store the real and complex parts of " << arg1
394 << " eigenvectors in real-valued vectors, their size (currently set to "
395 << arg2 << ") should be greater than or equal to " << arg1 + 1);
396
398 int,
399 int,
400 << "Number of wanted eigenvalues " << arg1
401 << " is larger that the size of eigenvalues " << arg2);
402
404 int,
405 int,
406 << "Number of Arnoldi vectors " << arg1
407 << " is larger that the size of the matrix " << arg2);
408
410 int,
411 int,
412 << "Number of Arnoldi vectors " << arg1
413 << " is too small to obtain " << arg2 << " eigenvalues");
414
416 int,
417 << "This ido " << arg1
418 << " is not supported. Check documentation of ARPACK");
419
421 int,
422 << "This mode " << arg1
423 << " is not supported. Check documentation of ARPACK");
424
426 int,
427 << "Error with dsaupd, info " << arg1
428 << ". Check documentation of ARPACK");
429
431 int,
432 << "Error with dnaupd, info " << arg1
433 << ". Check documentation of ARPACK");
434
436 int,
437 << "Error with dseupd, info " << arg1
438 << ". Check documentation of ARPACK");
439
441 int,
442 << "Error with dneupd, info " << arg1
443 << ". Check documentation of ARPACK");
444
446 int,
447 << "Maximum number " << arg1 << " of iterations reached.");
448
450 "No shifts could be applied during implicit"
451 " Arnoldi update, try increasing the number of"
452 " Arnoldi vectors.");
453};
454
455
457 const unsigned int number_of_arnoldi_vectors,
458 const WhichEigenvalues eigenvalue_of_interest,
459 const bool symmetric)
460 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
461 , eigenvalue_of_interest(eigenvalue_of_interest)
462 , symmetric(symmetric)
463{
464 // Check for possible options for symmetric problems
465 if (symmetric)
466 {
467 Assert(
470 "'largest real part' can only be used for non-symmetric problems!"));
471 Assert(
474 "'smallest real part' can only be used for non-symmetric problems!"));
475 Assert(
478 "'largest imaginary part' can only be used for non-symmetric problems!"));
479 Assert(
482 "'smallest imaginary part' can only be used for non-symmetric problems!"));
483 }
484 // Check for possible options for asymmetric problems
485 else
486 {
487 Assert(
490 "'largest algebraic part' can only be used for symmetric problems!"));
491 Assert(
494 "'smallest algebraic part' can only be used for symmetric problems!"));
497 "'both ends' can only be used for symmetric problems!"));
498 }
499}
500
501
503 const AdditionalData &data)
505 , additional_data(data)
507 , sigmar(0.0)
508 , sigmai(0.0)
509{}
510
511
512
513inline void
514ArpackSolver::set_shift(const std::complex<double> sigma)
515{
516 sigmar = sigma.real();
517 sigmai = sigma.imag();
518}
519
520
521
522template <typename VectorType>
523inline void
525{
527 resid.resize(vec.size());
528 for (size_type i = 0; i < vec.size(); ++i)
529 resid[i] = vec[i];
530}
531
532
533template <typename VectorType,
534 typename MatrixType1,
535 typename MatrixType2,
536 typename INVERSE>
537inline void
538ArpackSolver::solve(const MatrixType1 & /*system_matrix*/,
539 const MatrixType2 & mass_matrix,
540 const INVERSE & inverse,
541 std::vector<std::complex<double>> &eigenvalues,
542 std::vector<VectorType> & eigenvectors,
543 const unsigned int n_eigenvalues)
544{
545 // Problem size
546 unsigned int n = eigenvectors[0].size();
547
548 // Number of eigenvalues
549 const unsigned int nev_const =
550 (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
551 // nev for arpack, which might change by plus one during dneupd
552 unsigned int nev = nev_const;
553
554 // check input sizes
556 {
557 Assert(nev <= eigenvectors.size(),
559 }
560 else
561 Assert(nev + 1 <= eigenvectors.size(),
563 eigenvectors.size()));
564
565 Assert(nev <= eigenvalues.size(),
567
568 // check large enough problem size
570
574
575 // check whether we have enough Arnoldi vectors
579
580 // ARPACK mode for dsaupd/dnaupd, here only mode 3, i.e. shift-invert mode
581 int mode = 3;
582
583 // reverse communication parameter
584 int ido = 0;
585
586 // 'G' generalized eigenvalue problem 'I' standard eigenvalue problem
587 char bmat[2] = "G";
588
589 // Specify the eigenvalues of interest, possible parameters "LA" algebraically
590 // largest "SA" algebraically smallest "LM" largest magnitude "SM" smallest
591 // magnitude "LR" largest real part "SR" smallest real part "LI" largest
592 // imaginary part "SI" smallest imaginary part "BE" both ends of spectrum
593 // simultaneous.
594 char which[3];
596 {
598 std::strcpy(which, "LA");
599 break;
601 std::strcpy(which, "SA");
602 break;
604 std::strcpy(which, "LM");
605 break;
607 std::strcpy(which, "SM");
608 break;
610 std::strcpy(which, "LR");
611 break;
613 std::strcpy(which, "SR");
614 break;
616 std::strcpy(which, "LI");
617 break;
619 std::strcpy(which, "SI");
620 break;
621 case both_ends:
622 std::strcpy(which, "BE");
623 break;
624 }
625
626 // tolerance for ARPACK
627 double tol = control().tolerance();
628
629 // if the starting vector is used it has to be in resid
630 if (!initial_vector_provided || resid.size() != n)
631 resid.resize(n, 1.);
632
633 // number of Arnoldi basis vectors specified
634 // in additional_data
636
637 int ldv = n;
638 std::vector<double> v(ldv * ncv, 0.0);
639
640 // information to the routines
641 std::vector<int> iparam(11, 0);
642
643 iparam[0] = 1; // shift strategy
644
645 // maximum number of iterations
646 iparam[2] = control().max_steps();
647
648 // Set the mode of dsaupd. 1 is exact shifting, 2 is user-supplied shifts,
649 // 3 is shift-invert mode, 4 is buckling mode, 5 is Cayley mode.
650
651 iparam[6] = mode;
652 std::vector<int> ipntr(14, 0);
653
654 // work arrays for ARPACK
655 std::vector<double> workd(3 * n, 0.);
656 int lworkl =
657 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
658 std::vector<double> workl(lworkl, 0.);
659
660 // information out of the iteration
661 int info = 1;
662
663 while (ido != 99)
664 {
665 // call of ARPACK dsaupd/dnaupd routine
667 dsaupd_(&ido,
668 bmat,
669 &n,
670 which,
671 &nev,
672 &tol,
673 resid.data(),
674 &ncv,
675 v.data(),
676 &ldv,
677 iparam.data(),
678 ipntr.data(),
679 workd.data(),
680 workl.data(),
681 &lworkl,
682 &info);
683 else
684 dnaupd_(&ido,
685 bmat,
686 &n,
687 which,
688 &nev,
689 &tol,
690 resid.data(),
691 &ncv,
692 v.data(),
693 &ldv,
694 iparam.data(),
695 ipntr.data(),
696 workd.data(),
697 workl.data(),
698 &lworkl,
699 &info);
700
701 if (ido == 99)
702 break;
703
704 switch (mode)
705 {
706 case 3:
707 {
708 switch (ido)
709 {
710 case -1:
711 {
712 VectorType src, dst, tmp;
713 src.reinit(eigenvectors[0]);
714 dst.reinit(src);
715 tmp.reinit(src);
716
717
718 for (size_type i = 0; i < src.size(); ++i)
719 src(i) = workd[ipntr[0] - 1 + i];
720
721 // multiplication with mass matrix M
722 mass_matrix.vmult(tmp, src);
723 // solving linear system
724 inverse.vmult(dst, tmp);
725
726 for (size_type i = 0; i < dst.size(); ++i)
727 workd[ipntr[1] - 1 + i] = dst(i);
728 }
729 break;
730
731 case 1:
732 {
733 VectorType src, dst, tmp, tmp2;
734 src.reinit(eigenvectors[0]);
735 dst.reinit(src);
736 tmp.reinit(src);
737 tmp2.reinit(src);
738
739 for (size_type i = 0; i < src.size(); ++i)
740 {
741 src(i) = workd[ipntr[2] - 1 + i];
742 tmp(i) = workd[ipntr[0] - 1 + i];
743 }
744 // solving linear system
745 inverse.vmult(dst, src);
746
747 for (size_type i = 0; i < dst.size(); ++i)
748 workd[ipntr[1] - 1 + i] = dst(i);
749 }
750 break;
751
752 case 2:
753 {
754 VectorType src, dst;
755 src.reinit(eigenvectors[0]);
756 dst.reinit(src);
757
758 for (size_type i = 0; i < src.size(); ++i)
759 src(i) = workd[ipntr[0] - 1 + i];
760
761 // Multiplication with mass matrix M
762 mass_matrix.vmult(dst, src);
763
764 for (size_type i = 0; i < dst.size(); ++i)
765 workd[ipntr[1] - 1 + i] = dst(i);
766 }
767 break;
768
769 default:
770 Assert(false, ArpackExcArpackIdo(ido));
771 break;
772 }
773 }
774 break;
775 default:
776 Assert(false, ArpackExcArpackMode(mode));
777 break;
778 }
779 }
780
781 // Set number of used iterations in SolverControl
782 control().check(iparam[2], 0.);
783
784 if (info < 0)
785 {
787 {
788 Assert(false, ArpackExcArpackInfodsaupd(info));
789 }
790 else
791 Assert(false, ArpackExcArpackInfodnaupd(info));
792 }
793 else
794 {
795 // 1 - compute eigenvectors, 0 - only eigenvalues
796 int rvec = 1;
797
798 // which eigenvectors
799 char howmany = 'A';
800
801 std::vector<int> select(ncv, 1);
802
803 int ldz = n;
804
805 std::vector<double> eigenvalues_real(nev + 1, 0.);
806 std::vector<double> eigenvalues_im(nev + 1, 0.);
807
808 // call of ARPACK dseupd/dneupd routine
810 {
811 std::vector<double> z(ldz * nev, 0.);
812 dseupd_(&rvec,
813 &howmany,
814 select.data(),
815 eigenvalues_real.data(),
816 z.data(),
817 &ldz,
818 &sigmar,
819 bmat,
820 &n,
821 which,
822 &nev,
823 &tol,
824 resid.data(),
825 &ncv,
826 v.data(),
827 &ldv,
828 iparam.data(),
829 ipntr.data(),
830 workd.data(),
831 workl.data(),
832 &lworkl,
833 &info);
834 }
835 else
836 {
837 std::vector<double> workev(3 * ncv, 0.);
838 dneupd_(&rvec,
839 &howmany,
840 select.data(),
841 eigenvalues_real.data(),
842 eigenvalues_im.data(),
843 v.data(),
844 &ldz,
845 &sigmar,
846 &sigmai,
847 workev.data(),
848 bmat,
849 &n,
850 which,
851 &nev,
852 &tol,
853 resid.data(),
854 &ncv,
855 v.data(),
856 &ldv,
857 iparam.data(),
858 ipntr.data(),
859 workd.data(),
860 workl.data(),
861 &lworkl,
862 &info);
863 }
864
865 if (info == 1)
866 {
867 Assert(false, ArpackExcArpackInfoMaxIt(control().max_steps()));
868 }
869 else if (info == 3)
870 {
872 }
873 else if (info != 0)
874 {
876 {
877 Assert(false, ArpackExcArpackInfodseupd(info));
878 }
879 else
880 Assert(false, ArpackExcArpackInfodneupd(info));
881 }
882
883 for (unsigned int i = 0; i < nev; ++i)
884 for (unsigned int j = 0; j < n; ++j)
885 eigenvectors[i](j) = v[i * n + j];
886
887 for (unsigned int i = 0; i < nev_const; ++i)
888 eigenvalues[i] =
889 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
890 }
891}
892
893
894inline SolverControl &
896{
897 return solver_control;
898}
899
901
902
903#endif
904#endif
void dseupd_(int *rvec, char *howmany, int *select, double *d, double *z, int *ldz, double *sigmar, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dsaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dnaupd_(int *ido, char *bmat, unsigned int *n, char *which, unsigned int *nev, const double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
void dneupd_(int *rvec, char *howmany, int *select, double *d, double *di, double *z, int *ldz, double *sigmar, double *sigmai, double *workev, char *bmat, unsigned int *n, char *which, unsigned int *nev, double *tol, double *resid, int *ncv, double *v, int *ldv, int *iparam, int *ipntr, double *workd, double *workl, int *lworkl, int *info)
SolverControl & control() const
void set_initial_vector(const VectorType &vec)
bool initial_vector_provided
ArpackSolver(SolverControl &control, const AdditionalData &data=AdditionalData())
std::vector< double > resid
SolverControl & solver_control
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=0)
void set_shift(const std::complex< double > sigma)
const AdditionalData additional_data
virtual State check(const unsigned int step, const double check_value)
unsigned int max_steps() const
double tolerance() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ArpackExcInvalidEigenvalueSize(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackInfoMaxIt(int arg1)
static ::ExceptionBase & ArpackExcInvalidNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ArpackExcInvalidEigenvectorSizeNonsymmetric(int arg1, int arg2)
static ::ExceptionBase & ArpackExcSmallNumberofArnoldiVectors(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackNoShifts()
static ::ExceptionBase & ArpackExcInvalidEigenvectorSize(int arg1, int arg2)
static ::ExceptionBase & ArpackExcArpackMode(int arg1)
static ::ExceptionBase & ArpackExcInvalidNumberofEigenvalues(int arg1, int arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ArpackExcArpackInfodneupd(int arg1)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
static ::ExceptionBase & ArpackExcArpackInfodsaupd(int arg1)
static ::ExceptionBase & ArpackExcArpackInfodnaupd(int arg1)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:488
static ::ExceptionBase & ArpackExcArpackInfodseupd(int arg1)
static ::ExceptionBase & ArpackExcArpackIdo(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int global_dof_index
Definition types.h:82
const WhichEigenvalues eigenvalue_of_interest
const unsigned int number_of_arnoldi_vectors
AdditionalData(const unsigned int number_of_arnoldi_vectors=15, const WhichEigenvalues eigenvalue_of_interest=largest_magnitude, const bool symmetric=false)
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)