#include "Epetra_CrsMatrix.h"
#include "Epetra_InvOperator.h"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Ifpack.h"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
#include "ModeLaplace2DQ2.h"
int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
#endif
#ifdef HAVE_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
int nev = 10;
int blockSize = 10;
int numBlocks = 4;
int maxRestarts = 100;
double tol = 1.0e-8;
int numElements = 10;
bool verbose = true;
std::string which("SM");
bool usePrec = true;
double prec_dropTol = 1e-4;
int prec_lofill = 0;
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("nev",&nev,"Number of eigenpairs to compted.");
cmdp.setOption("blockSize",&blockSize,"Block size.");
cmdp.setOption("numBlocks",&numBlocks,"Number of blocks in basis.");
cmdp.setOption("maxRestarts",&maxRestarts,"Maximum number of restarts.");
cmdp.setOption("tol",&tol,"Relative convergence tolerance.");
cmdp.setOption("numElements",&numElements,"Number of elements in the discretization.");
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
cmdp.setOption("usePrec","noPrec",&usePrec,"Use Ifpack for preconditioning.");
cmdp.setOption("prec_dropTol",&prec_dropTol,"Preconditioner: drop tolerance.");
cmdp.setOption("prec_lofill",&prec_lofill,"Preconditioner: level of fill.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return -1;
}
if (verbose) {
}
printer.
stream(Errors) << Anasazi_Version() << std::endl << std::endl;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
printer.
stream(Errors) <<
"Generating problem matrices..." << std::flush;
const 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] = numElements;
elements[1] = numElements;
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 );
printer.
stream(Errors) <<
" done." << std::endl;
Teuchos::RCP<Ifpack_Preconditioner> prec;
Teuchos::RCP<Epetra_Operator> PrecOp;
if (usePrec) {
printer.
stream(Errors) <<
"Constructing Incomplete Cholesky preconditioner..." << std::flush;
Ifpack precFactory;
std::string precType = "IC stand-alone";
int overlapLevel = 0;
prec = Teuchos::rcp( precFactory.Create(precType,K.get(),overlapLevel) );
Teuchos::ParameterList precParams;
precParams.set("fact: drop tolerance",prec_dropTol);
precParams.set("fact: level-of-fill",prec_lofill);
IFPACK_CHK_ERR(prec->SetParameters(precParams));
IFPACK_CHK_ERR(prec->Initialize());
IFPACK_CHK_ERR(prec->Compute());
<< " done." << std::endl;
PrecOp = Teuchos::rcp( new Epetra_InvOperator(&*prec) );
}
Teuchos::RCP<Epetra_MultiVector> ivec
= Teuchos::rcp( new Epetra_MultiVector(K->OperatorDomainMap(), blockSize) );
ivec->Random();
Teuchos::RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
MyProblem->setHermitian(true);
if (usePrec) {
MyProblem->setPrec(PrecOp);
}
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (boolret != true) {
printer.
print(Errors,
"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return -1;
}
Teuchos::ParameterList MyPL;
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 );
printer.
stream(Errors) <<
"Solving eigenvalue problem..." << std::endl;
if (usePrec) {
printer.
stream(FinalSummary) << *prec << std::endl;
}
std::vector<Value<double> > evals = sol.
Evals;
Teuchos::RCP<MV> evecs = sol.
Evecs;
std::vector<double> normR(sol.
numVecs);
Epetra_MultiVector Kvec( K->OperatorDomainMap(), evecs->NumVectors() );
Epetra_MultiVector Mvec( M->OperatorDomainMap(), evecs->NumVectors() );
T.putScalar(0.0);
for (
int i=0; i<sol.
numVecs; i++) {
T(i,i) = evals[i].realpart;
}
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 == 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;
for (
int i=0; i<sol.
numVecs; i++) {
os<<std::setw(16)<<evals[i].realpart
<<std::setw(18)<<normR[i]/evals[i].realpart
<<std::endl;
}
os<<"------------------------------------------------------"<<std::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...
The Anasazi::BlockDavidsonSolMgr provides a solver manager for the BlockDavidson eigensolver.
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...
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 ...
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson 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.
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.
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.