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
slepc_solver.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2009 - 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
17#ifndef dealii_slepc_solver_h
18# define dealii_slepc_solver_h
19
20# include <deal.II/base/config.h>
21
22# ifdef DEAL_II_WITH_SLEPC
23
27
28# include <petscconf.h>
29# include <petscksp.h>
30
31# include <slepceps.h>
32
33# include <memory>
34
36
132{
147 {
148 public:
154
158 virtual ~SolverBase();
159
178 template <typename OutputVector>
179 void
181 std::vector<PetscScalar> & eigenvalues,
182 std::vector<OutputVector> & eigenvectors,
183 const unsigned int n_eigenpairs = 1);
184
190 template <typename OutputVector>
191 void
194 std::vector<PetscScalar> & eigenvalues,
195 std::vector<OutputVector> & eigenvectors,
196 const unsigned int n_eigenpairs = 1);
197
203 template <typename OutputVector>
204 void
207 std::vector<double> & real_eigenvalues,
208 std::vector<double> & imag_eigenvalues,
209 std::vector<OutputVector> & real_eigenvectors,
210 std::vector<OutputVector> & imag_eigenvectors,
211 const unsigned int n_eigenpairs = 1);
212
219 template <typename Vector>
220 void
221 set_initial_space(const std::vector<Vector> &initial_space);
222
226 void
228
233 void
234 set_target_eigenvalue(const PetscScalar &this_target);
235
242 void
243 set_which_eigenpairs(EPSWhich set_which);
244
253 void
254 set_problem_type(EPSProblemType set_problem);
255
261 void
263
268
273 int,
274 << " An error with error number " << arg1
275 << " occurred while calling a SLEPc function");
276
281 int,
282 int,
283 << " The number of converged eigenvectors is " << arg1
284 << " but " << arg2 << " were requested. ");
285
290 control() const;
291
292 protected:
298
303
310 void
311 solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
312
318 void
319 get_eigenpair(const unsigned int index,
320 PetscScalar & eigenvalues,
322
328 void
329 get_eigenpair(const unsigned int index,
330 double & real_eigenvalues,
331 double & imag_eigenvalues,
332 PETScWrappers::VectorBase &real_eigenvectors,
333 PETScWrappers::VectorBase &imag_eigenvectors);
334
339 void
341
346 void
349
350 protected:
354 EPS eps;
355
356 private:
360 EPSConvergedReason reason;
361
362
369 static int
371 PetscScalar real_eigenvalue,
372 PetscScalar imag_eigenvalue,
373 PetscReal residual_norm,
374 PetscReal * estimated_error,
375 void * solver_control);
376 };
377
378
379
393 {
394 public:
400 {};
401
408 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
409 const AdditionalData &data = AdditionalData());
410
411 protected:
416 };
417
418
419
433 {
434 public:
440 {
446
451 };
452
459 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
460 const AdditionalData &data = AdditionalData());
461
462 protected:
467 };
468
469
470
484 {
485 public:
491 {
495 EPSLanczosReorthogType reorthog;
496
502 const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
503 };
504
511 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
512 const AdditionalData &data = AdditionalData());
513
514 protected:
519 };
520
521
522
535 class SolverPower : public SolverBase
536 {
537 public:
543 {};
544
551 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
552 const AdditionalData &data = AdditionalData());
553
554 protected:
559 };
560
561
562
576 {
577 public:
583 {
588
592 AdditionalData(bool double_expansion = false);
593 };
594
601 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
602 const AdditionalData &data = AdditionalData());
603
604 protected:
609 };
610
611
612
626 {
627 public:
633 {};
634
641 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
642 const AdditionalData &data = AdditionalData());
643
644 protected:
649 };
650
651
652
666 {
667 public:
673 {};
674
681 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
682 const AdditionalData &data = AdditionalData());
683
684 protected:
689 };
690
691
692
693 // --------------------------- inline and template functions -----------
698 // todo: The logic of these functions can be simplified without breaking
699 // backward compatibility...
700
701 template <typename OutputVector>
702 void
704 std::vector<PetscScalar> & eigenvalues,
705 std::vector<OutputVector> & eigenvectors,
706 const unsigned int n_eigenpairs)
707 {
708 // Panic if the number of eigenpairs wanted is out of bounds.
709 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
711
712 // Set the matrices of the problem
713 set_matrices(A);
714
715 // and solve
716 unsigned int n_converged = 0;
717 solve(n_eigenpairs, &n_converged);
718
719 if (n_converged > n_eigenpairs)
720 n_converged = n_eigenpairs;
721 AssertThrow(n_converged == n_eigenpairs,
723 n_eigenpairs));
724
726 eigenvectors.resize(n_converged, eigenvectors.front());
727 eigenvalues.resize(n_converged);
728
729 for (unsigned int index = 0; index < n_converged; ++index)
730 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
731 }
732
733 template <typename OutputVector>
734 void
737 std::vector<PetscScalar> & eigenvalues,
738 std::vector<OutputVector> & eigenvectors,
739 const unsigned int n_eigenpairs)
740 {
741 // Guard against incompatible matrix sizes:
742 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
743 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
744
745 // Panic if the number of eigenpairs wanted is out of bounds.
746 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
748
749 // Set the matrices of the problem
750 set_matrices(A, B);
751
752 // and solve
753 unsigned int n_converged = 0;
754 solve(n_eigenpairs, &n_converged);
755
756 if (n_converged >= n_eigenpairs)
757 n_converged = n_eigenpairs;
758
759 AssertThrow(n_converged == n_eigenpairs,
761 n_eigenpairs));
763
764 eigenvectors.resize(n_converged, eigenvectors.front());
765 eigenvalues.resize(n_converged);
766
767 for (unsigned int index = 0; index < n_converged; ++index)
768 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
769 }
770
771 template <typename OutputVector>
772 void
775 std::vector<double> & real_eigenvalues,
776 std::vector<double> & imag_eigenvalues,
777 std::vector<OutputVector> & real_eigenvectors,
778 std::vector<OutputVector> & imag_eigenvectors,
779 const unsigned int n_eigenpairs)
780 {
781 // Guard against incompatible matrix sizes:
782 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
783 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
784
785 // and incompatible eigenvalue/eigenvector sizes
786 AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
787 ExcDimensionMismatch(real_eigenvalues.size(),
788 imag_eigenvalues.size()));
789 AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
790 ExcDimensionMismatch(real_eigenvectors.size(),
791 imag_eigenvectors.size()));
792
793 // Panic if the number of eigenpairs wanted is out of bounds.
794 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
796
797 // Set the matrices of the problem
798 set_matrices(A, B);
799
800 // and solve
801 unsigned int n_converged = 0;
802 solve(n_eigenpairs, &n_converged);
803
804 if (n_converged >= n_eigenpairs)
805 n_converged = n_eigenpairs;
806
807 AssertThrow(n_converged == n_eigenpairs,
809 n_eigenpairs));
810 AssertThrow((real_eigenvectors.size() != 0) &&
811 (imag_eigenvectors.size() != 0),
813
814 real_eigenvectors.resize(n_converged, real_eigenvectors.front());
815 imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
816 real_eigenvalues.resize(n_converged);
817 imag_eigenvalues.resize(n_converged);
818
819 for (unsigned int index = 0; index < n_converged; ++index)
820 get_eigenpair(index,
821 real_eigenvalues[index],
822 imag_eigenvalues[index],
823 real_eigenvectors[index],
824 imag_eigenvectors[index]);
825 }
826
827 template <typename Vector>
828 void
829 SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
830 {
831 std::vector<Vec> vecs(this_initial_space.size());
832
833 for (unsigned int i = 0; i < this_initial_space.size(); ++i)
834 {
835 Assert(this_initial_space[i].l2_norm() > 0.0,
836 ExcMessage("Initial vectors should be nonzero."));
837 vecs[i] = this_initial_space[i];
838 }
839
840 // if the eigensolver supports only a single initial vector, but several
841 // guesses are provided, then all except the first one will be discarded.
842 // One could still build a vector that is rich in the directions of all
843 // guesses, by taking a linear combination of them. (TODO: make function
844 // virtual?)
845
846 const PetscErrorCode ierr =
847 EPSSetInitialSpace(eps, vecs.size(), vecs.data());
848 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
849 }
850
851} // namespace SLEPcWrappers
852
854
855# endif // DEAL_II_WITH_SLEPC
856
857/*---------------------------- slepc_solver.h ---------------------------*/
858
859#endif
860
861/*---------------------------- slepc_solver.h ---------------------------*/
const AdditionalData additional_data
SolverControl & control() const
EPSConvergedReason reason
void set_target_eigenvalue(const PetscScalar &this_target)
const MPI_Comm mpi_communicator
void set_matrices(const PETScWrappers::MatrixBase &A)
void get_solver_state(const SolverControl::State state)
void set_initial_space(const std::vector< Vector > &initial_space)
void set_problem_type(EPSProblemType set_problem)
SolverControl & solver_control
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
static int convergence_test(EPS eps, PetscScalar real_eigenvalue, PetscScalar imag_eigenvalue, PetscReal residual_norm, PetscReal *estimated_error, void *solver_control)
void set_which_eigenpairs(EPSWhich set_which)
const AdditionalData additional_data
const AdditionalData additional_data
const AdditionalData additional_data
const AdditionalData additional_data
const AdditionalData additional_data
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
static ::ExceptionBase & ExcSLEPcError(int arg1)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:510
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
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)