Use "Simple LOBPCG" with Epetra test problem (computed here).
Use "Simple LOBPCG" with Epetra test problem (computed here).This example computes the eigenvalues of largest magnitude of an eigenvalue problem $A x = \lambda x$, using Anasazi's "simple" implementation of the LOBPCG method, with Epetra linear algebra. It constructs the test problem within the example itself.
#include "Epetra_CrsMatrix.h"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_Assert.hpp"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using Teuchos::rcp;
using std::endl;
#ifdef HAVE_MPI
MPI_Init (&argc, &argv);
#endif
#ifdef HAVE_MPI
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
BasicOutputManager<double> printer;
printer.stream(Errors) << Anasazi_Version() << std::endl << std::endl;
std::string which ("SM");
Teuchos::CommandLineProcessor cmdp (false, true);
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
if (cmdp.parse (argc, argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
MPI_Finalize ();
#endif
return -1;
}
const int nx = 10;
const int NumGlobalElements = nx*nx;
Epetra_Map Map (NumGlobalElements, 0, Comm);
int NumMyElements = Map.NumMyElements ();
std::vector<int> MyGlobalElements (NumMyElements);
Map.MyGlobalElements (&MyGlobalElements[0]);
std::vector<int> NumNz(NumMyElements);
for (int i=0; i<NumMyElements; ++i) {
if (MyGlobalElements[i] == 0 || MyGlobalElements[i] == NumGlobalElements-1 ||
MyGlobalElements[i] == nx-1 || MyGlobalElements[i] == nx*(nx-1) ) {
NumNz[i] = 3;
}
else if (MyGlobalElements[i] < nx || MyGlobalElements[i] > nx*(nx-1) ||
MyGlobalElements[i]%nx == 0 || (MyGlobalElements[i]+1)%nx == 0) {
NumNz[i] = 4;
}
else {
NumNz[i] = 5;
}
}
RCP<Epetra_CrsMatrix> A = rcp (new Epetra_CrsMatrix (Epetra_DataAccess::Copy, Map, &NumNz[0]));
const double one = 1.0;
std::vector<double> Values(4);
std::vector<int> Indices(4);
double rho = 0.0;
double h = one /(nx+1);
double h2 = h*h;
double c = 5.0e-01*rho/ h;
Values[0] = -one/h2 - c; Values[1] = -one/h2 + c; Values[2] = -one/h2; Values[3]= -one/h2;
double diag = 4.0 / h2;
int NumEntries;
for (int i=0; i<NumMyElements; ++i) {
if (MyGlobalElements[i]==0) {
Indices[0] = 1;
Indices[1] = nx;
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i] == nx*(nx-1)) {
Indices[0] = nx*(nx-1)+1;
Indices[1] = nx*(nx-2);
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i] == nx-1) {
Indices[0] = nx-2;
NumEntries = 1;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
Indices[0] = 2*nx-1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i] == NumGlobalElements-1) {
Indices[0] = NumGlobalElements-2;
NumEntries = 1;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
Indices[0] = nx*(nx-1)-1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i] < nx) {
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]+nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i] > nx*(nx-1)) {
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]-nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if (MyGlobalElements[i]%nx == 0) {
Indices[0] = MyGlobalElements[i]+1;
Indices[1] = MyGlobalElements[i]-nx;
Indices[2] = MyGlobalElements[i]+nx;
NumEntries = 3;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else if ((MyGlobalElements[i]+1)%nx == 0) {
Indices[0] = MyGlobalElements[i]-nx;
Indices[1] = MyGlobalElements[i]+nx;
NumEntries = 2;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
Indices[0] = MyGlobalElements[i]-1;
NumEntries = 1;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
else {
Indices[0] = MyGlobalElements[i]-1;
Indices[1] = MyGlobalElements[i]+1;
Indices[2] = MyGlobalElements[i]-nx;
Indices[3] = MyGlobalElements[i]+nx;
NumEntries = 4;
int info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
int info = A->InsertGlobalValues(MyGlobalElements[i], 1, &diag, &MyGlobalElements[i]);
TEUCHOS_ASSERT( info==0 );
}
int info = A->FillComplete ();
TEUCHOS_ASSERT( info==0 );
A->SetTracebackMode (1);
RCP<Epetra_CrsMatrix> M = rcp (new Epetra_CrsMatrix (Epetra_DataAccess::Copy, Map, 1));
for (int i=0; i<NumMyElements; i++) {
Values[0] = one;
Indices[0] = i;
NumEntries = 1;
info = M->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
info = M->FillComplete ();
TEUCHOS_ASSERT( info==0 );
M->SetTracebackMode (1);
const int nev = 10;
const int blockSize = 5;
const int maxIters = 500;
const double tol = 1.0e-8;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
typedef MultiVecTraits<double, Epetra_MultiVector> MVT;
RCP<Epetra_MultiVector> ivec = rcp (new Epetra_MultiVector (Map, blockSize));
ivec->Random ();
RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
rcp (new BasicEigenproblem<double, MV, OP> (A, ivec));
MyProblem->setHermitian (true);
MyProblem->setNEV (nev);
const bool success = MyProblem->setProblem ();
if (! success) {
printer.print (Errors, "Anasazi::BasicEigenproblem::setProblem() reported an error.\n");
#ifdef HAVE_MPI
MPI_Finalize ();
#endif
return -1;
}
Teuchos::ParameterList MyPL;
MyPL.set ("Which", which);
MyPL.set ("Block Size", blockSize);
MyPL.set ("Maximum Iterations", maxIters);
MyPL.set ("Convergence Tolerance", tol);
MyPL.set ("Verbosity", verbosity);
SimpleLOBPCGSolMgr<double, MV, OP> MySolverMan (MyProblem, MyPL);
Eigensolution<double, MV> sol = MyProblem->getSolution ();
std::vector<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);
Epetra_MultiVector 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);
}
std::ostringstream os;
os.setf (std::ios_base::right, std::ios_base::adjustfield);
os << "Solver manager returned "
<< (returnCode == Converged ? "converged." : "unconverged.") << endl;
os << endl;
os << "------------------------------------------------------" << endl;
os << std::setw(16) << "Eigenvalue"
<< std::setw(18) << "Direct Residual"
<< endl;
os << "------------------------------------------------------" << endl;
for (int i = 0; i < sol.numVecs; ++i) {
os << std::setw(16) << evals[i].realpart
<< std::setw(18) << normR[i] / evals[i].realpart
<< endl;
}
os << "------------------------------------------------------" << endl;
printer.print (Errors, os.str ());
#ifdef HAVE_MPI
MPI_Finalize ();
#endif
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
Basic output manager for sending information of select verbosity levels to the appropriate output str...
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
The Anasazi::SimpleLOBPCGSolMgr provides a simple solver manager over the LOBPCG eigensolver.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.