diff --git a/Makefile b/Makefile
index 47a07f46cb23371ccecbe5ef49379adbb5842565..ba82842577390602a55904c2cb8c71a8863909cc 100644
--- a/Makefile
+++ b/Makefile
@@ -18,7 +18,7 @@ lj_ns.cu:
# Targets
lj_cpu: lj_ns.cpp
# g++ -o lj_cpu lj_ns.cpp -DDEBUG
- mpic++ -O3 -o lj_cpu lj_ns.cpp -DDEBUG
+ mpic++ -O3 -o lj_cpu lj_ns.cpp runtime/pairs.cpp runtime/devices/dummy.cpp -DDEBUG
lj_gpu: lj_ns.cu
nvcc -o lj_gpu lj_ns.cu
diff --git a/runtime/array.hpp b/runtime/array.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a755c4b07a70e4837ec02cbf1786071a53e3ea1
--- /dev/null
+++ b/runtime/array.hpp
@@ -0,0 +1,37 @@
+#include "pairs_common.hpp"
+
+#pragma once
+
+namespace pairs {
+
+class Array {
+protected:
+ array_t id;
+ std::string name;
+ void *h_ptr, *d_ptr;
+ size_t size;
+ bool is_static;
+
+public:
+ Array(array_t id_, std::string name_, void *h_ptr_, void *d_ptr_, size_t size_, bool is_static_ = false) :
+ id(id_),
+ name(name_),
+ h_ptr(h_ptr_),
+ d_ptr(d_ptr_),
+ size(size_),
+ is_static(is_static_) {
+
+ PAIRS_ASSERT(size_ > 0);
+ }
+
+ property_t getId() { return id; }
+ std::string getName() { return name; }
+ void *getHostPointer() { return h_ptr; }
+ void *getDevicePointer() { return d_ptr; }
+ void setPointers(void *h_ptr_, void *d_ptr_) { h_ptr = h_ptr_, d_ptr = d_ptr_; }
+ void setSize(size_t size_) { size = size_; }
+ size_t getSize() { return size; };
+ bool isStatic() { return is_static; }
+};
+
+}
diff --git a/runtime/device_flags.hpp b/runtime/device_flags.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4b5085fb6b61f8df1c3e1e129541266765745fdf
--- /dev/null
+++ b/runtime/device_flags.hpp
@@ -0,0 +1,95 @@
+#include <cmath>
+#include <iostream>
+
+#pragma once
+
+namespace pairs {
+
+class DeviceFlags {
+private:
+ unsigned long long int *hflags;
+ unsigned long long int *dflags;
+ int narrays;
+ int nflags;
+ static const int narrays_per_flag = 64;
+public:
+ DeviceFlags(int narrays_) : narrays(narrays_) {
+ nflags = std::ceil((double) narrays_ / (double) narrays_per_flag);
+ hflags = new unsigned long long int[nflags];
+ dflags = new unsigned long long int[nflags];
+
+ for(int i = 0; i < nflags; ++i) {
+ hflags[i] = 0xfffffffffffffffful;
+ dflags[i] = 0x0;
+ }
+ }
+
+ inline bool isHostFlagSet(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ return (hflags[flag_index] & (1 << bit)) != 0;
+ }
+
+ inline void setHostFlag(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ hflags[flag_index] |= 1 << bit;
+ }
+
+ inline void clearHostFlag(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ hflags[flag_index] &= ~(1 << bit);
+ }
+
+ inline bool isDeviceFlagSet(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ return (dflags[flag_index] & (1 << bit)) != 0;
+ }
+
+ inline void setDeviceFlag(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ dflags[flag_index] |= 1 << bit;
+ }
+
+ inline void clearDeviceFlag(int array_id) {
+ int flag_index = array_id / narrays_per_flag;
+ unsigned long long int bit = array_id % narrays_per_flag;
+ dflags[flag_index] &= ~(1 << bit);
+ }
+
+ void printFlags() {
+ std::cerr << "hflags = ";
+ for(int i = 0; i < nflags; i++) {
+ for(int b = 63; b >= 0; b--) {
+ if(hflags[i] & (1 << b)) {
+ std::cerr << "1";
+ } else {
+ std::cerr << "0";
+ }
+ }
+ }
+
+ std::cerr << std::endl << "dflags = ";
+ for(int i = 0; i < nflags; i++) {
+ for(int b = 63; b >= 0; b--) {
+ if(dflags[i] & (1 << b)) {
+ std::cerr << "1";
+ } else {
+ std::cerr << "0";
+ }
+ }
+ }
+
+ std::cerr << std::endl;
+ }
+
+ ~DeviceFlags() {
+ delete[] hflags;
+ delete[] dflags;
+ }
+};
+
+}
diff --git a/runtime/devices/cuda.hpp b/runtime/devices/cuda.cu
similarity index 88%
rename from runtime/devices/cuda.hpp
rename to runtime/devices/cuda.cu
index 60c5ceba8b5a08d1c344fd9b232466d17861e39a..50b35006f2bbf70be3879abd539eaa882e66f4c2 100644
--- a/runtime/devices/cuda.hpp
+++ b/runtime/devices/cuda.cu
@@ -43,8 +43,8 @@ __host__ void copy_static_symbol_to_host(void *d_ptr, const void *h_ptr, size_t
//CUDA_ASSERT(cudaMemcpyFromSymbol(h_ptr, d_ptr, count));
}
-inline __device__ int atomic_add(int *addr, int val) { return atomicAdd(addr, val); }
-inline __device__ int atomic_add_resize_check(int *addr, int val, int *resize, int capacity) {
+__device__ int atomic_add(int *addr, int val) { return atomicAdd(addr, val); }
+__device__ int atomic_add_resize_check(int *addr, int val, int *resize, int capacity) {
const int add_res = *addr + val;
if(add_res >= capacity) {
*resize = add_res;
diff --git a/runtime/devices/device.hpp b/runtime/devices/device.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..30b1d1e44e60b4bbfa30171f5d50fc510eae5891
--- /dev/null
+++ b/runtime/devices/device.hpp
@@ -0,0 +1,24 @@
+#include <iostream>
+#include <cstddef>
+
+#pragma once
+
+#ifndef PAIRS_TARGET_CUDA
+# define __host__
+# define __device__
+typedef int cudaError_t;
+#endif
+
+namespace pairs {
+
+void cuda_assert(cudaError_t err, const char *file, int line);
+__host__ void *device_alloc(size_t size);
+__host__ void *device_realloc(void *ptr, size_t size);
+__host__ void copy_to_device(const void *h_ptr, void *d_ptr, size_t count);
+__host__ void copy_to_host(const void *d_ptr, void *h_ptr, size_t count);
+__host__ void copy_static_symbol_to_device(void *h_ptr, const void *d_ptr, size_t count);
+__host__ void copy_static_symbol_to_host(void *d_ptr, const void *h_ptr, size_t count);
+__device__ int atomic_add(int *addr, int val);
+__device__ int atomic_add_resize_check(int *addr, int val, int *resize, int capacity);
+
+}
diff --git a/runtime/devices/dummy.hpp b/runtime/devices/dummy.cpp
similarity index 65%
rename from runtime/devices/dummy.hpp
rename to runtime/devices/dummy.cpp
index 3daebcf5804046b21c12e700e88f27b34dabd13a..8b8c80a85027c2c4a861dfce9e0cb57d66175b65 100644
--- a/runtime/devices/dummy.hpp
+++ b/runtime/devices/dummy.cpp
@@ -1,11 +1,11 @@
-#pragma once
+#include "device.hpp"
namespace pairs {
-void *device_alloc(size_t size) { return NULL; }
-void *device_realloc(void *ptr, size_t size) { return NULL; }
-void copy_to_device(void *h_ptr, const void *d_ptr, size_t count) {}
-void copy_to_host(void *d_ptr, const void *h_ptr, size_t count) {}
+void *device_alloc(size_t size) { return nullptr; }
+void *device_realloc(void *ptr, size_t size) { return nullptr; }
+void copy_to_device(void const *h_ptr, void *d_ptr, size_t count) {}
+void copy_to_host(void const *d_ptr, void *h_ptr, size_t count) {}
void copy_static_symbol_to_device(void *h_ptr, const void *d_ptr, size_t count) {}
void copy_static_symbol_to_host(void *d_ptr, const void *h_ptr, size_t count) {}
int atomic_add(int *addr, int val) {
diff --git a/runtime/domain/domain_partitioning.hpp b/runtime/domain/domain_partitioning.hpp
index 5e628f867cd3de109ca70397e44fbbbb6e1c7a69..c4a231b3e14c5ab6ea911f1b830a348399f55677 100644
--- a/runtime/domain/domain_partitioning.hpp
+++ b/runtime/domain/domain_partitioning.hpp
@@ -1,37 +1,45 @@
#include <mpi.h>
+//---
+#include "../pairs_common.hpp"
#pragma once
-typedef double real_t;
-
namespace pairs {
-template<int ndims> class Regular6DStencil;
+class Regular6DStencil;
-template<int ndims> class DomainPartitioner {
- friend class Regular6DStencil<ndims>;
+class DomainPartitioner {
+ friend class Regular6DStencil;
protected:
- real_t grid_min[ndims];
- real_t grid_max[ndims];
-
+ real_t *grid_min;
+ real_t *grid_max;
+ int ndims;
+
public:
DomainPartitioner(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
- static_assert(ndims == 3, "DomainPartitioner(): number of dimensions mismatching!");
- //PAIRS_ASSERT(xmax > xmin);
- //PAIRS_ASSERT(ymax > ymin);
- //PAIRS_ASSERT(zmax > zmin);
+ PAIRS_ASSERT(xmax > xmin);
+ PAIRS_ASSERT(ymax > ymin);
+ PAIRS_ASSERT(zmax > zmin);
+ ndims = 3;
+ grid_min = new real_t[ndims];
+ grid_max = new real_t[ndims];
grid_min[0] = xmin;
- grid_max[0] = xmax;
grid_min[1] = ymin;
- grid_max[1] = ymax;
grid_min[2] = zmin;
+ grid_max[0] = xmax;
+ grid_max[1] = ymax;
grid_max[2] = zmax;
}
+ ~DomainPartitioner() {
+ delete[] grid_min;
+ delete[] grid_max;
+ }
+
virtual void initialize(int *argc, char ***argv) = 0;
- virtual void fillArrays(int neighbor_ranks[], int pbc[], real_t subdom[]) = 0;
+ virtual void fillArrays(int *neighbor_ranks, int *pbc, real_t *subdom) = 0;
virtual void communicateSizes(int dim, const int *nsend, int *nrecv) = 0;
virtual void communicateData(
int dim, int elem_size,
@@ -40,6 +48,4 @@ public:
virtual void finalize() = 0;
};
-//class DomainPartitioner<3>;
-
}
diff --git a/runtime/domain/no_partition.hpp b/runtime/domain/no_partition.hpp
index 51eaa2f3737424a91932216d41c8c213f1db0e0c..b9ac35748f5a4c9387c0853f722c141f687ecadc 100644
--- a/runtime/domain/no_partition.hpp
+++ b/runtime/domain/no_partition.hpp
@@ -6,8 +6,7 @@
namespace pairs {
-template <int ndims>
-class NoPartitioning : DimensionRanges<ndims> {
+class NoPartitioning : DimensionRanges {
public:
void initialize(int *argc, const char **argv) {
for(int d = 0; d < ndims; d++) {
diff --git a/runtime/domain/regular_6d_stencil.hpp b/runtime/domain/regular_6d_stencil.hpp
index 4542f2fe80fc32469fd171c23456524be8ac2a98..631ded149ca38869db7f9c797dbd61cc6aa80075 100644
--- a/runtime/domain/regular_6d_stencil.hpp
+++ b/runtime/domain/regular_6d_stencil.hpp
@@ -1,32 +1,49 @@
+#include "../pairs_common.hpp"
#include "domain_partitioning.hpp"
#pragma once
#define SMALL 0.00001
-typedef double real_t;
-
namespace pairs {
-template <int ndims>
-class Regular6DStencil : public DomainPartitioner<ndims> {
+class Regular6DStencil : public DomainPartitioner {
private:
int world_size, rank;
- int nranks[ndims];
- int prev[ndims];
- int next[ndims];
- int pbc_prev[ndims];
- int pbc_next[ndims];
- real_t subdom_min[ndims];
- real_t subdom_max[ndims];
+ int *nranks;
+ int *prev;
+ int *next;
+ int *pbc_prev;
+ int *pbc_next;
+ real_t *subdom_min;
+ real_t *subdom_max;
public:
Regular6DStencil(real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) :
- DomainPartitioner<ndims>(xmin, xmax, ymin, ymax, zmin, zmax) {}
+ DomainPartitioner(xmin, xmax, ymin, ymax, zmin, zmax) {
+
+ nranks = new int[ndims];
+ prev = new int[ndims];
+ next = new int[ndims];
+ pbc_prev = new int[ndims];
+ pbc_next = new int[ndims];
+ subdom_min = new real_t[ndims];
+ subdom_max = new real_t[ndims];
+ }
+
+ ~Regular6DStencil() {
+ delete[] nranks;
+ delete[] prev;
+ delete[] next;
+ delete[] pbc_prev;
+ delete[] pbc_next;
+ delete[] subdom_min;
+ delete[] subdom_max;
+ }
void setConfig() {
- static_assert(ndims == 3, "setConfig() only implemented for three dimensions!");
- real_t area[ndims];
+ PAIRS_ASSERT(ndims == 3);
+ real_t area[3];
real_t best_surf = 0.0;
int d = 0;
@@ -60,9 +77,9 @@ public:
void setBoundingBox() {
MPI_Comm cartesian;
- int myloc[ndims];
- int periods[ndims];
- real_t rank_length[ndims];
+ int *myloc = new int[ndims];
+ int *periods = new int[ndims];
+ real_t *rank_length = new real_t[ndims];
int reorder = 0;
for(int d = 0; d < ndims; d++) {
@@ -80,6 +97,9 @@ public:
subdom_max[d] = subdom_min[d] + rank_length[d];
}
+ delete[] myloc;
+ delete[] periods;
+ delete[] rank_length;
MPI_Comm_free(&cartesian);
}
@@ -103,7 +123,7 @@ public:
z >= subdom_min[2] && z < subdom_max[2] - SMALL;
}
- void fillArrays(int neighbor_ranks[], int pbc[], real_t subdom[]) {
+ void fillArrays(int *neighbor_ranks, int *pbc, real_t *subdom) {
for(int d = 0; d < ndims; d++) {
neighbor_ranks[d * 2 + 0] = prev[d];
neighbor_ranks[d * 2 + 1] = next[d];
diff --git a/runtime/pairs.cpp b/runtime/pairs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..2fcb2d60cadeffe7f722a75e7a86474ef9a759a9
--- /dev/null
+++ b/runtime/pairs.cpp
@@ -0,0 +1,172 @@
+#include <algorithm>
+#include <cmath>
+#include <iostream>
+#include <memory>
+#include <vector>
+//---
+#include "pairs.hpp"
+#include "pairs_common.hpp"
+#include "devices/device.hpp"
+#include "domain/regular_6d_stencil.hpp"
+
+namespace pairs {
+
+void PairsSimulation::initDomain(int *argc, char ***argv, real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
+ if(dom_part_type == DimRanges) {
+ dom_part = new Regular6DStencil(xmin, xmax, ymin, ymax, zmin, zmax);
+ } else {
+ PAIRS_EXCEPTION("Domain partitioning type not implemented!\n");
+ }
+
+ dom_part->initialize(argc, argv);
+}
+
+void PairsSimulation::addArray(Array array) {
+ int id = array.getId();
+ auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+ PAIRS_ASSERT(a == std::end(arrays));
+ arrays.push_back(array);
+}
+
+Array &PairsSimulation::getArray(array_t id) {
+ auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+ PAIRS_ASSERT(a != std::end(arrays));
+ return *a;
+}
+
+Array &PairsSimulation::getArrayByName(std::string name) {
+ auto a = std::find_if(arrays.begin(), arrays.end(), [name](Array a) { return a.getName() == name; });
+ PAIRS_ASSERT(a != std::end(arrays));
+ return *a;
+}
+
+void PairsSimulation::addProperty(Property prop) {
+ int id = prop.getId();
+ auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
+ PAIRS_ASSERT(p == std::end(properties));
+ properties.push_back(prop);
+}
+
+Property &PairsSimulation::getProperty(property_t id) {
+ auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
+ PAIRS_ASSERT(p != std::end(properties));
+ return *p;
+}
+
+Property &PairsSimulation::getPropertyByName(std::string name) {
+ auto p = std::find_if(properties.begin(), properties.end(), [name](Property p) { return p.getName() == name; });
+ PAIRS_ASSERT(p != std::end(properties));
+ return *p;
+}
+
+void PairsSimulation::copyArrayToDevice(Array &array) {
+ int array_id = array.getId();
+ if(!array_flags->isDeviceFlagSet(array_id)) {
+ if(array.isStatic()) {
+ PAIRS_DEBUG("Copying static array %s to device\n", array.getName().c_str());
+ pairs::copy_static_symbol_to_device(array.getHostPointer(), array.getDevicePointer(), array.getSize());
+ } else {
+ PAIRS_DEBUG("Copying array %s to device\n", array.getName().c_str());
+ pairs::copy_to_device(array.getHostPointer(), array.getDevicePointer(), array.getSize());
+ }
+ }
+}
+
+void PairsSimulation::copyArrayToHost(Array &array) {
+ int array_id = array.getId();
+ if(!array_flags->isHostFlagSet(array_id)) {
+ if(array.isStatic()) {
+ PAIRS_DEBUG("Copying static array %s to host\n", array.getName().c_str());
+ pairs::copy_static_symbol_to_host(array.getDevicePointer(), array.getHostPointer(), array.getSize());
+ } else {
+ PAIRS_DEBUG("Copying array %s to host\n", array.getName().c_str());
+ pairs::copy_to_host(array.getDevicePointer(), array.getHostPointer(), array.getSize());
+ }
+ }
+}
+
+void PairsSimulation::copyPropertyToDevice(Property &prop) {
+ if(!prop_flags->isDeviceFlagSet(prop.getId())) {
+ PAIRS_DEBUG("Copying property %s to device\n", prop.getName().c_str());
+ pairs::copy_to_device(prop.getHostPointer(), prop.getDevicePointer(), prop.getTotalSize());
+ }
+}
+
+void PairsSimulation::copyPropertyToHost(Property &prop) {
+ if(!prop_flags->isHostFlagSet(prop.getId())) {
+ PAIRS_DEBUG("Copying property %s to host\n", prop.getName().c_str());
+ pairs::copy_to_host(prop.getDevicePointer(), prop.getHostPointer(), prop.getTotalSize());
+ }
+}
+
+void PairsSimulation::communicateSizes(int dim, const int *send_sizes, int *recv_sizes) {
+ this->getDomainPartitioner()->communicateSizes(dim, send_sizes, recv_sizes);
+ PAIRS_DEBUG("send_sizes=[%d, %d], recv_sizes=[%d, %d]\n", send_sizes[dim * 2 + 0], send_sizes[dim * 2 + 1], recv_sizes[dim * 2 + 0], recv_sizes[dim * 2 + 1]);
+}
+
+void PairsSimulation::communicateData(
+ int dim, int elem_size,
+ const real_t *send_buf, const int *send_offsets, const int *nsend,
+ real_t *recv_buf, const int *recv_offsets, const int *nrecv) {
+
+ this->getDomainPartitioner()->communicateData(dim, elem_size, send_buf, send_offsets, nsend, recv_buf, recv_offsets, nrecv);
+
+ /*
+ // Debug messages
+ const int elems_to_print = 5;
+
+ // Send buffer debug
+ for(int i = 0; i < 2; i++) {
+ int nsnd = nsend[dim * 2 + i];
+
+ PAIRS_DEBUG("send_buf=[");
+ for(int j = 0; j < MIN(elems_to_print, nsnd); j++) {
+ for(int k = 0; k < elem_size; k++) {
+ PAIRS_DEBUG("%f,", send_buf[(send_offsets[dim * 2 + i] + j) * elem_size + k]);
+ }
+ }
+
+ if(elems_to_print * 2 < nsnd) {
+ PAIRS_DEBUG("\b ... ");
+ }
+
+ for(int j = MAX(elems_to_print, nsnd - elems_to_print); j < nsnd; j++) {
+ for(int k = 0; k < elem_size; k++) {
+ PAIRS_DEBUG("%f,", send_buf[(send_offsets[dim * 2 + i] + j) * elem_size + k]);
+ }
+ }
+
+ PAIRS_DEBUG("\b]\n");
+ }
+
+ // Receive buffer debug
+ for(int i = 0; i < 2; i++) {
+ int nrec = nrecv[dim * 2 + i];
+
+ PAIRS_DEBUG("recv_buf=[");
+ for(int j = 0; j < MIN(elems_to_print, nrec); j++) {
+ for(int k = 0; k < elem_size; k++) {
+ PAIRS_DEBUG("%f,", recv_buf[(recv_offsets[dim * 2 + i] + j) * elem_size + k]);
+ }
+ }
+
+ if(elems_to_print * 2 < nrec) {
+ PAIRS_DEBUG("\b ... ");
+ }
+
+ for(int j = MAX(elems_to_print, nrec - elems_to_print); j < nrec; j++) {
+ for(int k = 0; k < elem_size; k++) {
+ PAIRS_DEBUG("%f,", recv_buf[(recv_offsets[dim * 2 + i] + j) * elem_size + k]);
+ }
+ }
+
+ PAIRS_DEBUG("\b]\n");
+ }
+ */
+}
+
+void PairsSimulation::fillCommunicationArrays(int *neighbor_ranks, int *pbc, real_t *subdom) {
+ this->getDomainPartitioner()->fillArrays(neighbor_ranks, pbc, subdom);
+}
+
+}
diff --git a/runtime/pairs.hpp b/runtime/pairs.hpp
index 3c37795dc327780f1b43c0a6366d8d52fc10c28f..7c5c80d65503466a8c2f937e18a62e81d64fdd9a 100644
--- a/runtime/pairs.hpp
+++ b/runtime/pairs.hpp
@@ -1,264 +1,24 @@
#include <algorithm>
#include <cmath>
-#include <iostream>
#include <memory>
#include <vector>
//---
-#ifdef PAIRS_TARGET_CUDA
-# include "devices/cuda.hpp"
-#else
-# include "devices/dummy.hpp"
-#endif
-
+#include "array.hpp"
+#include "device_flags.hpp"
+#include "pairs_common.hpp"
+#include "property.hpp"
+#include "vector3.hpp"
+#include "devices/device.hpp"
#include "domain/regular_6d_stencil.hpp"
#pragma once
-#ifdef DEBUG
-# include <assert.h>
-# define PAIRS_DEBUG(...) fprintf(stderr, __VA_ARGS__)
-# define PAIRS_ASSERT(a) assert(a)
-# define PAIRS_EXCEPTION(a)
-#else
-# define PAIRS_DEBUG(...)
-# define PAIRS_ASSERT(a)
-# define PAIRS_EXCEPTION(a)
-#endif
-
-#define PAIRS_ERROR(...) fprintf(stderr, __VA_ARGS__)
-#define MIN(a,b) ((a) < (b) ? (a) : (b))
-#define MAX(a,b) ((a) > (b) ? (a) : (b))
-
namespace pairs {
-typedef int array_t;
-typedef int property_t;
-typedef int layout_t;
-
-enum PropertyType {
- Prop_Invalid = -1,
- Prop_Integer = 0,
- Prop_Float,
- Prop_Vector
-};
-
-enum DataLayout {
- Invalid = -1,
- AoS = 0,
- SoA
-};
-
-enum DomainPartitioning {
- DimRanges = 0,
- BoxList,
-};
-
-template<typename real_t>
-class Vector3 {
-private:
- real_t x, y, z;
-
-public:
- Vector3() : x((real_t) 0.0), y((real_t) 0.0), z((real_t) 0.0) {}
- Vector3(real_t v) : x(v), y(v), z(v) {}
- Vector3(real_t x_, real_t y_, real_t z_) : x(x_), y(y_), z(z_) {}
-};
-
-class Array {
-protected:
- array_t id;
- std::string name;
- void *h_ptr, *d_ptr;
- size_t size;
- bool is_static;
-
-public:
- Array(array_t id_, std::string name_, void *h_ptr_, void *d_ptr_, size_t size_, bool is_static_ = false) :
- id(id_),
- name(name_),
- h_ptr(h_ptr_),
- d_ptr(d_ptr_),
- size(size_),
- is_static(is_static_) {
-
- PAIRS_ASSERT(size_ > 0);
- }
-
- property_t getId() { return id; }
- std::string getName() { return name; }
- void *getHostPointer() { return h_ptr; }
- void *getDevicePointer() { return d_ptr; }
- void setPointers(void *h_ptr_, void *d_ptr_) { h_ptr = h_ptr_, d_ptr = d_ptr_; }
- void setSize(size_t size_) { size = size_; }
- size_t getSize() { return size; };
- bool isStatic() { return is_static; }
-};
-
-class Property {
-protected:
- property_t id;
- std::string name;
- void *h_ptr, *d_ptr;
- PropertyType type;
- layout_t layout;
- size_t sx, sy;
-
-public:
- /*Property(property_t id_, std::string name_, void *h_ptr_, void *d_ptr_, PropertyType type_, layout_t layout, size_t sx_) :
- id(id_),
- name(name_),
- h_ptr(h_ptr_),
- d_ptr(d_ptr_),
- type(type_),
- layout(Invalid),
- sx(sx_), sy(1) {
-
- PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0);
- }*/
-
- Property(property_t id_, std::string name_, void *h_ptr_, void *d_ptr_, PropertyType type_, layout_t layout_, size_t sx_, size_t sy_=1) :
- id(id_),
- name(name_),
- h_ptr(h_ptr_),
- d_ptr(d_ptr_),
- type(type_),
- layout(layout_),
- sx(sx_), sy(sy_) {
-
- PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
- }
-
- property_t getId() { return id; }
- std::string getName() { return name; }
- void *getHostPointer() { return h_ptr; }
- void *getDevicePointer() { return d_ptr; }
- void setPointers(void *h_ptr_, void *d_ptr_) { h_ptr = h_ptr_, d_ptr = d_ptr_; }
- void setSizes(size_t sx_, size_t sy_) { sx = sx_, sy = sy_; }
- size_t getTotalSize() { return sx * sy * getElemSize(); };
- PropertyType getType() { return type; }
- layout_t getLayout() { return layout; }
- size_t getElemSize() {
- return (type == Prop_Integer) ? sizeof(int) :
- (type == Prop_Float) ? sizeof(double) :
- (type == Prop_Vector) ? sizeof(double) : 0;
- }
-};
-
-class IntProperty : public Property {
-public:
- inline int &operator()(int i) { return static_cast<int *>(h_ptr)[i]; }
-};
-
-class FloatProperty : public Property {
-public:
- inline double &operator()(int i) { return static_cast<double *>(h_ptr)[i]; }
-};
-
-class VectorProperty : public Property {
-public:
- double &operator()(int i, int j) {
- PAIRS_ASSERT(type != Prop_Invalid && layout != Invalid && sx > 0 && sy > 0);
- double *dptr = static_cast<double *>(h_ptr);
- if(layout == AoS) { return dptr[i * sy + j]; }
- if(layout == SoA) { return dptr[j * sx + i]; }
- PAIRS_ERROR("VectorProperty::operator[]: Invalid data layout!");
- return dptr[0];
- }
-};
-
-class DeviceFlags {
-private:
- unsigned long long int *hflags;
- unsigned long long int *dflags;
- int narrays;
- int nflags;
- static const int narrays_per_flag = 64;
-public:
- DeviceFlags(int narrays_) : narrays(narrays_) {
- nflags = std::ceil((double) narrays_ / (double) narrays_per_flag);
- hflags = new unsigned long long int[nflags];
- dflags = new unsigned long long int[nflags];
-
- for(int i = 0; i < nflags; ++i) {
- hflags[i] = 0xfffffffffffffffful;
- dflags[i] = 0x0;
- }
- }
-
- inline bool isHostFlagSet(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- return (hflags[flag_index] & (1 << bit)) != 0;
- }
-
- inline void setHostFlag(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- hflags[flag_index] |= 1 << bit;
- }
-
- inline void clearHostFlag(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- hflags[flag_index] &= ~(1 << bit);
- }
-
- inline bool isDeviceFlagSet(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- return (dflags[flag_index] & (1 << bit)) != 0;
- }
-
- inline void setDeviceFlag(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- dflags[flag_index] |= 1 << bit;
- }
-
- inline void clearDeviceFlag(int array_id) {
- int flag_index = array_id / narrays_per_flag;
- unsigned long long int bit = array_id % narrays_per_flag;
- dflags[flag_index] &= ~(1 << bit);
- }
-
- void printFlags() {
- fprintf(stderr, "hflags = ");
- for(int i = 0; i < nflags; i++) {
- for(int b = 63; b >= 0; b--) {
- if(hflags[i] & (1 << b)) {
- fprintf(stderr, "1");
- } else {
- fprintf(stderr, "0");
- }
- }
- }
-
- fprintf(stderr, "\n");
- fprintf(stderr, "dflags = ");
- for(int i = 0; i < nflags; i++) {
- for(int b = 63; b >= 0; b--) {
- if(dflags[i] & (1 << b)) {
- fprintf(stderr, "1");
- } else {
- fprintf(stderr, "0");
- }
- }
- }
-
- fprintf(stderr, "\n");
- }
-
- ~DeviceFlags() {
- delete[] hflags;
- delete[] dflags;
- }
-};
-
-template<int ndims>
class PairsSimulation {
private:
- Regular6DStencil<ndims> *dom_part;
- //DomainPartitioner<ndims> *dom_part;
+ Regular6DStencil *dom_part;
+ //DomainPartitioner *dom_part;
std::vector<Property> properties;
std::vector<Array> arrays;
DeviceFlags *prop_flags, *array_flags;
@@ -277,179 +37,32 @@ public:
delete array_flags;
}
- void initDomain(int *argc, char ***argv, real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax) {
- if(dom_part_type == DimRanges) {
- dom_part = new Regular6DStencil<ndims>(xmin, xmax, ymin, ymax, zmin, zmax);
- } else {
- PAIRS_EXCEPTION("Domain partitioning type not implemented!\n");
- }
-
- dom_part->initialize(argc, argv);
- }
-
- Regular6DStencil<ndims> *getDomainPartitioner() { return dom_part; }
-
- template<typename T_ptr>
- void addArray(array_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, size_t size) {
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) malloc(size);
- PAIRS_ASSERT(*h_ptr != nullptr);
- addArray(Array(id, name, *h_ptr, nullptr, size, false));
- }
-
- template<typename T_ptr>
- void addArray(array_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, size_t size) {
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) malloc(size);
- *d_ptr = (T_ptr *) pairs::device_alloc(size);
- PAIRS_ASSERT(*h_ptr != nullptr && *d_ptr != nullptr);
- addArray(Array(id, name, *h_ptr, *d_ptr, size, false));
- }
-
- template<typename T_ptr>
- void addStaticArray(array_t id, std::string name, T_ptr *h_ptr, std::nullptr_t, size_t size) {
- addArray(Array(id, name, h_ptr, nullptr, size, true));
- }
-
- template<typename T_ptr>
- void addStaticArray(array_t id, std::string name, T_ptr *h_ptr, T_ptr *d_ptr, size_t size) {
- addArray(Array(id, name, h_ptr, d_ptr, size, true));
- }
-
- void addArray(Array array) {
- int id = array.getId();
- auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
- PAIRS_ASSERT(a == std::end(arrays));
- arrays.push_back(array);
- }
-
- template<typename T_ptr>
- void reallocArray(array_t id, T_ptr **h_ptr, std::nullptr_t, size_t size) {
- // This should be a pointer (and not a reference) in order to be modified
- auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
- PAIRS_ASSERT(a != std::end(arrays));
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) realloc(*h_ptr, size);
- PAIRS_ASSERT(*h_ptr != nullptr);
-
- a->setPointers(*h_ptr, nullptr);
- a->setSize(size);
- }
-
- template<typename T_ptr>
- void reallocArray(array_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t size) {
- // This should be a pointer (and not a reference) in order to be modified
- auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
- PAIRS_ASSERT(a != std::end(arrays));
- PAIRS_ASSERT(size > 0);
-
- void *new_h_ptr = realloc(*h_ptr, size);
- void *new_d_ptr = pairs::device_realloc(*d_ptr, size);
- PAIRS_ASSERT(new_h_ptr != nullptr && new_d_ptr != nullptr);
-
- a->setPointers(new_h_ptr, new_d_ptr);
- a->setSize(size);
-
- *h_ptr = (T_ptr *) new_h_ptr;
- *d_ptr = (T_ptr *) new_d_ptr;
- if(array_flags->isDeviceFlagSet(id)) {
- copyArrayToDevice(id);
- }
- }
+ void initDomain(int *argc, char ***argv, real_t xmin, real_t xmax, real_t ymin, real_t ymax, real_t zmin, real_t zmax);
+ Regular6DStencil *getDomainPartitioner() { return dom_part; }
- Array &getArray(array_t id) {
- auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
- PAIRS_ASSERT(a != std::end(arrays));
- return *a;
- }
+ template<typename T_ptr> void addArray(array_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, size_t size);
+ template<typename T_ptr> void addArray(array_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, size_t size);
+ template<typename T_ptr> void addStaticArray(array_t id, std::string name, T_ptr *h_ptr, std::nullptr_t, size_t size);
+ template<typename T_ptr> void addStaticArray(array_t id, std::string name, T_ptr *h_ptr, T_ptr *d_ptr, size_t size);
+ void addArray(Array array);
- Array &getArrayByName(std::string name) {
- auto a = std::find_if(arrays.begin(), arrays.end(), [name](Array a) { return a.getName() == name; });
- PAIRS_ASSERT(a != std::end(arrays));
- return *a;
- }
+ template<typename T_ptr> void reallocArray(array_t id, T_ptr **h_ptr, std::nullptr_t, size_t size);
+ template<typename T_ptr> void reallocArray(array_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t size);
- template<typename T_ptr>
- void addProperty(property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, PropertyType type, layout_t layout, size_t sx, size_t sy=1) {
- size_t size = sx * sy * sizeof(T_ptr);
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) malloc(size);
- PAIRS_ASSERT(*h_ptr != nullptr);
- addProperty(Property(id, name, *h_ptr, nullptr, type, layout, sx, sy));
- }
+ Array &getArray(array_t id);
+ Array &getArrayByName(std::string name);
- template<typename T_ptr>
- void addProperty(property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, PropertyType type, layout_t layout, size_t sx, size_t sy=1) {
- size_t size = sx * sy * sizeof(T_ptr);
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) malloc(size);
- *d_ptr = (T_ptr *) pairs::device_alloc(size);
- PAIRS_ASSERT(*h_ptr != nullptr && *d_ptr != nullptr);
- addProperty(Property(id, name, *h_ptr, *d_ptr, type, layout, sx, sy));
- }
-
- void addProperty(Property prop) {
- int id = prop.getId();
- auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
- PAIRS_ASSERT(p == std::end(properties));
- properties.push_back(prop);
- }
-
- template<typename T_ptr>
- void reallocProperty(property_t id, T_ptr **h_ptr, std::nullptr_t, size_t sx = 1, size_t sy = 1) {
- // This should be a pointer (and not a reference) in order to be modified
- auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
- PAIRS_ASSERT(p != std::end(properties));
-
- size_t size = sx * sy * p->getElemSize();
- PAIRS_ASSERT(size > 0);
-
- *h_ptr = (T_ptr *) realloc(*h_ptr, size);
- PAIRS_ASSERT(*h_ptr != nullptr);
-
- p->setPointers(*h_ptr, nullptr);
- p->setSizes(sx, sy);
- }
+ template<typename T_ptr> void addProperty(
+ property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, PropertyType type, layout_t layout, size_t sx, size_t sy = 1);
+ template<typename T_ptr> void addProperty(
+ property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, PropertyType type, layout_t layout, size_t sx, size_t sy = 1);
+ void addProperty(Property prop);
- template<typename T_ptr>
- void reallocProperty(property_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t sx = 1, size_t sy = 1) {
- // This should be a pointer (and not a reference) in order to be modified
- auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
- PAIRS_ASSERT(p != std::end(properties));
+ template<typename T_ptr> void reallocProperty(property_t id, T_ptr **h_ptr, std::nullptr_t, size_t sx = 1, size_t sy = 1);
+ template<typename T_ptr> void reallocProperty(property_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t sx = 1, size_t sy = 1);
- size_t size = sx * sy * p->getElemSize();
- PAIRS_ASSERT(size > 0);
-
- void *new_h_ptr = realloc(*h_ptr, size);
- void *new_d_ptr = pairs::device_realloc(*d_ptr, size);
- PAIRS_ASSERT(new_h_ptr != nullptr && new_d_ptr != nullptr);
-
- p->setPointers(new_h_ptr, new_d_ptr);
- p->setSizes(sx, sy);
-
- *h_ptr = (T_ptr *) new_h_ptr;
- *d_ptr = (T_ptr *) new_d_ptr;
- if(prop_flags->isDeviceFlagSet(id)) {
- copyPropertyToDevice(id);
- }
- }
-
- Property &getProperty(property_t id) {
- auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
- PAIRS_ASSERT(p != std::end(properties));
- return *p;
- }
-
- Property &getPropertyByName(std::string name) {
- auto p = std::find_if(properties.begin(), properties.end(), [name](Property p) { return p.getName() == name; });
- PAIRS_ASSERT(p != std::end(properties));
- return *p;
- }
+ Property &getProperty(property_t id);
+ Property &getPropertyByName(std::string name);
inline IntProperty &getAsIntegerProperty(Property &prop) { return static_cast<IntProperty&>(prop); }
inline FloatProperty &getAsFloatProperty(Property &prop) { return static_cast<FloatProperty&>(prop); }
@@ -463,130 +76,164 @@ public:
void clearArrayDeviceFlag(array_t id) { clearArrayDeviceFlag(getArray(id)); }
void clearArrayDeviceFlag(Array &array) { array_flags->clearDeviceFlag(array.getId()); }
void copyArrayToDevice(array_t id) { copyArrayToDevice(getArray(id)); }
- void copyArrayToDevice(Array &array) {
- int array_id = array.getId();
- if(!array_flags->isDeviceFlagSet(array_id)) {
- if(array.isStatic()) {
- PAIRS_DEBUG("Copying static array %s to device\n", array.getName().c_str());
- pairs::copy_static_symbol_to_device(array.getHostPointer(), array.getDevicePointer(), array.getSize());
- } else {
- PAIRS_DEBUG("Copying array %s to device\n", array.getName().c_str());
- pairs::copy_to_device(array.getHostPointer(), array.getDevicePointer(), array.getSize());
- }
- }
- }
+ void copyArrayToDevice(Array &array);
void setArrayHostFlag(array_t id) { setArrayHostFlag(getArray(id)); }
void setArrayHostFlag(Array &array) { array_flags->setHostFlag(array.getId()); }
void clearArrayHostFlag(array_t id) { clearArrayHostFlag(getArray(id)); }
void clearArrayHostFlag(Array &array) { array_flags->clearHostFlag(array.getId()); }
void copyArrayToHost(array_t id) { copyArrayToHost(getArray(id)); }
- void copyArrayToHost(Array &array) {
- int array_id = array.getId();
- if(!array_flags->isHostFlagSet(array_id)) {
- if(array.isStatic()) {
- PAIRS_DEBUG("Copying static array %s to host\n", array.getName().c_str());
- pairs::copy_static_symbol_to_host(array.getDevicePointer(), array.getHostPointer(), array.getSize());
- } else {
- PAIRS_DEBUG("Copying array %s to host\n", array.getName().c_str());
- pairs::copy_to_host(array.getDevicePointer(), array.getHostPointer(), array.getSize());
- }
- }
- }
+ void copyArrayToHost(Array &array);
void setPropertyDeviceFlag(property_t id) { setPropertyDeviceFlag(getProperty(id)); }
void setPropertyDeviceFlag(Property &prop) { prop_flags->setDeviceFlag(prop.getId()); }
void clearPropertyDeviceFlag(property_t id) { clearPropertyDeviceFlag(getProperty(id)); }
void clearPropertyDeviceFlag(Property &prop) { prop_flags->clearDeviceFlag(prop.getId()); }
void copyPropertyToDevice(property_t id) { copyPropertyToDevice(getProperty(id)); }
- void copyPropertyToDevice(Property &prop) {
- if(!prop_flags->isDeviceFlagSet(prop.getId())) {
- PAIRS_DEBUG("Copying property %s to device\n", prop.getName().c_str());
- pairs::copy_to_device(prop.getHostPointer(), prop.getDevicePointer(), prop.getTotalSize());
- }
- }
+ void copyPropertyToDevice(Property &prop);
void setPropertyHostFlag(property_t id) { setPropertyHostFlag(getProperty(id)); }
void setPropertyHostFlag(Property &prop) { prop_flags->setHostFlag(prop.getId()); }
void clearPropertyHostFlag(property_t id) { clearPropertyHostFlag(getProperty(id)); }
void clearPropertyHostFlag(Property &prop) { prop_flags->clearHostFlag(prop.getId()); }
void copyPropertyToHost(property_t id) { copyPropertyToHost(getProperty(id)); }
- void copyPropertyToHost(Property &prop) {
- if(!prop_flags->isHostFlagSet(prop.getId())) {
- PAIRS_DEBUG("Copying property %s to host\n", prop.getName().c_str());
- pairs::copy_to_host(prop.getDevicePointer(), prop.getHostPointer(), prop.getTotalSize());
- }
- }
-
- void communicateSizes(int dim, const int *send_sizes, int *recv_sizes) {
- this->getDomainPartitioner()->communicateSizes(dim, send_sizes, recv_sizes);
- PAIRS_DEBUG("send_sizes=[%d, %d], recv_sizes=[%d, %d]\n", send_sizes[dim * 2 + 0], send_sizes[dim * 2 + 1], recv_sizes[dim * 2 + 0], recv_sizes[dim * 2 + 1]);
- }
+ void copyPropertyToHost(Property &prop);
+ void communicateSizes(int dim, const int *send_sizes, int *recv_sizes);
void communicateData(
int dim, int elem_size,
const real_t *send_buf, const int *send_offsets, const int *nsend,
- real_t *recv_buf, const int *recv_offsets, const int *nrecv) {
-
- this->getDomainPartitioner()->communicateData(dim, elem_size, send_buf, send_offsets, nsend, recv_buf, recv_offsets, nrecv);
-
- /*
- // Debug messages
- const int elems_to_print = 5;
-
- // Send buffer debug
- for(int i = 0; i < 2; i++) {
- int nsnd = nsend[dim * 2 + i];
-
- PAIRS_DEBUG("send_buf=[");
- for(int j = 0; j < MIN(elems_to_print, nsnd); j++) {
- for(int k = 0; k < elem_size; k++) {
- PAIRS_DEBUG("%f,", send_buf[(send_offsets[dim * 2 + i] + j) * elem_size + k]);
- }
- }
-
- if(elems_to_print * 2 < nsnd) {
- PAIRS_DEBUG("\b ... ");
- }
-
- for(int j = MAX(elems_to_print, nsnd - elems_to_print); j < nsnd; j++) {
- for(int k = 0; k < elem_size; k++) {
- PAIRS_DEBUG("%f,", send_buf[(send_offsets[dim * 2 + i] + j) * elem_size + k]);
- }
- }
-
- PAIRS_DEBUG("\b]\n");
- }
-
- // Receive buffer debug
- for(int i = 0; i < 2; i++) {
- int nrec = nrecv[dim * 2 + i];
-
- PAIRS_DEBUG("recv_buf=[");
- for(int j = 0; j < MIN(elems_to_print, nrec); j++) {
- for(int k = 0; k < elem_size; k++) {
- PAIRS_DEBUG("%f,", recv_buf[(recv_offsets[dim * 2 + i] + j) * elem_size + k]);
- }
- }
-
- if(elems_to_print * 2 < nrec) {
- PAIRS_DEBUG("\b ... ");
- }
-
- for(int j = MAX(elems_to_print, nrec - elems_to_print); j < nrec; j++) {
- for(int k = 0; k < elem_size; k++) {
- PAIRS_DEBUG("%f,", recv_buf[(recv_offsets[dim * 2 + i] + j) * elem_size + k]);
- }
- }
-
- PAIRS_DEBUG("\b]\n");
- }
- */
+ real_t *recv_buf, const int *recv_offsets, const int *nrecv);
+
+ void fillCommunicationArrays(int neighbor_ranks[], int pbc[], real_t subdom[]);
+};
+
+template<typename T_ptr>
+void PairsSimulation::addArray(array_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, size_t size) {
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) malloc(size);
+ PAIRS_ASSERT(*h_ptr != nullptr);
+ addArray(Array(id, name, *h_ptr, nullptr, size, false));
+}
+
+template<typename T_ptr>
+void PairsSimulation::addArray(array_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, size_t size) {
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) malloc(size);
+ *d_ptr = (T_ptr *) pairs::device_alloc(size);
+ PAIRS_ASSERT(*h_ptr != nullptr && *d_ptr != nullptr);
+ addArray(Array(id, name, *h_ptr, *d_ptr, size, false));
+}
+
+template<typename T_ptr>
+void PairsSimulation::addStaticArray(array_t id, std::string name, T_ptr *h_ptr, std::nullptr_t, size_t size) {
+ addArray(Array(id, name, h_ptr, nullptr, size, true));
+}
+
+template<typename T_ptr>
+void PairsSimulation::addStaticArray(array_t id, std::string name, T_ptr *h_ptr, T_ptr *d_ptr, size_t size) {
+ addArray(Array(id, name, h_ptr, d_ptr, size, true));
+}
+
+template<typename T_ptr>
+void PairsSimulation::reallocArray(array_t id, T_ptr **h_ptr, std::nullptr_t, size_t size) {
+ // This should be a pointer (and not a reference) in order to be modified
+ auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+ PAIRS_ASSERT(a != std::end(arrays));
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) realloc(*h_ptr, size);
+ PAIRS_ASSERT(*h_ptr != nullptr);
+
+ a->setPointers(*h_ptr, nullptr);
+ a->setSize(size);
+}
+
+template<typename T_ptr>
+void PairsSimulation::reallocArray(array_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t size) {
+ // This should be a pointer (and not a reference) in order to be modified
+ auto a = std::find_if(arrays.begin(), arrays.end(), [id](Array a) { return a.getId() == id; });
+ PAIRS_ASSERT(a != std::end(arrays));
+ PAIRS_ASSERT(size > 0);
+
+ void *new_h_ptr = realloc(*h_ptr, size);
+ void *new_d_ptr = pairs::device_realloc(*d_ptr, size);
+ PAIRS_ASSERT(new_h_ptr != nullptr && new_d_ptr != nullptr);
+
+ a->setPointers(new_h_ptr, new_d_ptr);
+ a->setSize(size);
+
+ *h_ptr = (T_ptr *) new_h_ptr;
+ *d_ptr = (T_ptr *) new_d_ptr;
+ if(array_flags->isDeviceFlagSet(id)) {
+ copyArrayToDevice(id);
}
+}
+
+template<typename T_ptr>
+void PairsSimulation::addProperty(
+ property_t id, std::string name, T_ptr **h_ptr, std::nullptr_t, PropertyType type, layout_t layout, size_t sx, size_t sy) {
- void fillCommunicationArrays(int neighbor_ranks[], int pbc[], real_t subdom[]) {
- this->getDomainPartitioner()->fillArrays(neighbor_ranks, pbc, subdom);
+ size_t size = sx * sy * sizeof(T_ptr);
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) malloc(size);
+ PAIRS_ASSERT(*h_ptr != nullptr);
+ addProperty(Property(id, name, *h_ptr, nullptr, type, layout, sx, sy));
+}
+
+template<typename T_ptr>
+void PairsSimulation::addProperty(
+ property_t id, std::string name, T_ptr **h_ptr, T_ptr **d_ptr, PropertyType type, layout_t layout, size_t sx, size_t sy) {
+
+ size_t size = sx * sy * sizeof(T_ptr);
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) malloc(size);
+ *d_ptr = (T_ptr *) pairs::device_alloc(size);
+ PAIRS_ASSERT(*h_ptr != nullptr && *d_ptr != nullptr);
+ addProperty(Property(id, name, *h_ptr, *d_ptr, type, layout, sx, sy));
+}
+
+template<typename T_ptr>
+void PairsSimulation::reallocProperty(property_t id, T_ptr **h_ptr, std::nullptr_t, size_t sx, size_t sy) {
+ // This should be a pointer (and not a reference) in order to be modified
+ auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
+ PAIRS_ASSERT(p != std::end(properties));
+
+ size_t size = sx * sy * p->getElemSize();
+ PAIRS_ASSERT(size > 0);
+
+ *h_ptr = (T_ptr *) realloc(*h_ptr, size);
+ PAIRS_ASSERT(*h_ptr != nullptr);
+
+ p->setPointers(*h_ptr, nullptr);
+ p->setSizes(sx, sy);
+}
+
+template<typename T_ptr>
+void PairsSimulation::reallocProperty(property_t id, T_ptr **h_ptr, T_ptr **d_ptr, size_t sx, size_t sy) {
+ // This should be a pointer (and not a reference) in order to be modified
+ auto p = std::find_if(properties.begin(), properties.end(), [id](Property p) { return p.getId() == id; });
+ PAIRS_ASSERT(p != std::end(properties));
+
+ size_t size = sx * sy * p->getElemSize();
+ PAIRS_ASSERT(size > 0);
+
+ void *new_h_ptr = realloc(*h_ptr, size);
+ void *new_d_ptr = pairs::device_realloc(*d_ptr, size);
+ PAIRS_ASSERT(new_h_ptr != nullptr && new_d_ptr != nullptr);
+
+ p->setPointers(new_h_ptr, new_d_ptr);
+ p->setSizes(sx, sy);
+
+ *h_ptr = (T_ptr *) new_h_ptr;
+ *d_ptr = (T_ptr *) new_d_ptr;
+ if(prop_flags->isDeviceFlagSet(id)) {
+ copyPropertyToDevice(id);
}
-};
+}
}
diff --git a/runtime/pairs_common.hpp b/runtime/pairs_common.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e74453c60692b204650f41333124a3ca693a75e7
--- /dev/null
+++ b/runtime/pairs_common.hpp
@@ -0,0 +1,41 @@
+#include <iostream>
+
+#pragma once
+
+typedef double real_t;
+typedef int array_t;
+typedef int property_t;
+typedef int layout_t;
+
+enum PropertyType {
+ Prop_Invalid = -1,
+ Prop_Integer = 0,
+ Prop_Float,
+ Prop_Vector
+};
+
+enum DataLayout {
+ Invalid = -1,
+ AoS = 0,
+ SoA
+};
+
+enum DomainPartitioning {
+ DimRanges = 0,
+ BoxList,
+};
+
+#ifdef DEBUG
+# include <assert.h>
+# define PAIRS_DEBUG(...) fprintf(stderr, __VA_ARGS__)
+# define PAIRS_ASSERT(a) assert(a)
+# define PAIRS_EXCEPTION(a)
+#else
+# define PAIRS_DEBUG(...)
+# define PAIRS_ASSERT(a)
+# define PAIRS_EXCEPTION(a)
+#endif
+
+#define PAIRS_ERROR(...) fprintf(stderr, __VA_ARGS__)
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
diff --git a/runtime/property.hpp b/runtime/property.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4690a8f14a435045867a9baa104640ccd2e87f1f
--- /dev/null
+++ b/runtime/property.hpp
@@ -0,0 +1,67 @@
+#include "pairs_common.hpp"
+
+#pragma once
+
+namespace pairs {
+
+class Property {
+protected:
+ property_t id;
+ std::string name;
+ void *h_ptr, *d_ptr;
+ PropertyType type;
+ layout_t layout;
+ size_t sx, sy;
+
+public:
+ Property(property_t id_, std::string name_, void *h_ptr_, void *d_ptr_, PropertyType type_, layout_t layout_, size_t sx_, size_t sy_=1) :
+ id(id_),
+ name(name_),
+ h_ptr(h_ptr_),
+ d_ptr(d_ptr_),
+ type(type_),
+ layout(layout_),
+ sx(sx_), sy(sy_) {
+
+ PAIRS_ASSERT(type != Prop_Invalid && layout_ != Invalid && sx_ > 0 && sy_ > 0);
+ }
+
+ property_t getId() { return id; }
+ std::string getName() { return name; }
+ void *getHostPointer() { return h_ptr; }
+ void *getDevicePointer() { return d_ptr; }
+ void setPointers(void *h_ptr_, void *d_ptr_) { h_ptr = h_ptr_, d_ptr = d_ptr_; }
+ void setSizes(size_t sx_, size_t sy_) { sx = sx_, sy = sy_; }
+ size_t getTotalSize() { return sx * sy * getElemSize(); };
+ PropertyType getType() { return type; }
+ layout_t getLayout() { return layout; }
+ size_t getElemSize() {
+ return (type == Prop_Integer) ? sizeof(int) :
+ (type == Prop_Float) ? sizeof(double) :
+ (type == Prop_Vector) ? sizeof(double) : 0;
+ }
+};
+
+class IntProperty : public Property {
+public:
+ inline int &operator()(int i) { return static_cast<int *>(h_ptr)[i]; }
+};
+
+class FloatProperty : public Property {
+public:
+ inline double &operator()(int i) { return static_cast<double *>(h_ptr)[i]; }
+};
+
+class VectorProperty : public Property {
+public:
+ double &operator()(int i, int j) {
+ PAIRS_ASSERT(type != Prop_Invalid && layout != Invalid && sx > 0 && sy > 0);
+ double *dptr = static_cast<double *>(h_ptr);
+ if(layout == AoS) { return dptr[i * sy + j]; }
+ if(layout == SoA) { return dptr[j * sx + i]; }
+ PAIRS_ERROR("VectorProperty::operator[]: Invalid data layout!");
+ return dptr[0];
+ }
+};
+
+}
diff --git a/runtime/read_from_file.hpp b/runtime/read_from_file.hpp
index 7a681dc0b564bd005a357e6b7a6694155a015990..ba80da72897937e8db58fea068b6268fdc0cddf0 100644
--- a/runtime/read_from_file.hpp
+++ b/runtime/read_from_file.hpp
@@ -4,13 +4,13 @@
#include <sstream>
//---
#include "pairs.hpp"
+#include "pairs_common.hpp"
#pragma once
namespace pairs {
-template<int ndims>
-void read_grid_data(PairsSimulation<ndims> *ps, const char *filename, double *grid_buffer) {
+void read_grid_data(PairsSimulation *ps, const char *filename, double *grid_buffer) {
std::ifstream in_file(filename, std::ifstream::in);
std::string line;
@@ -21,7 +21,7 @@ void read_grid_data(PairsSimulation<ndims> *ps, const char *filename, double *gr
int i = 0;
while(std::getline(line_stream, in0, ',')) {
- PAIRS_ASSERT(i < ndims * 2);
+ //PAIRS_ASSERT(i < ndims * 2);
grid_buffer[i] = std::stod(in0);
i++;
}
@@ -30,8 +30,7 @@ void read_grid_data(PairsSimulation<ndims> *ps, const char *filename, double *gr
}
}
-template<int ndims>
-size_t read_particle_data(PairsSimulation<ndims> *ps, const char *filename, const property_t properties[], size_t nprops) {
+size_t read_particle_data(PairsSimulation *ps, const char *filename, const property_t properties[], size_t nprops) {
std::ifstream in_file(filename, std::ifstream::in);
std::string line;
size_t n = 0;
diff --git a/runtime/vector3.hpp b/runtime/vector3.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..15277183eb15835e0e625c6f92ef65cccfd6620e
--- /dev/null
+++ b/runtime/vector3.hpp
@@ -0,0 +1,18 @@
+#include "pairs_common.hpp"
+
+#pragma once
+
+namespace pairs {
+
+template<typename real_t>
+class Vector3 {
+private:
+ real_t x, y, z;
+
+public:
+ Vector3() : x((real_t) 0.0), y((real_t) 0.0), z((real_t) 0.0) {}
+ Vector3(real_t v) : x(v), y(v), z(v) {}
+ Vector3(real_t x_, real_t y_, real_t z_) : x(x_), y(y_), z(z_) {}
+};
+
+}
diff --git a/runtime/vtk.hpp b/runtime/vtk.hpp
index eb10bce12a4c7902078629062371323c3d12200e..e2041a787dd0468e4a37995721f1e731ebc52c05 100644
--- a/runtime/vtk.hpp
+++ b/runtime/vtk.hpp
@@ -8,8 +8,7 @@
namespace pairs {
-template <int ndims>
-void vtk_write_data(PairsSimulation<ndims> *ps, const char *filename, int start, int end, int timestep) {
+void vtk_write_data(PairsSimulation *ps, const char *filename, int start, int end, int timestep) {
std::string output_filename(filename);
auto masses = ps->getAsFloatProperty(ps->getPropertyByName("mass"));
auto positions = ps->getAsVectorProperty(ps->getPropertyByName("position"));
diff --git a/src/pairs/code_gen/cgen.py b/src/pairs/code_gen/cgen.py
index bd41f7ac9a5fbd20fedb132f5fe67d8d2fbfee07..4efca1aae8ce954ee517f502fdc90537db8b683a 100644
--- a/src/pairs/code_gen/cgen.py
+++ b/src/pairs/code_gen/cgen.py
@@ -58,10 +58,10 @@ class CGen:
self.print("#include \"runtime/read_from_file.hpp\"")
self.print("#include \"runtime/vtk.hpp\"")
- if self.target.is_gpu():
- self.print("#include \"runtime/devices/cuda.hpp\"")
- else:
- self.print("#include \"runtime/devices/dummy.hpp\"")
+ #if self.target.is_gpu():
+ # self.print("#include \"runtime/devices/cuda.hpp\"")
+ #else:
+ # self.print("#include \"runtime/devices/dummy.hpp\"")
self.print("")
self.print("using namespace pairs;")
@@ -91,7 +91,7 @@ class CGen:
nprops = module.sim.properties.nprops()
narrays = module.sim.arrays.narrays()
self.print("int main(int argc, char **argv) {")
- self.print(f" PairsSimulation<{ndims}> *pairs = new PairsSimulation<{ndims}>({nprops}, {narrays}, DimRanges);")
+ self.print(f" PairsSimulation *pairs = new PairsSimulation({nprops}, {narrays}, DimRanges);")
self.generate_statement(module.block)
self.print(" return 0;")
self.print("}")