Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Amesos2 : Use Teuchos to gather matrix #2 #13740

Merged
merged 2 commits into from
Jan 23, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions packages/amesos2/cmake/Amesos2_config.h.in
Original file line number Diff line number Diff line change
Expand Up @@ -116,3 +116,4 @@
#ifdef HAVE_AMESOS2_ZOLTAN2CORE
# define HAVE_AMESOS2_ZOLTAN2
#endif
#cmakedefine HAVE_AMESOS2_GALERI
2 changes: 1 addition & 1 deletion packages/amesos2/cmake/Dependencies.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
SET(LIB_REQUIRED_DEP_PACKAGES Teuchos Tpetra TrilinosSS Kokkos)
SET(LIB_OPTIONAL_DEP_PACKAGES Epetra EpetraExt ShyLU_NodeBasker ShyLU_NodeTacho)
SET(TEST_REQUIRED_DEP_PACKAGES)
SET(TEST_OPTIONAL_DEP_PACKAGES ShyLU_NodeBasker ShyLU_NodeTacho Kokkos TrilinosSS Xpetra Zoltan2Core)
SET(TEST_OPTIONAL_DEP_PACKAGES ShyLU_NodeBasker ShyLU_NodeTacho Kokkos TrilinosSS Xpetra Zoltan2Core Galeri)
# SET(LIB_REQUIRED_DEP_TPLS SuperLU)
SET(LIB_REQUIRED_DEP_TPLS )
SET(LIB_OPTIONAL_DEP_TPLS MPI SuperLU SuperLUMT SuperLUDist LAPACK UMFPACK PARDISO_MKL CSS_MKL ParMETIS METIS Cholmod MUMPS STRUMPACK CUSPARSE CUSOLVER)
Expand Down
199 changes: 127 additions & 72 deletions packages/amesos2/example/SimpleSolve.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@
#include <Teuchos_oblackholestream.hpp>
#include <Teuchos_Tuple.hpp>
#include <Teuchos_VerboseObject.hpp>
#include <Teuchos_StackedTimer.hpp>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_ParameterXMLFileReader.hpp>

#include <Tpetra_Core.hpp>
#include <Tpetra_Map.hpp>
Expand All @@ -31,7 +34,10 @@

#include "Amesos2.hpp"
#include "Amesos2_Version.hpp"

#if defined(HAVE_AMESOS2_XPETRA) && defined(HAVE_AMESOS2_GALERI)
#include "Galeri_XpetraMaps.hpp"
#include "Galeri_XpetraProblemFactory.hpp"
#endif

int main(int argc, char *argv[]) {
Tpetra::ScopeGuard tpetraScope(&argc,&argv);
Expand All @@ -40,6 +46,7 @@ int main(int argc, char *argv[]) {
typedef Tpetra::Map<>::local_ordinal_type LO;
typedef Tpetra::Map<>::global_ordinal_type GO;

typedef Tpetra::Map<LO,GO> MAP;
typedef Tpetra::CrsMatrix<Scalar,LO,GO> MAT;
typedef Tpetra::MultiVector<Scalar,LO,GO> MV;

Expand All @@ -48,9 +55,24 @@ int main(int argc, char *argv[]) {
using Teuchos::RCP;
using Teuchos::rcp;

bool verbose = false;
GO nx = 1;
std::string GaleriName3D {"Laplace3D"};
std::string solvername("Superlu");
std::string xml_filename("");
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
cmdp.setOption("solvername",&solvername,"Name of solver.");
cmdp.setOption("xml_filename",&xml_filename,"XML Filename for Solver parameters.");
cmdp.setOption("nx",&nx,"Dimension of 3D problem.");
cmdp.setOption ("galeriMatrixName", &GaleriName3D, "Name of 3D Galeri Matrix");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}

// Before we do anything, check that SuperLU is enabled
if( !Amesos2::query("SuperLU") ){
std::cerr << "SuperLU not enabled. Exiting..." << std::endl;
if( !Amesos2::query(solvername) ){
std::cerr << solvername << " not enabled. Exiting..." << std::endl;
return EXIT_SUCCESS; // Otherwise CTest will pick it up as
// failure, which it isn't really
}
Expand All @@ -66,82 +88,115 @@ int main(int argc, char *argv[]) {

const size_t numVectors = 1;

// create a Map
global_size_t nrows = 6;
RCP<Tpetra::Map<LO,GO> > map
= rcp( new Tpetra::Map<LO,GO>(nrows,0,comm) );

RCP<MAT> A = rcp( new MAT(map,3) ); // max of three entries in a row

/*
* We will solve a system with a known solution, for which we will be using
* the following matrix:
*
* [ [ 7, 0, -3, 0, -1, 0 ]
* [ 2, 8, 0, 0, 0, 0 ]
* [ 0, 0, 1, 0, 0, 0 ]
* [ -3, 0, 0, 5, 0, 0 ]
* [ 0, -1, 0, 0, 4, 0 ]
* [ 0, 0, 0, -2, 0, 6 ] ]
*
*/
// Construct matrix
if( myRank == 0 ){
A->insertGlobalValues(0,tuple<GO>(0,2,4),tuple<Scalar>(7,-3,-1));
A->insertGlobalValues(1,tuple<GO>(0,1),tuple<Scalar>(2,8));
A->insertGlobalValues(2,tuple<GO>(2),tuple<Scalar>(1));
A->insertGlobalValues(3,tuple<GO>(0,3),tuple<Scalar>(-3,5));
A->insertGlobalValues(4,tuple<GO>(1,4),tuple<Scalar>(-1,4));
A->insertGlobalValues(5,tuple<GO>(3,5),tuple<Scalar>(-2,6));
RCP<MAT> A;
RCP<const MAP> map;
if (nx > 0) {
#if defined(HAVE_AMESOS2_XPETRA) && defined(HAVE_AMESOS2_GALERI)
typedef Galeri::Xpetra::Problem<MAP, MAT, MV> Galeri_t;

Teuchos::ParameterList galeriList;
Tpetra::global_size_t nGlobalElements = nx * nx * nx;
galeriList.set("nx", nx);
galeriList.set("ny", nx);
galeriList.set("nz", nx);
if (GaleriName3D == "Elasticity3D") {
GO mx = 1;
galeriList.set("mx", mx);
galeriList.set("my", mx);
galeriList.set("mz", mx);
nGlobalElements *= 3;
}
map = rcp(new MAP(nGlobalElements, 0, comm));
RCP<Galeri_t> galeriProblem =
Galeri::Xpetra::BuildProblem<Scalar, LO, GO, MAP, MAT, MV>
(GaleriName3D, map, galeriList);
A = galeriProblem->BuildMatrix();
#else
std::cerr << "Galeri or Xpetra not enabled. Exiting..." << std::endl;
return EXIT_SUCCESS; // Otherwise CTest will pick it up as
#endif
} else {
// create a Map
global_size_t nrows = 6;
map = rcp( new MAP(nrows,0,comm) );

RCP<MAT> A = rcp( new MAT(map,3) ); // max of three entries in a row

/*
* We will solve a system with a known solution, for which we will be using
* the following matrix:
*
* [ [ 7, 0, -3, 0, -1, 0 ]
* [ 2, 8, 0, 0, 0, 0 ]
* [ 0, 0, 1, 0, 0, 0 ]
* [ -3, 0, 0, 5, 0, 0 ]
* [ 0, -1, 0, 0, 4, 0 ]
* [ 0, 0, 0, -2, 0, 6 ] ]
*
*/
// Construct matrix
if( myRank == 0 ){
A->insertGlobalValues(0,tuple<GO>(0,2,4),tuple<Scalar>(7,-3,-1));
A->insertGlobalValues(1,tuple<GO>(0,1),tuple<Scalar>(2,8));
A->insertGlobalValues(2,tuple<GO>(2),tuple<Scalar>(1));
A->insertGlobalValues(3,tuple<GO>(0,3),tuple<Scalar>(-3,5));
A->insertGlobalValues(4,tuple<GO>(1,4),tuple<Scalar>(-1,4));
A->insertGlobalValues(5,tuple<GO>(3,5),tuple<Scalar>(-2,6));
}
A->fillComplete();
}
A->fillComplete();

// Create random X
// Create X
RCP<MV> X = rcp(new MV(map,numVectors));
X->randomize();
X->putScalar(1);

/* Create B
*
* Use RHS:
*
* [[-7]
* [18]
* [ 3]
* [17]
* [18]
* [28]]
*/
/* Create B */
RCP<MV> B = rcp(new MV(map,numVectors));
int data[6] = {-7,18,3,17,18,28};
for( int i = 0; i < 6; ++i ){
if( B->getMap()->isNodeGlobalElement(i) ){
B->replaceGlobalValue(i,0,data[i]);
}
}

// Create solver interface to Superlu with Amesos2 factory method
RCP<Amesos2::Solver<MAT,MV> > solver = Amesos2::create<MAT,MV>("Superlu", A, X, B);

solver->symbolicFactorization().numericFactorization().solve();


/* Print the solution
*
* Should be:
*
* [[1]
* [2]
* [3]
* [4]
* [5]
* [6]]
*/
RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
A->apply(*X, *B);
X->randomize();

*fos << "Solution :" << std::endl;
X->describe(*fos,Teuchos::VERB_EXTREME);
*fos << std::endl;
// Create solver interface with Amesos2 factory method
RCP<Amesos2::Solver<MAT,MV> > solver = Amesos2::create<MAT,MV>(solvername, A, X, B);
if (xml_filename != "") {
Teuchos::ParameterList test_params =
Teuchos::ParameterXMLFileReader(xml_filename).getParameters();
Teuchos::ParameterList& amesos2_params = test_params.sublist("Amesos2");
solver->setParameters( Teuchos::rcpFromRef(amesos2_params) );
}

RCP<Teuchos::StackedTimer> stackedTimer;
stackedTimer = rcp(new Teuchos::StackedTimer("Amesos2 SimpleSolve-File"));
Teuchos::TimeMonitor::setStackedTimer(stackedTimer);
{
solver->symbolicFactorization().numericFactorization().solve();
}
stackedTimer->stopBaseTimer();
{
Teuchos::StackedTimer::OutputOptions options;
options.num_histogram=3;
options.print_warnings = false;
options.output_histogram = true;
options.output_fraction=true;
options.output_minmax = true;
stackedTimer->report(std::cout, comm, options);
}
if (verbose) {
/* Print the solution
*
* Should be:
*
* [[1]
* [1]
* [1]
* [1]
* [1]
* [1]]
*/
RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(out));
*fos << "Solution :" << std::endl;
X->describe(*fos,Teuchos::VERB_EXTREME);
*fos << std::endl;
}
// We are done.
return 0;
}
4 changes: 2 additions & 2 deletions packages/amesos2/src/Amesos2_CssMKL_def.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -264,7 +264,7 @@ namespace Amesos2 {
}

return( 0 );
}
}


template <class Matrix, class Vector>
Expand Down Expand Up @@ -470,7 +470,7 @@ CssMKL<Matrix,Vector>::loadA_impl(EPhase current_phase)
#if 1
// Only reinex GIDs
css_rowmap_ = this->matrixA_->getRowMap(); // use original map to redistribute vectors in solve
Teuchos::RCP<const MatrixAdapter<Matrix> > contig_mat = this->matrixA_->reindex(css_contig_rowmap_, css_contig_colmap_);
Teuchos::RCP<const MatrixAdapter<Matrix> > contig_mat = this->matrixA_->reindex(css_contig_rowmap_, css_contig_colmap_, current_phase);
#else
// Redistribued matrixA into contiguous GIDs
Teuchos::RCP<const MatrixAdapter<Matrix> > contig_mat = this->matrixA_->get(ptrInArg(*css_rowmap_));
Expand Down
Loading
Loading