Line data    Source code

       1              : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
       2              : // SPDX-License-Identifier: Apache-2.0
       3              : 
       4              : #ifndef PALACE_LINALG_EPS_HPP
       5              : #define PALACE_LINALG_EPS_HPP
       6              : 
       7              : #include <complex>
       8              : #include "linalg/ksp.hpp"
       9              : #include "linalg/operator.hpp"
      10              : #include "linalg/vector.hpp"
      11              : 
      12              : namespace palace
      13              : {
      14              : 
      15              : template <typename VecType>
      16              : class DivFreeSolver;
      17              : 
      18              : //
      19              : // Pure abstract base class for solving generalized linear eigenvalue problems problems or
      20              : // quadratic polynomial eigenvalue problems.
      21              : //
      22              : class EigenvalueSolver
      23              : {
      24              : public:
      25              :   enum class ScaleType
      26              :   {
      27              :     NONE,
      28              :     NORM_2
      29              :   };
      30              : 
      31              :   enum class WhichType
      32              :   {
      33              :     LARGEST_MAGNITUDE,
      34              :     SMALLEST_MAGNITUDE,
      35              :     LARGEST_REAL,
      36              :     SMALLEST_REAL,
      37              :     LARGEST_IMAGINARY,
      38              :     SMALLEST_IMAGINARY,
      39              :     TARGET_MAGNITUDE,
      40              :     TARGET_REAL,
      41              :     TARGET_IMAGINARY
      42              :   };
      43              : 
      44              :   enum class ErrorType
      45              :   {
      46              :     ABSOLUTE,
      47              :     RELATIVE,
      48              :     BACKWARD
      49              :   };
      50              : 
      51              : public:
      52              :   EigenvalueSolver() = default;
      53              :   virtual ~EigenvalueSolver() = default;
      54              : 
      55              :   // Set operators for the generalized eigenvalue problem, quadratic polynomial
      56              :   // eigenvalue problem, or nonlinear eigenvalue problem.
      57            0 :   virtual void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
      58              :                             ScaleType type)
      59              :   {
      60            0 :     MFEM_ABORT("SetOperators not defined!");
      61              :   }
      62              : 
      63            0 :   virtual void SetOperators(const ComplexOperator &K, const ComplexOperator &C,
      64              :                             const ComplexOperator &M, ScaleType type)
      65              :   {
      66            0 :     MFEM_ABORT("SetOperators not defined!");
      67              :   }
      68              : 
      69            0 :   virtual void SetOperators(const ComplexOperator &K, const ComplexOperator &M,
      70              :                             std::function<const ComplexOperator &(std::complex<double>)> A2,
      71              :                             ScaleType type)
      72              :   {
      73            0 :     MFEM_ABORT("SetOperators not defined!");
      74              :   }
      75              : 
      76            0 :   virtual void SetExtraSystemMatrix(std::function<std::unique_ptr<ComplexOperator>(double)>)
      77              :   {
      78            0 :     MFEM_ABORT("SetExtraSystemMatrix not defined!");
      79              :   }
      80              : 
      81            0 :   virtual void SetPreconditionerUpdate(
      82              :       std::function<std::unique_ptr<ComplexOperator>(
      83              :           std::complex<double>, std::complex<double>, std::complex<double>, double)>)
      84              :   {
      85            0 :     MFEM_ABORT("SetPreconditionerUpdate not defined!");
      86              :   }
      87              : 
      88              :   // For the linear generalized case, the linear solver should be configured to compute the
      89              :   // action of M⁻¹ (with no spectral transformation) or (K - σ M)⁻¹. For the quadratic
      90              :   // case, the linear solver should be configured to compute the action of M⁻¹ (with no
      91              :   // spectral transformation) or P(σ)⁻¹.
      92              :   virtual void SetLinearSolver(ComplexKspSolver &ksp) = 0;
      93              : 
      94              :   // Set the projection operator for enforcing the divergence-free constraint.
      95              :   virtual void SetDivFreeProjector(const DivFreeSolver<ComplexVector> &divfree) = 0;
      96              : 
      97              :   // Set optional B matrix used for weighted inner products. This must be set explicitly
      98              :   // even for generalized problems, otherwise the identity will be used.
      99              :   virtual void SetBMat(const Operator &B) = 0;
     100              : 
     101              :   // Get scaling factors used by the solver.
     102              :   virtual double GetScalingGamma() const = 0;
     103              :   virtual double GetScalingDelta() const = 0;
     104              : 
     105              :   // Set the number of required eigenmodes.
     106              :   virtual void SetNumModes(int num_eig, int num_vec = 0) = 0;
     107              : 
     108              :   // Set solver tolerance.
     109              :   virtual void SetTol(double tol) = 0;
     110              : 
     111              :   // Set maximum number of Arnoldi update iterations.
     112              :   virtual void SetMaxIter(int max_it) = 0;
     113              : 
     114              :   // Set target spectrum for the eigensolver. When a spectral transformation is used, this
     115              :   // applies to the spectrum of the shifted operator.
     116              :   virtual void SetWhichEigenpairs(WhichType type) = 0;
     117              : 
     118              :   // Set shift-and-invert spectral transformation.
     119              :   virtual void SetShiftInvert(std::complex<double> s, bool precond = false) = 0;
     120              : 
     121              :   // Set an initial vector for the solution subspace.
     122              :   virtual void SetInitialSpace(const ComplexVector &v) = 0;
     123              : 
     124              :   // Solve the eigenvalue problem. Returns the number of converged eigenvalues.
     125              :   virtual int Solve() = 0;
     126              : 
     127              :   // Get the corresponding eigenvalue.
     128              :   virtual std::complex<double> GetEigenvalue(int i) const = 0;
     129              : 
     130              :   // Get the corresponding eigenvector. Eigenvectors are normalized such that ||x||₂ = 1,
     131              :   // unless the B-matrix is set for weighted inner products.
     132              :   virtual void GetEigenvector(int i, ComplexVector &x) const = 0;
     133              : 
     134              :   // Get the corresponding eigenpair error.
     135              :   virtual double GetError(int i, ErrorType type) const = 0;
     136              : 
     137              :   // Re-normalize the given number of eigenvectors, for example if the matrix B for weighted
     138              :   // inner products has changed. This does not perform re-orthogonalization with respect to
     139              :   // the new matrix, only normalization.
     140              :   virtual void RescaleEigenvectors(int num_eig) = 0;
     141              : };
     142              : 
     143              : }  // namespace palace
     144              : 
     145              : #endif  // PALACE_LINALG_EPS_HPP