Skip to content
Permalink
Browse files

Merge #3075

3075: Store non-bonded interaction parameters only once per pair r=KaiSzuttor a=fweik

Description of changes:
 - Index function for upper triangular matrices
 - Store non-bonded interaction parameters only once per pair


Co-authored-by: Florian Weik <fweik@icp.uni-stuttgart.de>
  • Loading branch information...
bors and fweik committed Aug 12, 2019
2 parents b94daa3 + d284e1f commit 3c0d86badac7b61f07abc2ca6f53982b27ddfe82
@@ -340,8 +340,6 @@ void mpi_bcast_ia_params(int i, int j) {
if (j >= 0) {
/* non-bonded interaction parameters */
boost::mpi::broadcast(comm_cart, *get_ia_param(i, j), 0);

*get_ia_param(j, i) = *get_ia_param(i, j);
} else {
/* bonded interaction parameters */
MPI_Bcast(&(bonded_ia_params[i]), sizeof(Bonded_ia_parameters), MPI_BYTE, 0,
@@ -361,9 +359,6 @@ void mpi_bcast_ia_params_slave(int i, int j) {
if (j >= 0) { /* non-bonded interaction parameters */

boost::mpi::broadcast(comm_cart, *get_ia_param(i, j), 0);

*get_ia_param(j, i) = *get_ia_param(i, j);

} else { /* bonded interaction parameters */
make_bond_type_exist(i); /* realloc bonded_ia_params on slave nodes! */
MPI_Bcast(&(bonded_ia_params[i]), sizeof(Bonded_ia_parameters), MPI_BYTE, 0,
@@ -22,38 +22,11 @@
* Implementation of nonbonded_interaction_data.hpp
*/
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "actor/DipolarDirectSum.hpp"
#include "bonded_interactions/bonded_interaction_data.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include "dpd.hpp"
#include "electrostatics_magnetostatics/magnetic_non_p3m_methods.hpp"
#include "electrostatics_magnetostatics/mdlc_correction.hpp"
#include "errorhandling.hpp"
#include "event.hpp"
#include "grid.hpp"
#include "nonbonded_interaction_data.hpp"
#include "nonbonded_interactions/buckingham.hpp"
#include "nonbonded_interactions/gaussian.hpp"
#include "nonbonded_interactions/gay_berne.hpp"
#include "nonbonded_interactions/hat.hpp"
#include "nonbonded_interactions/hertzian.hpp"
#include "nonbonded_interactions/lj.hpp"
#include "nonbonded_interactions/ljcos.hpp"
#include "nonbonded_interactions/ljcos2.hpp"
#include "nonbonded_interactions/ljgen.hpp"
#include "nonbonded_interactions/morse.hpp"
#include "nonbonded_interactions/nonbonded_tab.hpp"
#include "nonbonded_interactions/smooth_step.hpp"
#include "nonbonded_interactions/soft_sphere.hpp"
#ifdef DIPOLAR_BARNES_HUT
#include "actor/DipolarBarnesHut.hpp"
#endif
#include "layered.hpp"
#include "object-in-fluid/affinity.hpp"
#include "object-in-fluid/membrane_collision.hpp"
#include "pressure.hpp"
#include "rattle.hpp"

#include "serialization/IA_parameters.hpp"

#include <boost/archive/binary_iarchive.hpp>
@@ -64,13 +37,8 @@
#include <boost/serialization/string.hpp>
#include <boost/serialization/vector.hpp>

#include <cstdlib>
#include <cstring>

#include "electrostatics_magnetostatics/coulomb.hpp"
#include "electrostatics_magnetostatics/debye_hueckel.hpp"
#include "electrostatics_magnetostatics/dipole.hpp"
#include "electrostatics_magnetostatics/p3m.hpp"

/****************************************
* variables
@@ -252,13 +220,13 @@ void realloc_ia_params(int nsize) {
if (nsize <= max_seen_particle_type)
return;

auto new_params = std::vector<IA_parameters>(nsize * nsize);
auto new_params = std::vector<IA_parameters>(nsize * (nsize + 1) / 2);

/* if there is an old field, move entries */
for (int i = 0; i < max_seen_particle_type; i++)
for (int j = 0; j < max_seen_particle_type; j++) {
new_params[i * nsize + j] =
std::move(ia_params[i * max_seen_particle_type + j]);
for (int j = i; j < max_seen_particle_type; j++) {
new_params.at(Utils::upper_triangular(i, j, nsize)) =
std::move(*get_ia_param(i, j));
}

max_seen_particle_type = nsize;
@@ -28,6 +28,7 @@
#include "dpd.hpp"
#include "particle_data.hpp"

#include <utils/index.hpp>
#include <utils/math/sqr.hpp>

/** Cutoff for deactivated interactions. Must be negative, so that even
@@ -322,11 +323,23 @@ extern double min_global_cut;
* exported functions
************************************************/

/** Get interaction parameters between particle sorts i and j */
/**
* @brief Get interaction parameters between particle types i and j
*
* This is symmetric, e.g. it holds that get_ia_param(i, j) and
* get_ia_param(j, i) point to the same data.
*
* @param i First type, has to be be smaller than @ref max_seen_particle_type.
* @param j Second type, has to be be smaller than @ref max_seen_particle_type.
*
* @return Pointer to interaction parameters for the type pair.
* */
inline IA_parameters *get_ia_param(int i, int j) {
extern std::vector<IA_parameters> ia_params;
extern int max_seen_particle_type;
return &ia_params[i * max_seen_particle_type + j];
assert(i >= 0 && i < max_seen_particle_type);
assert(j >= 0 && j < max_seen_particle_type);

return &ia_params[Utils::upper_triangular(std::min(i, j), std::max(i, j),
max_seen_particle_type)];
}

/** Get interaction parameters between particle sorts i and j.
@@ -96,6 +96,27 @@ get_linear_index(const Vector3i &ind, const Vector3i &adim,
return get_linear_index(ind[0], ind[1], ind[2], adim, memory_order);
}

/**
* @brief Linear index into a upper triangular matrix.
*
* This is row-major.
*
* @tparam T Integral
* @param i row index
* @param j column index
* @param n matrix size
* @return linear index
*/
template <class T> T upper_triangular(T i, T j, T n) {
/* n is a valid size */
assert(n >= 0);
/* i is a valid row index */
assert((i >= 0) && (i < n));
/* j is in the upper triangle */
assert((j >= i) && (j < n));
return (n * (n - 1)) / 2 - ((n - i) * (n - i - 1)) / 2 + j;
}

/*@}*/

} // namespace Utils
@@ -21,11 +21,13 @@

#define BOOST_TEST_MODULE Utils::ravel_index test
#define BOOST_TEST_DYN_LINK

#include <utils/Vector.hpp>
#include <utils/index.hpp>

#include <array>
#include <boost/test/unit_test.hpp>

#include "utils/index.hpp"

BOOST_AUTO_TEST_CASE(ravel_index_test) {
const std::array<std::size_t, 4> unravelled_indices{{12, 23, 5, 51}};
const std::array<std::size_t, 4> dimensions{{15, 35, 6, 52}};
@@ -68,4 +70,24 @@ BOOST_AUTO_TEST_CASE(get_linear_index) {
BOOST_CHECK_EQUAL(get_linear_index(index[0], index[1], index[2], grid_size),
get_linear_index(index, grid_size));
}
}

BOOST_AUTO_TEST_CASE(upper_triangular_test) {
// clang-format off
const Utils::VectorXi<2> A[] =
{
{0, 0}, {0, 1}, {0, 2},
{1, 1}, {1, 2},
{2, 2}
};
// clang-format on

for (int i = 0; i < 3; i++) {
for (int j = i; j < 3; j++) {
auto const expected = Utils::VectorXi<2>{i, j};
auto const result = A[Utils::upper_triangular(i, j, 3)];

BOOST_CHECK(expected == result);
}
}
}

0 comments on commit 3c0d86b

Please sign in to comment.
You can’t perform that action at this time.