This is an example of how to use the TraceMinDavidsonSolMgr solver manager to solve a standard eigenvalue problem where the matrix is not explicitly available, using Tpetra data stuctures.
#include "Tpetra_Core.hpp"
#include "Tpetra_Version.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Operator.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ArrayViewDecl.hpp"
#include "Teuchos_ParameterList.hpp"
using Teuchos::rcp;
using Teuchos::RCP;
using std::cout;
using std::endl;
typedef double Scalar;
typedef Tpetra::MultiVector<Scalar> TMV;
typedef Tpetra::Operator<Scalar> TOP;
typedef Tpetra::Map<> Map;
typedef Tpetra::Import<> Import;
typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
typedef TMV::local_ordinal_type LocalOrdinal;
typedef TMV::global_ordinal_type GlobalOrdinal;
typedef TMV::node_type Node;
class MyOp : public virtual TOP {
public:
MyOp (const GlobalOrdinal n,
const Teuchos::RCP<const Teuchos::Comm<int> > comm)
{
TEUCHOS_TEST_FOR_EXCEPTION(
comm.is_null (), std::invalid_argument,
"MyOp constructor: The input Comm object must be nonnull.");
const GlobalOrdinal indexBase = 0;
opMap_ = rcp (new Map (n, indexBase, comm));
myRank_ = comm->getRank ();
numProcs_ = comm->getSize ();
LocalOrdinal nlocal = opMap_->getLocalNumElements ();
if (myRank_ > 0) {
++nlocal;
}
if (myRank_ < numProcs_ - 1) {
++nlocal;
}
std::vector<GlobalOrdinal> indices;
indices.reserve (nlocal);
if (myRank_ > 0) {
indices.push_back (opMap_->getMinGlobalIndex () - 1);
}
for (GlobalOrdinal i = opMap_->getMinGlobalIndex (); i <= opMap_->getMaxGlobalIndex (); ++i) {
indices.push_back (i);
}
if (myRank_ < numProcs_ - 1) {
indices.push_back (opMap_->getMaxGlobalIndex () + 1);
}
Teuchos::ArrayView<const GlobalOrdinal> elementList (indices);
const GlobalOrdinal numGlobalElements = n + 2*(numProcs_ - 1);
redistMap_ = rcp (new Map (numGlobalElements, elementList, indexBase, comm));
importer_= rcp (new Import (opMap_, redistMap_));
};
virtual ~MyOp() {};
RCP<const Map> getDomainMap() const { return opMap_; };
RCP<const Map> getRangeMap() const { return opMap_; };
void
apply (const TMV& X,
TMV& Y,
Teuchos::ETransp mode = Teuchos::NO_TRANS,
Scalar alpha = Teuchos::ScalarTraits<Scalar>::one (),
Scalar beta = Teuchos::ScalarTraits<Scalar>::zero ()) const
{
typedef Teuchos::ScalarTraits<Scalar> STS;
TEUCHOS_TEST_FOR_EXCEPTION(
alpha != STS::one() || beta != STS::zero(), std::logic_error,
"MyOp::apply was given alpha != 1 or beta != 0. "
"These cases are not implemented.");
const LocalOrdinal nlocRows = X.getLocalLength ();
int numVecs = static_cast<int> (X.getNumVectors ());
RCP<TMV> redistData = rcp (new TMV (redistMap_, numVecs));
redistData->doImport (X, *importer_, Tpetra::INSERT);
for (int c = 0; c < numVecs; ++c) {
Teuchos::ArrayRCP<Scalar> colView = redistData->getDataNonConst (c);
LocalOrdinal offset;
if (myRank_ > 0) {
Y.replaceLocalValue (0, c, -colView[0] + 2*colView[1] - colView[2]);
offset = 0;
}
else {
Y.replaceLocalValue (0, c, 2*colView[0] - colView[1]);
offset = 1;
}
for (LocalOrdinal r = 1; r < nlocRows - 1; ++r) {
const Scalar newVal =
-colView[r-offset] + 2*colView[r+1-offset] - colView[r+2-offset];
Y.replaceLocalValue (r, c, newVal);
}
if (myRank_ < numProcs_ - 1) {
const Scalar newVal =
-colView[nlocRows-1-offset] + 2*colView[nlocRows-offset]
- colView[nlocRows+1-offset];
Y.replaceLocalValue (nlocRows-1, c, newVal);
}
else {
const Scalar newVal =
-colView[nlocRows-1-offset] + 2*colView[nlocRows-offset];
Y.replaceLocalValue (nlocRows-1, c, newVal);
}
}
}
private:
RCP<const Map> opMap_, redistMap_;
RCP<const Import> importer_;
int myRank_, numProcs_;
};
int
main (int argc, char *argv[])
{
Tpetra::ScopeGuard tpetraScope (&argc, &argv);
RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm ();
const int myRank = comm->getRank ();
Scalar tol = 1e-5;
GlobalOrdinal n = 100;
int nev = 1;
int blockSize = 1;
bool verbose = true;
std::string whenToShift = "Always";
Teuchos::CommandLineProcessor cmdp (false, true);
cmdp.setOption ("n", &n, "Number of rows of our operator.");
cmdp.setOption ("tolerance", &tol, "Relative residual used for solver.");
cmdp.setOption ("nev", &nev, "Number of desired eigenpairs.");
cmdp.setOption ("blocksize", &blockSize, "Number of vectors to add to the "
"subspace at each iteration.");
cmdp.setOption ("verbose","quiet", &verbose,
"Whether to print a lot of info or a little bit.");
cmdp.setOption ("whenToShift", &whenToShift, "When to perform Ritz shifts. "
"Options: Never, After Trace Levels, Always.");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
RCP<MyOp> K = rcp (new MyOp (n, comm));
int verbosity;
int numRestartBlocks = 2*nev/blockSize;
int numBlocks = 10*nev/blockSize;
if (verbose) {
} else {
}
Teuchos::ParameterList MyPL;
MyPL.set ("Verbosity", verbosity);
MyPL.set ("Saddle Solver Type", "Projected Krylov");
MyPL.set ("Block Size", blockSize);
MyPL.set ("Convergence Tolerance", tol);
MyPL.set ("Num Restart Blocks", numRestartBlocks);
MyPL.set ("Num Blocks", numBlocks);
MyPL.set ("When To Shift", whenToShift);
RCP<TMV> ivec = rcp (new TMV (K->getDomainMap (), numRestartBlocks*blockSize));
TMVT::MvRandom (*ivec);
RCP<Anasazi::BasicEigenproblem<Scalar,TMV,TOP> > MyProblem =
MyProblem->setHermitian (true);
MyProblem->setNEV (nev);
const bool setProblemCorrectly = MyProblem->setProblem ();
if (! setProblemCorrectly) {
if (myRank == 0) {
cout << "Anasazi::BasicEigenproblem::setProblem() failed to set the "
"problem correctly." << endl;
}
return -1;
}
solver_type MySolverMgr (MyProblem, MyPL);
if (myRank == 0) {
cout << "Anasazi::EigensolverMgr::solve() returned ";
cout << "unconverged.";
} else {
cout << "converged.";
}
cout << endl;
}
std::vector<Anasazi::Value<Scalar> > evals = sol.
Evals;
RCP<TMV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
Teuchos::SerialDenseMatrix<int,Scalar> T (numev, numev);
TMV tempvec (K->getDomainMap (), TMVT::GetNumberVecs (*evecs));
std::vector<Scalar> normR (sol.numVecs);
TMV Kvec (K->getRangeMap (), TMVT::GetNumberVecs (*evecs));
TOPT::Apply (*K, *evecs, Kvec );
TMVT::MvTransMv (1.0, Kvec, *evecs, T);
TMVT::MvTimesMatAddMv (-1.0, *evecs, T, 1.0, Kvec);
TMVT::MvNorm (Kvec, normR);
if (myRank == 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) << "Error" << endl;
cout<<"------------------------------------------------------" << endl;
for (int i = 0; i < numev; ++i) {
cout << std::setw(16) << T(i,i) << std::setw(16) << normR[i]/T(i,i)
<< endl;
}
cout << "------------------------------------------------------" << endl;
}
}
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for creating operators that can interface with the Anasazi::OperatorTraits cl...
Partial specialization of Anasazi::MultiVecTraits and Anasazi::OperatorTraits for Tpetra objects.
The Anasazi::TraceMinDavidsonSolMgr provides a solver manager for the TraceMinDavidson eigensolver wi...
This provides a basic implementation for defining standard or generalized eigenvalue problems.
The Anasazi::TraceMinDavidsonSolMgr provides a flexible solver manager over the TraceMinDavidson eige...
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.