This is an example of how to use the TraceMinDavidsonSolMgr solver manager to compute the Fiedler vector, using Tpetra data stuctures and an Ifpack2 preconditioner.
#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_ArrayViewDecl.hpp"
#include "Teuchos_ParameterList.hpp"
#ifdef HAVE_ANASAZI_IFPACK2
#include "Ifpack2_Factory.hpp"
#include "Ifpack2_Preconditioner.hpp"
#endif
using Teuchos::RCP;
using std::cout;
using std::cin;
using Scalar = double;
using CrsMatrix = Tpetra::CrsMatrix<Scalar>;
using Vector = Tpetra::Vector<Scalar>;
using TMV = Tpetra::MultiVector<Scalar>;
using TOP = Tpetra::Operator<Scalar>;
void formLaplacian(const RCP<const CrsMatrix>& A, const bool weighted, const bool normalized, RCP<CrsMatrix>& L, RCP<Vector>& auxVec);
int main(int argc, char *argv[]) {
Tpetra::ScopeGuard tpetraScope(&argc, &argv);
RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
const int myRank = comm->getRank();
std::string inputFilename("/home/amklinv/matrices/mesh1em6_Laplacian.mtx");
std::string outputFilename("/home/amklinv/matrices/mesh1em6_Fiedler.mtx");
Scalar tol = 1e-6;
int nev = 1;
int blockSize = 1;
bool usePrec = false;
bool useNormalizedLaplacian = false;
bool useWeightedLaplacian = false;
bool verbose = true;
std::string whenToShift = "Always";
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("fin",&inputFilename, "Filename for Matrix-Market test matrix.");
cmdp.setOption("fout",&outputFilename, "Filename for Fiedler vector.");
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("precondition","no-precondition",&usePrec, "Whether to use a diagonal preconditioner.");
cmdp.setOption("normalization","no-normalization",&useNormalizedLaplacian, "Whether to normalize the laplacian.");
cmdp.setOption("weighted","unweighted",&useWeightedLaplacian, "Whether to normalize the laplacian.");
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> fileMat =
Tpetra::MatrixMarket::Reader<CrsMatrix>::readSparseFile(inputFilename, comm);
RCP<CrsMatrix> L;
RCP<Vector> auxVec;
formLaplacian(fileMat, useWeightedLaplacian, useNormalizedLaplacian, L, auxVec);
RCP<const CrsMatrix> K = L;
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("Saddle Solver Type", "Block Diagonal Preconditioned Minres");
RCP<TMV> ivec = Teuchos::rcp( new TMV(K->getRowMap(), blockSize) );
TMVT::MvRandom( *ivec );
RCP<Anasazi::BasicEigenproblem<Scalar,TMV,TOP> > MyProblem =
MyProblem->setHermitian(true);
MyProblem->setNEV( nev );
if(usePrec)
{
#ifdef HAVE_ANASAZI_IFPACK2
Ifpack2::Factory factory;
const std::string precType = "RELAXATION";
Teuchos::ParameterList PrecPL;
PrecPL.set( "relaxation: type", "Jacobi");
RCP<Ifpack2::Preconditioner<Scalar> > Prec = factory.create(precType, K);
assert(Prec != Teuchos::null);
Prec->setParameters(PrecPL);
Prec->initialize();
Prec->compute();
MyProblem->setPrec(Prec);
#else
if(myRank == 0)
cout << "You did not build Trilinos with Ifpack2 preconditioning enabled. Please either\n1. Reinstall Trilinos with Ifpack2 enabled\n2. Try running this driver again without preconditioning enabled\n";
return -1;
#endif
}
MyProblem->setAuxVecs(auxVec);
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<TMV> evecs = sol.Evecs;
int numev = sol.numVecs;
if (numev > 0) {
Teuchos::SerialDenseMatrix<int,Scalar> T(numev,numev);
TMV tempvec(K->getRowMap(), TMVT::GetNumberVecs( *evecs ));
std::vector<Scalar> normR(sol.numVecs);
TMV Kvec( K->getRowMap(), 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) : "<<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)<<T(i,i)
<<std::setw(16)<<normR[i]/mat_norm
<<std::endl;
}
cout<<"------------------------------------------------------"<<std::endl;
}
}
if (numev > 0) {
Tpetra::MatrixMarket::Writer<CrsMatrix>::writeDenseFile(outputFilename,evecs,"","Fiedler vector of "+inputFilename);
}
return 0;
}
void formLaplacian(const RCP<const CrsMatrix>& A, const bool weighted, const bool normalized, RCP<CrsMatrix>& L, RCP<Vector>& auxVec)
{
using SCT = Teuchos::ScalarTraits<Scalar>;
using LO = Tpetra::Map<>::local_ordinal_type;
using GO = Tpetra::Map<>::global_ordinal_type;
Scalar ONE = SCT::one();
Scalar ZERO = SCT::zero();
std::vector<GO> diagIndex(static_cast<GO>(1));
std::vector<LO> diagIndexLcl(static_cast<LO>(1));
std::vector<Scalar> value(1,ONE);
Teuchos::ArrayView<const GO> cols(diagIndex);
Teuchos::ArrayView<const LO> colsLcl(diagIndexLcl);
Teuchos::ArrayView<const Scalar> vals(value);
const GO n = A->getGlobalNumRows();
L = Tpetra::MatrixMatrix::add(ONE,false,*A,ONE,true,*A);
RCP<const Tpetra::Map<> > rowMap = L->getRowMap();
L->resumeFill();
if(weighted)
{
using indices_view = typename CrsMatrix::nonconst_global_inds_host_view_type;
using values_view = typename CrsMatrix::nonconst_values_host_view_type;
indices_view colIndices("colIndices");
values_view values("values");
RCP<Vector> diagonal = Teuchos::rcp(new Vector(rowMap));
for(GO i=0; i<n; i++)
{
if(rowMap->isNodeGlobalElement(i))
{
size_t numentries = L->getNumEntriesInGlobalRow(i);
Kokkos::resize(colIndices,numentries);
Kokkos::resize(values,numentries);
L->getGlobalRowCopy(i,colIndices,values,numentries);
for(size_t j=0; j<colIndices.size(); j++)
{
if(i == rowMap->getGlobalElement(colIndices[j]))
values[j] = ZERO;
else
values[j] = -abs(values[j]);
diagonal->sumIntoGlobalValue(i,-values[j]);
}
L->replaceGlobalValues(i, colIndices, values);
}
}
Teuchos::ArrayRCP<const Scalar> diagView = diagonal->getData();
using size_type = Teuchos::ArrayRCP<const Scalar>::size_type;
auxVec = Teuchos::rcp(new Vector(rowMap,false));
if (normalized) {
for (size_type i = 0; i < diagView.size (); ++i) {
auxVec->replaceLocalValue (static_cast<LO> (i), sqrt (diagView[i]));
}
}
else {
auxVec->putScalar(ONE);
}
Scalar vecNorm = auxVec->norm2();
auxVec->scale(ONE/vecNorm);
if (normalized) {
Vector scaleVec (rowMap, false);
for(size_type i=0; i<diagView.size(); i++)
{
scaleVec.replaceLocalValue(static_cast<LO>(i),ONE/sqrt(diagView[i]));
}
L->fillComplete();
L->leftScale(scaleVec);
L->rightScale(scaleVec);
L->resumeFill();
for(GO i=0; i<n; i++)
{
diagIndex[0] = i;
if(rowMap->isNodeGlobalElement(i)) L->replaceGlobalValues(i,cols,vals);
}
}
else
{
for(size_type i=0; i<diagView.size(); i++)
{
diagIndex[0] = rowMap->getGlobalElement(i);
value[0] = diagView[i];
L->replaceLocalValues(i,colsLcl,vals);
}
}
L->fillComplete();
}
else {
for (GO i = 0; i < n; ++i) {
diagIndex[0] = i;
if(rowMap->isNodeGlobalElement(i)) {
L->replaceGlobalValues(i,cols,vals);
}
}
L->setAllToScalar(-ONE);
L->fillComplete();
auxVec = Teuchos::rcp( new Vector (rowMap,false) );
if (normalized) {
for (GO i = 0; i < n; ++i) {
if (rowMap->isNodeGlobalElement(i)) {
Scalar temp;
size_t nnzInRow = L->getNumEntriesInGlobalRow(i) - static_cast<size_t> (1);
temp = sqrt(nnzInRow);
auxVec->replaceGlobalValue(i,temp);
}
}
}
else {
auxVec->putScalar(ONE);
}
Scalar vecNorm = auxVec->norm2();
auxVec->scale(ONE/vecNorm);
if (normalized) {
Vector scalars(rowMap,false);
for (GO i = 0; i < n; ++i) {
if(rowMap->isNodeGlobalElement(i)) {
Scalar temp;
size_t nnzInRow = L->getNumEntriesInGlobalRow(i) - static_cast<size_t> (1);
temp = ONE/sqrt(nnzInRow);
scalars.replaceGlobalValue(i,temp);
}
}
L->leftScale(scalars);
L->rightScale(scalars);
L->resumeFill();
for (GO i = 0; i < n; ++i) {
diagIndex[0] = i;
if(rowMap->isNodeGlobalElement(i)) L->replaceGlobalValues(i,cols,vals);
}
}
else {
L->resumeFill();
for (GO i = 0; i < n; ++i) {
if(rowMap->isNodeGlobalElement(i)) {
size_t nnzInRow = L->getNumEntriesInGlobalRow(i) - static_cast<size_t> (1);
diagIndex[0] = i;
value[0] = nnzInRow;
L->replaceGlobalValues(i,cols,vals);
}
}
}
L->fillComplete();
}
}
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.