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;
 
  //=======================================================================
  //=======================================================================
  class EigenSolver : public DistributableLinearAlgebraObject
  {
  protected:
    double Sigma_real;
 
  public:
    EigenSolver() : Sigma_real(0.0) {}
 
    EigenSolver(const EigenSolver&) {}
 
    virtual ~EigenSolver() {}
 
    virtual void solve_eigenproblem(Problem* const& problem_pt,
                                    const int& n_eval,
                                    Vector<std::complex<double>>& eigenvalue,
                                    Vector<DoubleVector>& eigenvector) = 0;
 
 
    void set_shift(const double& shift_value)
    {
      Sigma_real = shift_value;
    }
 
    const double& get_shift() const
    {
      return Sigma_real;
    }
  };
 
 
  //=====================================================================
  //=====================================================================
  class ARPACK : public EigenSolver
  {
  private:
    LinearSolver* Linear_solver_pt;
 
    LinearSolver* Default_linear_solver_pt;
 
    int Spectrum;
 
 
    int NArnoldi;
 
 
    bool Small;
 
    bool Compute_eigenvectors;
 
 
  public:
    ARPACK();
 
    ARPACK(const ARPACK&) {}
 
    virtual ~ARPACK();
 
    int& narnoldi()
    {
      return NArnoldi;
    }
 
    const int& narnoldi() const
    {
      return NArnoldi;
    }
 
    void enable_compute_eigenvectors()
    {
      Compute_eigenvectors = true;
    }
 
    void disable_compute_eigenvectors()
    {
      Compute_eigenvectors = false;
    }
 
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector);
 
    // void find_eigenvalues(const DoubleMatrixBase &A, const int &n_eval,
    //                      Vector<std::complex<double> > &eigenvalue,
    //                      Vector<Vector<double> > &eigenvector);
 
 
    void get_eigenvalues_left_of_shift()
    {
      Small = true;
    }
 
    void get_eigenvalues_right_of_shift()
    {
      Small = false;
    }
 
    void track_eigenvalue_real_part()
    {
      Spectrum = 1;
    }
 
    void track_eigenvalue_imaginary_part()
    {
      Spectrum = -1;
    }
 
    void track_eigenvalue_magnitude()
    {
      Spectrum = 0;
    }
 
    LinearSolver*& linear_solver_pt()
    {
      return Linear_solver_pt;
    }
 
    LinearSolver* const& linear_solver_pt() const
    {
      return Linear_solver_pt;
    }
  };
 
 
  //=====================================================================
  //=====================================================================
  class LAPACK_QZ : public EigenSolver
  {
  public:
    LAPACK_QZ() : EigenSolver() {}
 
    LAPACK_QZ(const LAPACK_QZ&) {}
 
    virtual ~LAPACK_QZ() {}
 
    void solve_eigenproblem(Problem* const& problem_pt,
                            const int& n_eval,
                            Vector<std::complex<double>>& eigenvalue,
                            Vector<DoubleVector>& eigenvector);
 
    void find_eigenvalues(const ComplexMatrixBase& A,
                          const ComplexMatrixBase& M,
                          Vector<std::complex<double>>& eigenvalue,
                          Vector<Vector<std::complex<double>>>& eigenvector);
 
    void get_eigenvalues_right_of_shift() {}
  };
 
} // namespace oomph
 
#endif