Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
BlockKrylovSchur/BlockKrylovSchurEpetraExGenAmesos.cpp

Compute smallest eigenvalues of a generalized eigenvalue problem, using block Krylov-Schur with Epetra and an Amesos direct solver.This example computes the eigenvalues of smallest magnitude of a generalized eigenvalue problem $K x = \lambda M x$, using Anasazi's implementation of the block Krylov-Schur method, with Epetra linear algebra and a direct solver from the Amesos package.

Anasazi computes the smallest-magnitude eigenvalues using a shift-and-invert strategy. For simplicity, this example uses a shift of zero. It illustrates the general pattern for using Anasazi for this problem:

  1. Construct an "operator" A such that $Az = K^{-1} M z$.
  2. Use Anasazi to solve $Az = \sigma z$, which is a spectral transformation of the original problem $K x = \lambda M x$.
  3. The eigenvalues $\lambda$ of the original problem are the inverses of the eigenvalues $\sigma$ of the transformed problem.

In the example, the "operator A such that \f$A z = K^{-1} M z\f$" is a subclass of Epetra_Operator. The Apply method of that operator takes the vector b, and computes $x = K^{-1} M b$. It does so first by applying the matrix M, and then by solving the linear system $K x = M b$ for x. Trilinos implements many different ways to solve linear systems. The example uses the sparse direct solver KLU to do so. Trilinos' Amesos package has an interface to KLU.

// @HEADER
// *****************************************************************************
// Anasazi: Block Eigensolvers Package
//
// Copyright 2004 NTESS and the Anasazi contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
// This example computes the eigenvalues of smallest magnitude of a
// generalized eigenvalue problem $K x = \lambda M x$, using Anasazi's
// implementation of the block Krylov-Schur method. It implements
// inverse iteration using the KLU sparse direct linear solver, and
// accesses KLU through Trilinos' Amesos package.
//
// Anasazi computes the smallest-magnitude eigenvalues using a
// shift-and-invert strategy. For simplicity, the code below uses a
// shift of zero. It illustrates the general pattern for using
// Anasazi for this problem:
//
// 1. Construct an "operator" A such that $Az = K^{-1} M z$.
// 2. Use Anasazi to solve $Az = \sigma z$, which is a spectral
// transformation of the original problem $K x = \lambda M x$.
// 3. The eigenvalues $\lambda$ of the original problem are the
// inverses of the eigenvalues $\sigma$ of the transformed
// problem.
//
// The "operator A such that $A z = K^{-1} M z$" is a subclass of
// Epetra_Operator. The Apply method of that operator takes the
// vector b, and computes $x = K^{-1} M b$. It does so first by
// applying the matrix M, and then by solving the linear system $K x =
// M b$ for x. Trilinos implements many different ways to solve
// linear systems. The example below uses the sparse direct solver
// KLU to do so. Trilinos' Amesos package has an interface to KLU.
// Include header for block Krylov-Schur solver
// Include header to define basic eigenproblem Ax = \lambda*Bx
// Include header to provide Anasazi with Epetra adapters. If you
// plan to use Tpetra objects instead of Epetra objects, include
// AnasaziTpetraAdapter.hpp instead; do analogously if you plan to use
// Thyra objects instead of Epetra objects.
// Include header for Epetra sparse matrix, Map (representation of
// parallel distributions), and linear problem. Amesos uses the
// latter to encapsulate linear systems to solve.
#include "Epetra_Map.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
// Include header for Amesos solver.
#include "Amesos.h"
// Include header for Teuchos serial dense matrix, which we use to
// compute the eigenvectors.
// Include header for the problem definition
#include "ModeLaplace2DQ2.h"
// Include selected communicator class required by Epetra objects
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif
// \class AmesosGenOp
// \brief Operator that computes \f$Y = K^{-1} M X\f$ or \f$Y = M
// K^{-T} X\f$, where solves with the Epetra_CrsMatrix K use an
// Amesos solver, and M is an Epetra_Operator.
// \author Heidi K. Thornquist
// \date Jan 9, 2006
//
// This operator calls an inverse method in Apply(). It may use a
// Belos solve, Amesos solve, or AztecOO solve. This is used to
// implement shift and invert, with a shift to zero. Shift before you
// invert (i.e., solve). This operator goes into Block Krylov-Schur.
//
// Anasazi will ask for the largest eigenvalues of the inverse, thus,
// the smallest eigenvalues of the original (shifted) matrix. If you
// wanted a shift, just apply the shift first, and send that to this
// operator.
class AmesosGenOp : public virtual Epetra_Operator {
public:
// \brief Constructor
//
// \param solver [in] The Amesos solver to use for solving linear
// systems with K. It must already have its Epetra_LinearProblem
// set, and that Epetra_LinearProblem must have K.
// \param massMtx [in] The "mass matrix" M.
// \param useTranspose [in] If false, apply \f$K^{-1} M\f$; else,
// apply \f$M K^{-T}\f$.
AmesosGenOp (const Teuchos::RCP<Amesos_BaseSolver>& solver,
const bool useTranspose = false );
// Destructor
virtual ~AmesosGenOp () {}
// \brief Compute Y such that \f$K^{-1} M X = Y\f$ or \f$M K^{-T} X
// = Y\f$, where solves with the Epetra_CrsMatrix K use an Amesos
// solver.
int Apply (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
// The Operator's (human-readable) label.
const char* Label() const {
return "Operator that applies K^{-1} M or M K^{-T}";
}
// Whether to apply \f$K^{-1} M\f$ (false) or \f$K^{-T} M\f$ (true).
bool UseTranspose() const { return useTranspose_; };
// \brief Set whether to apply \f$K^{-1} M\f$ (false) or
// \f$K^{-T} M\f$ (true).
//
// \param useTranspose [in] If true, tell this Operator to apply
// \f$K^{-T} M\f$, from now until the next call to
// SetUseTranspose(). If false, tell this Operator to apply
// \f$K^{-1} M\f$, until the next call to SetUseTranspose().
int SetUseTranspose (bool useTranspose);
// The Operator's communicator.
const Epetra_Comm& Comm () const {
return solver_->Comm ();
}
// The Operator's domain Map.
const Epetra_Map& OperatorDomainMap () const {
return massMtx_->OperatorDomainMap ();
}
// The Operator's range Map.
const Epetra_Map& OperatorRangeMap () const {
return massMtx_->OperatorRangeMap ();
}
// \brief NOT IMPLEMENTED: Apply \f$( K^{-1} M )^{-1}\f$ or
// \f$( K^{-1} M )^{-T}\f$.
//
// Returning a nonzero value means that this Operator doesn't know
// how to apply its inverse.
return -1;
};
// NOT IMPLEMENTED: Whether this Operator can compute its infinity norm.
bool HasNormInf() const { return false; }
// \brief NOT IMPLEMENTED: Infinity norm of \f$( K^{-1} M )^{-1}\f$
// or \f$( K^{-1} M )^{-T}\f$.
//
// Returning -1.0 means that the Operator does not know how to
// compute its infinity norm.
double NormInf () const { return -1.0; }
private:
// Default constructor: You may not call this.
AmesosGenOp () {};
// Copy constructor: You may not call this.
AmesosGenOp (const AmesosGenOp& genOp);
bool useTranspose_;
};
// ****************************************************************************
// BEGIN MAIN ROUTINE
// ****************************************************************************
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using Teuchos::rcp;
using std::cerr;
using std::cout;
using std::endl;
// Anasazi solvers have the following template parameters:
//
// - Scalar: The type of dot product results.
// - MV: The type of (multi)vectors.
// - OP: The type of operators (functions from multivector to
// multivector). A matrix (like Epetra_CrsMatrix) is an example
// of an operator; an Ifpack preconditioner is another example.
//
// Here, Scalar is double, MV is Epetra_MultiVector, and OP is
// Epetra_Operator.
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
#ifdef EPETRA_MPI
// Initialize MPI
MPI_Init (&argc, &argv);
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
#endif // EPETRA_MPI
const int MyPID = Comm.MyPID ();
//
// Set up the test problem
//
// Dimensionality of the spatial domain to discretize
const int space_dim = 2;
// Size of each of the dimensions of the (discrete) domain
std::vector<double> brick_dim (space_dim);
brick_dim[0] = 1.0;
brick_dim[1] = 1.0;
// Number of elements in each of the dimensions of the domain
std::vector<int> elements (space_dim);
elements[0] = 10;
elements[1] = 10;
// Create the test problem.
RCP<ModalProblem> testCase =
rcp (new ModeLaplace2DQ2 (Comm, brick_dim[0], elements[0],
brick_dim[1], elements[1]));
// Get the stiffness and mass matrices.
//
// rcp (T*, false) returns a nonowning (doesn't deallocate) RCP.
RCP<Epetra_CrsMatrix> K =
rcp (const_cast<Epetra_CrsMatrix* > (testCase->getStiffness ()), false);
RCP<Epetra_CrsMatrix> M =
rcp (const_cast<Epetra_CrsMatrix* > (testCase->getMass ()), false);
//
// Create linear solver for linear systems with K
//
// Anasazi uses shift and invert, with a "shift" of zero, to find
// the eigenvalues of least magnitude. In this example, we
// implement the "invert" part of shift and invert by using an
// Amesos direct solver.
//
// Create Epetra linear problem class for solving linear systems
// with K. This implements the inverse operator for shift and
// invert.
Epetra_LinearProblem AmesosProblem;
// Tell the linear problem about the matrix K. Epetra_LinearProblem
// doesn't know about RCP, so we have to give it a raw pointer.
AmesosProblem.SetOperator (K.get ());
// Create Amesos factory and solver for solving linear systems with
// K. The solver uses the KLU library to do a sparse LU
// factorization.
//
// Note that the AmesosProblem object "absorbs" K. Anasazi doesn't
// see K, just the operator that implements K^{-1} M.
Amesos amesosFactory;
RCP<Amesos_BaseSolver> AmesosSolver;
// Amesos can interface to many different solvers. The following
// loop picks a solver that Amesos supports. The loop order
// reflects solver preference, only in the sense that using LAPACK
// here is a suboptimal fall-back. (With the LAPACK option, Amesos
// makes a dense version of the sparse matrix and uses LAPACK to
// compute the factorization. The other options are true sparse
// direct factorizations.)
const int numSolverNames = 9;
const char* solverNames[9] = {
"Klu", "Umfpack", "Superlu", "Superludist", "Mumps",
"Paradiso", "Taucs", "CSparse", "Lapack"
};
for (int k = 0; k < numSolverNames; ++k) {
const char* const solverName = solverNames[k];
if (amesosFactory.Query (solverName)) {
AmesosSolver = rcp (amesosFactory.Create (solverName, AmesosProblem));
if (MyPID == 0) {
cout << "Amesos solver: \"" << solverName << "\"" << endl;
}
}
}
if (AmesosSolver.is_null ()) {
throw std::runtime_error ("Amesos appears not to have any solvers enabled.");
}
// The AmesosGenOp class assumes that the symbolic and numeric
// factorizations have already been performed on the linear problem.
AmesosSolver->SymbolicFactorization ();
AmesosSolver->NumericFactorization ();
//
// Set parameters for the block Krylov-Schur eigensolver
//
double tol = 1.0e-8; // convergence tolerance
int nev = 10; // number of eigenvalues for which to solve
int blockSize = 3; // block size (number of eigenvectors processed at once)
int numBlocks = 3 * nev / blockSize; // restart length
int maxRestarts = 5; // maximum number of restart cycles
// We're looking for the largest-magnitude eigenvalues of the
// _inverse_ operator, thus, the smallest-magnitude eigenvalues of
// the original operator.
std::string which = "LM";
// Create ParameterList to pass into eigensolver
MyPL.set ("Verbosity", verbosity);
MyPL.set ("Which", which);
MyPL.set ("Block Size", blockSize);
MyPL.set ("Num Blocks", numBlocks);
MyPL.set ("Maximum Restarts", maxRestarts);
MyPL.set ("Convergence Tolerance", tol);
// Create an initial set of vectors to start the eigensolver. Note:
// This needs to have the same number of columns as the block size.
RCP<MV> ivec = rcp (new MV (K->Map (), blockSize));
ivec->Random ();
// Create the Epetra_Operator for the spectral transformation using
// the Amesos direct solver.
//
// The AmesosGenOp object is the operator we give to Anasazi. Thus,
// Anasazi just sees an operator that computes y = K^{-1} M x. The
// matrix K got absorbed into AmesosProblem (the
// Epetra_LinearProblem object). Later, when we set up the Anasazi
// eigensolver, we will need to tell it about M, so that it can
// orthogonalize basis vectors with respect to the inner product
// defined by M (since it is symmetric positive definite).
RCP<AmesosGenOp> Aop = rcp (new AmesosGenOp (AmesosSolver, M));
// Create the eigenproblem. This object holds all the stuff about
// your problem that Anasazi will see.
//
// Anasazi only needs M so that it can orthogonalize basis vectors
// with respect to the M inner product. Wouldn't it be nice if
// Anasazi didn't require M in two different places? Alas, this is
// not currently the case.
RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
// Tell the eigenproblem that the matrix pencil (K,M) is symmetric.
MyProblem->setHermitian (true);
// Set the number of eigenvalues requested
MyProblem->setNEV (nev);
// Tell the eigenproblem that you are finished passing it information.
const bool boolret = MyProblem->setProblem ();
if (boolret != true) {
if (MyPID == 0) {
cerr << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return -1;
}
// Create the Block Krylov-Schur eigensolver.
Anasazi::BlockKrylovSchurSolMgr<double, MV, OP> MySolverMgr (MyProblem, MyPL);
// Solve the eigenvalue problem.
//
// Note that creating the eigensolver is separate from solving it.
// After creating the eigensolver, you may call solve() multiple
// times with different parameters or initial vectors. This lets
// you reuse intermediate state, like allocated basis vectors.
Anasazi::ReturnType returnCode = MySolverMgr.solve ();
if (returnCode != Anasazi::Converged && MyPID == 0) {
cout << "Anasazi eigensolver did not converge." << endl;
}
// Get the eigenvalues and eigenvectors from the eigenproblem.
Anasazi::Eigensolution<double,MV> sol = MyProblem->getSolution ();
// Anasazi returns eigenvalues as Anasazi::Value, so that if
// Anasazi's Scalar type is real-valued (as it is in this case), but
// some eigenvalues are complex, you can still access the
// eigenvalues correctly. In this case, there are no complex
// eigenvalues, since the matrix pencil is symmetric.
std::vector<Anasazi::Value<double> > evals = sol.Evals;
RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
// Reconstruct the eigenvalues. The ones that Anasazi gave back
// are the inverses of the original eigenvalues. Reconstruct the
// eigenvectors too.
MV tempvec (K->Map (), MVT::GetNumberVecs (*evecs));
K->Apply (*evecs, tempvec);
MVT::MvTransMv (1.0, tempvec, *evecs, dmatr);
if (MyPID == 0) {
double compeval = 0.0;
cout.setf (std::ios_base::right, std::ios_base::adjustfield);
cout << "Actual Eigenvalues (obtained by Rayleigh quotient) : " << endl;
cout << "------------------------------------------------------" << endl;
cout << std::setw(16) << "Real Part"
<< std::setw(16) << "Rayleigh Error" << endl;
cout << "------------------------------------------------------" << endl;
for (int i = 0; i < numev; ++i) {
compeval = dmatr(i,i);
cout << std::setw(16) << compeval
<< std::setw(16)
<< std::fabs (compeval - 1.0/evals[i].realpart)
<< endl;
}
cout << "------------------------------------------------------" << endl;
}
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return 0;
}
AmesosGenOp::
AmesosGenOp (const Teuchos::RCP<Amesos_BaseSolver>& solver,
const bool useTranspose)
: solver_ (solver),
massMtx_ (massMtx),
problem_ (NULL),
useTranspose_ (useTranspose)
{
if (solver.is_null ()) {
throw std::invalid_argument ("AmesosGenOp constructor: The 'solver' "
"input argument is null.");
}
if (massMtx.is_null ()) {
throw std::invalid_argument ("AmesosGenOp constructor: The 'massMtx' "
"input argument is null.");
}
Epetra_LinearProblem* problem = const_cast<Epetra_LinearProblem*> (solver->GetProblem ());
if (problem == NULL) {
throw std::invalid_argument ("The solver's GetProblem() method returned "
"NULL. This probably means that its "
"LinearProblem has not yet been set.");
}
problem_ = problem;
if (solver_->UseTranspose ()) {
solver_->SetUseTranspose (! useTranspose);
} else {
solver_->SetUseTranspose (useTranspose);
}
if (massMtx_->UseTranspose ()) {
massMtx_->SetUseTranspose (! useTranspose);
} else {
massMtx_->SetUseTranspose (useTranspose);
}
}
int
AmesosGenOp::SetUseTranspose (bool useTranspose)
{
int err = 0;
if (problem_ == NULL) {
throw std::logic_error ("AmesosGenOp::SetUseTranspose: problem_ is NULL");
}
if (massMtx_.is_null ()) {
throw std::logic_error ("AmesosGenOp::SetUseTranspose: massMtx_ is null");
}
if (solver_.is_null ()) {
throw std::logic_error ("AmesosGenOp::SetUseTranspose: solver_ is null");
}
const bool solverUsesTranspose = solver_->UseTranspose ();
if (solverUsesTranspose) {
err = solver_->SetUseTranspose (! useTranspose);
} else {
err = solver_->SetUseTranspose (useTranspose);
}
// If SetUseTranspose returned zero above, then the Amesos solver
// doesn't know how to change the transpose state.
if (err != 0) {
return err;
}
if (massMtx_->UseTranspose ()) {
err = massMtx_->SetUseTranspose (! useTranspose);
} else {
err = massMtx_->SetUseTranspose (useTranspose);
}
// If SetUseTranspose returned zero above, then the mass matrix
// doesn't know how to change the transpose state.
if (err != 0) {
// Put the solver back like we found it.
(void) solver_->SetUseTranspose (solverUsesTranspose);
return err;
}
useTranspose_ = useTranspose;
return 0; // the method completed correctly
}
int
AmesosGenOp::Apply (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
if (problem_ == NULL) {
throw std::logic_error ("AmesosGenOp::Apply: problem_ is NULL");
}
if (massMtx_.is_null ()) {
throw std::logic_error ("AmesosGenOp::Apply: massMtx_ is null");
}
if (solver_.is_null ()) {
throw std::logic_error ("AmesosGenOp::Apply: solver_ is null");
}
if (! useTranspose_) {
// Storage for M*X
Epetra_MultiVector MX (X.Map (), X.NumVectors ());
// Apply M*X
massMtx_->Apply (X, MX);
Y.PutScalar (0.0);
// Set the LHS and RHS
problem_->SetRHS (&MX);
problem_->SetLHS (&Y);
// Solve the linear system A*Y = MX
solver_->Solve ();
}
else { // apply the transposed operator
// Storage for A^{-T}*X
Epetra_MultiVector ATX (X.Map (), X.NumVectors ());
Epetra_MultiVector tmpX = const_cast<Epetra_MultiVector&> (X);
// Set the LHS and RHS
problem_->SetRHS (&tmpX);
problem_->SetLHS (&ATX);
// Solve the linear system A^T*Y = X
solver_->Solve ();
// Apply M*ATX
massMtx_->Apply (ATX, Y);
}
return 0; // the method completed correctly
}