#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;
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.
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 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.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.