#include "Didasko_ConfigDefs.h"
#if defined(HAVE_DIDASKO_EPETRA) && defined(HAVE_DIDASKO_ANASAZI) && defined(HAVE_DIDASKO_TEUCHOS) && defined(HAVE_DIDASKO_TRIUTILS)
#ifdef HAVE_MPI
#include "mpi.h"
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_MultiVector.h"
#include "Trilinos_Util_CrsMatrixGallery.h"
#include "AnasaziBasicEigenproblem.hpp"
#include "AnasaziEpetraAdapter.hpp"
#include "AnasaziBlockKrylovSchurSolMgr.hpp"
int main(int argc, char *argv[]) {
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
typedef Anasazi::MultiVecTraits<double, MV> MVT;
typedef Anasazi::OperatorTraits<double, MV, OP> OPT;
bool ierr;
int MyPID = Comm.MyPID();
bool verbose = (MyPID==0);
Trilinos_Util::CrsMatrixGallery Gallery("recirc_2d", Comm);
Gallery.Set("problem_size", 100);
if (verbose) {
cout << "Anasazi Example: Block Kyrlov Schur" << endl;
cout << "Problem info:" << endl;
cout << Gallery;
cout << endl;
}
int nev = 4;
int blocksize = 4;
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( Gallery.GetMatrix(), false );
const Epetra_Map * Map = &(A->RowMap());
Teuchos::RCP<Epetra_MultiVector> ivec =
Teuchos::rcp( new Epetra_MultiVector(*Map,blocksize) );
ivec->Random();
Teuchos::RCP< Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(A, ivec) );
MyProblem->setHermitian(false);
MyProblem->setNEV( nev );
ierr = MyProblem->setProblem();
assert(ierr);
int verbosity = Anasazi::Warnings + Anasazi::Errors + Anasazi::FinalSummary;
std::string which("SM");
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blocksize );
MyPL.set( "Num Blocks", 5 );
MyPL.set( "Maximum Restarts", 300 );
MyPL.set( "Convergence Tolerance", 1.0e-8 );
Anasazi::BlockKrylovSchurSolMgr<double,MV,OP>
MyBlockKrylovSchur(MyProblem, MyPL );
Anasazi::ReturnType solverreturn = MyBlockKrylovSchur.solve();
switch (solverreturn) {
case Anasazi::Unconverged:
if (verbose)
cout << "Anasazi::BlockKrylovSchur::solve() did not converge!" << endl;
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return 0;
break;
case Anasazi::Converged:
if (verbose)
cout << "Anasazi::BlockKrylovSchur::solve() converged!" << endl;
break;
}
Anasazi::Eigensolution<double, Epetra_MultiVector> sol = MyProblem->getSolution();
int numev = sol.numVecs;
std::vector<int> index = sol.index;
Teuchos::RCP<Epetra_MultiVector> evecs = sol.Evecs;
std::vector<Anasazi::Value<double> > evals = sol.Evals;
Teuchos::LAPACK<int,double> lapack;
std::vector<double> normA(numev);
int i=0;
std::vector<int> curind(1);
std::vector<double> resnorm(1), tempnrm(1);
Teuchos::RCP<MV> evecr, eveci, tempAevec;
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<double>::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[i], resnorm[i] ) /
lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
normA[i+1] = normA[i];
i=i+2;
}
}
if(verbose) {
cout << scientific << showpos << setprecision(6) << showpoint;
cout << "******************************************************\n"
<< " Results (outside of eigensolver) \n"
<< "------------------------------------------------------\n"
<< " real(evals)\t\timag(evals)\tnormR\n"
<< "------------------------------------------------------\n";
for(int i=0; i<numev; ++i) {
cout << " " << evals[i].realpart << "\t\t" << evals[i].imagpart
<< "\t" << normA[i] << endl;
}
cout << "******************************************************\n" << endl;
}
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return(EXIT_SUCCESS);
}
#else
#include <stdlib.h>
#include <stdio.h>
#ifdef HAVE_MPI
#include "mpi.h"
#endif
int main(int argc, char *argv[])
{
#ifdef HAVE_MPI
MPI_Init(&argc, &argv);
#endif
puts("Please configure Didasko with:");
puts("--enable-epetra");
puts("--enable-teuchos");
puts("--enable-triutils");
puts("--enable-anasazi");
#ifdef HAVE_MPI
MPI_Finalize();
#endif
return 0;
}
#endif