oomph::LAPACK_QZ Class Reference

Class for the LAPACK eigensolver. More...

#include <eigen_solver.h>

Inheritance diagram for oomph::LAPACK_QZ:

Public Member Functions
	LAPACK_QZ ()
	Empty constructor. More...
	LAPACK_QZ (const LAPACK_QZ &)
	Empty copy constructor. More...
virtual	~LAPACK_QZ ()
	Empty desctructor. More...
void	solve_eigenproblem (Problem *const &problem_pt, const int &n_eval, Vector< std::complex< double >> &eigenvalue, Vector< DoubleVector > &eigenvector)
	Solve the eigen problem. More...
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 ()
Public Member Functions inherited from oomph::EigenSolver
	EigenSolver ()
	Empty constructor. More...
	EigenSolver (const EigenSolver &)
	Empty copy constructor. More...
virtual	~EigenSolver ()
	Empty destructor. More...
void	set_shift (const double &shift_value)
	Set the value of the shift. More...
const double &	get_shift () const
	Return the value of the shift (const version) More...
Public Member Functions inherited from oomph::DistributableLinearAlgebraObject
	DistributableLinearAlgebraObject ()
	Default constructor - create a distribution. More...
	DistributableLinearAlgebraObject (const DistributableLinearAlgebraObject &matrix)=delete
	Broken copy constructor. More...
void	operator= (const DistributableLinearAlgebraObject &)=delete
	Broken assignment operator. More...
virtual	~DistributableLinearAlgebraObject ()
	Destructor. More...
LinearAlgebraDistribution *	distribution_pt () const
	access to the LinearAlgebraDistribution More...
unsigned	nrow () const
	access function to the number of global rows. More...
unsigned	nrow_local () const
	access function for the num of local rows on this processor. More...
unsigned	nrow_local (const unsigned &p) const
	access function for the num of local rows on this processor. More...
unsigned	first_row () const
	access function for the first row on this processor More...
unsigned	first_row (const unsigned &p) const
	access function for the first row on this processor More...
bool	distributed () const
	distribution is serial or distributed More...
bool	distribution_built () const
void	build_distribution (const LinearAlgebraDistribution *const dist_pt)
void	build_distribution (const LinearAlgebraDistribution &dist)

Additional Inherited Members
Protected Member Functions inherited from oomph::DistributableLinearAlgebraObject
void	clear_distribution ()
Protected Attributes inherited from oomph::EigenSolver
double	Sigma_real

Detailed Description

Class for the LAPACK eigensolver.

Constructor & Destructor Documentation

LAPACK_QZ() [1/2]

oomph::LAPACK_QZ::LAPACK_QZ ( )

inline

Empty constructor.

: EigenSolver() {}

LAPACK_QZ() [2/2]

oomph::LAPACK_QZ::LAPACK_QZ ( const LAPACK_QZ &  )

inline

Empty copy constructor.

{}

~LAPACK_QZ()

virtual oomph::LAPACK_QZ::~LAPACK_QZ ( )

inlinevirtual

Empty desctructor.

{}

Member Function Documentation

find_eigenvalues()

void oomph::LAPACK_QZ::find_eigenvalues	(	const ComplexMatrixBase &	A,
		const ComplexMatrixBase &	M,
		Vector< std::complex< double >> &	eigenvalue,
		Vector< Vector< std::complex< double >>> &	eigenvector
	)

Find the eigenvalues of a generalised eigenvalue problem specified by \( Ax = \lambda Mx \)

Use LAPACK to solve a complex eigen problem specified by the given matrices.

  {
    // Some character identifiers for use in the LAPACK routine
    // Do not calculate the left eigenvectors
    char no_eigvecs[2] = "N";
    // Do caculate the eigenvectors
    char eigvecs[2] = "V";
 
    // Get the dimension of the matrix
    int n = A.nrow(); // Total size of matrix
 
    // Storage for the matrices in the packed form required by the LAPACK
    // routine
    double* M_linear = new double[2 * n * n];
    double* A_linear = new double[2 * n * n];
 
    // Now convert the matrices into the appropriate packed form
    unsigned index = 0;
    for (int i = 0; i < n; ++i)
    {
      for (int j = 0; j < n; ++j)
      {
        M_linear[index] = M(j, i).real();
        A_linear[index] = A(j, i).real();
        ++index;
        M_linear[index] = M(j, i).imag();
        A_linear[index] = A(j, i).imag();
        ++index;
      }
    }
 
    // Temporary eigenvalue storage
    double* alpha = new double[2 * n];
    double* beta = new double[2 * n];
    // Temporary eigenvector storage
    double* vec_left = new double[2];
    double* vec_right = new double[2 * n * n];
 
    // Workspace for the LAPACK routine
    std::vector<double> work(2, 0.0);
    std::vector<double> rwork(8 * n, 0.0);
 
    // Info flag for the LAPACK routine
    int info = 0;
 
    // Call FORTRAN LAPACK to get the required workspace
    // Note the use of the padding flag
    LAPACK_ZGGEV(no_eigvecs,
                 eigvecs,
                 n,
                 &A_linear[0],
                 n,
                 &M_linear[0],
                 n,
                 alpha,
                 beta,
                 vec_left,
                 1,
                 vec_right,
                 n,
                 &work[0],
                 -1,
                 &rwork[0],
                 info);
 
    // Get the amount of requires workspace
    int required_workspace = (int)work[0];
    // If we need it
    work.resize(2 * required_workspace);
 
    // call FORTRAN LAPACK again with the optimum workspace
    LAPACK_ZGGEV(no_eigvecs,
                 eigvecs,
                 n,
                 &A_linear[0],
                 n,
                 &M_linear[0],
                 n,
                 alpha,
                 beta,
                 vec_left,
                 1,
                 vec_right,
                 n,
                 &work[0],
                 required_workspace,
                 &rwork[0],
                 info);
 
    // Now resize storage for the eigenvalues and eigenvectors
    // We get them all!
    eigenvalue.resize(n);
    eigenvector.resize(n);
 
    // Now convert the output into our format
    for (int i = 0; i < n; ++i)
    {
      // We have temporary complex numbers
      std::complex<double> num(alpha[2 * i], alpha[2 * i + 1]);
      std::complex<double> den(beta[2 * i], beta[2 * i + 1]);
 
      // N.B. This is naive and dangerous according to the documentation
      // beta could be very small giving over or under flow
      eigenvalue[i] = num / den;
 
      // Resize the eigenvector storage
      eigenvector[i].resize(n);
      // Load up the eigenvector (assume that it's real)
      for (int k = 0; k < n; ++k)
      {
        eigenvector[i][k] = std::complex<double>(
          vec_right[2 * i * n + 2 * k], vec_right[2 * i * n + 2 * k + 1]);
      }
    }
 
    // Delete all allocated storage
    delete[] vec_right;
    delete[] vec_left;
    delete[] beta;
    delete[] alpha;
    delete[] A_linear;
    delete[] M_linear;
  }

References alpha, beta, i, info, int(), j, k, LAPACK_ZGGEV, and n.

get_eigenvalues_right_of_shift()

void oomph::LAPACK_QZ::get_eigenvalues_right_of_shift ( )

inline

Set the desired eigenvalues to be right of the shift Dummy at the moment

{}

solve_eigenproblem()

void oomph::LAPACK_QZ::solve_eigenproblem ( Problem *const &  problem_pt, const int &  n_eval, Vector< std::complex< double >> &  eigenvalue, Vector< DoubleVector > &  eigenvector  )

virtual

Solve the eigen problem.

Use LAPACK to solve an eigen problem that is assembled by elements in a mesh in a Problem object.

Implements oomph::EigenSolver.

  {
    // Some character identifiers for use in the LAPACK routine
    // Do not calculate the left eigenvectors
    char no_eigvecs[2] = "N";
    // Do caculate the eigenvectors
    char eigvecs[2] = "V";
 
    // Get the dimension of the matrix
    int n = problem_pt->ndof(); // Total size of matrix
 
    // Initialise a padding integer
    int padding = 0;
    // If the dimension of the matrix is even, then pad the arrays to
    // make the size odd. This somehow sorts out a strange run-time behaviour
    // identified by Rich Hewitt.
    if (n % 2 == 0)
    {
      padding = 1;
    }
 
    // Get the real and imaginary parts of the shift
    double sigmar = Sigma_real; // sigmai = 0.0;
 
    // Actual size of matrix that will be allocated
    int padded_n = n + padding;
 
    // Storage for the matrices in the packed form required by the LAPACK
    // routine
    double* M = new double[padded_n * padded_n];
    double* A = new double[padded_n * padded_n];
 
    // TEMPORARY
    // only use non-distributed matrices and vectors
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n, false);
    this->build_distribution(dist);
 
    // Enclose in a separate scope so that memory is cleaned after assembly
    {
      // Allocated Row compressed matrices for the mass matrix and shifted main
      // matrix
      CRDoubleMatrix temp_M(this->distribution_pt()),
        temp_AsigmaM(this->distribution_pt());
 
      // Assemble the matrices
      problem_pt->get_eigenproblem_matrices(temp_M, temp_AsigmaM, sigmar);
 
      // Now convert these matrices into the appropriate packed form
      unsigned index = 0;
      for (int i = 0; i < n; ++i)
      {
        for (int j = 0; j < n; ++j)
        {
          M[index] = temp_M(j, i);
          A[index] = temp_AsigmaM(j, i);
          ++index;
        }
        // If necessary pad the columns with zeros
        if (padding)
        {
          M[index] = 0.0;
          A[index] = 0.0;
          ++index;
        }
      }
      // No need to pad the final row because it is never accessed by the
      // routine.
    }
 
    // Temporary eigenvalue storage
    double* alpha_r = new double[n];
    double* alpha_i = new double[n];
    double* beta = new double[n];
    // Temporary eigenvector storage
    double* vec_left = new double[1];
    double* vec_right = new double[n * n];
 
    // Workspace for the LAPACK routine
    std::vector<double> work(1, 0.0);
    // Info flag for the LAPACK routine
    int info = 0;
 
    // Call FORTRAN LAPACK to get the required workspace
    // Note the use of the padding flag
    LAPACK_DGGEV(no_eigvecs,
                 eigvecs,
                 n,
                 &A[0],
                 padded_n,
                 &M[0],
                 padded_n,
                 alpha_r,
                 alpha_i,
                 beta,
                 vec_left,
                 1,
                 vec_right,
                 n,
                 &work[0],
                 -1,
                 info);
 
    // Get the amount of requires workspace
    int required_workspace = (int)work[0];
    // If we need it
    work.resize(required_workspace);
 
    // call FORTRAN LAPACK again with the optimum workspace
    LAPACK_DGGEV(no_eigvecs,
                 eigvecs,
                 n,
                 &A[0],
                 padded_n,
                 &M[0],
                 padded_n,
                 alpha_r,
                 alpha_i,
                 beta,
                 vec_left,
                 1,
                 vec_right,
                 n,
                 &work[0],
                 required_workspace,
                 info);
 
    // Now resize storage for the eigenvalues and eigenvectors
    // We get them all!
    eigenvalue.resize(n);
    eigenvector.resize(n);
 
    // Now convert the output into our format
    for (int i = 0; i < n; ++i)
    {
      // N.B. This is naive and dangerous according to the documentation
      // beta could be very small giving over or under flow
      // Remember to fix the shift back again
      eigenvalue[i] = std::complex<double>(sigmar + alpha_r[i] / beta[i],
                                           alpha_i[i] / beta[i]);
 
      // Resize the eigenvector  storage
      eigenvector[i].build(this->distribution_pt(), 0.0);
      // Load up the eigenvector (assume that it's real)
      for (int k = 0; k < n; ++k)
      {
        eigenvector[i][k] = vec_right[i * n + k];
      }
    }
 
    // Delete all allocated storage
    delete[] vec_right;
    delete[] vec_left;
    delete[] beta;
    delete[] alpha_r;
    delete[] alpha_i;
    delete[] A;
    delete[] M;
  }

References beta, oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::Problem::communicator_pt(), oomph::DistributableLinearAlgebraObject::distribution_pt(), oomph::Problem::get_eigenproblem_matrices(), i, info, int(), j, k, LAPACK_DGGEV, n, oomph::Problem::ndof(), and oomph::EigenSolver::Sigma_real.

---

The documentation for this class was generated from the following files:

• eigen_solver.h
• eigen_solver.cc