#include "Epetra_CrsMatrix.h"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
#include "Teuchos_Assert.hpp"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#include <mpi.h>
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Map.h"
int main(int argc, char *argv[]) {
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
#endif
bool success = false;
try {
#ifdef EPETRA_MPI
Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
Epetra_SerialComm Comm;
#endif
int i, j, info;
const double one = 1.0;
const double zero = 0.0;
Teuchos::LAPACK<int,double> lapack;
int MyPID = Comm.MyPID();
int m = 500;
int n = 100;
Epetra_Map RowMap(m, 0, Comm);
Epetra_Map ColMap(n, 0, Comm);
int NumMyRowElements = RowMap.NumMyElements();
std::vector<int> MyGlobalRowElements(NumMyRowElements);
RowMap.MyGlobalElements(&MyGlobalRowElements[0]);
Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::Copy, RowMap, n) );
std::vector<double> Values(n);
std::vector<int> Indices(n);
double inv_mp1 = one/(m+1);
double inv_np1 = one/(n+1);
for (i=0; i<n; i++) { Indices[i] = i; }
for (i=0; i<NumMyRowElements; i++) {
for (j=0; j<n; j++) {
if ( MyGlobalRowElements[i] <= j ) {
Values[j] = inv_np1 * ( (MyGlobalRowElements[i]+one)*inv_mp1 ) * ( (j+one)*inv_np1 - one );
}
else {
Values[j] = inv_np1 * ( (j+one)*inv_np1 ) * ( (MyGlobalRowElements[i]+one)*inv_mp1 - one );
}
}
info = A->InsertGlobalValues(MyGlobalRowElements[i], n, &Values[0], &Indices[0]);
TEUCHOS_ASSERT( info==0 );
}
info = A->FillComplete(ColMap, RowMap);
TEUCHOS_ASSERT( info==0 );
info = A->OptimizeStorage();
TEUCHOS_ASSERT( info==0 );
A->SetTracebackMode(1);
int nev = 4;
int blockSize = 1;
int numBlocks = 10;
int maxRestarts = 20;
double tol = lapack.LAMCH('E');
std::string which = "LM";
Teuchos::ParameterList MyPL;
MyPL.set( "Verbosity", verbosity );
MyPL.set( "Which", which );
MyPL.set( "Block Size", blockSize );
MyPL.set( "Num Blocks", numBlocks );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Convergence Tolerance", tol );
ivec->MvRandom();
Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> > MyProblem =
MyProblem->setHermitian(true);
MyProblem->setNEV( nev );
bool boolret = MyProblem->setProblem();
if (!boolret) {
throw "Anasazi::BasicEigenproblem::setProblem() returned with error.";
}
std::cout << "Anasazi::EigensolverMgr::solve() returned unconverged." << std::endl;
}
std::vector<Anasazi::Value<double> > evals = sol.
Evals;
if (numev > 0) {
if (MyPID==0) {
std::cout<<"------------------------------------------------------"<<std::endl;
std::cout<<"Computed Singular Values: "<<std::endl;
std::cout<<"------------------------------------------------------"<<std::endl;
}
for (i=0; i<numev; i++) { evals[i].realpart = Teuchos::ScalarTraits<double>::squareroot( evals[i].realpart ); }
std::vector<double> tempnrm(numev), directnrm(numev);
std::vector<int> index(numev);
for (i=0; i<numev; i++) { index[i] = i; }
Teuchos::SerialDenseMatrix<int,double> S(numev,numev);
A->Apply( *evecs, Av );
Av.MvNorm( tempnrm );
for (i=0; i<numev; i++) { S(i,i) = one/tempnrm[i]; };
for (i=0; i<numev; i++) { S(i,i) = evals[i].realpart; }
Av.MvTimesMatAddMv( -one, u, S, one );
Av.MvNorm( directnrm );
if (MyPID==0) {
std::cout.setf(std::ios_base::right, std::ios_base::adjustfield);
std::cout<<std::setw(16)<<"Singular Value"
<<std::setw(20)<<"Direct Residual"
<<std::endl;
std::cout<<"------------------------------------------------------"<<std::endl;
for (i=0; i<numev; i++) {
std::cout<<std::setw(16)<<evals[i].realpart
<<std::setw(20)<<directnrm[i]
<<std::endl;
}
std::cout<<"------------------------------------------------------"<<std::endl;
}
}
success = true;
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(true, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize() ;
#endif
return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
}
Basic implementation of the Anasazi::Eigenproblem class.
The Anasazi::BlockKrylovSchurSolMgr class provides a user interface for the block Krylov-Schur eigens...
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Declarations of Anasazi multi-vector and operator classes using Epetra_MultiVector and Epetra_Operato...
This provides a basic implementation for defining standard or generalized eigenvalue problems.
The Anasazi::BlockKrylovSchurSolMgr provides a flexible solver manager over the BlockKrylovSchur eige...
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
Basic adapter class for Anasazi::MultiVec that uses Epetra_MultiVector.
void MvTimesMatAddMv(double alpha, const MultiVec< double > &A, const Teuchos::SerialDenseMatrix< int, double > &B, double beta)
Update *this with .
Adapter class for creating a symmetric operator from an Epetra_Operator.
Interface for multivectors used by Anasazi's linear solvers.
Anasazi's templated virtual class for constructing an operator that can interface with the OperatorTr...
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.