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-2022 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 generated by autoconfig
#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() {}
 
    /// Actual eigensolver. This takes a pointer to a problem and returns
    /// a vector of complex numbers representing the eigenvalues
    /// and a corresponding vector of eigenvectors
    virtual void solve_eigenproblem(Problem* const& problem_pt,
                                    const int& n_eval,
                                    Vector<std::complex<double>>& eigenvalue,
                                    Vector<DoubleVector>& eigenvector) = 0;
 
 
    /// Set the value of the shift
    void set_shift(const double& shift_value)
    {
      Sigma_real = shift_value;
    }
 
    /// Return the value of the shift (const version)
    const double& get_shift() const
    {
      return Sigma_real;
    }
  };
 
 
  //=====================================================================
  /// Class for the ARPACK eigensolver
  //=====================================================================
  class ARPACK : public EigenSolver
  {
  private:
    /// Pointer to a linear solver
    LinearSolver* Linear_solver_pt;
 
    /// Pointer to a default linear solver
    LinearSolver* Default_linear_solver_pt;
 
    /// Integer to set whether the real, imaginary or magnitude is
    /// required
    /// to be small or large.
    int Spectrum;
 
 
    /// Number of Arnoldi vectors to compute
    int NArnoldi;
 
 
    /// Boolean to set which part of the spectrum left (default) or right
    /// of the shifted value.
    bool Small;
 
    /// Boolean to indicate whether or not to compute the eigenvectors
    bool Compute_eigenvectors;
 
 
  public:
    /// Constructor
    ARPACK();
 
    /// Empty copy constructor
    ARPACK(const ARPACK&) {}
 
    /// Destructor, delete the linear solver
    virtual ~ARPACK();
 
    /// Access function for the number of Arnoldi vectors
    int& narnoldi()
    {
      return NArnoldi;
    }
 
    /// Access function for the number of Arnoldi vectors (const version)
    const int& narnoldi() const
    {
      return NArnoldi;
    }
 
    /// Set to enable the computation of the eigenvectors (default)
    void enable_compute_eigenvectors()
    {
      Compute_eigenvectors = true;
    }
 
    /// Set to disable the computation of the eigenvectors
    void disable_compute_eigenvectors()
    {
      Compute_eigenvectors = false;
    }
 
    /// Solve the eigen problem
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector);
 
    /// Use the eigensolver to find the eigenvalues of a given matrix
    // void find_eigenvalues(const DoubleMatrixBase &A, const int &n_eval,
    //                      Vector<std::complex<double> > &eigenvalue,
    //                      Vector<Vector<double> > &eigenvector);
 
 
    /// Set the desired eigenvalues to be left of the shift
    void get_eigenvalues_left_of_shift()
    {
      Small = true;
    }
 
    /// Set the desired eigenvalues to be right of the shift
    void get_eigenvalues_right_of_shift()
    {
      Small = false;
    }
 
    /// Set the real part to be the quantity of interest (default)
    void track_eigenvalue_real_part()
    {
      Spectrum = 1;
    }
 
    /// Set the imaginary part fo the quantity of interest
    void track_eigenvalue_imaginary_part()
    {
      Spectrum = -1;
    }
 
    /// Set the magnitude to be the quantity of interest
    void track_eigenvalue_magnitude()
    {
      Spectrum = 0;
    }
 
    /// Return a pointer to the linear solver object
    LinearSolver*& linear_solver_pt()
    {
      return Linear_solver_pt;
    }
 
    /// Return a pointer to the linear solver object (const version)
    LinearSolver* const& linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
 
  //=====================================================================
  /// Class for the LAPACK eigensolver
  //=====================================================================
  class LAPACK_QZ : public EigenSolver
  {
  public:
    /// Empty constructor
    LAPACK_QZ() : EigenSolver() {}
 
    /// Empty copy constructor
    LAPACK_QZ(const LAPACK_QZ&) {}
 
    /// Empty desctructor
    virtual ~LAPACK_QZ() {}
 
    /// Solve the eigen problem
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector);
 
    /// Find the eigenvalues of a generalised eigenvalue problem
    /// specified by \f$ Ax = \lambda  Mx \f$
    void find_eigenvalues(const ComplexMatrixBase& A,
                          const ComplexMatrixBase& M,
                          Vector<std::complex<double>>& eigenvalue,
                          Vector<Vector<std::complex<double>>>& eigenvector);
 
    /// Set the desired eigenvalues to be right of the shift
    /// Dummy at the moment
    void get_eigenvalues_right_of_shift() {}
  };
 
} // namespace oomph
 
#endif