Skip to content

Commit 38ed0a4

Browse files
author
Dominic Marcello
committed
Efficieny improved
1 parent 8268bf9 commit 38ed0a4

5 files changed

Lines changed: 170 additions & 117 deletions

File tree

‎src/defs.hpp‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -23,14 +23,15 @@
2323

2424
#include <array>
2525

26-
#define CORRECTION_ON
27-
#define CORRECTION_OPTIMIZE
26+
//#define CORRECTION_ON
27+
//#define CORRECTION_OPTIMIZE
2828

29-
using real = double;
3029
using integer = long long;
3130

31+
#include "real.hpp"
32+
3233
constexpr integer NDIM = 3;
33-
constexpr integer ncrit = 8;
34+
constexpr integer ncrit = 16;
3435
constexpr integer nparts = 100000;
3536
constexpr integer nchild = 8;
3637
constexpr integer LP = 20;

‎src/expansion.hpp‎

Lines changed: 40 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -24,11 +24,7 @@
2424

2525
template<class Type>
2626
class expansion: public std::array<Type, LP> {
27-
private:
28-
const size_t map2[3][3] = { { 0, 1, 2 }, { 1, 3, 4 }, { 2, 4, 5 } };
29-
const size_t map3[3][3][3] = { { { 0, 1, 2 }, { 1, 3, 4 }, { 2, 4, 5 } }, {
30-
{ 1, 3, 4 }, { 3, 6, 7 }, { 4, 7, 8 } }, { { 2, 4, 5 }, { 4, 7, 8 },
31-
{ 5, 8, 9 } } };
27+
3228
public:
3329
expansion<Type>& operator*=(Type r) {
3430
for (integer i = 0; i != LP; ++i) {
@@ -84,19 +80,31 @@ inline Type& expansion<Type>::operator ()(int i) {
8480

8581
template<class Type>
8682
inline Type expansion<Type>::operator ()(int i, int j) const {
83+
static constexpr size_t
84+
map2[3][3] = { {0, 1, 2}, {1, 3, 4}, {2, 4, 5}};
8785
return (*this)[4 + map2[i][j]];
8886
}
8987
template<class Type>
9088
inline Type& expansion<Type>::operator ()(int i, int j) {
89+
static constexpr size_t
90+
map2[3][3] = { {0, 1, 2}, {1, 3, 4}, {2, 4, 5}};
9191
return (*this)[4 + map2[i][j]];
9292
}
9393

9494
template<class Type>
9595
inline Type expansion<Type>::operator ()(int i, int j, int k) const {
96+
static constexpr size_t
97+
map3[3][3][3] = { { {0, 1, 2}, {1, 3, 4}, {2, 4, 5}}, { {1, 3, 4}, {3, 6, 7}, {4, 7, 8}},
98+
{ { 2, 4, 5}, {4, 7, 8}, {5, 8, 9}}};
99+
96100
return (*this)[10 + map3[i][j][k]];
97101
}
98102
template<class Type>
99103
inline Type& expansion<Type>::operator ()(int i, int j, int k) {
104+
static constexpr size_t
105+
map3[3][3][3] = { { {0, 1, 2}, {1, 3, 4}, {2, 4, 5}}, { {1, 3, 4}, {3, 6, 7}, {4, 7, 8}},
106+
{ { 2, 4, 5}, {4, 7, 8}, {5, 8, 9}}};
107+
100108
return (*this)[10 + map3[i][j][k]];
101109
}
102110

@@ -134,7 +142,7 @@ inline expansion<Type>& expansion<Type>::operator<<=(
134142
}
135143
for (integer a = 0; a < 3; a++) {
136144
for (integer b = 0; b < 3; b++) {
137-
me() += me(a, b) * dX[a] * dX[b] * 0.5;
145+
me() += me(a, b) * dX[a] * dX[b] * real(0.5);
138146
}
139147
}
140148
for (integer a = 0; a < 3; a++) {
@@ -152,7 +160,7 @@ inline expansion<Type>& expansion<Type>::operator<<=(
152160
for (integer a = 0; a < 3; a++) {
153161
for (integer b = 0; b < 3; b++) {
154162
for (integer c = 0; c < 3; c++) {
155-
me(a) += me(a, b, c) * dX[b] * dX[c] * 0.5;
163+
me(a) += me(a, b, c) * dX[b] * dX[c] * real(0.5);
156164
}
157165
}
158166
}
@@ -176,7 +184,7 @@ inline Type expansion<Type>::translate_to_particle(
176184
}
177185
for (integer a = 0; a < 3; a++) {
178186
for (integer b = 0; b < 3; b++) {
179-
this_phi += L(a, b) * dX[a] * dX[b] * 0.5;
187+
this_phi += L(a, b) * dX[a] * dX[b] * real(0.5);
180188
}
181189
}
182190
for (integer a = 0; a < 3; a++) {
@@ -229,6 +237,23 @@ template<class Type>
229237
inline expansion<Type>::~expansion() {
230238
}
231239

240+
static expansion<real> factor;
241+
242+
__attribute((constructor))
243+
static void init_factors() {
244+
factor = 0.0;
245+
factor() += 1.0;
246+
for (integer a = 0; a < NDIM; ++a) {
247+
factor(a) += 1.0;
248+
for (integer b = 0; b < NDIM; ++b) {
249+
factor(a, b) += 1.0;
250+
for (integer c = 0; c < NDIM; ++c) {
251+
factor(a, b, c) += 1.0;
252+
}
253+
}
254+
}
255+
}
256+
232257
template<class Type, class Type2>
233258
inline void multipole_interaction(expansion<Type2>& L1,
234259
const multipole<Type>& M1, const multipole<Type2>& M2,
@@ -309,8 +334,8 @@ inline void multipole_interaction(expansion<Type2>& L1,
309334

310335
L1() += M2() * D();
311336
for (integer a = 0; a < 3; a++) {
312-
for (integer b = 0; b < 3; b++) {
313-
L1() += M2(a, b) * D(a, b) * 0.5;
337+
for (integer b = a; b < 3; b++) {
338+
L1() += M2(a, b) * D(a, b) * real(0.5) * factor(a, b);
314339
}
315340
}
316341

@@ -319,8 +344,8 @@ inline void multipole_interaction(expansion<Type2>& L1,
319344
}
320345
for (integer a = 0; a < 3; a++) {
321346
for (integer b = 0; b < 3; b++) {
322-
for (integer c = 0; c < 3; c++) {
323-
L1(a) += M2(c, b) * D(a, b, c) * 0.5;
347+
for (integer c = b; c < 3; c++) {
348+
L1(a) += M2(c, b) * D(a, b, c) * real(0.5) * factor(c, b);
324349
}
325350
}
326351
}
@@ -342,11 +367,11 @@ inline void multipole_interaction(expansion<Type2>& L1,
342367
#ifdef CORRECTION_ON
343368
for (integer i = 0; i != NDIM; ++i) {
344369
for (integer j = 0; j != NDIM; ++j) {
345-
for (integer k = 0; k != NDIM; ++k) {
346-
for (integer l = 0; l != NDIM; ++l) {
370+
for (integer k = j; k != NDIM; ++k) {
371+
for (integer l = k; l != NDIM; ++l) {
347372
L1(i) -= D4[i](j, k, l)
348373
* (M2(j, k, l) - M1(j, k, l) * M2() / M1())
349-
* (Type(1) / Type(6));
374+
* (Type(1) / Type(6)) * factor(j, k, l);
350375
}
351376
}
352377
}

0 commit comments

Comments
 (0)