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;
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;
}