Use LOBPCG with Tpetra, with custom StatusTest.
Use LOBPCG with Tpetra, with custom StatusTest.
#include "Teuchos_GlobalMPISession.hpp"
#include "Tpetra_Core.hpp"
#include "MatrixMarket_Tpetra.hpp"
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::ArrayView;
using std::cout;
using std::endl;
typedef double Scalar;
typedef Tpetra::MultiVector<Scalar> TMV;
typedef Tpetra::Vector<Scalar> Vector;
typedef Tpetra::Operator<Scalar> TOP;
namespace {
class FoldOp : public TOP {
public:
typedef Tpetra::Map<> map_type;
FoldOp (const RCP<const TOP> A) { A_ = A; };
int SetUseTranspose (bool UseTranspose) { return -1; };
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;
RCP<const map_type> getDomainMap () const { return A_->getDomainMap (); };
RCP<const map_type> getRangeMap () const { return A_->getRangeMap (); };
private:
RCP<const TOP> A_;
};
void
FoldOp::apply (const TMV &X, TMV &Y, Teuchos::ETransp mode,
Scalar alpha, Scalar beta) const
{
TMV Y1 (X.getMap (), X.getNumVectors (), false);
A_->apply (X, Y1, mode, alpha, beta);
A_->apply (Y1, Y, mode, alpha, beta);
}
public:
StatusTestFolding (Scalar tol, int quorum = -1,
bool scaled = true,
bool throwExceptionOnNan = true,
const RCP<const TOP>& A = Teuchos::null);
virtual ~StatusTestFolding() {};
std::vector<int>
whichVecs ()
const {
return ind_; }
int howMany ()
const {
return ind_.size (); }
ind_.resize (0);
}
std::ostream&
print (std::ostream &os,
int indent=0)
const;
private:
Scalar tol_;
std::vector<int> ind_;
int quorum_;
bool scaled_;
RCP<const TOP> A_;
const Scalar ONE;
};
StatusTestFolding::
StatusTestFolding (Scalar tol, int quorum, bool scaled,
bool ,
const RCP<const TOP>& A)
tol_ (tol),
quorum_ (quorum),
scaled_ (scaled),
A_ (A),
ONE (Teuchos::ScalarTraits<Scalar>::one ())
{}
{
const int numev = X->getNumVectors ();
std::vector<Scalar> res (numev);
TMV AX (X->getMap (), numev, false);
Teuchos::SerialDenseMatrix<int, Scalar> T (numev, numev);
A_->apply (*X, AX);
TMVT::MvTransMv (1.0, AX, *X, T);
TMVT::MvTimesMatAddMv (-1.0, *X, T, 1.0, AX);
TMVT::MvNorm (AX, res);
if (scaled_) {
for (int i = 0; i < numev; ++i) {
res[i] /= std::abs (T(i,i));
}
}
ind_.resize (0);
for (int i = 0; i < numev; ++i) {
if (res[i] < tol_) {
ind_.push_back (i);
}
}
const int have = ind_.size ();
const int need = (quorum_ == -1) ? numev : quorum_;
return state_;
}
std::ostream&
StatusTestFolding::print (std::ostream& os, int indent) const
{
std::string ind (indent, ' ');
os << ind << "- StatusTestFolding: ";
switch (state_) {
os << "Passed\n";
break;
os << "Failed\n";
break;
os << "Undefined\n";
break;
}
os << ind << " (Tolerance, WhichNorm,Scaled,Quorum): "
<< "(" << tol_
<< ",RES_2NORM"
<< "," << (scaled_ ? "true" : "false")
<< "," << quorum_
<< ")\n";
os << ind << " Which vectors: ";
if (ind_.size () > 0) {
for (size_t i = 0; i < ind_.size (); ++i) {
os << ind_[i] << " ";
}
os << std::endl;
}
else {
os << "[empty]\n";
}
}
return os;
}
}
int
main (int argc, char* argv[])
{
typedef Tpetra::CrsMatrix<> CrsMatrix;
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 fileA ("/u/slotnick_s2/aklinvex/matrices/anderson4.mtx");
Teuchos::CommandLineProcessor cmdp (false, true);
cmdp.setOption ("fileA", &fileA, "Filename for the Matrix-Market stiffness matrix.");
if (cmdp.parse (argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
RCP<const CrsMatrix> A = Reader::readSparseFile (fileA, comm);
RCP<FoldOp> K = rcp (new FoldOp (A));
int blockSize = 4;
double tol = 1e-5;
bool scaled = true;
int nev = 4;
Teuchos::ParameterList MyPL;
MyPL.set ("Which", "SR");
MyPL.set ("Maximum Restarts", 10000);
MyPL.set ("Maximum Iterations", 1000000);
MyPL.set ("Block Size", blockSize);
MyPL.set ("Convergence Tolerance", tol );
MyPL.set ("Relative Convergence Tolerance", scaled);
MyPL.set ("Relative Locking Tolerance", scaled);
RCP<TMV> ivec = rcp (new TMV (A->getRowMap (), blockSize));
TMVT::MvRandom (*ivec);
RCP<Problem> MyProblem = rcp (new Problem (K, ivec));
MyProblem->setHermitian (true);
MyProblem->setNEV (nev);
MyProblem->setProblem ();
RCP<StatusTestFolding> convTest =
rcp (new StatusTestFolding (tol, nev, scaled, true, A));
RCP<StatusTestFolding> lockTest =
rcp (new StatusTestFolding (tol/10., 1, scaled, true, A));
solver.setGlobalStatusTest (convTest);
solver.setLockingStatusTest (lockTest);
cout << "The solve did NOT converge." << endl;
} else if (myRank == 0) {
cout << "The solve converged." << endl;
}
std::vector<Anasazi::Value<Scalar> > evals = sol.
Evals;
RCP<TMV> evecs = sol.
Evecs;
if (numev > 0) {
std::vector<Scalar> normR (sol.
numVecs);
TMV Avec (A->getRowMap (), TMVT::GetNumberVecs (*evecs));
TOPT::Apply (*A, *evecs, Avec);
Teuchos::SerialDenseMatrix<int,Scalar> T (numev, numev);
TMVT::MvTransMv (1.0, Avec, *evecs, T);
TMVT::MvTimesMatAddMv (-1.0, *evecs, T, 1.0, Avec);
TMVT::MvNorm (Avec, 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)<<T(i,i)
<<std::setw(16)<<normR[i]/std::abs(T(i,i))
<<std::endl;
}
cout<<"------------------------------------------------------"<<std::endl;
}
}
return 0;
}
Basic implementation of the Anasazi::Eigenproblem class.
The Anasazi::LOBPCGSolMgr provides a powerful solver manager for the LOBPCG eigensolver.
Forward declaration of pure virtual base class Anasazi::StatusTest.
Partial specialization of Anasazi::MultiVecTraits and Anasazi::OperatorTraits for Tpetra objects.
Types and exceptions used within Anasazi solvers and interfaces.
This provides a basic implementation for defining standard or generalized eigenvalue problems.
The Eigensolver is a templated virtual base class that defines the basic interface that any eigensolv...
virtual Teuchos::RCP< const MV > getRitzVectors()=0
Get the Ritz vectors from the previous iteration. These are indexed using getRitzIndex().
User interface for the LOBPCG eigensolver.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Common interface of stopping criteria for Anasazi's solvers.
virtual std::vector< int > whichVecs() const =0
Get the indices for the vectors that passed the test.
virtual int howMany() const =0
Get the number of vectors that passed the test.
virtual void reset()=0
Informs the status test that it should reset its internal configuration to the uninitialized state.
virtual TestStatus checkStatus(Eigensolver< ScalarType, MV, OP > *solver)=0
virtual std::ostream & print(std::ostream &os, int indent=0) const =0
Output formatted description of stopping test to output stream.
virtual void clearStatus()=0
Clears the results of the last status test.
virtual TestStatus getStatus() const =0
Return the result of the most recent checkStatus call, or undefined if it has not been run.
ReturnType
Enumerated type used to pass back information from a solver manager.
TestStatus
Enumerated type used to pass back information from a StatusTest.
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.