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
139
140
|
//
// Notes:
//
// 1) strategy: one thread per node in the 2D block;
// after initialisation it marches in the k-direction
// working with 3 planes of data at a time
//
// 2) each thread also loads in data for at most one halo node;
// assumes the number of halo nodes is not more than the
// number of interior nodes
//
// 3) corner halo nodes are included because they are needed
// for more general applications with cross-derivatives
//
// 4) could try double-buffering in the future fairly easily
//
// definition to use efficient __mul24 intrinsic
#define INDEX(i,j,j_off) (i +__mul24(j,j_off))
// device code
__global__ void GPU_laplace3d(int NX, int NY, int NZ, int pitch,
float *d_u1, float *d_u2)
{
int indg, indg_h, indg0;
int i, j, k, ind, ind_h, halo, active;
float u2, sixth=1.0f/6.0f;
int NXM1 = NX-1;
int NYM1 = NY-1;
int NZM1 = NZ-1;
//
// define local array offsets
//
#define IOFF 1
#define JOFF (BLOCK_X+2)
#define KOFF (BLOCK_X+2)*(BLOCK_Y+2)
__shared__ float u1[3*KOFF];
//
// first set up indices for halos
//
k = threadIdx.x + threadIdx.y*BLOCK_X;
halo = k < 2*(BLOCK_X+BLOCK_Y+2);
if (halo) {
if (threadIdx.y<2) { // y-halos (coalesced)
i = threadIdx.x;
j = threadIdx.y*(BLOCK_Y+1) - 1;
}
else { // x-halos (not coalesced)
i = (k%2)*(BLOCK_X+1) - 1;
j = k/2 - BLOCK_X - 1;
}
ind_h = INDEX(i+1,j+1,JOFF) + KOFF;
i = INDEX(i,blockIdx.x,BLOCK_X); // global indices
j = INDEX(j,blockIdx.y,BLOCK_Y);
indg_h = INDEX(i,j,pitch);
halo = (i>=0) && (i<NX) && (j>=0) && (j<NY);
}
//
// then set up indices for main block
//
i = threadIdx.x;
j = threadIdx.y;
ind = INDEX(i+1,j+1,JOFF) + KOFF;
i = INDEX(i,blockIdx.x,BLOCK_X); // global indices
j = INDEX(j,blockIdx.y,BLOCK_Y);
indg = INDEX(i,j,pitch);
active = (i<NX) && (j<NY);
//
// read initial plane of u1 array
//
if (active) u1[ind+KOFF] = d_u1[indg];
if (halo) u1[ind_h+KOFF] = d_u1[indg_h];
//
// loop over k-planes
//
for (k=0; k<NZ; k++) {
// move two planes down and read in new plane k+1
if (active) {
indg0 = indg;
indg = INDEX(indg,NY,pitch);
u1[ind-KOFF] = u1[ind];
u1[ind] = u1[ind+KOFF];
if (k<NZM1)
u1[ind+KOFF] = d_u1[indg];
}
if (halo) {
indg_h = INDEX(indg_h,NY,pitch);
u1[ind_h-KOFF] = u1[ind_h];
u1[ind_h] = u1[ind_h+KOFF];
if (k<NZM1)
u1[ind_h+KOFF] = d_u1[indg_h];
}
__syncthreads();
//
// perform Jacobi iteration to set values in u2
//
if (active) {
if (i==0 || i==NXM1 || j==0 || j==NYM1 || k==0 || k==NZM1) {
u2 = u1[ind]; // Dirichlet b.c.'s
}
else {
u2 = ( u1[ind-IOFF] + u1[ind+IOFF]
+ u1[ind-JOFF] + u1[ind+JOFF]
+ u1[ind-KOFF] + u1[ind+KOFF] ) * sixth;
}
d_u2[indg0] = u2;
}
__syncthreads();
}
}
|
