Skip to content

Commit

Permalink
Amesos2 : add Galeri to SimpleSolve example
Browse files Browse the repository at this point in the history 
Signed-off-by: iyamazaki <[email protected]>
  • Loading branch information
@iyamazaki
iyamazaki committed Jan 17, 2025
1 parent c517a5f commit 65a3e9b
Show file tree
Hide file tree
Showing 3 changed files with 129 additions and 73 deletions.
1 change: 1 addition & 0 deletions packages/amesos2/cmake/Amesos2_config.h.in
View file
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
View file
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
View file
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;
}

0 comments on commit 65a3e9b

Please sign in to comment.