Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08: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
slepc_solver.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 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
16#ifndef dealii_slepc_solver_h
17# define dealii_slepc_solver_h
18
19# include <deal.II/base/config.h>
20
21# ifdef DEAL_II_WITH_SLEPC
22
26
27# include <petscconf.h>
28# include <petscksp.h>
29
30# include <slepceps.h>
31
32# include <memory>
33
35
131{
146 {
147 public:
153
157 virtual ~SolverBase();
158
177 template <typename OutputVector>
178 void
180 std::vector<PetscScalar> &eigenvalues,
181 std::vector<OutputVector> &eigenvectors,
182 const unsigned int n_eigenpairs = 1);
183
189 template <typename OutputVector>
190 void
193 std::vector<PetscScalar> &eigenvalues,
194 std::vector<OutputVector> &eigenvectors,
195 const unsigned int n_eigenpairs = 1);
196
202 template <typename OutputVector>
203 void
206 std::vector<double> &real_eigenvalues,
207 std::vector<double> &imag_eigenvalues,
208 std::vector<OutputVector> &real_eigenvectors,
209 std::vector<OutputVector> &imag_eigenvectors,
210 const unsigned int n_eigenpairs = 1);
211
218 template <typename Vector>
219 void
220 set_initial_space(const std::vector<Vector> &initial_space);
221
225 void
227
232 void
233 set_target_eigenvalue(const PetscScalar &this_target);
234
241 void
242 set_which_eigenpairs(EPSWhich set_which);
243
252 void
253 set_problem_type(EPSProblemType set_problem);
254
260 void
262
267
272 int,
273 << " An error with error number " << arg1
274 << " occurred while calling a SLEPc function");
275
280 int,
281 int,
282 << " The number of converged eigenvectors is " << arg1
283 << " but " << arg2 << " were requested. ");
284
289 control() const;
290
291 protected:
297
302
309 void
310 solve(const unsigned int n_eigenpairs, unsigned int *n_converged);
311
317 void
318 get_eigenpair(const unsigned int index,
319 PetscScalar &eigenvalues,
321
327 void
328 get_eigenpair(const unsigned int index,
329 double &real_eigenvalues,
330 double &imag_eigenvalues,
331 PETScWrappers::VectorBase &real_eigenvectors,
332 PETScWrappers::VectorBase &imag_eigenvectors);
333
338 void
340
345 void
348
349 protected:
353 EPS eps;
354
355 private:
359 EPSConvergedReason reason;
360
361
368 static int
370 PetscScalar real_eigenvalue,
371 PetscScalar imag_eigenvalue,
372 PetscReal residual_norm,
373 PetscReal *estimated_error,
374 void *solver_control);
375 };
376
377
378
392 {
393 public:
399 {};
400
406 explicit SolverKrylovSchur(
407 SolverControl &cn,
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 {
445 explicit AdditionalData(const bool delayed_reorthogonalization = false);
446
451 };
452
458 explicit SolverArnoldi(SolverControl &cn,
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
501 explicit AdditionalData(
502 const EPSLanczosReorthogType r = EPS_LANCZOS_REORTHOG_FULL);
503 };
504
510 explicit SolverLanczos(SolverControl &cn,
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
550 explicit SolverPower(SolverControl &cn,
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 explicit AdditionalData(bool double_expansion = false);
593 };
594
601 SolverControl &cn,
602 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
603 const AdditionalData &data = AdditionalData());
604
605 protected:
610 };
611
612
613
627 {
628 public:
634 {};
635
641 explicit SolverJacobiDavidson(
642 SolverControl &cn,
643 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
644 const AdditionalData &data = AdditionalData());
645
646 protected:
651 };
652
653
654
668 {
669 public:
675 {};
676
682 explicit SolverLAPACK(SolverControl &cn,
683 const MPI_Comm mpi_communicator = PETSC_COMM_SELF,
684 const AdditionalData &data = AdditionalData());
685
686 protected:
691 };
692
693
694
695 // --------------------------- inline and template functions -----------
700 // todo: The logic of these functions can be simplified without breaking
701 // backward compatibility...
702
703 template <typename OutputVector>
704 void
706 std::vector<PetscScalar> &eigenvalues,
707 std::vector<OutputVector> &eigenvectors,
708 const unsigned int n_eigenpairs)
709 {
710 // Panic if the number of eigenpairs wanted is out of bounds.
711 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
713
714 // Set the matrices of the problem
715 set_matrices(A);
716
717 // and solve
718 unsigned int n_converged = 0;
719 solve(n_eigenpairs, &n_converged);
720
721 if (n_converged > n_eigenpairs)
722 n_converged = n_eigenpairs;
723 AssertThrow(n_converged == n_eigenpairs,
725 n_eigenpairs));
726
728 eigenvectors.resize(n_converged, eigenvectors.front());
729 eigenvalues.resize(n_converged);
730
731 for (unsigned int index = 0; index < n_converged; ++index)
732 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
733 }
734
735 template <typename OutputVector>
736 void
739 std::vector<PetscScalar> &eigenvalues,
740 std::vector<OutputVector> &eigenvectors,
741 const unsigned int n_eigenpairs)
742 {
743 // Guard against incompatible matrix sizes:
744 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
745 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
746
747 // Panic if the number of eigenpairs wanted is out of bounds.
748 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
750
751 // Set the matrices of the problem
752 set_matrices(A, B);
753
754 // and solve
755 unsigned int n_converged = 0;
756 solve(n_eigenpairs, &n_converged);
757
758 if (n_converged >= n_eigenpairs)
759 n_converged = n_eigenpairs;
760
761 AssertThrow(n_converged == n_eigenpairs,
763 n_eigenpairs));
765
766 eigenvectors.resize(n_converged, eigenvectors.front());
767 eigenvalues.resize(n_converged);
768
769 for (unsigned int index = 0; index < n_converged; ++index)
770 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
771 }
772
773 template <typename OutputVector>
774 void
777 std::vector<double> &real_eigenvalues,
778 std::vector<double> &imag_eigenvalues,
779 std::vector<OutputVector> &real_eigenvectors,
780 std::vector<OutputVector> &imag_eigenvectors,
781 const unsigned int n_eigenpairs)
782 {
783 // Guard against incompatible matrix sizes:
784 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
785 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
786
787 // and incompatible eigenvalue/eigenvector sizes
788 AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
789 ExcDimensionMismatch(real_eigenvalues.size(),
790 imag_eigenvalues.size()));
791 AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
792 ExcDimensionMismatch(real_eigenvectors.size(),
793 imag_eigenvectors.size()));
794
795 // Panic if the number of eigenpairs wanted is out of bounds.
796 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
798
799 // Set the matrices of the problem
800 set_matrices(A, B);
801
802 // and solve
803 unsigned int n_converged = 0;
804 solve(n_eigenpairs, &n_converged);
805
806 if (n_converged >= n_eigenpairs)
807 n_converged = n_eigenpairs;
808
809 AssertThrow(n_converged == n_eigenpairs,
811 n_eigenpairs));
812 AssertThrow((real_eigenvectors.size() != 0) &&
813 (imag_eigenvectors.size() != 0),
815
816 real_eigenvectors.resize(n_converged, real_eigenvectors.front());
817 imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
818 real_eigenvalues.resize(n_converged);
819 imag_eigenvalues.resize(n_converged);
820
821 for (unsigned int index = 0; index < n_converged; ++index)
822 get_eigenpair(index,
823 real_eigenvalues[index],
824 imag_eigenvalues[index],
825 real_eigenvectors[index],
826 imag_eigenvectors[index]);
827 }
828
829 template <typename Vector>
830 void
831 SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
832 {
833 std::vector<Vec> vecs(this_initial_space.size());
834
835 for (unsigned int i = 0; i < this_initial_space.size(); ++i)
836 {
837 Assert(this_initial_space[i].l2_norm() > 0.0,
838 ExcMessage("Initial vectors should be nonzero."));
839 vecs[i] = this_initial_space[i];
840 }
841
842 // if the eigensolver supports only a single initial vector, but several
843 // guesses are provided, then all except the first one will be discarded.
844 // One could still build a vector that is rich in the directions of all
845 // guesses, by taking a linear combination of them. (TODO: make function
846 // virtual?)
847
848 const PetscErrorCode ierr =
849 EPSSetInitialSpace(eps, vecs.size(), vecs.data());
850 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
851 }
852
853} // namespace SLEPcWrappers
854
856
857# endif // DEAL_II_WITH_SLEPC
858
859/*---------------------------- slepc_solver.h ---------------------------*/
860
861#endif
862
863/*---------------------------- 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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DeclException0(Exception0)
Definition exceptions.h:471
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:539
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition exceptions.h:516
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)