Skip to content

Commit

Permalink
Merge pull request #12166 from NexGenAnalytics/nga-fy23-pr-candidate-8
Browse files Browse the repository at this point in the history 
Anasazi: Add global MPI communicator that can be used instead of hardcoded `MPI_COMM_WORLD`
  • Loading branch information
@hkthorn
hkthorn authored Aug 25, 2023
2 parents 9f2a806 + 8d7e5d7 commit 8db46da
Show file tree
Hide file tree
Showing 11 changed files with 238 additions and 160 deletions.
28 changes: 14 additions & 14 deletions packages/anasazi/epetra/src-rbgen/RBGenDriver_EpetraMV.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_Assert.hpp"

int main( int argc, char* argv[] )
{

Expand Down Expand Up @@ -111,13 +111,13 @@ int main( int argc, char* argv[] )
// ---------------------------------------------------------------
//
Teuchos::RCP<Teuchos::ParameterList> BasisParams = RBGen::createParams( xml_file );
if (verbose && Comm.MyPID() == 0)
if (verbose && Comm.MyPID() == 0)
{
std::cout<<"-------------------------------------------------------"<<std::endl;
std::cout<<"Input Parameter List: "<<std::endl;
std::cout<<"-------------------------------------------------------"<<std::endl;
BasisParams->print();
}
}
//
// ---------------------------------------------------------------
// CREATE THE FILE I/O HANDLER
Expand All @@ -134,17 +134,17 @@ int main( int argc, char* argv[] )
//
Teuchos::RCP< RBGen::FileIOHandler<Epetra_MultiVector> > mvFileIO;
Teuchos::RCP< RBGen::FileIOHandler<Epetra_Operator> > opFileIO =
Teuchos::rcp( new RBGen::EpetraCrsMatrixFileIOHandler() );
Teuchos::rcp( new RBGen::EpetraCrsMatrixFileIOHandler() );
{
Teuchos::TimeMonitor lcltimer( *timerFileIO );
mvFileIO = fio_factory.create( *BasisParams );
//
//
// Initialize file IO handlers
//
mvFileIO->Initialize( BasisParams );
opFileIO->Initialize( BasisParams );
}
if (verbose && Comm.MyPID() == 0)
}
if (verbose && Comm.MyPID() == 0)
{
std::cout<<"-------------------------------------------------------"<<std::endl;
std::cout<<"File I/O Handlers Generated"<<std::endl;
Expand All @@ -164,7 +164,7 @@ int main( int argc, char* argv[] )
{
Teuchos::TimeMonitor lcltimer( *timerSnapshotIn );
testMV = mvFileIO->Read( *filenames );
}
}

RBGen::EpetraMVPreprocessorFactory preprocess_factory;

Expand All @@ -180,14 +180,14 @@ int main( int argc, char* argv[] )
prep->Initialize( BasisParams, mvFileIO );
}

Teuchos::RCP<Teuchos::Time> timerPreprocess = Teuchos::rcp( new Teuchos::Time("Preprocess Snapshot Set") );
Teuchos::RCP<Teuchos::Time> timerPreprocess = Teuchos::rcp( new Teuchos::Time("Preprocess Snapshot Set") );
timersRBGen.push_back( timerPreprocess );
{
Teuchos::TimeMonitor lcltimer( *timerPreprocess );
prep->Preprocess( testMV );
}

if (verbose && Comm.MyPID() == 0)
if (verbose && Comm.MyPID() == 0)
{
std::cout<<"-------------------------------------------------------"<<std::endl;
std::cout<<"Snapshot Set Imported and Preprocessed"<<std::endl;
Expand All @@ -208,7 +208,7 @@ int main( int argc, char* argv[] )
timersRBGen.push_back( timerCreateMethod );
Teuchos::RCP<RBGen::Method<Epetra_MultiVector,Epetra_Operator> > method;
{
Teuchos::TimeMonitor lcltimer( *timerCreateMethod );
Teuchos::TimeMonitor lcltimer( *timerCreateMethod );
method = mthd_factory.create( *BasisParams );
//
// Initialize reduced basis method.
Expand All @@ -221,7 +221,7 @@ int main( int argc, char* argv[] )
Teuchos::RCP<Teuchos::Time> timerComputeBasis = Teuchos::rcp( new Teuchos::Time("Reduced Basis Computation") );
timersRBGen.push_back( timerComputeBasis );
{
Teuchos::TimeMonitor lcltimer( *timerComputeBasis );
Teuchos::TimeMonitor lcltimer( *timerComputeBasis );
method->computeBasis();
}
//
Expand All @@ -239,8 +239,8 @@ int main( int argc, char* argv[] )
std::cout<<"Computed Singular Values : "<<std::endl;
std::cout<<"-------------------------------------------------------"<<std::endl;
for (unsigned int i=0; i<sv.size(); ++i) { std::cout << sv[i] << std::endl; }
}
}

if (Comm.MyPID() == 0) {
std::cout<<"-------------------------------------------------------"<<std::endl;
std::cout<<"RBGen Computation Time Breakdown (seconds) : "<<std::endl;
Expand Down
67 changes: 34 additions & 33 deletions packages/anasazi/epetra/src-rbgen/RBGen_AnasaziPOD.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,17 +52,18 @@

#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#include "AnasaziGlobalComm.hpp"
#else
#include "Epetra_SerialComm.h"
#endif

namespace RBGen {

AnasaziPOD::AnasaziPOD() :
isInitialized_( false ),
isInner_( true ),
basis_size_( 16 ),
comp_time_( 0.0 )
comp_time_( 0.0 )
{
}

Expand All @@ -76,8 +77,8 @@ namespace RBGen {

if ( rbmethod_params.get("Basis Size", 16) < ss->NumVectors() ) {
basis_size_ = rbmethod_params.get("Basis Size", 16);
}
else {
}
else {
basis_size_ = ss->NumVectors();
}
// Get the inner / outer product form of the operator
Expand All @@ -97,20 +98,20 @@ namespace RBGen {
inner_prod_op_ = fileio->Read( filename );
}

// Resize the singular value vector
sv_.resize( basis_size_ );
// Resize the singular value vector
sv_.resize( basis_size_ );

// Set the snapshot set
ss_ = ss;
}

void AnasaziPOD::computeBasis()
{
{
#ifdef EPETRA_MPI
Epetra_MpiComm comm( MPI_COMM_WORLD );
Epetra_MpiComm comm( get_global_comm() );
#else
Epetra_SerialComm comm;
#endif
#endif

int MyPID = comm.MyPID();
//
Expand Down Expand Up @@ -156,7 +157,7 @@ namespace RBGen {
MyPL.set( "Block Size", blockSize );
MyPL.set( "Num Blocks", maxBlocks );
MyPL.set( "Maximum Restarts", maxRestarts );
MyPL.set( "Which", which );
MyPL.set( "Which", which );
MyPL.set( "Step Size", step );
MyPL.set( "Convergence Tolerance", tol );
//
Expand Down Expand Up @@ -203,24 +204,24 @@ namespace RBGen {
// Create the eigenproblem
Teuchos::RCP<Anasazi::BasicEigenproblem<double,MV,OP> > MyProblem =
Teuchos::rcp( new Anasazi::BasicEigenproblem<double,MV,OP>(Amat, ivec) );

// Inform the eigenproblem that the operator A is symmetric
MyProblem->setHermitian( true );
MyProblem->setHermitian( true );

// Set the number of eigenvalues requested
MyProblem->setNEV( nev );

// Inform the eigenproblem that you are finishing passing it information
bool boolret = MyProblem->setProblem();
if (boolret != true) {
if (MyPID == 0) {
std::cout << "Anasazi::BasicEigenproblem::setProblem() returned with error." << std::endl;
}
}

// Initialize the Block Arnoldi solver
Anasazi::BlockKrylovSchurSolMgr<double,MV,OP> MySolverMgr(MyProblem, MyPL);

timer.ResetStartTime();

// Solve the problem to the specified tolerances or length
Expand All @@ -239,43 +240,43 @@ namespace RBGen {
if (numev < nev && MyPID==0) {
std::cout << "RBGen::AnasaziPOD will return " << numev << " singular vectors." << std::endl;
}

// Resize the singular value vector to the number of computed singular pairs.
sv_.resize( numev );
sv_.resize( numev );

if (numev > 0) {
// Retrieve eigenvectors
std::vector<int> index(numev);
for (i=0; i<numev; i++) { index[i] = i; }
Teuchos::RCP<const Anasazi::EpetraMultiVec> evecs = Teuchos::rcp_dynamic_cast<const Anasazi::EpetraMultiVec>(
Teuchos::RCP<const Anasazi::EpetraMultiVec> evecs = Teuchos::rcp_dynamic_cast<const Anasazi::EpetraMultiVec>(
Anasazi::MultiVecTraits<double,MV>::CloneView(*sol.Evecs, index) );
// const Anasazi::EpetraMultiVec* evecs = dynamic_cast<const Anasazi::EpetraMultiVec *>(sol.Evecs->CloneView( index ));
// const Anasazi::EpetraMultiVec* evecs = dynamic_cast<const Anasazi::EpetraMultiVec *>(sol.Evecs->CloneView( index ));
//
// Compute singular values which are the square root of the eigenvalues
//
for (i=0; i<numev; i++) {
sv_[i] = Teuchos::ScalarTraits<double>::squareroot( Teuchos::ScalarTraits<double>::magnitude(evals[i].realpart) );
for (i=0; i<numev; i++) {
sv_[i] = Teuchos::ScalarTraits<double>::squareroot( Teuchos::ScalarTraits<double>::magnitude(evals[i].realpart) );
}
//
// If we are using inner product formulation,
// If we are using inner product formulation,
// then we must compute left singular vectors:
// u = Av/sigma
//
std::vector<double> tempnrm( numev );
if (isInner_) {
basis_ = Teuchos::rcp( new Epetra_MultiVector(ss_->Map(), numev) );
Epetra_MultiVector AV( ss_->Map(), numev );
/* A*V */

/* A*V */
AV.Multiply( 'N', 'N', 1.0, *ss_, *evecs, 0.0 );
AV.Norm2( &tempnrm[0] );

/* U = A*V(i)/S(i) */
Epetra_LocalMap localMap2( numev, 0, ss_->Map().Comm() );
Epetra_MultiVector S( localMap2, numev );
for( i=0; i<numev; i++ ) { S[i][i] = 1.0/tempnrm[i]; }
basis_->Multiply( 'N', 'N', 1.0, AV, S, 0.0 );

/* Compute direct residuals : || Av - sigma*u ||
for (i=0; i<nev; i++) { S[i][i] = sv_[i]; }
info = AV.Multiply( 'N', 'N', -1.0, *basis_, S, 1.0 );
Expand All @@ -292,19 +293,19 @@ namespace RBGen {
} else {
basis_ = Teuchos::rcp( new Epetra_MultiVector( Copy, *evecs, 0, numev ) );
Epetra_MultiVector ATU( localMap, numev );
Epetra_MultiVector V( localMap, numev );
/* A^T*U */
Epetra_MultiVector V( localMap, numev );

/* A^T*U */
ATU.Multiply( 'T', 'N', 1.0, *ss_, *evecs, 0.0 );
ATU.Norm2( &tempnrm[0] );

/* V = A^T*U(i)/S(i) */
Epetra_LocalMap localMap2( numev, 0, ss_->Map().Comm() );
Epetra_MultiVector S( localMap2, numev );
for( i=0; i<numev; i++ ) { S[i][i] = 1.0/tempnrm[i]; }
V.Multiply( 'N', 'N', 1.0, ATU, S, 0.0 );
/* Compute direct residuals : || (A^T)u - sigma*v ||

/* Compute direct residuals : || (A^T)u - sigma*v ||
for (i=0; i<nev; i++) { S[i][i] = sv_[i]; }
info = ATU.Multiply( 'N', 'N', -1.0, V, S, 1.0 );
ATU.Norm2( tempnrm );
Expand Down
Loading

0 comments on commit 8db46da

Please sign in to comment.