This is an example of how to use the TraceMinDavidsonSolMgr solver manager to compute the largest eigenpairs of a matrix via an invisible spectral transformation, using Tpetra data stuctures.
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Core.hpp"
#include "Tpetra_Version.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_MultiVector.hpp"
#include "Tpetra_Operator.hpp"
#include "Tpetra_Vector.hpp"
#include <MatrixMarket_Tpetra.hpp>
#include <TpetraExt_MatrixMatrix_def.hpp>
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ParameterList.hpp"
int
main (int argc, char *argv[])
{
using Teuchos::RCP;
using std::cout;
using std::cin;
typedef Tpetra::CrsMatrix<>::scalar_type Scalar;
typedef Tpetra::CrsMatrix<Scalar> CrsMatrix;
typedef Tpetra::MultiVector<Scalar> MV;
typedef Tpetra::Operator<Scalar> OP;
Tpetra::ScopeGuard tpetraScope(&argc,&argv);
RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
const int myRank = comm->getRank();
std::string filenameA("/home/amklinv/matrices/bcsstk06.mtx");
Scalar tol = 1e-6;
int nev = 4;
int blockSize = 1;
bool verbose = true;
std::string whenToShift = "Always";
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("fileA",&filenameA, "Filename for the Matrix-Market stiffness matrix.");
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<const CrsMatrix> K = Tpetra::MatrixMarket::Reader<CrsMatrix>::readSparseFile(filenameA, comm);
Scalar mat_norm = K->getFrobeniusNorm();
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*mat_norm );
MyPL.set( "Relative Convergence Tolerance", false);
MyPL.set( "Use Locking", true);
MyPL.set( "Relative Locking Tolerance", false);
MyPL.set("Num Restart Blocks", numRestartBlocks);
MyPL.set("Num Blocks", numBlocks);
MyPL.set("When To Shift", whenToShift);
MyPL.set("Which", "LM");
RCP<MV> ivec = Teuchos::rcp( new MV(K->getRowMap(), numRestartBlocks*blockSize) );
MVT::MvRandom( *ivec );
RCP<Anasazi::BasicEigenproblem<Scalar,MV,OP> > MyProblem =
MyProblem->setHermitian(true);
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (boolret != true) {
if (myRank == 0) {
cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << std::endl;
}
return -1;
}
cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
else if (myRank == 0)
cout << "Anasazi::EigensolverMgr::solve() returned converged." << std::endl;
std::vector<Anasazi::Value<Scalar> > evals = sol.
Evals;
RCP<MV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
Teuchos::SerialDenseMatrix<int,Scalar> T(numev,numev);
for(int i=0; i < numev; i++)
T(i,i) = evals[i].realpart;
std::vector<Scalar> normR(sol.numVecs);
MV Kvec( K->getRowMap(), MVT::GetNumberVecs( *evecs ) );
OPT::Apply( *K, *evecs, Kvec );
MVT::MvTimesMatAddMv( -1.0, *evecs, T, 1.0, Kvec );
MVT::MvNorm( Kvec, normR );
if (myRank == 0) {
cout.setf(std::ios_base::right, std::ios_base::adjustfield);
cout<<"Actual Eigenvalues (obtained by Rayleigh quotient) : "<<std::endl;
cout<<"------------------------------------------------------"<<std::endl;
cout<<std::setw(16)<<"Real Part"
<<std::setw(16)<<"Error"<<std::endl;
cout<<"------------------------------------------------------"<<std::endl;
for (int i=0; i<numev; i++) {
cout<<std::setw(16)<<evals[i].realpart
<<std::setw(16)<<normR[i]/mat_norm
<<std::endl;
}
cout<<"------------------------------------------------------"<<std::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.