deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+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
arpack_solver.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2010 - 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_arpack_solver_h
16#define dealii_arpack_solver_h
17
18#include <deal.II/base/config.h>
19
23#include <deal.II/base/types.h>
24
26
27#include <cstring>
28
29
30#ifdef DEAL_II_WITH_ARPACK
31
33
34
35extern "C" void
37 char *bmat,
38 unsigned int *n,
39 char *which,
40 unsigned int *nev,
41 const double *tol,
42 double *resid,
43 int *ncv,
44 double *v,
45 int *ldv,
46 int *iparam,
47 int *ipntr,
48 double *workd,
49 double *workl,
50 int *lworkl,
51 int *info);
52
53extern "C" void
55 char *bmat,
56 unsigned int *n,
57 char *which,
58 unsigned int *nev,
59 double *tol,
60 double *resid,
61 int *ncv,
62 double *v,
63 int *ldv,
64 int *iparam,
65 int *ipntr,
66 double *workd,
67 double *workl,
68 int *lworkl,
69 int *info);
70
71extern "C" void
73 char *howmany,
74 int *select,
75 double *d,
76 double *di,
77 double *z,
78 int *ldz,
79 double *sigmar,
80 double *sigmai,
81 double *workev,
82 char *bmat,
83 unsigned int *n,
84 char *which,
85 unsigned int *nev,
86 double *tol,
87 double *resid,
88 int *ncv,
89 double *v,
90 int *ldv,
91 int *iparam,
92 int *ipntr,
93 double *workd,
94 double *workl,
95 int *lworkl,
96 int *info);
97
98extern "C" void
100 char *howmany,
101 int *select,
102 double *d,
103 double *z,
104 int *ldz,
105 double *sigmar,
106 char *bmat,
107 unsigned int *n,
108 char *which,
109 unsigned int *nev,
110 double *tol,
111 double *resid,
112 int *ncv,
113 double *v,
114 int *ldv,
115 int *iparam,
116 int *ipntr,
117 double *workd,
118 double *workl,
119 int *lworkl,
120 int *info);
121
170{
171public:
176
177
224
229 {
235 explicit AdditionalData(
236 const unsigned int number_of_arnoldi_vectors = 15,
238 const bool symmetric = false);
239
245 const unsigned int number_of_arnoldi_vectors;
246
251
255 const bool symmetric;
256 };
257
262 control() const;
263
269
273 template <typename VectorType>
274 void
275 set_initial_vector(const VectorType &vec);
276
282 void
283 set_shift(const std::complex<double> sigma);
284
334 template <typename VectorType,
335 typename MatrixType1,
336 typename MatrixType2,
337 typename INVERSE>
338 void
339 solve(const MatrixType1 &A,
340 const MatrixType2 &B,
341 const INVERSE &inverse,
342 std::vector<std::complex<double>> &eigenvalues,
343 std::vector<VectorType> &eigenvectors,
344 const unsigned int n_eigenvalues = 0);
345
346protected:
352
357
362 std::vector<double> resid;
363
367 double sigmar;
368
372 double sigmai;
373
374
375private:
380 int,
381 int,
382 << "Number of wanted eigenvalues " << arg1
383 << " is larger that the size of the matrix " << arg2);
384
386 int,
387 int,
388 << "Number of wanted eigenvalues " << arg1
389 << " is larger that the size of eigenvectors " << arg2);
390
393 int,
394 int,
395 << "To store the real and complex parts of " << arg1
396 << " eigenvectors in real-valued vectors, their size (currently set to "
397 << arg2 << ") should be greater than or equal to " << arg1 + 1);
398
400 int,
401 int,
402 << "Number of wanted eigenvalues " << arg1
403 << " is larger that the size of eigenvalues " << arg2);
404
406 int,
407 int,
408 << "Number of Arnoldi vectors " << arg1
409 << " is larger that the size of the matrix " << arg2);
410
412 int,
413 int,
414 << "Number of Arnoldi vectors " << arg1
415 << " is too small to obtain " << arg2 << " eigenvalues");
416
418 int,
419 << "This ido " << arg1
420 << " is not supported. Check documentation of ARPACK");
421
423 int,
424 << "This mode " << arg1
425 << " is not supported. Check documentation of ARPACK");
426
428 int,
429 << "Error with dsaupd, info " << arg1
430 << ". Check documentation of ARPACK");
431
433 int,
434 << "Error with dnaupd, info " << arg1
435 << ". Check documentation of ARPACK");
436
438 int,
439 << "Error with dseupd, info " << arg1
440 << ". Check documentation of ARPACK");
441
443 int,
444 << "Error with dneupd, info " << arg1
445 << ". Check documentation of ARPACK");
446
448 int,
449 << "Maximum number " << arg1 << " of iterations reached.");
450
452 "No shifts could be applied during implicit"
453 " Arnoldi update, try increasing the number of"
454 " Arnoldi vectors.");
455};
456
457
459 const unsigned int number_of_arnoldi_vectors,
460 const WhichEigenvalues eigenvalue_of_interest,
461 const bool symmetric)
462 : number_of_arnoldi_vectors(number_of_arnoldi_vectors)
463 , eigenvalue_of_interest(eigenvalue_of_interest)
464 , symmetric(symmetric)
465{
466 // Check for possible options for symmetric problems
467 if (symmetric)
468 {
469 Assert(
472 "'largest real part' can only be used for non-symmetric problems!"));
473 Assert(
476 "'smallest real part' can only be used for non-symmetric problems!"));
477 Assert(
480 "'largest imaginary part' can only be used for non-symmetric problems!"));
481 Assert(
484 "'smallest imaginary part' can only be used for non-symmetric problems!"));
485 }
486 // Check for possible options for asymmetric problems
487 else
488 {
489 Assert(
492 "'largest algebraic part' can only be used for symmetric problems!"));
493 Assert(
496 "'smallest algebraic part' can only be used for symmetric problems!"));
499 "'both ends' can only be used for symmetric problems!"));
500 }
501}
502
503
512
513
514
515inline void
516ArpackSolver::set_shift(const std::complex<double> sigma)
517{
518 sigmar = sigma.real();
519 sigmai = sigma.imag();
520}
521
522
523
524template <typename VectorType>
525inline void
527{
529 resid.resize(vec.size());
530 for (size_type i = 0; i < vec.size(); ++i)
531 resid[i] = vec[i];
532}
533
534
535template <typename VectorType,
536 typename MatrixType1,
537 typename MatrixType2,
538 typename INVERSE>
539inline void
540ArpackSolver::solve(const MatrixType1 & /*system_matrix*/,
541 const MatrixType2 &mass_matrix,
542 const INVERSE &inverse,
543 std::vector<std::complex<double>> &eigenvalues,
544 std::vector<VectorType> &eigenvectors,
545 const unsigned int n_eigenvalues)
546{
547 // Problem size
548 unsigned int n = eigenvectors[0].size();
549
550 // Number of eigenvalues
551 const unsigned int nev_const =
552 (n_eigenvalues == 0) ? eigenvalues.size() : n_eigenvalues;
553 // nev for arpack, which might change by plus one during dneupd
554 unsigned int nev = nev_const;
555
556 // check input sizes
557 if (additional_data.symmetric)
558 {
559 Assert(nev <= eigenvectors.size(),
561 }
562 else
563 Assert(nev + 1 <= eigenvectors.size(),
565 eigenvectors.size()));
566
567 Assert(nev <= eigenvalues.size(),
569
570 // check large enough problem size
572
573 Assert(additional_data.number_of_arnoldi_vectors < n,
575 additional_data.number_of_arnoldi_vectors, n));
576
577 // check whether we have enough Arnoldi vectors
578 Assert(additional_data.number_of_arnoldi_vectors > 2 * nev + 1,
580 additional_data.number_of_arnoldi_vectors, nev));
581
582 // ARPACK mode for dsaupd/dnaupd, here only mode 3, i.e. shift-invert mode
583 int mode = 3;
584
585 // reverse communication parameter
586 int ido = 0;
587
588 // 'G' generalized eigenvalue problem 'I' standard eigenvalue problem
589 char bmat[2] = "G";
590
591 // Specify the eigenvalues of interest, possible parameters "LA" algebraically
592 // largest "SA" algebraically smallest "LM" largest magnitude "SM" smallest
593 // magnitude "LR" largest real part "SR" smallest real part "LI" largest
594 // imaginary part "SI" smallest imaginary part "BE" both ends of spectrum
595 // simultaneous.
596 char which[3];
597 switch (additional_data.eigenvalue_of_interest)
598 {
600 std::strcpy(which, "LA");
601 break;
603 std::strcpy(which, "SA");
604 break;
606 std::strcpy(which, "LM");
607 break;
609 std::strcpy(which, "SM");
610 break;
612 std::strcpy(which, "LR");
613 break;
615 std::strcpy(which, "SR");
616 break;
618 std::strcpy(which, "LI");
619 break;
621 std::strcpy(which, "SI");
622 break;
623 case both_ends:
624 std::strcpy(which, "BE");
625 break;
626 }
627
628 // tolerance for ARPACK
629 double tol = control().tolerance();
630
631 // if the starting vector is used it has to be in resid
632 if (!initial_vector_provided || resid.size() != n)
633 resid.resize(n, 1.);
634
635 // number of Arnoldi basis vectors specified
636 // in additional_data
637 int ncv = additional_data.number_of_arnoldi_vectors;
638
639 int ldv = n;
640 std::vector<double> v(ldv * ncv, 0.0);
641
642 // information to the routines
643 std::vector<int> iparam(11, 0);
644
645 iparam[0] = 1; // shift strategy
646
647 // maximum number of iterations
648 iparam[2] = control().max_steps();
649
650 // Set the mode of dsaupd. 1 is exact shifting, 2 is user-supplied shifts,
651 // 3 is shift-invert mode, 4 is buckling mode, 5 is Cayley mode.
652
653 iparam[6] = mode;
654 std::vector<int> ipntr(14, 0);
655
656 // work arrays for ARPACK
657 std::vector<double> workd(3 * n, 0.);
658 int lworkl =
659 additional_data.symmetric ? ncv * ncv + 8 * ncv : 3 * ncv * ncv + 6 * ncv;
660 std::vector<double> workl(lworkl, 0.);
661
662 // information out of the iteration
663 int info = 1;
664
665 while (ido != 99)
666 {
667 // call of ARPACK dsaupd/dnaupd routine
668 if (additional_data.symmetric)
669 dsaupd_(&ido,
670 bmat,
671 &n,
672 which,
673 &nev,
674 &tol,
675 resid.data(),
676 &ncv,
677 v.data(),
678 &ldv,
679 iparam.data(),
680 ipntr.data(),
681 workd.data(),
682 workl.data(),
683 &lworkl,
684 &info);
685 else
686 dnaupd_(&ido,
687 bmat,
688 &n,
689 which,
690 &nev,
691 &tol,
692 resid.data(),
693 &ncv,
694 v.data(),
695 &ldv,
696 iparam.data(),
697 ipntr.data(),
698 workd.data(),
699 workl.data(),
700 &lworkl,
701 &info);
702
703 if (ido == 99)
704 break;
705
706 switch (mode)
707 {
708 case 3:
709 {
710 switch (ido)
711 {
712 case -1:
713 {
714 VectorType src, dst, tmp;
715 src.reinit(eigenvectors[0]);
716 dst.reinit(src);
717 tmp.reinit(src);
718
719
720 for (size_type i = 0; i < src.size(); ++i)
721 src(i) = workd[ipntr[0] - 1 + i];
722
723 // multiplication with mass matrix M
724 mass_matrix.vmult(tmp, src);
725 // solving linear system
726 inverse.vmult(dst, tmp);
727
728 for (size_type i = 0; i < dst.size(); ++i)
729 workd[ipntr[1] - 1 + i] = dst(i);
730 }
731 break;
732
733 case 1:
734 {
735 VectorType src, dst, tmp, tmp2;
736 src.reinit(eigenvectors[0]);
737 dst.reinit(src);
738 tmp.reinit(src);
739 tmp2.reinit(src);
740
741 for (size_type i = 0; i < src.size(); ++i)
742 {
743 src(i) = workd[ipntr[2] - 1 + i];
744 tmp(i) = workd[ipntr[0] - 1 + i];
745 }
746 // solving linear system
747 inverse.vmult(dst, src);
748
749 for (size_type i = 0; i < dst.size(); ++i)
750 workd[ipntr[1] - 1 + i] = dst(i);
751 }
752 break;
753
754 case 2:
755 {
756 VectorType src, dst;
757 src.reinit(eigenvectors[0]);
758 dst.reinit(src);
759
760 for (size_type i = 0; i < src.size(); ++i)
761 src(i) = workd[ipntr[0] - 1 + i];
762
763 // Multiplication with mass matrix M
764 mass_matrix.vmult(dst, src);
765
766 for (size_type i = 0; i < dst.size(); ++i)
767 workd[ipntr[1] - 1 + i] = dst(i);
768 }
769 break;
770
771 default:
773 break;
774 }
775 }
776 break;
777 default:
778 Assert(false, ArpackExcArpackMode(mode));
779 break;
780 }
781 }
782
783 // Set number of used iterations in SolverControl
784 control().check(iparam[2], 0.);
785
786 if (info < 0)
787 {
788 if (additional_data.symmetric)
789 {
791 }
792 else
794 }
795 else
796 {
797 // 1 - compute eigenvectors, 0 - only eigenvalues
798 int rvec = 1;
799
800 // which eigenvectors
801 char howmany = 'A';
802
803 std::vector<int> select(ncv, 1);
804
805 int ldz = n;
806
807 std::vector<double> eigenvalues_real(nev + 1, 0.);
808 std::vector<double> eigenvalues_im(nev + 1, 0.);
809
810 // call of ARPACK dseupd/dneupd routine
811 if (additional_data.symmetric)
812 {
813 std::vector<double> z(ldz * nev, 0.);
814 dseupd_(&rvec,
815 &howmany,
816 select.data(),
817 eigenvalues_real.data(),
818 z.data(),
819 &ldz,
820 &sigmar,
821 bmat,
822 &n,
823 which,
824 &nev,
825 &tol,
826 resid.data(),
827 &ncv,
828 v.data(),
829 &ldv,
830 iparam.data(),
831 ipntr.data(),
832 workd.data(),
833 workl.data(),
834 &lworkl,
835 &info);
836 }
837 else
838 {
839 std::vector<double> workev(3 * ncv, 0.);
840 dneupd_(&rvec,
841 &howmany,
842 select.data(),
843 eigenvalues_real.data(),
844 eigenvalues_im.data(),
845 v.data(),
846 &ldz,
847 &sigmar,
848 &sigmai,
849 workev.data(),
850 bmat,
851 &n,
852 which,
853 &nev,
854 &tol,
855 resid.data(),
856 &ncv,
857 v.data(),
858 &ldv,
859 iparam.data(),
860 ipntr.data(),
861 workd.data(),
862 workl.data(),
863 &lworkl,
864 &info);
865 }
866
867 if (info == 1)
868 {
869 Assert(false, ArpackExcArpackInfoMaxIt(control().max_steps()));
870 }
871 else if (info == 3)
872 {
874 }
875 else if (info != 0)
876 {
877 if (additional_data.symmetric)
878 {
880 }
881 else
883 }
884
885 for (unsigned int i = 0; i < nev; ++i)
886 for (unsigned int j = 0; j < n; ++j)
887 eigenvectors[i](j) = v[i * n + j];
888
889 for (unsigned int i = 0; i < nev_const; ++i)
890 eigenvalues[i] =
891 std::complex<double>(eigenvalues_real[i], eigenvalues_im[i]);
892 }
893}
894
895
896inline SolverControl &
898{
899 return solver_control;
900}
901
903
904
905#endif
906#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:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
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:536
static ::ExceptionBase & ArpackExcArpackInfodsaupd(int arg1)
static ::ExceptionBase & ArpackExcArpackInfodnaupd(int arg1)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:491
static ::ExceptionBase & ArpackExcArpackInfodseupd(int arg1)
static ::ExceptionBase & ArpackExcArpackIdo(int arg1)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:513
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
Definition mpi.cc:740
unsigned int global_dof_index
Definition types.h:90
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 > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)