#include "Epetra_CrsMatrix.h"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_StandardCatchMacros.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(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]);
assert( info==0 );
}
info = A->FillComplete(ColMap, RowMap);
assert( info==0 );
info = A->OptimizeStorage();
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 );
}