Skip to content

[core] Test GPGPU (OpenCL/CUDA) backend with AMGCL solver #3831

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.

Sign up for GitHub

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

Merged
ddemidov merged 6 commits into from
Feb 13, 2019
Merged

[core] Test GPGPU (OpenCL/CUDA) backend with AMGCL solver #3831

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
4b907c5
[wip] [skip ci] Test OpenCL backend with amgcl solver
ddemidov Jan 15, 2019
1a576a9
See if cotire works this way
ddemidov Jan 23, 2019
9b83d1a
Patch COTIRE to support cmake generator expressions
ddemidov Jan 23, 2019
271886c
Add vexcl source to external_libraries
ddemidov Feb 8, 2019
a034674
Patch vexcl cmake to resolve cotire issue
ddemidov Feb 8, 2019
22f0b38
Get rid of explicit template instantiation
ddemidov Feb 13, 2019
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 applications/trilinos_application/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ set( KRATOS_TRILINOS_APPLICATION_SOURCES
${CMAKE_CURRENT_SOURCE_DIR}/custom_factories/trilinos_linear_solver_factory.cpp
${CMAKE_SOURCE_DIR}/applications/FluidDynamicsApplication/fluid_dynamics_application_variables.cpp #TODO: this should REALLY NOT BE HERE
${CMAKE_CURRENT_SOURCE_DIR}/custom_utilities/mpi_normal_calculation_utilities.cpp;
${CMAKE_CURRENT_SOURCE_DIR}/amgcl_mpi_solver_impl.cpp;
)

## Kratos tests sources. Enabled by default
Expand Down
216 changes: 216 additions & 0 deletions applications/trilinos_application/amgcl_mpi_solver_impl.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,216 @@
#ifndef KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
#define KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H

#include <boost/range/iterator_range.hpp>
#include <boost/property_tree/ptree.hpp>

#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/adapter/epetra.hpp>
#include <amgcl/adapter/ublas.hpp>
#include <amgcl/adapter/zero_copy.hpp>
#include <amgcl/adapter/block_matrix.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/value_type/static_matrix.hpp>
#include <amgcl/solver/runtime.hpp>

#include <amgcl/mpi/util.hpp>
#include <amgcl/mpi/make_solver.hpp>
#include <amgcl/mpi/preconditioner.hpp>

#ifdef AMGCL_GPGPU
# include <amgcl/backend/vexcl.hpp>
# include <amgcl/backend/vexcl_static_matrix.hpp>
#endif

#include "Epetra_FECrsMatrix.h"
#include "Epetra_FEVector.h"
#include "trilinos_space.h"

#include "external_includes/amgcl_mpi_solve_functions.h"

namespace Kratos
{

#ifdef AMGCL_GPGPU
vex::Context& vexcl_context();

template <int TBlockSize>
void register_vexcl_static_matrix_type();
#endif

// Spacialization of AMGCLScalarSolve for distribued systems.
void AMGCLScalarSolve(
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType& rA,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber,
double& rResidual,
const boost::property_tree::ptree &amgclParams,
int verbosity_level,
bool use_gpgpu
)
{
#ifdef AMGCL_GPGPU
if (use_gpgpu && vexcl_context()) {
auto &ctx = vexcl_context();

typedef amgcl::backend::vexcl<double> Backend;

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

Backend::params bprm;
bprm.q = ctx;

Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams, bprm);

std::size_t n = rA.NumMyRows();

vex::vector<double> b(ctx, n, rB.Values());
vex::vector<double> x(ctx, n, rX.Values());

std::tie(rIterationNumber, rResidual) = solve(b, x);

vex::copy(x.begin(), x.end(), rX.Values());
} else
#endif
{
typedef amgcl::backend::builtin<double> Backend;

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams);

std::size_t n = rA.NumMyRows();

auto b_range = boost::make_iterator_range(rB.Values(), rB.Values() + n);
auto x_range = boost::make_iterator_range(rX.Values(), rX.Values() + n);

std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
}
}

// Spacialization of AMGCLBlockSolve for distribued systems.
template <int TBlockSize>
void AMGCLBlockSolve(
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType & rA,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber,
double& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level,
bool use_gpgpu
)
{
if(amgclParams.get<std::string>("precond.class") != "amg")
amgclParams.erase("precond.coarsening");
else
amgclParams.put("precond.coarsening.aggr.block_size",1);

typedef amgcl::static_matrix<double, TBlockSize, TBlockSize> val_type;
typedef amgcl::static_matrix<double, TBlockSize, 1> rhs_type;

std::size_t n = rA.RowMap().NumMyElements();
std::size_t nb = n / TBlockSize;

#ifdef AMGCL_GPGPU
if (use_gpgpu && vexcl_context()) {
auto &ctx = vexcl_context();
register_vexcl_static_matrix_type<TBlockSize>();

typedef amgcl::backend::vexcl<val_type> Backend;

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

typename Backend::params bprm;
bprm.q = ctx;

Solver solve(
MPI_COMM_WORLD,
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
amgclParams, bprm
);

auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());

vex::vector<rhs_type> x(ctx, nb, x_begin);
vex::vector<rhs_type> b(ctx, nb, b_begin);

std::tie(rIterationNumber, rResidual) = solve(b, x);

vex::copy(x.begin(), x.end(), x_begin);
} else
#endif
{
typedef amgcl::backend::builtin<val_type> Backend;

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

Solver solve(
MPI_COMM_WORLD,
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
amgclParams
);

auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());

auto b_range = boost::make_iterator_range(b_begin, b_begin + nb);
auto x_range = boost::make_iterator_range(x_begin, x_begin + nb);

std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
}
}

void AMGCLSolve(
int block_size,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType& rA,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber,
double& rResidual,
const boost::property_tree::ptree &amgclParams,
int verbosity_level,
bool use_gpgpu
)
{
switch (block_size) {
case 2:
AMGCLBlockSolve<2>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
case 3:
AMGCLBlockSolve<3>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
case 4:
AMGCLBlockSolve<4>(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
default:
AMGCLScalarSolve(rA, rX, rB, rIterationNumber, rResidual, amgclParams, verbosity_level, use_gpgpu);
return;
}
}

} // namespace Kratos

#endif
99 changes: 10 additions & 89 deletions applications/trilinos_application/external_includes/amgcl_mpi_solve_functions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,101 +1,22 @@
#ifndef KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H
#define KRATOS_AMGCL_MPI_SOLVE_FUNCTIONS_H

#include <boost/range/iterator_range.hpp>

#include <amgcl/adapter/crs_tuple.hpp>
#include <amgcl/adapter/epetra.hpp>
#include <amgcl/adapter/ublas.hpp>
#include <amgcl/adapter/zero_copy.hpp>
#include <amgcl/adapter/block_matrix.hpp>
#include <amgcl/backend/builtin.hpp>
#include <amgcl/value_type/static_matrix.hpp>
#include <amgcl/solver/runtime.hpp>

#include <amgcl/mpi/util.hpp>
#include <amgcl/mpi/make_solver.hpp>
#include <amgcl/mpi/preconditioner.hpp>
#include <boost/property_tree/ptree.hpp>

namespace Kratos
{

// Spacialization of AMGCLScalarSolve for distribued systems.
template <class TSparseSpaceType>
typename std::enable_if<TSparseSpaceType::IsDistributed(), void>::type
AMGCLScalarSolve(
typename TSparseSpaceType::MatrixType& rA,
typename TSparseSpaceType::VectorType& rX,
typename TSparseSpaceType::VectorType& rB,
typename TSparseSpaceType::IndexType& rIterationNumber,
void AMGCLSolve(
int block_size,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::MatrixType& rA,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rX,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::VectorType& rB,
TrilinosSpace<Epetra_FECrsMatrix, Epetra_FEVector>::IndexType& rIterationNumber,
double& rResidual,
const boost::property_tree::ptree &amgclParams,
int verbosity_level
)
{
typedef amgcl::backend::builtin<double> Backend;

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

Solver solve(MPI_COMM_WORLD, amgcl::adapter::map(rA), amgclParams);

std::size_t n = rA.NumMyRows();

auto b_range = boost::make_iterator_range(rB.Values(), rB.Values() + n);
auto x_range = boost::make_iterator_range(rX.Values(), rX.Values() + n);

std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
}

// Spacialization of AMGCLBlockSolve for distribued systems.
template <int TBlockSize, class TSparseSpaceType>
typename std::enable_if<TSparseSpaceType::IsDistributed(), void>::type
AMGCLBlockSolve(
typename TSparseSpaceType::MatrixType & rA,
typename TSparseSpaceType::VectorType& rX,
typename TSparseSpaceType::VectorType& rB,
typename TSparseSpaceType::IndexType& rIterationNumber,
double& rResidual,
boost::property_tree::ptree amgclParams,
int verbosity_level
)
{
if(amgclParams.get<std::string>("precond.class") != "amg")
amgclParams.erase("precond.coarsening");
else
amgclParams.put("precond.coarsening.aggr.block_size",1);

typedef amgcl::static_matrix<double, TBlockSize, TBlockSize> val_type;
typedef amgcl::static_matrix<double, TBlockSize, 1> rhs_type;
typedef amgcl::backend::builtin<val_type> Backend;

std::size_t n = rA.RowMap().NumMyElements();

typedef
amgcl::mpi::make_solver<
amgcl::runtime::mpi::preconditioner<Backend>,
amgcl::runtime::solver::wrapper
>
Solver;

Solver solve(
MPI_COMM_WORLD,
amgcl::adapter::block_matrix<val_type>(amgcl::adapter::map(rA)),
amgclParams
);

auto b_begin = reinterpret_cast<const rhs_type*>(rB.Values());
auto x_begin = reinterpret_cast<rhs_type*>(rX.Values());

auto b_range = boost::make_iterator_range(b_begin, b_begin + n / TBlockSize);
auto x_range = boost::make_iterator_range(x_begin, x_begin + n / TBlockSize);

std::tie(rIterationNumber, rResidual) = solve(b_range, x_range);
}
int verbosity_level,
bool use_gpgpu
);

} // namespace Kratos

Expand Down
Loading