Reference documentation for deal.II version 9.4.1
\(\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 - 2021 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 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
642 const MPI_Comm &mpi_communicator = PETSC_COMM_SELF,
643 const AdditionalData &data = AdditionalData());
644
645 protected:
650 };
651
652
653
667 {
668 public:
674 {};
675
682 const MPI_Comm & mpi_communicator = PETSC_COMM_SELF,
683 const AdditionalData &data = AdditionalData());
684
685 protected:
690 };
691
692
693
694 // --------------------------- inline and template functions -----------
699 // todo: The logic of these functions can be simplified without breaking
700 // backward compatibility...
701
702 template <typename OutputVector>
703 void
705 std::vector<PetscScalar> & eigenvalues,
706 std::vector<OutputVector> & eigenvectors,
707 const unsigned int n_eigenpairs)
708 {
709 // Panic if the number of eigenpairs wanted is out of bounds.
710 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
712
713 // Set the matrices of the problem
714 set_matrices(A);
715
716 // and solve
717 unsigned int n_converged = 0;
718 solve(n_eigenpairs, &n_converged);
719
720 if (n_converged > n_eigenpairs)
721 n_converged = n_eigenpairs;
722 AssertThrow(n_converged == n_eigenpairs,
724 n_eigenpairs));
725
727 eigenvectors.resize(n_converged, eigenvectors.front());
728 eigenvalues.resize(n_converged);
729
730 for (unsigned int index = 0; index < n_converged; ++index)
731 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
732 }
733
734 template <typename OutputVector>
735 void
738 std::vector<PetscScalar> & eigenvalues,
739 std::vector<OutputVector> & eigenvectors,
740 const unsigned int n_eigenpairs)
741 {
742 // Guard against incompatible matrix sizes:
743 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
744 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
745
746 // Panic if the number of eigenpairs wanted is out of bounds.
747 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
749
750 // Set the matrices of the problem
751 set_matrices(A, B);
752
753 // and solve
754 unsigned int n_converged = 0;
755 solve(n_eigenpairs, &n_converged);
756
757 if (n_converged >= n_eigenpairs)
758 n_converged = n_eigenpairs;
759
760 AssertThrow(n_converged == n_eigenpairs,
762 n_eigenpairs));
764
765 eigenvectors.resize(n_converged, eigenvectors.front());
766 eigenvalues.resize(n_converged);
767
768 for (unsigned int index = 0; index < n_converged; ++index)
769 get_eigenpair(index, eigenvalues[index], eigenvectors[index]);
770 }
771
772 template <typename OutputVector>
773 void
776 std::vector<double> & real_eigenvalues,
777 std::vector<double> & imag_eigenvalues,
778 std::vector<OutputVector> & real_eigenvectors,
779 std::vector<OutputVector> & imag_eigenvectors,
780 const unsigned int n_eigenpairs)
781 {
782 // Guard against incompatible matrix sizes:
783 AssertThrow(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
784 AssertThrow(A.n() == B.n(), ExcDimensionMismatch(A.n(), B.n()));
785
786 // and incompatible eigenvalue/eigenvector sizes
787 AssertThrow(real_eigenvalues.size() == imag_eigenvalues.size(),
788 ExcDimensionMismatch(real_eigenvalues.size(),
789 imag_eigenvalues.size()));
790 AssertThrow(real_eigenvectors.size() == imag_eigenvectors.size(),
791 ExcDimensionMismatch(real_eigenvectors.size(),
792 imag_eigenvectors.size()));
793
794 // Panic if the number of eigenpairs wanted is out of bounds.
795 AssertThrow((n_eigenpairs > 0) && (n_eigenpairs <= A.m()),
797
798 // Set the matrices of the problem
799 set_matrices(A, B);
800
801 // and solve
802 unsigned int n_converged = 0;
803 solve(n_eigenpairs, &n_converged);
804
805 if (n_converged >= n_eigenpairs)
806 n_converged = n_eigenpairs;
807
808 AssertThrow(n_converged == n_eigenpairs,
810 n_eigenpairs));
811 AssertThrow((real_eigenvectors.size() != 0) &&
812 (imag_eigenvectors.size() != 0),
814
815 real_eigenvectors.resize(n_converged, real_eigenvectors.front());
816 imag_eigenvectors.resize(n_converged, imag_eigenvectors.front());
817 real_eigenvalues.resize(n_converged);
818 imag_eigenvalues.resize(n_converged);
819
820 for (unsigned int index = 0; index < n_converged; ++index)
821 get_eigenpair(index,
822 real_eigenvalues[index],
823 imag_eigenvalues[index],
824 real_eigenvectors[index],
825 imag_eigenvectors[index]);
826 }
827
828 template <typename Vector>
829 void
830 SolverBase::set_initial_space(const std::vector<Vector> &this_initial_space)
831 {
832 std::vector<Vec> vecs(this_initial_space.size());
833
834 for (unsigned int i = 0; i < this_initial_space.size(); ++i)
835 {
836 Assert(this_initial_space[i].l2_norm() > 0.0,
837 ExcMessage("Initial vectors should be nonzero."));
838 vecs[i] = this_initial_space[i];
839 }
840
841 // if the eigensolver supports only a single initial vector, but several
842 // guesses are provided, then all except the first one will be discarded.
843 // One could still build a vector that is rich in the directions of all
844 // guesses, by taking a linear combination of them. (TODO: make function
845 // virtual?)
846
847 const PetscErrorCode ierr =
848 EPSSetInitialSpace(eps, vecs.size(), vecs.data());
849 AssertThrow(ierr == 0, ExcSLEPcError(ierr));
850 }
851
852} // namespace SLEPcWrappers
853
855
856# endif // DEAL_II_WITH_SLEPC
857
858/*---------------------------- slepc_solver.h ---------------------------*/
859
860#endif
861
862/*---------------------------- slepc_solver.h ---------------------------*/
const AdditionalData additional_data
Definition: slepc_solver.h:466
SolverControl & control() const
EPSConvergedReason reason
Definition: slepc_solver.h:360
void set_target_eigenvalue(const PetscScalar &this_target)
const MPI_Comm mpi_communicator
Definition: slepc_solver.h:302
void set_matrices(const PETScWrappers::MatrixBase &A)
Definition: slepc_solver.cc:77
void get_solver_state(const SolverControl::State state)
void set_initial_space(const std::vector< Vector > &initial_space)
Definition: slepc_solver.h:830
void set_problem_type(EPSProblemType set_problem)
SolverControl & solver_control
Definition: slepc_solver.h:297
void solve(const PETScWrappers::MatrixBase &A, std::vector< PetscScalar > &eigenvalues, std::vector< OutputVector > &eigenvectors, const unsigned int n_eigenpairs=1)
Definition: slepc_solver.h:704
void get_eigenpair(const unsigned int index, PetscScalar &eigenvalues, PETScWrappers::VectorBase &eigenvectors)
void set_transformation(SLEPcWrappers::TransformationBase &this_transformation)
Definition: slepc_solver.cc:98
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
Definition: slepc_solver.h:649
const AdditionalData additional_data
Definition: slepc_solver.h:415
const AdditionalData additional_data
Definition: slepc_solver.h:689
const AdditionalData additional_data
Definition: slepc_solver.h:518
const AdditionalData additional_data
Definition: slepc_solver.h:558
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DeclException0(Exception0)
Definition: exceptions.h:464
static ::ExceptionBase & ExcSLEPcEigenvectorConvergenceMismatchError(int arg1, int arg2)
static ::ExceptionBase & ExcSLEPcError(int arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:532
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
static ::ExceptionBase & ExcSLEPcWrappersUsageError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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)