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

Use LOBPCG with Epetra test problem from Galeri.This example computes the eigenvalues of largest magnitude of an eigenvalue problem $A x = x$, using Anasazi's implementation of the LOBPCG method, with Epetra linear algebra. It uses the Galeri package to construct the test problem.

// @HEADER
// *****************************************************************************
// Anasazi: Block Eigensolvers Package
//
// Copyright 2004 NTESS and the Anasazi contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
// Include header for LOBPCG eigensolver
// Include header to define eigenproblem Ax = \lambda*x
// 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 and multivector.
#include "Epetra_CrsMatrix.h"
#include "Epetra_MultiVector.h"
// The Trilinos package Galeri has many example problems.
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
// Include selected communicator class required by Epetra objects
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif // EPETRA_MPI
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;
// Anasazi's interface to the linear algebra implementation's
// vectors and multivectors is MultiVecTraits. This is a C++ traits
// interface. It is templated on Scalar and MV (see above).
#ifdef EPETRA_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.
//
// We use Trilinos' Galeri package to construct a test problem.
// Here, we use a discretization of the 2-D Laplacian operator.
// The global mesh size is nx * nx.
//
const int nx = 30;
GaleriList.set ("n", nx * nx);
GaleriList.set ("nx", nx);
GaleriList.set ("ny", nx);
RCP<Epetra_Map> Map = rcp (Galeri::CreateMap ("Linear", Comm, GaleriList));
RCP<Epetra_RowMatrix> A =
rcp (Galeri::CreateCrsMatrix ("Laplace2D", &*Map, GaleriList));
// Set eigensolver parameters.
const double tol = 1.0e-8; // convergence tolerance
const int nev = 10; // number of eigenvalues for which to solve
const int blockSize = 5; // block size (number of eigenvectors processed at once)
const int maxIters = 500; // maximum number of iterations
// Create a set of initial vectors to start the eigensolver.
// This needs to have the same number of columns as the block size.
RCP<MV> ivec = rcp (new MV (*Map, blockSize));
ivec->Random ();
// Create the eigenproblem. This object holds all the stuff about
// your problem that Anasazi will see. In this case, it knows about
// the matrix A and the inital vectors.
RCP<Anasazi::BasicEigenproblem<double, MV, OP> > problem =
// Tell the eigenproblem that the operator A is symmetric.
problem->setHermitian (true);
// Set the number of eigenvalues requested
problem->setNEV (nev);
// Tell the eigenproblem that you are finishing passing it information.
const bool boolret = problem->setProblem();
if (boolret != true) {
if (MyPID == 0) {
cerr << "Anasazi::BasicEigenproblem::setProblem() returned an error." << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize ();
#endif // EPETRA_MPI
return -1;
}
// Create a ParameterList, to pass parameters into the LOBPCG
// eigensolver.
anasaziPL.set ("Which", "SM");
anasaziPL.set ("Block Size", blockSize);
anasaziPL.set ("Maximum Iterations", maxIters);
anasaziPL.set ("Convergence Tolerance", tol);
anasaziPL.set ("Full Ortho", true);
anasaziPL.set ("Use Locking", true);
anasaziPL.set ("Verbosity", Anasazi::Errors + Anasazi::Warnings +
// Create the LOBPCG eigensolver.
Anasazi::LOBPCGSolMgr<double, MV, OP> anasaziSolver (problem, anasaziPL);
// 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 = anasaziSolver.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 = problem->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;
// Compute residuals.
std::vector<double> normR (sol.numVecs);
if (sol.numVecs > 0) {
MV tempAevec (*Map, sol.numVecs);
T.putScalar (0.0);
for (int i=0; i<sol.numVecs; ++i) {
T(i,i) = evals[i].realpart;
}
A->Apply (*evecs, tempAevec);
MVT::MvTimesMatAddMv (-1.0, *evecs, T, 1.0, tempAevec);
MVT::MvNorm (tempAevec, normR);
}
// Print the results on MPI process 0.
if (MyPID == 0) {
cout << endl << "Solver manager returned "
<< (returnCode == Anasazi::Converged ? "converged." : "unconverged.")
<< endl << endl
<< "------------------------------------------------------" << endl
<< std::setw(16) << "Eigenvalue"
<< std::setw(18) << "Direct Residual"
<< endl
<< "------------------------------------------------------" << endl;
for (int i=0; i<sol.numVecs; ++i) {
cout << std::setw(16) << evals[i].realpart
<< std::setw(18) << normR[i] / evals[i].realpart
<< endl;
}
cout << "------------------------------------------------------" << endl;
}
#ifdef EPETRA_MPI
MPI_Finalize () ;
#endif // EPETRA_MPI
return 0;
}