-
Notifications
You must be signed in to change notification settings - Fork 76
/
ilp.h
138 lines (113 loc) · 3.68 KB
/
ilp.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#ifndef ILP_H
#define ILP_H
#include <stdint.h>
#ifndef ILP_PART_PER_TILE
#define ILP_PART_PER_TILE 4096 /* 4096*24 ~ 100k */
#endif
#define ILP_ALIGN_BITS 2
#define ILP_ALIGN_SIZE (1<<ILP_ALIGN_BITS)
#define ILP_ALIGN_MASK (ILP_ALIGN_SIZE-1)
typedef struct ilPart {
double m,x,y,z;
double fourh2;
uint64_t iOrder;
} ILP;
#if 0
/*
** We use a union here so that the compiler can properly align the values.
*/
typedef union {
float f[ILP_PART_PER_TILE];
} ilpFloat;
typedef union {
int64_t i[ILP_PART_PER_TILE]; /* Signed because negative marks softened cells */
} ilpInt64;
typedef struct ilpTile {
struct ilpTile *next;
struct ilpTile *prev;
uint32_t nMaxPart; /* Maximum number of particles in this tile */
uint32_t nPart; /* Current number of particles */
ilpFloat dx, dy, dz; /* Offset from ilp->cx, cy, cz */
ilpFloat d2; /* Distance squared: calculated */
ilpFloat m; /* Mass */
ilpFloat fourh2; /* Softening: calculated */
/* #ifdef HERMITE */
ilpFloat vx, vy, vz;
/* #endif */
/* #if defined(SYMBA) || defined(PLANETS) */
ilpInt64 iOrder;
/* #endif */
} *ILPTILE;
typedef struct ilpContext {
ILPTILE first; /* first tile in the chain */
ILPTILE tile; /* Current tile in the chain */
double cx, cy, cz; /* Center coordinates */
uint32_t nPrevious; /* Particles in tiles prior to "tile" */
} *ILP;
typedef struct {
ILPTILE tile;
uint32_t nPart;
uint32_t nPrevious;
} ILPCHECKPT;
ILPTILE ilpExtend(ILP ilp); /* Add tile and return new tile */
ILPTILE ilpClear(ILP ilp); /* Go back to, and return first tile (empty) */
void ilpInitialize(ILP *ilp);
void ilpFinish(ILP ilp);
size_t ilpMemory(ILP ilp);
static inline void ilpCheckPt(ILP ilp,ILPCHECKPT *cp) {
cp->tile = ilp->tile;
cp->nPart = ilp->tile->nPart;
cp->nPrevious = ilp->nPrevious;
}
static inline void ilpRestore(ILP ilp,ILPCHECKPT *cp) {
ilp->tile = cp->tile;
ilp->nPrevious = cp->nPrevious;
ilp->tile->nPart = cp->nPart;
}
static inline uint32_t ilpCount(ILP ilp) {
return ilp->nPrevious + ilp->tile->nPart;
}
/* #if defined(SYMBA) || defined(PLANETS) */
#define ilpAppend_1(ilp,I) tile->iOrder.i[ILP_APPEND_i] = (I);
/* #else */
/* #define ilpAppend_1(ilp,I) */
/* #endif */
/* #if defined(HERMITE) */
#define ilpAppend_2(ilp,VX,VY,VZ) \
tile->vx.f[ILP_APPEND_i] = (VX); \
tile->vy.f[ILP_APPEND_i] = (VY); \
tile->vz.f[ILP_APPEND_i] = (VZ);
/* #else */
/* #define ilpAppend_2(ilp,VX,VY,VZ) */
/* #endif */
#define ilpAppend(ilp,X,Y,Z,M,S,I) \
{ \
ILPTILE tile = (ilp)->tile; \
uint_fast32_t ILP_APPEND_i; \
if ( tile->nPart == tile->nMaxPart ) tile = ilpExtend((ilp)); \
ILP_APPEND_i = tile->nPart; \
tile->dx.f[ILP_APPEND_i] = (ilp)->cx - (X); \
tile->dy.f[ILP_APPEND_i] = (ilp)->cy - (Y); \
tile->dz.f[ILP_APPEND_i] = (ilp)->cz - (Z); \
assert( (M) > 0.0 ); \
tile->m.f[ILP_APPEND_i] = (M); \
tile->fourh2.f[ILP_APPEND_i] = (S); \
ilpAppend_1((ilp),I); \
++tile->nPart; \
}
#define ILP_LOOP(ilp,ptile) for( ptile=(ilp)->first; ptile!=(ilp)->tile->next; ptile=ptile->next )
static inline void ilpCompute(ILP ilp, float fx, float fy, float fz ) {
ILPTILE tile;
uint32_t j;
ILP_LOOP(ilp,tile) {
for (j=0;j<tile->nPart;++j) {
tile->dx.f[j] += fx;
tile->dy.f[j] += fy;
tile->dz.f[j] += fz;
tile->d2.f[j] = tile->dx.f[j]*tile->dx.f[j]
+ tile->dy.f[j]*tile->dy.f[j] + tile->dz.f[j]*tile->dz.f[j];
}
}
}
#endif
#endif