oomph::ARPACK Class Reference

Class for the ARPACK eigensolver. More...

#include <eigen_solver.h>

Inheritance diagram for oomph::ARPACK:

Public Member Functions
	ARPACK ()
	Constructor. More...
	ARPACK (const ARPACK &)
	Empty copy constructor. More...
virtual	~ARPACK ()
	Destructor, delete the linear solver. More...
int &	narnoldi ()
	Access function for the number of Arnoldi vectors. More...
const int &	narnoldi () const
	Access function for the number of Arnoldi vectors (const version) More...
void	enable_compute_eigenvectors ()
	Set to enable the computation of the eigenvectors (default) More...
void	disable_compute_eigenvectors ()
	Set to disable the computation of the eigenvectors. 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	get_eigenvalues_left_of_shift ()
	Use the eigensolver to find the eigenvalues of a given matrix. More...
void	get_eigenvalues_right_of_shift ()
	Set the desired eigenvalues to be right of the shift. More...
void	track_eigenvalue_real_part ()
	Set the real part to be the quantity of interest (default) More...
void	track_eigenvalue_imaginary_part ()
	Set the imaginary part fo the quantity of interest. More...
void	track_eigenvalue_magnitude ()
	Set the magnitude to be the quantity of interest. More...
LinearSolver *&	linear_solver_pt ()
	Return a pointer to the linear solver object. More...
LinearSolver *const &	linear_solver_pt () const
	Return a pointer to the linear solver object (const version) More...
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)

Private Attributes
LinearSolver *	Linear_solver_pt
	Pointer to a linear solver. More...
LinearSolver *	Default_linear_solver_pt
	Pointer to a default linear solver. More...
int	Spectrum
int	NArnoldi
	Number of Arnoldi vectors to compute. More...
bool	Small
bool	Compute_eigenvectors
	Boolean to indicate whether or not to compute the eigenvectors. More...

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 ARPACK eigensolver.

Constructor & Destructor Documentation

ARPACK() [1/2]

oomph::ARPACK::ARPACK

(

)

Constructor.

Constructor, set default values and set the initial linear solver to be superlu

    : EigenSolver(),
      Spectrum(1),
      NArnoldi(30),
      Small(true),
      Compute_eigenvectors(true)
  {
    Default_linear_solver_pt = Linear_solver_pt = new SuperLUSolver;
  }

References Default_linear_solver_pt, and Linear_solver_pt.

ARPACK() [2/2]

oomph::ARPACK::ARPACK ( const ARPACK &  )

inline

Empty copy constructor.

{}

~ARPACK()

oomph::ARPACK::~ARPACK ( )

virtual

Destructor, delete the linear solver.

Destructor, delete the default linear solver.

  {
    delete Default_linear_solver_pt;
  }

References Default_linear_solver_pt.

Member Function Documentation

disable_compute_eigenvectors()

void oomph::ARPACK::disable_compute_eigenvectors ( )

inline

Set to disable the computation of the eigenvectors.

    {
      Compute_eigenvectors = false;
    }

References Compute_eigenvectors.

enable_compute_eigenvectors()

void oomph::ARPACK::enable_compute_eigenvectors ( )

inline

Set to enable the computation of the eigenvectors (default)

    {
      Compute_eigenvectors = true;
    }

References Compute_eigenvectors.

get_eigenvalues_left_of_shift()

void oomph::ARPACK::get_eigenvalues_left_of_shift ( )

inline

Use the eigensolver to find the eigenvalues of a given matrix.

Set the desired eigenvalues to be left of the shift

    {
      Small = true;
    }

References Small.

get_eigenvalues_right_of_shift()

void oomph::ARPACK::get_eigenvalues_right_of_shift ( )

inline

Set the desired eigenvalues to be right of the shift.

    {
      Small = false;
    }

References Small.

linear_solver_pt() [1/2]

LinearSolver*& oomph::ARPACK::linear_solver_pt ( )

inline

Return a pointer to the linear solver object.

    {
      return Linear_solver_pt;
    }

References Linear_solver_pt.

linear_solver_pt() [2/2]

LinearSolver* const& oomph::ARPACK::linear_solver_pt ( ) const

inline

Return a pointer to the linear solver object (const version)

    {
      return Linear_solver_pt;
    }

References Linear_solver_pt.

narnoldi() [1/2]

int& oomph::ARPACK::narnoldi ( )

inline

Access function for the number of Arnoldi vectors.

    {
      return NArnoldi;
    }

References NArnoldi.

narnoldi() [2/2]

const int& oomph::ARPACK::narnoldi ( ) const

inline

Access function for the number of Arnoldi vectors (const version)

    {
      return NArnoldi;
    }

References NArnoldi.

solve_eigenproblem()

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

virtual

Solve the eigen problem.

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

Implements oomph::EigenSolver.

  {
    bool Verbose = false;
 
    // Control parameters
    int ido = 0; // Reverse communication flag
    char bmat[2]; // Specifies the type of matrix on the RHS
    // Must be 'I' for standard eigenvalue problem
    //        'G' for generalised eigenvalue problem
    int n; // Dimension of the problem
    char which[3]; // Set which eigenvalues are required.
    int nev; // Number of eigenvalues desired
    double tol = 0.0; // Stopping criteria
    int ncv; // Number of columns of the matrix V (Dimension of Arnoldi
             // subspace)
    int iparam[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // Control parameters
    // Pointers to starting locations in the workspace vectors
    int ipntr[14] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    int info =
      0; // Setting to zero gives random initial guess for vector to start
    // Arnoldi iteration
 
    // Set up the sizes of the matrix
    n = problem_pt->ndof(); // Total size of matrix
    nev = n_eval; // Number of desired eigenvalues
    ncv = NArnoldi < n ? NArnoldi : n; // Number of Arnoldi vectors allowed
                                       // Maximum possible value is max
                                       // dimension of matrix
 
    // If we don't have enough Arnoldi vectors to compute the desired number
    // of eigenvalues then complain
    if (nev > ncv)
    {
      std::ostringstream warning_stream;
      warning_stream << "Number of requested eigenvalues " << nev << "\n"
                     << "is greater than the number of Arnoldi vectors:" << ncv
                     << "\n";
      // Increase the number of Arnoldi vectors
      ncv = nev + 10;
      if (ncv > n)
      {
        ncv = n;
      }
 
      warning_stream << "Increasing number of Arnoldi vectors to " << ncv
                     << "\n but you may want to increase further using\n"
                     << "ARPACK::narnoldi()\n"
                     << "which will also get rid of this warning.\n";
 
      OomphLibWarning(
        warning_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    // Allocate some workspace
    int lworkl = 3 * ncv * ncv + 6 * ncv;
 
    // Array that holds the final set of Arnoldi basis vectors
    double** v = new double*;
    *v = new double[ncv * n];
    // Leading dimension of v (n)
    int ldv = n;
 
    // Residual vector
    double* resid = new double[n];
    // Workspace vector
    double* workd = new double[3 * n];
    // Workspace vector
    double* workl = new double[lworkl];
 
    bmat[0] = 'G';
    // If we require eigenvalues to the left of the shift
    if (Small)
    {
      which[0] = 'S';
    }
    // Otherwise we require eigenvalues to the right of the shift
    else
    {
      which[0] = 'L';
    }
 
    // Which part of the eigenvalue interests us
    switch (Spectrum)
    {
        // Imaginary part
      case -1:
        which[1] = 'I';
        break;
 
        // Magnitude
      case 0:
        which[1] = 'M';
        break;
 
        // Real part
      case 1:
        which[1] = 'R';
        break;
 
      default:
        std::ostringstream error_stream;
        error_stream
          << "Spectrum is set to an invalid value. It must be 0, 1 or -1\n"
          << "Not " << Spectrum << std::endl;
 
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    //     %---------------------------------------------------%
    //     | This program uses exact shifts with respect to    |
    //     | the current Hessenberg matrix (IPARAM(1) = 1).    |
    //     | IPARAM(3) specifies the maximum number of Arnoldi |
    //     | iterations allowed.  Mode 3 of DNAUPD is used     |
    //     | (IPARAM(7) = 3). All these options can be changed |
    //     | by the user. For details see the documentation in |
    //     | DNAUPD.                                           |
    //     %---------------------------------------------------%
 
    int ishifts = 1;
    int maxitr = 300;
    int mode = 3; // M is symetric and semi-definite
 
    iparam[0] = ishifts; // Exact shifts wrt Hessenberg matrix H
    iparam[2] = maxitr; // Maximum number of allowed iterations
    iparam[3] = 1; // Bloacksize to be used in the recurrence
    iparam[6] = mode; // Mode is shift and invert in real arithmetic
 
    // Real and imaginary parts of the shift
    double sigmar = Sigma_real, sigmai = 0.0;
 
    // TEMPORARY
    // only use non distributed matrice and vectors
    LinearAlgebraDistribution dist(problem_pt->communicator_pt(), n, false);
    this->build_distribution(dist);
 
    // Before we get into the Arnoldi loop solve the shifted matrix problem
    // Allocated Row compressed matrices for the mass matrix and shifted main
    // matrix
    CRDoubleMatrix M(this->distribution_pt()), AsigmaM(this->distribution_pt());
 
    // Assemble the matrices
    problem_pt->get_eigenproblem_matrices(M, AsigmaM, sigmar);
 
    // Allocate storage for the vectors to be used in matrix vector products
    DoubleVector rhs(this->distribution_pt(), 0.0);
    DoubleVector x(this->distribution_pt(), 0.0);
 
 
    bool LOOP_FLAG = true;
    bool First = true;
 
    // Enable resolves for the linear solver
    Linear_solver_pt->enable_resolve();
 
    // Do not report the time taken
    Linear_solver_pt->disable_doc_time();
 
    do
    {
      //        %---------------------------------------------%
      //        | Repeatedly call the routine DNAUPD and take |
      //        | actions indicated by parameter IDO until    |
      //        | either convergence is indicated or maxitr   |
      //        | has been exceeded.                          |
      //        %---------------------------------------------%
      //
 
      DNAUPD(ido,
             bmat,
             n,
             which,
             nev,
             tol,
             resid,
             ncv,
             v,
             ldv,
             iparam,
             ipntr,
             workd,
             workl,
             lworkl,
             info);
 
      if (ido == -1)
      {
        //           %-------------------------------------------%
        //           | Perform matrix vector multiplication      |
        //           |                y <--- OP*x                |
        //           | The user should supply his/her own        |
        //           | matrix vector multiplication routine here |
        //           | that takes workd(ipntr(1)) as the input   |
        //           | vector, and return the matrix vector      |
        //           | product to workd(ipntr(2)).               |
        //           %-------------------------------------------%
        //
 
        // Firstly multiply by mx
        for (int i = 0; i < n; i++)
        {
          x[i] = workd[ipntr[0] - 1 + i];
        }
        M.multiply(x, rhs);
 
        // Now solve the system
        DoubleVector temp(rhs);
        if (First)
        {
          Linear_solver_pt->solve(&AsigmaM, temp, rhs);
          First = false;
        }
        else
        {
          Linear_solver_pt->resolve(temp, rhs);
        }
        temp.clear();
        // AsigmaM.lubksub(rhs);
        // Put the solution into the workspace
        for (int i = 0; i < n; i++)
        {
          workd[ipntr[1] - 1 + i] = rhs[i];
        }
 
        continue;
      }
      else if (ido == 1)
      {
        // Already done multiplication by Mx
        // Need to load the rhs vector
        for (int i = 0; i < n; i++)
        {
          rhs[i] = workd[ipntr[2] - 1 + i];
        }
        // Now solve the system
        // AsigmaM.lubksub(rhs);
        DoubleVector temp(rhs);
        if (First)
        {
          Linear_solver_pt->solve(&AsigmaM, temp, rhs);
          First = false;
        }
        else
        {
          Linear_solver_pt->resolve(temp, rhs);
        }
        // Put the solution into the workspace
        for (int i = 0; i < n; i++)
        {
          workd[ipntr[1] - 1 + i] = rhs[i];
        }
        continue;
      }
      else if (ido == 2)
      {
        // Need to multiply by mass matrix
        // Vector<double> x(n);
        for (int i = 0; i < n; i++)
        {
          x[i] = workd[ipntr[0] - 1 + i];
        }
        M.multiply(x, rhs);
 
        for (int i = 0; i < n; i++)
        {
          workd[ipntr[1] - 1 + i] = rhs[i];
        }
        continue;
      }
 
      if (info < 0)
      {
        oomph_info << "ERROR" << std::endl;
        oomph_info << "Error with _naupd, info = '" << info << std::endl;
        if (info == -7)
        {
          oomph_info << "Not enough memory " << std::endl;
        }
      }
 
      //        %-------------------------------------------%
      //        | No fatal errors occurred.                 |
      //        | Post-Process using DNEUPD.                |
      //        |                                           |
      //        | Computed eigenvalues may be extracted.    |
      //        |                                           |
      //        | Eigenvectors may also be computed now if  |
      //        | desired.  (indicated by rvec = .true.)    |
      //        %-------------------------------------------%
      //
      LOOP_FLAG = false;
    } while (LOOP_FLAG);
 
    // Report any problems
    if (info == 1)
    {
      oomph_info << "Maximum number of iterations reached." << std::endl;
    }
    else if (info == 3)
    {
      oomph_info << "No shifts could be applied during implicit Arnoldi "
                 << "update, try increasing NCV. " << std::endl;
    }
 
    // Disable resolves for the linear solver
    Linear_solver_pt->disable_resolve();
 
    // Compute Ritz or Schur vectors, if desired
    int rvec = Compute_eigenvectors;
    // Specify the form of the basis computed
    char howmany[2];
    howmany[0] = 'A';
    // Find out the number of converged eigenvalues
    int nconv = iparam[4];
 
    // Note that there is an error (feature) in ARPACK that
    // means in certain cases, if we request eigenvectors,
    // consequent Schur factorization of the matrix spanned
    // by the Arnoldi vectors leads to more converged eigenvalues
    // than reported by DNAPUD. This is a pain because it
    // causes a segmentation fault and in every case that I've found
    // the eigenvalues/vectors are not those desired.
    // At the moment, I'm just going to leave it, but I note
    // the problem here to remind myself about it.
 
    // Array to select which Ritz vectors are computed
    int select[ncv];
    // Array that holds the real part of the eigenvalues
    double* eval_real = new double[nconv + 1];
    // Array that holds the imaginary part of the eigenvalues
    double* eval_imag = new double[nconv + 1];
 
    // Workspace
    double* workev = new double[3 * ncv];
    // Array to hold the eigenvectors
    double** z = new double*;
    *z = new double[(nconv + 1) * n];
 
    // Error flag
    int ierr;
 
    DNEUPD(rvec,
           howmany,
           select,
           eval_real,
           eval_imag,
           z,
           ldv,
           sigmar,
           sigmai,
           workev,
           bmat,
           n,
           which,
           nev,
           tol,
           resid,
           ncv,
           v,
           ldv,
           iparam,
           ipntr,
           workd,
           workl,
           lworkl,
           ierr);
    //
    //        %-----------------------------------------------%
    //        | The real part of the eigenvalue is returned   |
    //        | in the first column of the two dimensional    |
    //        | array D, and the imaginary part is returned   |
    //        | in the second column of D.  The corresponding |
    //        | eigenvectors are returned in the first NEV    |
    //        | columns of the two dimensional array V if     |
    //        | requested.  Otherwise, an orthogonal basis    |
    //        | for the invariant subspace corresponding to   |
    //        | the eigenvalues in D is returned in V.        |
    //        %-----------------------------------------------%
 
    if (ierr != 0)
    {
      oomph_info << "Error with _neupd, info = " << ierr << std::endl;
    }
    else
    // Print the output anyway
    {
      // Resize the eigenvalue output vector
      eigenvalue.resize(nconv);
      for (int j = 0; j < nconv; j++)
      {
        // Turn the output from ARPACK into a complex number
        std::complex<double> eigen(eval_real[j], eval_imag[j]);
        // Add the eigenvalues to the output vector
        eigenvalue[j] = eigen;
      }
 
      if (Compute_eigenvectors)
      {
        // Load the eigenvectors
        eigenvector.resize(nconv);
        for (int j = 0; j < nconv; j++)
        {
          eigenvector[j].build(this->distribution_pt(), 0.0);
          for (int i = 0; i < n; i++)
          {
            eigenvector[j][i] = z[0][j * n + i];
          }
        }
      }
      else
      {
        eigenvector.resize(0);
      }
 
      // Report the information
      if (Verbose)
      {
        oomph_info << " ARPACK " << std::endl;
        oomph_info << " ====== " << std::endl;
        oomph_info << " Size of the matrix is " << n << std::endl;
        oomph_info << " The number of Ritz values requested is " << nev
                   << std::endl;
        oomph_info << " The number of Arnoldi vectors generated (NCV) is "
                   << ncv << std::endl;
        oomph_info << " What portion of the spectrum: " << which << std::endl;
        oomph_info << " The number of converged Ritz values is " << nconv
                   << std::endl;
        oomph_info
          << " The number of Implicit Arnoldi update iterations taken is "
          << iparam[2] << std::endl;
        oomph_info << " The number of OP*x is " << iparam[8] << std::endl;
        oomph_info << " The convergence criterion is " << tol << std::endl;
      }
    }
 
    // Clean up the allocated memory
    delete[] * z;
    *z = 0;
    delete z;
    z = 0;
    delete[] workev;
    workev = 0;
    delete[] eval_imag;
    eval_imag = 0;
    delete[] eval_real;
    eval_real = 0;
 
    // Clean up the memory
    delete[] workl;
    workl = 0;
    delete[] workd;
    workd = 0;
    delete[] resid;
    resid = 0;
    delete[] * v;
    *v = 0;
    delete v;
    v = 0;
  }

References oomph::DistributableLinearAlgebraObject::build_distribution(), oomph::DoubleVector::clear(), oomph::Problem::communicator_pt(), Compute_eigenvectors, oomph::LinearSolver::disable_doc_time(), oomph::LinearSolver::disable_resolve(), oomph::DistributableLinearAlgebraObject::distribution_pt(), DNAUPD, DNEUPD, oomph::LinearSolver::enable_resolve(), oomph::Problem::get_eigenproblem_matrices(), i, info, j, Linear_solver_pt, n, NArnoldi, oomph::Problem::ndof(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::oomph_info, oomph::LinearSolver::resolve(), oomph::EigenSolver::Sigma_real, Small, oomph::LinearSolver::solve(), Spectrum, v, and plotDoE::x.

track_eigenvalue_imaginary_part()

void oomph::ARPACK::track_eigenvalue_imaginary_part ( )

inline

Set the imaginary part fo the quantity of interest.

    {
      Spectrum = -1;
    }

References Spectrum.

track_eigenvalue_magnitude()

void oomph::ARPACK::track_eigenvalue_magnitude ( )

inline

Set the magnitude to be the quantity of interest.

    {
      Spectrum = 0;
    }

References Spectrum.

track_eigenvalue_real_part()

void oomph::ARPACK::track_eigenvalue_real_part ( )

inline

Set the real part to be the quantity of interest (default)

    {
      Spectrum = 1;
    }

References Spectrum.

Member Data Documentation

Compute_eigenvectors

bool oomph::ARPACK::Compute_eigenvectors

private

Boolean to indicate whether or not to compute the eigenvectors.

Referenced by disable_compute_eigenvectors(), enable_compute_eigenvectors(), and solve_eigenproblem().

Default_linear_solver_pt

LinearSolver* oomph::ARPACK::Default_linear_solver_pt

private

Pointer to a default linear solver.

Referenced by ARPACK(), and ~ARPACK().

Linear_solver_pt

LinearSolver* oomph::ARPACK::Linear_solver_pt

private

Pointer to a linear solver.

Referenced by ARPACK(), linear_solver_pt(), and solve_eigenproblem().

NArnoldi

int oomph::ARPACK::NArnoldi

private

Number of Arnoldi vectors to compute.

Referenced by narnoldi(), and solve_eigenproblem().

Small

bool oomph::ARPACK::Small

private

Boolean to set which part of the spectrum left (default) or right of the shifted value.

Referenced by get_eigenvalues_left_of_shift(), get_eigenvalues_right_of_shift(), and solve_eigenproblem().

Spectrum

int oomph::ARPACK::Spectrum

private

Integer to set whether the real, imaginary or magnitude is required to be small or large.

Referenced by solve_eigenproblem(), track_eigenvalue_imaginary_part(), track_eigenvalue_magnitude(), and track_eigenvalue_real_part().

---

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

• eigen_solver.h
• eigen_solver.cc