Use LOBPCG with Epetra test problem from Galeri.
Use LOBPCG with Epetra test problem from Galeri.This example computes the eigenvalues of largest magnitude of an eigenvalue problem $A x = \lambda x$, using Anasazi's implementation of the LOBPCG method, with Epetra linear algebra. It uses the Galeri package to construct the test problem.
#include "Epetra_CrsMatrix.h"
#include "Epetra_MultiVector.h"
#include "Galeri_Maps.h"
#include "Galeri_CrsMatrices.h"
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using Teuchos::rcp;
using std::cerr;
using std::cout;
using std::endl;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
#ifdef EPETRA_MPI
MPI_Init (&argc, &argv);
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
const int MyPID = Comm.MyPID ();
Teuchos::ParameterList GaleriList;
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));
const double tol = 1.0e-8;
const int nev = 10;
const int blockSize = 5;
const int maxIters = 500;
RCP<MV> ivec = rcp (new MV (*Map, blockSize));
ivec->Random ();
RCP<Anasazi::BasicEigenproblem<double, MV, OP> > problem =
problem->setHermitian (true);
problem->setNEV (nev);
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
return -1;
}
Teuchos::ParameterList anasaziPL;
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);
cout << "Anasazi eigensolver did not converge." << endl;
}
std::vector<Anasazi::Value<double> > evals = sol.
Evals;
RCP<MV> evecs = sol.Evecs;
std::vector<double> normR (sol.numVecs);
if (sol.numVecs > 0) {
Teuchos::SerialDenseMatrix<int,double> T (sol.numVecs, sol.numVecs);
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);
}
if (MyPID == 0) {
cout << endl << "Solver manager returned "
<< 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
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
The Anasazi::LOBPCGSolMgr provides a powerful solver manager for the LOBPCG eigensolver.
This provides a basic implementation for defining standard or generalized eigenvalue problems.
User interface for the LOBPCG eigensolver.
Traits class which defines basic operations on multivectors.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.