Use LOBPCG with Epetra, with shifted eigenvalue problem.
This example computes the eigenvalues of largest magnitude of the discretized 2-D Laplacian operator, using Anasazi's implementation of the LOBPCG method. This problem constructs a shifted eigenproblem that targets the smallest eigenvalues around a certain value (sigma). This operator is discretized using linear finite elements and constructed as an Epetra matrix, then passed shifted using EpetraExt utilities.
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "ModeLaplace2DQ2.h"
#include "EpetraExt_MatrixMatrix.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "Teuchos_StandardCatchMacros.hpp"
int
main (int argc, char *argv[])
{
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
bool success = false;
try {
int i;
int MyPID = Comm.MyPID();
int space_dim = 2;
std::vector<double> brick_dim( space_dim );
brick_dim[0] = 1.0;
brick_dim[1] = 1.0;
std::vector<int> elements( space_dim );
elements[0] = 10;
elements[1] = 10;
Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );
Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
double sigma = 1.0;
Teuchos::RCP<Epetra_CrsMatrix> Kshift = Teuchos::rcp( new Epetra_CrsMatrix( *K ) );
int addErr = EpetraExt::MatrixMatrix::Add( *M, false, -sigma, *Kshift, 1.0 );
if (addErr != 0) {
throw -1;
}
const int nev = 10;
const int blockSize = 5;
const int maxIters = 500;
const double tol = 1.0e-8;
std::string which = "SM";
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Maximum Iterations", maxIters );
MyPL.set( "Convergence Tolerance", tol );
MyPL.set( "Full Ortho", true );
MyPL.set( "Use Locking", true );
MyPL.set( "Locking Tolerance", tol/10 );
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->Map(), blockSize) );
MVT::MvRandom( *ivec );
Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
MyProblem->setHermitian(true);
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (!boolret) {
throw -1;
}
std::cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
std::vector<Anasazi::Value<double> > evals = sol.
Evals;
Teuchos::RCP<MV> evecs = sol.
Evecs;
if (numev > 0) {
std::vector<double> compEvals(numev);
for (i=0; i<numev; ++i) {
compEvals[i] = evals[i].realpart + sigma;
}
std::vector<double> normR(sol.
numVecs);
Epetra_MultiVector Kvec( K->OperatorDomainMap(), evecs->NumVectors() );
Epetra_MultiVector Mvec( M->OperatorDomainMap(), evecs->NumVectors() );
T.putScalar(0.0);
T(i,i) = compEvals[i];
}
K->Apply( *evecs, Kvec );
M->Apply( *evecs, Mvec );
MVT::MvTimesMatAddMv( -1.0, Mvec, T, 1.0, Kvec );
MVT::MvNorm( Kvec, normR );
std::ostringstream os;
os.setf(std::ios_base::right, std::ios_base::adjustfield);
os<<
"Solver manager returned " << (returnCode ==
Anasazi::Converged ?
"converged." :
"unconverged.") << std::endl;
os<<std::endl;
os<<"------------------------------------------------------"<<std::endl;
os<<std::setw(16)<<"Eigenvalue"
<<std::setw(18)<<"Direct Residual"
<<std::endl;
os<<"------------------------------------------------------"<<std::endl;
os<<std::setw(16)<<compEvals[i]
<<std::setw(18)<<normR[i]/compEvals[i]
<<std::endl;
}
os<<"------------------------------------------------------"<<std::endl;
}
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize();
#endif
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}
Basic implementation of the Anasazi::Eigenproblem class.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
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.
Anasazi's basic output manager for sending information of select verbosity levels to the appropriate ...
User interface for the LOBPCG eigensolver.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Traits class which defines basic operations on multivectors.
virtual Teuchos::FancyOStream & stream(MsgType type)
Create a stream for outputting to.
virtual void print(MsgType type, const std::string output)
Send output to the output manager.
ReturnType
Enumerated type used to pass back information from a solver manager.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.