#include "Epetra_Map.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_LinearProblem.h"
#include "AztecOO.h"
#include "AztecOO_Operator.h"
#include "Ifpack.h"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "ModeLaplace2DQ2.h"
#ifdef EPETRA_MPI
# include "Epetra_MpiComm.h"
#else
# include "Epetra_SerialComm.h"
#endif
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using Teuchos::rcp;
using std::cerr;
using std::cout;
using std::endl;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
#ifdef EPETRA_MPI
MPI_Init (&argc, &argv);
Epetra_MpiComm Comm (MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
const int MyPID = Comm.MyPID ();
const 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;
RCP<ModalProblem> testCase =
rcp (new ModeLaplace2DQ2 (Comm, brick_dim[0], elements[0],
brick_dim[1], elements[1]));
RCP<Epetra_CrsMatrix> K =
rcp (const_cast<Epetra_CrsMatrix* > (testCase->getStiffness ()), false);
RCP<Epetra_CrsMatrix> M =
rcp (const_cast<Epetra_CrsMatrix*> (testCase->getMass ()), false);
Ifpack Factory;
std::string PrecType = "ICT";
const int OverlapLevel = 0;
RCP<Ifpack_Preconditioner> Prec =
rcp (Factory.Create (PrecType, &*K, OverlapLevel));
if (Prec.is_null ()) {
throw std::logic_error ("Ifpack's factory returned a NULL preconditioner!");
}
Teuchos::ParameterList ifpackList;
ifpackList.set ("fact: drop tolerance", 1e-4);
ifpackList.set ("fact: ict level-of-fill", 0.0);
ifpackList.set ("schwarz: combine mode", "Add");
IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));
IFPACK_CHK_ERR(Prec->Initialize());
IFPACK_CHK_ERR(Prec->Compute());
Epetra_LinearProblem precProblem;
precProblem.SetOperator (K.getRawPtr ());
AztecOO precSolver (precProblem);
precSolver.SetPrecOperator (Prec.get ());
precSolver.SetAztecOption (AZ_output, AZ_none);
precSolver.SetAztecOption (AZ_solver, AZ_cg);
RCP<AztecOO_Operator> precOperator =
rcp (new AztecOO_Operator (&precSolver, K->NumGlobalRows (), 1e-12));
double tol = 1.0e-8;
int nev = 10;
int blockSize = 3;
int numBlocks = 3 * nev / blockSize;
int maxRestarts = 10;
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);
RCP<MV> ivec = rcp (new MV (K->Map (), blockSize));
MVT::MvRandom (*ivec);
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 EPETRA_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;
int numev = sol.numVecs;
if (numev > 0) {
Teuchos::SerialDenseMatrix<int,double> dmatr(numev,numev);
MV tempvec (K->Map (), MVT::GetNumberVecs (*evecs));
K->Apply (*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 (int i = 0; i < numev; ++i) {
compeval = dmatr(i,i);
cout << std::setw(16) << compeval
<< std::setw(16)
<< std::fabs (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...
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...
Adapter class for creating an operators often used in solving generalized eigenproblems.
Traits class which defines basic operations on multivectors.
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.