This is an example of how to use the TraceMinDavidsonSolMgr solver manager to solve a generalized eigenvalue problem, 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 "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_ArrayViewDecl.hpp"
#include "Teuchos_ParameterList.hpp"
int main(int argc, char *argv[]) {
using Teuchos::RCP;
using Teuchos::rcp;
using std::cout;
typedef double Scalar;
typedef Tpetra::CrsMatrix<Scalar> CrsMatrix;
typedef Tpetra::MultiVector<Scalar> MV;
typedef Tpetra::Operator<Scalar> OP;
typedef Tpetra::MatrixMarket::Reader<CrsMatrix> Reader;
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");
std::string filenameB ("/home/amklinv/matrices/bcsstm06.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("fileB",&filenameB, "Filename for the Matrix-Market mass 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 = Reader::readSparseFile(filenameA, comm);
RCP<const CrsMatrix> M = Reader::readSparseFile(filenameB, comm);
Scalar mat_norm = std::max(K->getFrobeniusNorm(),M->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);
RCP<MV> ivec = rcp (new MV (K->getRowMap (), 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 ) );
MV Mvec( M->getRowMap(), MVT::GetNumberVecs( *evecs ) );
OPT::Apply( *K, *evecs, Kvec );
OPT::Apply( *M, *evecs, Mvec );
MVT::MvTimesMatAddMv( -1.0, Mvec, T, 1.0, Kvec );
MVT::MvNorm( Kvec, normR );
if (myRank == 0) {
cout.setf(std::ios_base::right, std::ios_base::adjustfield);
cout<<"Actual Eigenvalues: "<<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.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
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.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
int numVecs
The number of computed eigenpairs.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.