#include "Epetra_CrsMatrix.h"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "EpetraExt_readEpetraLinearSystem.h"
#include "Ifpack.h"
#include "Ifpack_Preconditioner.h"
#include "Epetra_InvOperator.h"
#ifdef EPETRA_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 std::cout;
using std::endl;
bool haveM = true;
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
Epetra_SerialComm Comm;
#endif
int MyPID = Comm.MyPID();
bool verbose=false;
bool isHermitian=false;
std::string k_filename = "bfw782a.mtx";
std::string m_filename = "bfw782b.mtx";
std::string which = "LR";
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,or LR).");
cmdp.setOption("K-filename",&k_filename,"Filename and path of the stiffness matrix.");
cmdp.setOption("M-filename",&m_filename,"Filename and path of the mass matrix.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return -1;
}
Teuchos::RCP<Epetra_Map> Map;
Teuchos::RCP<Epetra_CrsMatrix> K, M;
EpetraExt::readEpetraLinearSystem( k_filename, Comm, &K, &Map );
if (haveM) {
EpetraExt::readEpetraLinearSystem( m_filename, Comm, &M, &Map );
}
Ifpack factory;
std::string ifpack_type = "ILUT";
int overlap = 0;
Teuchos::RCP<Ifpack_Preconditioner> ifpack_prec = Teuchos::rcp(
factory.Create( ifpack_type, K.get(), overlap ) );
Teuchos::ParameterList ifpack_params;
double droptol = 1e-2;
double fill = 2.0;
ifpack_params.set("fact: drop tolerance",droptol);
ifpack_params.set("fact: ilut level-of-fill",fill);
ifpack_prec->SetParameters(ifpack_params);
ifpack_prec->Initialize();
ifpack_prec->Compute();
Teuchos::RCP<Epetra_Operator> prec = Teuchos::rcp(
new Epetra_InvOperator(ifpack_prec.get()) );
int nev = 5;
int blockSize = 5;
int maxDim = 40;
int maxRestarts = 10;
double tol = 1.0e-8;
if (verbose) {
}
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Maximum Subspace Dimension", maxDim );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Convergence Tolerance", tol );
MyPL.set( "Initial Guess", "User" );
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(*Map, blockSize) );
ivec->Random();
Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem;
if (haveM) {
MyProblem->setA(K);
MyProblem->setM(M);
MyProblem->setPrec(prec);
MyProblem->setInitVec(ivec);
}
else {
MyProblem->setA(K);
MyProblem->setPrec(prec);
MyProblem->setInitVec(ivec);
}
MyProblem->setHermitian( isHermitian );
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (boolret != true) {
if (verbose && 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;
std::vector<int> index = sol.index;
int numev = sol.numVecs;
if (numev > 0) {
Teuchos::LAPACK<int,double> lapack;
std::vector<double> normR(numev);
int i=0;
std::vector<int> curind(1);
std::vector<double> resnorm(1), tempnrm(1);
Teuchos::RCP<MV> tempKevec, Mevecs;
Teuchos::RCP<const MV> tempeveci, tempMevec;
Epetra_MultiVector Kevecs(*Map,numev);
OPT::Apply( *K, *evecs, Kevecs );
if (haveM) {
Mevecs = Teuchos::rcp( new Epetra_MultiVector(*Map,numev) );
OPT::Apply( *M, *evecs, *Mevecs );
}
else {
Mevecs = evecs;
}
Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
while (i<numev) {
if (index[i]==0) {
curind[0] = i;
tempMevec = MVT::CloneView( *Mevecs, curind );
tempKevec = MVT::CloneCopy( Kevecs, curind );
Breal(0,0) = evals[i].realpart;
MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );
MVT::MvNorm( *tempKevec, resnorm );
normR[i] = resnorm[0];
i++;
} else {
curind[0] = i;
tempMevec = MVT::CloneView( *Mevecs, curind );
tempKevec = MVT::CloneCopy( Kevecs, curind );
curind[0] = i+1;
tempeveci = MVT::CloneView( *Mevecs, curind );
Breal(0,0) = evals[i].realpart;
Bimag(0,0) = evals[i].imagpart;
MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );
MVT::MvTimesMatAddMv( 1.0, *tempeveci, Bimag, 1.0, *tempKevec );
MVT::MvNorm( *tempKevec, tempnrm );
tempKevec = MVT::CloneCopy( Kevecs, curind );
MVT::MvTimesMatAddMv( -1.0, *tempMevec, Bimag, 1.0, *tempKevec );
MVT::MvTimesMatAddMv( -1.0, *tempeveci, Breal, 1.0, *tempKevec );
MVT::MvNorm( *tempKevec, resnorm );
normR[i] = lapack.LAPY2( tempnrm[0], resnorm[0] );
normR[i+1] = normR[i];
i=i+2;
}
}
if (verbose && MyPID==0) {
cout.setf(std::ios_base::right, std::ios_base::adjustfield);
cout<<endl<< "Actual Residuals"<<endl;
cout<< std::setw(16) << "Real Part"
<< std::setw(16) << "Imag Part"
<< std::setw(20) << "Direct Residual"<< endl;
cout<<"-----------------------------------------------------------"<<endl;
for (int j=0; j<numev; j++) {
cout<< std::setw(16) << evals[j].realpart
<< std::setw(16) << evals[j].imagpart
<< std::setw(20) << normR[j] << endl;
}
cout<<"-----------------------------------------------------------"<<endl;
}
}
#ifdef EPETRA_MPI
MPI_Finalize() ;
#endif
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
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::GeneralizedDavidsonSolMgr provides a solver manager for the GeneralizedDavidson eigensol...
This provides a basic implementation for defining standard or generalized eigenvalue problems.
Solver Manager for GeneralizedDavidson.
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.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.