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