eigen_solver.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2025 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for a class that defines interfaces to Eigensolvers
 
// Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_EIGEN_SOLVER_HEADER
#define OOMPH_EIGEN_SOLVER_HEADER
 
// Include the header
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif
 
#include <complex>
#include "Vector.h"
#include "complex_matrices.h"
 
namespace oomph
{
  // Forward definition of problem class
  class Problem;
 
  // Forward definition of matrix class
  class DoubleMatrixBase;
 
  // Forward definition of linear solver class
  class LinearSolver;
 
  //=======================================================================
  /// Base class for all EigenProblem solves. This simply defines standard
  /// interfaces so that different solvers can be used easily.
  //=======================================================================
  class EigenSolver : public DistributableLinearAlgebraObject
  {
  protected:
    /// Double value that represents the real part of the shift in
    /// shifted eigensolvers
    double Sigma_real;
 
  public:
    /// Empty constructor
    EigenSolver() : Sigma_real(0.0) {}
 
    /// Empty copy constructor
    EigenSolver(const EigenSolver&) {}
 
    /// Empty destructor
    virtual ~EigenSolver() {}
 
    /// Solve the real eigenproblem that is assembled by elements in
    /// a mesh in a Problem object. Note that the assembled matrices include the
    /// shift and are real. The eigenvalues and eigenvectors are,
    /// in general, complex, and the eigenvalues can be NaNs or Infs.
    /// This function is therefore merely provided as a convenience, to be
    /// used if the user is confident that NaNs don't arise (e.g. in Arnoldi
    /// based solvers where typically only a small number of (finite)
    /// eigenvalues are computed), or if the users is happy to deal with NaNs in
    /// the subsequent post-processing.
    /// Function is virtual so it can be overloaded for Arnoldi type solvers
    /// that compute the (finite) eigenvalues directly
    /// At least n_eval eigenvalues are computed.
    virtual void solve_eigenproblem(Problem* const& problem_pt,
                                    const int& n_eval,
                                    Vector<std::complex<double>>& eigenvalue,
                                    Vector<DoubleVector>& eigenvector_real,
                                    Vector<DoubleVector>& eigenvector_imag,
                                    const bool& do_adjoint_problem = false)
    {
      Vector<std::complex<double>> alpha;
      Vector<double> beta;
 
      // Call the "safe" version
      solve_eigenproblem(problem_pt,
                         n_eval,
                         alpha,
                         beta,
                         eigenvector_real,
                         eigenvector_imag,
                         do_adjoint_problem);
 
      // Now do the brute force conversion, possibly creating NaNs and Infs...
      unsigned n = alpha.size();
      eigenvalue.resize(n);
      for (unsigned i = 0; i < n; i++)
      {
        eigenvalue[i] = alpha[i] / beta[i];
      }
    }
 
    /// Solve the real eigenproblem that is assembled by elements in
    /// a mesh in a Problem object. Note that the assembled matrices include the
    /// shift and are real. The eigenvalues and eigenvectors are,
    /// in general, complex. Eigenvalues may be infinite and are therefore
    /// returned as
    /// \f$ \lambda_i = \alpha_i / \beta_i \f$ where \f$ \alpha_i \f$ is complex
    /// while \f$ \beta_i \f$ is real. The actual eigenvalues may then be
    /// computed by doing the division, checking for zero betas to avoid NaNs.
    /// There's a convenience wrapper to this function that simply computes
    /// these eigenvalues regardless. That version may die in NaN checking is
    /// enabled (via the fenv.h header and the associated feenable function).
    /// At least n_eval eigenvalues are computed.
    virtual void solve_eigenproblem(Problem* const& problem_pt,
                                    const int& n_eval,
                                    Vector<std::complex<double>>& alpha,
                                    Vector<double>& beta,
                                    Vector<DoubleVector>& eigenvector_real,
                                    Vector<DoubleVector>& eigenvector_imag,
                                    const bool& do_adjoint_problem = false) = 0;
 
    /// Set the value of the (real) shift
    void set_shift(const double& shift_value)
    {
      Sigma_real = shift_value;
    }
 
    /// Return the value of the (real) shift (const version)
    const double& get_shift() const
    {
      return Sigma_real;
    }
  };
 
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////
 
  //=====================================================================
  /// Class for the LAPACK QZ eigensolver
  //=====================================================================
  class LAPACK_QZ : public EigenSolver
  {
  public:
    /// Empty constructor
    LAPACK_QZ() : EigenSolver(), Tolerance_for_ccness_check(1.0e-12) {}
 
    /// Broken copy constructor
    LAPACK_QZ(const LAPACK_QZ&) = delete;
 
    /// Broken assignment operator
    void operator=(const LAPACK_QZ&) = delete;
 
    /// Empty desctructor
    virtual ~LAPACK_QZ() {}
 
    /// Solve the real eigenproblem that is assembled by elements in
    /// a mesh in a Problem object. Note that the assembled matrices include the
    /// shift and are real. The eigenvalues and eigenvectors are,
    /// in general, complex. Eigenvalues may be infinite and are therefore
    /// returned as
    /// \f$ \lambda_i = \alpha_i / \beta_i \f$ where \f$ \alpha_i \f$ is complex
    /// while \f$ \beta_i \f$ is real. The actual eigenvalues may then be
    /// computed by doing the division, checking for zero betas to avoid NaNs.
    /// There's a convenience wrapper to this function that simply computes
    /// these eigenvalues regardless. That version may die in NaN checking is
    /// enabled (via the fenv.h header and the associated feenable function).
    /// At least n_eval eigenvalues are computed.
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& alpha,
                            Vector<double>& beta,
                            Vector<DoubleVector>& eigenvector_real,
                            Vector<DoubleVector>& eigenvector_imag,
                            const bool& do_adjoint_problem = false);
 
    /// Find the eigenvalues of a complex generalised eigenvalue problem
    /// specified by \f$ Ax = \lambda  Mx \f$. Note: the (real) shift
    /// that's specifiable via the EigenSolver base class is ignored here.
    /// A warning gets issued if it's set to a nonzero value.
    void find_eigenvalues(const ComplexMatrixBase& A,
                          const ComplexMatrixBase& M,
                          Vector<std::complex<double>>& eigenvalue,
                          Vector<Vector<std::complex<double>>>& eigenvector);
 
 
    /// Access to tolerance for checking complex conjugateness of eigenvalues
    /// (const version)
    double tolerance_for_ccness_check() const
    {
      return Tolerance_for_ccness_check;
    }
 
 
    /// Access to tolerance for checking complex conjugateness of eigenvalues
    double& tolerance_for_ccness_check()
    {
      return Tolerance_for_ccness_check;
    }
 
  private:
    /// Helper function called from "raw" lapack
    /// code
    void solve_eigenproblem_helper(Problem* const& problem_pt,
                                   const int& n_eval,
                                   Vector<std::complex<double>>& alpha,
                                   Vector<double>& beta,
                                   Vector<DoubleVector>& eigenvector);
 
 
    /// Provide diagonstic for DGGEV error return
    void DGGEV_error(const int& info, const int& n)
    {
      // Throw an error
      std::ostringstream error_stream;
      error_stream << "Failure in LAPACK_DGGEV(...).\n"
                   << "info = " << info << std::endl;
      error_stream
        << "Diagnostics below are from \n\n"
        << "http://www.netlib.org/lapack/explore-html/d9/d8e/"
           "group__double_g_eeigen_ga4f59e87e670a755b41cbdd7e97f36bea.html"
        << std::endl;
      if (info < 0)
      {
        error_stream << -info << "-th input arg had an illegal value\n";
      }
      else if (info <= n)
      {
        error_stream << "The QZ iteration failed.  No eigenvectors have been\n"
                     << "calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)\n"
                     << "should be correct for j=INFO+1,...,N, where \n"
                     << "info = " << info << " and N = " << n << std::endl;
      }
      else if (info == (n + 1))
      {
        error_stream << "QZ iteration failed in DHGEQZ.\n";
      }
      else if (info == (n + 2))
      {
        error_stream << "error return from DTGEVC.\n";
      }
      error_stream
        << "Aborting here; if you know how to proceed then\n"
        << "then implement ability to catch this error and continue\n";
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Provide diagonstic for ZGGEV error return
    void ZGGEV_error(const int& info, const int& n)
    {
      // Throw an error
      std::ostringstream error_stream;
      error_stream << "Failure in LAPACK_ZGGEV(...).\n"
                   << "info = " << info << std::endl;
      error_stream << "Diagnostics below are from \n\n"
                   << "http://www.netlib.org/lapack/explore-html/"
                   << "db/d55/group__complex16_g_eeigen_ga79fcce20c"
                   << "617429ccf985e6f123a6171.html" << std::endl;
      if (info < 0)
      {
        error_stream << -info << "-th input arg had an illegal value\n";
      }
      else if (info <= n)
      {
        error_stream << "The QZ iteration failed.  No eigenvectors have been\n"
                     << "calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)\n"
                     << "should be correct for j=INFO+1,...,N, where \n"
                     << "info = " << info << " and N = " << n << std::endl;
      }
      else if (info == (n + 1))
      {
        error_stream << "QZ iteration failed in ZHGEQZ.\n";
      }
      else if (info == (n + 2))
      {
        error_stream << "error return from ZTGEVC.\n";
      }
      error_stream
        << "Aborting here; if you know how to proceed then\n"
        << "then implement ability to catch this error and continue\n";
 
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Tolerance for checking complex conjugateness of eigenvalues
    double Tolerance_for_ccness_check;
  };
 
} // namespace oomph
 
#endif