#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
#include "Epetra_InvOperator.h"
#include "BelosEpetraOperator.h"
#include "BelosEpetraAdapter.hpp"
#include "Ifpack.h"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "ModeLaplace2DQ2.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
int main(int argc, char *argv[]) {
using std::cout;
using std::endl;
int i;
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
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 );
Teuchos::ParameterList ifpackList;
Ifpack Factory;
std::string PrecType = "ICT";
int OverlapLevel = 0;
Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*K, OverlapLevel) );
assert(Prec != Teuchos::null);
ifpackList.set("fact: drop tolerance", 1e-4);
ifpackList.set("fact: ict level-of-fill", 0.);
ifpackList.set("schwarz: combine mode", "Add");
IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));
IFPACK_CHK_ERR(Prec->Initialize());
IFPACK_CHK_ERR(Prec->Compute());
int blockSize = 3;
int maxits = K->NumGlobalRows();
Teuchos::RCP<Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> >
My_LP = Teuchos::rcp( new Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>() );
My_LP->setOperator( K );
Teuchos::RCP<Epetra_Operator> belosPrec = Teuchos::rcp( new Epetra_InvOperator( Prec.get() ) );
My_LP->setLeftPrec( belosPrec );
Teuchos::RCP<Teuchos::ParameterList> My_List = Teuchos::rcp( new Teuchos::ParameterList() );
My_List->set( "Solver", "BlockCG" );
My_List->set( "Maximum Iterations", maxits );
My_List->set( "Block Size", 1 );
My_List->set( "Convergence Tolerance", 1e-12 );
Teuchos::RCP<Belos::EpetraOperator> BelosOp =
Teuchos::rcp( new Belos::EpetraOperator( My_LP, My_List ));
double tol = 1.0e-8;
int nev = 10;
int numBlocks = 3*nev/blockSize;
int maxRestarts = 5;
std::string which = "LM";
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 );
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::EpetraGenOp> Aop = Teuchos::rcp(
new Anasazi::EpetraGenOp(BelosOp, M,
false) );
Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
MyProblem->setHermitian(true);
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (boolret != true) {
if (MyPID == 0) {
cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << endl;
}
#ifdef HAVE_MPI
MPI_Finalize() ;
#endif
return -1;
}
cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << endl;
}
std::vector<Anasazi::Value<double> > evals = sol.
Evals;
Teuchos::RCP<MV> evecs = sol.
Evecs;
if (numev > 0) {
Teuchos::SerialDenseMatrix<int,double> dmatr(numev,numev);
Epetra_MultiVector tempvec(K->Map(), MVT::GetNumberVecs( *evecs ));
OPT::Apply( *K, *evecs, tempvec );
MVT::MvTransMv( 1.0, tempvec, *evecs, dmatr );
if (MyPID==0) {
double compeval = 0.0;
cout.setf(std::ios_base::right, std::ios_base::adjustfield);
cout<<"Actual Eigenvalues (obtained by Rayleigh quotient) : "<<endl;
cout<<"------------------------------------------------------"<<endl;
cout<<std::setw(16)<<"Real Part"
<<std::setw(16)<<"Rayleigh Error"<<endl;
cout<<"------------------------------------------------------"<<endl;
for (i=0; i<numev; i++) {
compeval = dmatr(i,i);
cout<<std::setw(16)<<compeval
<<std::setw(16)<<Teuchos::ScalarTraits<double>::magnitude(compeval-1.0/evals[i].realpart)
<<endl;
}
cout<<"------------------------------------------------------"<<endl;
}
}
#ifdef EPETRA_MPI
MPI_Finalize();
#endif
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
The Anasazi::BlockKrylovSchurSolMgr class provides a user interface for the block Krylov-Schur eigens...
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.
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Adapter class for creating an operators often used in solving generalized eigenproblems.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
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.