#include "Epetra_CrsMatrix.h"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#include "Teuchos_Assert.hpp"
#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;
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
#endif
bool success = false;
bool verbose = true;
try {
#ifdef EPETRA_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
bool boolret;
int MyPID = Comm.MyPID();
bool debug = false;
bool dynXtraNev = false;
std::string which("SM");
int nx = 10;
int nev = 4;
int blockSize = 1;
int numBlocks = 20;
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,LR,SI,or LI).");
cmdp.setOption("nx",&nx,"Number of discretization points in each direction; n=nx*nx.");
cmdp.setOption("nev",&nev,"Number of eigenvalues to compute.");
cmdp.setOption("blocksize",&blockSize,"Block Size.");
cmdp.setOption("numblocks",&numBlocks,"Number of blocks for the Krylov-Schur form.");
cmdp.setOption("dynrestart","nodynrestart",&dynXtraNev,"Use dynamic restart boundary to accelerate convergence.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return -1;
}
typedef double ScalarType;
typedef Teuchos::ScalarTraits<ScalarType> SCT;
typedef SCT::magnitudeType MagnitudeType;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
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;
}
}
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, Map, &NumNz[0]) );
double rho = 0.0;
const double one = 1.0;
std::vector<double> Values(4);
std::vector<int> Indices(4);
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, info;
for (int i=0; i<NumMyElements; i++)
{
if (MyGlobalElements[i]==0)
{
Indices[0] = 1;
Indices[1] = nx;
NumEntries = 2;
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;
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;
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;
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;
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;
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;
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;
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;
info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
info = A->InsertGlobalValues(MyGlobalElements[i], 1, &diag, &MyGlobalElements[i]);
TEUCHOS_ASSERT( info==0 );
}
info = A->FillComplete();
TEUCHOS_ASSERT( info==0 );
A->SetTracebackMode(1);
int maxRestarts = 500;
double tol = 1e-8;
Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > MySort =
if (verbose) {
}
if (debug) {
}
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Sort Manager", MySort );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Num Blocks", numBlocks );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Convergence Tolerance", tol );
MyPL.set("Dynamic Extra NEV",dynXtraNev);
Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(Map, blockSize) );
ivec->Random();
Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
MyProblem->setHermitian(rho==0.0);
MyProblem->setNEV( nev );
boolret = MyProblem->setProblem();
if (boolret != true) {
if (verbose && MyPID == 0) {
std::cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << std::endl;
}
#ifdef HAVE_MPI
MPI_Finalize() ;
#endif
return -1;
}
std::cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
std::vector<Anasazi::Value<double> > ritzValues = MySolverMgr.
getRitzValues();
if (verbose && MyPID==0) {
int numritz = (int)ritzValues.size();
std::cout.setf(std::ios_base::right, std::ios_base::adjustfield);
std::cout<<std::endl<< "Computed Ritz Values"<< std::endl;
if (MyProblem->isHermitian()) {
std::cout<< std::setw(16) << "Real Part"
<< std::endl;
std::cout<<"-----------------------------------------------------------"<<std::endl;
for (int i=0; i<numritz; i++) {
std::cout<< std::setw(16) << ritzValues[i].realpart
<< std::endl;
}
std::cout<<"-----------------------------------------------------------"<<std::endl;
}
else {
std::cout<< std::setw(16) << "Real Part"
<< std::setw(16) << "Imag Part"
<< std::endl;
std::cout<<"-----------------------------------------------------------"<<std::endl;
for (int i=0; i<numritz; i++) {
std::cout<< std::setw(16) << ritzValues[i].realpart
<< std::setw(16) << ritzValues[i].imagpart
<< std::endl;
}
std::cout<<"-----------------------------------------------------------"<<std::endl;
}
}
std::vector<Anasazi::Value<ScalarType> > 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> normA(numev);
if (MyProblem->isHermitian()) {
Epetra_MultiVector Aevecs(Map,numev);
Teuchos::SerialDenseMatrix<int,double> B(numev,numev);
B.putScalar(0.0);
for (int i=0; i<numev; i++) {B(i,i) = evals[i].realpart;}
OPT::Apply( *A, *evecs, Aevecs );
MVT::MvTimesMatAddMv( -1.0, *evecs, B, 1.0, Aevecs );
MVT::MvNorm( Aevecs, normA );
for (int i=0; i<numev; i++) {
normA[i] /= Teuchos::ScalarTraits<double>::magnitude( evals[i].realpart );
}
} else {
int i=0;
std::vector<int> curind(1);
std::vector<double> resnorm(1), tempnrm(1);
Teuchos::RCP<MV> tempAevec;
Teuchos::RCP<const MV> evecr, eveci;
Epetra_MultiVector Aevec(Map,numev);
OPT::Apply( *A, *evecs, Aevec );
Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
while (i<numev) {
if (index[i]==0) {
curind[0] = i;
evecr = MVT::CloneView( *evecs, curind );
tempAevec = MVT::CloneCopy( Aevec, curind );
Breal(0,0) = evals[i].realpart;
MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
MVT::MvNorm( *tempAevec, resnorm );
normA[i] = resnorm[0]/Teuchos::ScalarTraits<MagnitudeType>::magnitude( evals[i].realpart );
i++;
} else {
curind[0] = i;
evecr = MVT::CloneView( *evecs, curind );
tempAevec = MVT::CloneCopy( Aevec, curind );
curind[0] = i+1;
eveci = MVT::CloneView( *evecs, curind );
Breal(0,0) = evals[i].realpart;
Bimag(0,0) = evals[i].imagpart;
MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
MVT::MvTimesMatAddMv( 1.0, *eveci, Bimag, 1.0, *tempAevec );
MVT::MvNorm( *tempAevec, tempnrm );
tempAevec = MVT::CloneCopy( Aevec, curind );
MVT::MvTimesMatAddMv( -1.0, *evecr, Bimag, 1.0, *tempAevec );
MVT::MvTimesMatAddMv( -1.0, *eveci, Breal, 1.0, *tempAevec );
MVT::MvNorm( *tempAevec, resnorm );
normA[i] = lapack.LAPY2( tempnrm[0], resnorm[0] ) /
lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
normA[i+1] = normA[i];
i=i+2;
}
}
}
if (verbose && MyPID==0) {
std::cout.setf(std::ios_base::right, std::ios_base::adjustfield);
std::cout<<std::endl<< "Actual Residuals"<<std::endl;
if (MyProblem->isHermitian()) {
std::cout<< std::setw(16) << "Real Part"
<< std::setw(20) << "Direct Residual"<< std::endl;
std::cout<<"-----------------------------------------------------------"<<std::endl;
for (int i=0; i<numev; i++) {
std::cout<< std::setw(16) << evals[i].realpart
<< std::setw(20) << normA[i] << std::endl;
}
std::cout<<"-----------------------------------------------------------"<<std::endl;
}
else {
std::cout<< std::setw(16) << "Real Part"
<< std::setw(16) << "Imag Part"
<< std::setw(20) << "Direct Residual"<< std::endl;
std::cout<<"-----------------------------------------------------------"<<std::endl;
for (int i=0; i<numev; i++) {
std::cout<< std::setw(16) << evals[i].realpart
<< std::setw(16) << evals[i].imagpart
<< std::setw(20) << normA[i] << std::endl;
}
std::cout<<"-----------------------------------------------------------"<<std::endl;
}
}
}
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize();
#endif
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}
Basic implementation of the Anasazi::Eigenproblem class.
Basic implementation of the Anasazi::SortManager 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.
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
std::vector< Value< ScalarType > > getRitzValues() const
Return the Ritz values from the most recent solve.
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.