9.2
general documentation
cs_cuda_reduce.h
Go to the documentation of this file.
1#ifndef CS_CUDA_REDUCE_H
2#define CS_CUDA_REDUCE_H
3
4/*============================================================================
5 * BLAS (Basic Linear Algebra Subroutine) functions
6 *============================================================================*/
7
8/*
9 This file is part of code_saturne, a general-purpose CFD tool.
10
11 Copyright (C) 1998-2026 EDF S.A.
12
13 This program is free software; you can redistribute it and/or modify it under
14 the terms of the GNU General Public License as published by the Free Software
15 Foundation; either version 2 of the License, or (at your option) any later
16 version.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
21 details.
22
23 You should have received a copy of the GNU General Public License along with
24 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
25 Street, Fifth Floor, Boston, MA 02110-1301, USA.
26*/
27
28/*----------------------------------------------------------------------------*/
29
30#include "base/cs_defs.h"
31
32/*----------------------------------------------------------------------------
33 * External library headers
34 *----------------------------------------------------------------------------*/
35
36/*----------------------------------------------------------------------------
37 * Local headers
38 *----------------------------------------------------------------------------*/
39
40#include "base/cs_base.h"
41#include "base/cs_base_cuda.h"
42
43/*============================================================================
44 * Macro definitions
45 *============================================================================*/
46
47/*============================================================================
48 * Type definitions
49 *============================================================================*/
50
51/*============================================================================
52 * Templated function definitions
53 *============================================================================*/
54
55#if defined(__CUDACC__)
56
57/*----------------------------------------------------------------------------*/
68/*----------------------------------------------------------------------------*/
69
70template <size_t blockSize, size_t stride, typename T>
71__device__ static void __forceinline__
72cs_cuda_reduce_warp_reduce_sum(volatile T *stmp,
73 size_t tid)
74{
75 if (stride == 1) {
76
77 if (blockSize >= 64) stmp[tid] += stmp[tid + 32];
78 if (blockSize >= 32) stmp[tid] += stmp[tid + 16];
79 if (blockSize >= 16) stmp[tid] += stmp[tid + 8];
80 if (blockSize >= 8) stmp[tid] += stmp[tid + 4];
81 if (blockSize >= 4) stmp[tid] += stmp[tid + 2];
82 if (blockSize >= 2) stmp[tid] += stmp[tid + 1];
83
84 }
85 else {
86
87 if (blockSize >= 64) {
88 #pragma unroll
89 for (size_t i = 0; i < stride; i++)
90 stmp[tid*stride + i] += stmp[(tid + 32)*stride + i];
91 }
92 if (blockSize >= 32) {
93 #pragma unroll
94 for (size_t i = 0; i < stride; i++)
95 stmp[tid*stride + i] += stmp[(tid + 16)*stride + i];
96 }
97 if (blockSize >= 16) {
98 #pragma unroll
99 for (size_t i = 0; i < stride; i++)
100 stmp[tid*stride + i] += stmp[(tid + 8)*stride + i];
101 }
102 if (blockSize >= 8) {
103 #pragma unroll
104 for (size_t i = 0; i < stride; i++)
105 stmp[tid*stride + i] += stmp[(tid + 4)*stride + i];
106 }
107 if (blockSize >= 4) {
108 #pragma unroll
109 for (size_t i = 0; i < stride; i++)
110 stmp[tid*stride + i] += stmp[(tid + 2)*stride + i];
111 }
112 if (blockSize >= 2) {
113 #pragma unroll
114 for (size_t i = 0; i < stride; i++)
115 stmp[tid*stride + i] += stmp[(tid + 1)*stride + i];
116 }
117 }
118}
119
120/*----------------------------------------------------------------------------*/
138/*----------------------------------------------------------------------------*/
139
140template <size_t blockSize, size_t stride, typename T>
141__device__ static void __forceinline__
142cs_cuda_reduce_block_reduce_sum(T *stmp,
143 size_t tid,
144 T *sum_block)
145{
146 __syncthreads();
147
148 /* Loop explicitely unrolled */
149
150 if (stride == 1) { /* Scalar data */
151
152 if (blockSize >= 1024) {
153 if (tid < 512) {
154 stmp[tid] += stmp[tid + 512];
155 }
156 __syncthreads();
157 }
158 if (blockSize >= 512) {
159 if (tid < 256) {
160 stmp[tid] += stmp[tid + 256];
161 }
162 __syncthreads();
163 }
164 if (blockSize >= 256) {
165 if (tid < 128) {
166 stmp[tid] += stmp[tid + 128];
167 } __syncthreads();
168 }
169 if (blockSize >= 128) {
170 if (tid < 64) {
171 stmp[tid] += stmp[tid + 64];
172 } __syncthreads();
173 }
174
175 if (tid < 32) {
176 cs_cuda_reduce_warp_reduce_sum<blockSize, stride>(stmp, tid);
177 }
178
179 // Output: b_res for this block
180
181 if (tid == 0) sum_block[blockIdx.x] = stmp[0];
182
183 }
184
185 else { /* Vector data */
186
187 if (blockSize >= 1024) {
188 if (tid < 512) {
189 #pragma unroll
190 for (size_t i = 0; i < stride; i++)
191 stmp[tid*stride + i] += stmp[(tid + 512)*stride + i];
192 }
193 __syncthreads();
194 }
195 if (blockSize >= 512) {
196 if (tid < 256) {
197 #pragma unroll
198 for (size_t i = 0; i < stride; i++)
199 stmp[tid*stride + i] += stmp[(tid + 256)*stride + i];
200 }
201 __syncthreads();
202 }
203 if (blockSize >= 256) {
204 if (tid < 128) {
205 #pragma unroll
206 for (size_t i = 0; i < stride; i++)
207 stmp[tid*stride + i] += stmp[(tid + 128)*stride + i];
208 } __syncthreads();
209 }
210 if (blockSize >= 128) {
211 if (tid < 64) {
212 #pragma unroll
213 for (size_t i = 0; i < stride; i++)
214 stmp[tid*stride + i] += stmp[(tid + 64)*stride + i];
215 } __syncthreads();
216 }
217
218 if (tid < 32)
219 cs_cuda_reduce_warp_reduce_sum<blockSize, stride>(stmp, tid);
220
221 // Output: b_res for this block
222
223 if (tid == 0) {
224 #pragma unroll
225 for (size_t i = 0; i < stride; i++)
226 sum_block[(blockIdx.x)*stride + i] = stmp[i];
227 }
228
229 }
230}
231
232/*----------------------------------------------------------------------------*/
242/*----------------------------------------------------------------------------*/
243
244template <size_t blockSize, size_t stride, typename T>
245__global__ static void
246cs_cuda_reduce_sum_single_block(size_t n,
247 T *g_idata,
248 T *g_odata)
249{
250 __shared__ T sdata[blockSize * stride];
251
252 size_t tid = threadIdx.x;
253 T r_s[stride];
254
255 if (stride == 1) {
256
257 r_s[0] = 0;
258 sdata[tid] = 0.;
259
260 for (size_t i = threadIdx.x; i < n; i+= blockSize)
261 r_s[0] += g_idata[i];
262
263 sdata[tid] = r_s[0];
264 __syncthreads();
265
266 for (size_t j = blockSize/2; j > CS_CUDA_WARP_SIZE; j /= 2) {
267 if (tid < j) {
268 sdata[tid] += sdata[tid + j];
269 }
270 __syncthreads();
271 }
272
273 if (tid < 32) cs_cuda_reduce_warp_reduce_sum<blockSize, stride>(sdata, tid);
274 if (tid == 0) *g_odata = sdata[0];
275
276 }
277 else {
278
279 #pragma unroll
280 for (size_t k = 0; k < stride; k++) {
281 r_s[k] = 0;
282 sdata[tid*stride + k] = 0.;
283 }
284
285 for (size_t i = threadIdx.x; i < n; i+= blockSize) {
286 #pragma unroll
287 for (size_t k = 0; k < stride; k++) {
288 r_s[k] += g_idata[i*stride + k];
289 }
290 }
291
292 #pragma unroll
293 for (size_t k = 0; k < stride; k++)
294 sdata[tid*stride + k] = r_s[k];
295 __syncthreads();
296
297 for (size_t j = blockSize/2; j > CS_CUDA_WARP_SIZE; j /= 2) {
298 if (tid < j) {
299 #pragma unroll
300 for (size_t k = 0; k < stride; k++)
301 sdata[tid*stride + k] += sdata[(tid + j)*stride + k];
302 }
303 __syncthreads();
304 }
305
306 if (tid < 32) cs_cuda_reduce_warp_reduce_sum<blockSize, stride>(sdata, tid);
307 if (tid == 0) {
308 #pragma unroll
309 for (size_t k = 0; k < stride; k++)
310 g_odata[k] = sdata[k];
311 }
312
313 }
314}
315
316/*----------------------------------------------------------------------------*/
327/*----------------------------------------------------------------------------*/
328
329template <size_t blockSize, typename R, typename T>
330__device__ static void __forceinline__
331cs_cuda_reduce_warp_reduce(volatile T *stmp,
332 size_t tid)
333{
334 R reducer;
335
336 if (blockSize >= 64) reducer.combine(stmp[tid], stmp[tid + 32]);
337 if (blockSize >= 32) reducer.combine(stmp[tid], stmp[tid + 16]);
338 if (blockSize >= 16) reducer.combine(stmp[tid], stmp[tid + 8]);
339 if (blockSize >= 8) reducer.combine(stmp[tid], stmp[tid + 4]);
340 if (blockSize >= 4) reducer.combine(stmp[tid], stmp[tid + 2]);
341 if (blockSize >= 2) reducer.combine(stmp[tid], stmp[tid + 1]);
342}
343
344/*----------------------------------------------------------------------------*/
359/*----------------------------------------------------------------------------*/
360
361template <size_t blockSize, typename R, typename T>
362__device__ static void __forceinline__
363cs_cuda_reduce_block_reduce(T *stmp,
364 size_t tid,
365 T *rd_block)
366{
367 R reducer;
368 __syncthreads();
369
370 /* Loop explicitely unrolled */
371
372 if (blockSize >= 1024) {
373 if (tid < 512) {
374 reducer.combine(stmp[tid], stmp[tid + 512]);
375 }
376 __syncthreads();
377 }
378 if (blockSize >= 512) {
379 if (tid < 256) {
380 reducer.combine(stmp[tid], stmp[tid + 256]);
381 }
382 __syncthreads();
383 }
384 if (blockSize >= 256) {
385 if (tid < 128) {
386 reducer.combine(stmp[tid], stmp[tid + 128]);
387 } __syncthreads();
388 }
389 if (blockSize >= 128) {
390 if (tid < 64) {
391 reducer.combine(stmp[tid], stmp[tid + 64]);
392 } __syncthreads();
393 }
394
395 if (tid < 32) {
396 cs_cuda_reduce_warp_reduce<blockSize, R>(stmp, tid);
397 }
398
399 // Output: rd_block for this block
400
401 if (tid == 0) rd_block[blockIdx.x] = stmp[0];
402}
403
404/*----------------------------------------------------------------------------*/
416/*----------------------------------------------------------------------------*/
417
418template <size_t blockSize, typename R, typename T>
419__global__ static void
420cs_cuda_reduce_single_block(size_t n,
421 T *g_idata,
422 T *g_odata)
423{
424 extern __shared__ int p_stmp[];
425 T *sdata = reinterpret_cast<T *>(p_stmp);
426 R reducer;
427
428 size_t tid = threadIdx.x;
429 T r_s[1];
430
431 reducer.identity(r_s[0]);
432 reducer.identity(sdata[tid]);
433
434 for (size_t i = threadIdx.x; i < n; i+= blockSize)
435 reducer.combine(r_s[0], g_idata[i]);
436
437 sdata[tid] = r_s[0];
438 __syncthreads();
439
440 for (size_t j = blockSize/2; j > CS_CUDA_WARP_SIZE; j /= 2) {
441 if (tid < j) {
442 reducer.combine(sdata[tid], sdata[tid + j]);
443 }
444 __syncthreads();
445 }
446
447 if (tid < 32) cs_cuda_reduce_warp_reduce<blockSize, R>(sdata, tid);
448 if (tid == 0) *g_odata = sdata[0];
449}
450
451#endif /* defined(__CUDACC__) */
452
453/*============================================================================
454 * Public function prototypes
455 *============================================================================*/
456
457/*----------------------------------------------------------------------------*/
458
459#endif /* CS_CUDA_REDUCE_H */
cs_base.h
cs_base_cuda.h
cs_defs.h
k
@ k
Definition: cs_field_pointer.h:68