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//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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:

 /// \short 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() {}

 /// \short 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:

 /// \short Pointer to a linear solver
 LinearSolver* Linear_solver_pt;

 /// \short Pointer to a default linear solver
 LinearSolver* Default_linear_solver_pt;

 /// \short Integer to set whether the real, imaginary or magnitude is required
 ///to be small or large.
 int Spectrum;


 /// \short Number of Arnoldi vectors to compute
 int NArnoldi;


 ///\short Boolean to set which part of the spectrum left (default) or right
 /// of the shifted value.
 bool Small;

 /// \short 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;}

 /// \short Set to enable the computation of the eigenvectors (default)
 void enable_compute_eigenvectors() {Compute_eigenvectors=true;}

 /// \short 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() { }

};

}

#endif