9.2
general documentation
cs_blas_cuda.h
Go to the documentation of this file.
1#ifndef CS_BLAS_CUDA_H
2#define CS_BLAS_CUDA_H
3
4/*============================================================================
5 * BLAS (Basic Linear Algebra Subroutine) functions using CUDA.
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_blas_cuda_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_blas_cuda_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_blas_cuda_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_blas_cuda_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_blas_cuda_reduce_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_blas_cuda_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_blas_cuda_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#endif /* defined(__CUDACC__) */
317
318/*============================================================================
319 * Public function prototypes
320 *============================================================================*/
321
322/*----------------------------------------------------------------------------*/
323/*
324 * \brief Finalize CUDA BLAS API.
325 *
326 * This frees resources such as the cuBLAS handle, if used.
327 */
328/*----------------------------------------------------------------------------*/
329
330extern "C" void
331cs_blas_cuda_finalize(void);
332
333#if defined(__CUDACC__)
334
335/*----------------------------------------------------------------------------*/
336/*
337 * \brief Return the absolute sum of vector values using CUDA.
338 *
339 * \param[in] stream associated CUDA stream
340 * \param[in] n size of array x
341 * \param[in] x array of floating-point values (on device)
342 *
343 * \return sum of absolute array values
344 */
345/*----------------------------------------------------------------------------*/
346
347double
348cs_blas_cuda_asum(cudaStream_t stream,
349 cs_lnum_t n,
350 const cs_real_t x[]);
351
352/*----------------------------------------------------------------------------*/
353/*
354 * \brief Return the dot product of 2 vectors: x.y using CUDA.
355 *
356 * \param[in] stream associated CUDA stream
357 * \param[in] n size of arrays x and y
358 * \param[in] x array of floating-point values (on device)
359 * \param[in] y array of floating-point values (on device)
360 *
361 * \return dot product
362 */
363/*----------------------------------------------------------------------------*/
364
365double
366cs_blas_cuda_dot(cudaStream_t stream,
367 cs_lnum_t n,
368 const cs_real_t x[],
369 const cs_real_t y[]);
370
371#if defined(HAVE_CUBLAS)
372
373/*----------------------------------------------------------------------------*/
374/*
375 * \brief Return the absolute sum of vector values using cuBLAS.
376 *
377 * \param[in] stream associated CUDA stream
378 * \param[in] n size of arrays x and y
379 * \param[in] x array of floating-point values (on device)
380 *
381 * \return sum of absolute array values
382 */
383/*----------------------------------------------------------------------------*/
384
385double
386cs_blas_cublas_asum(cudaStream_t stream,
387 cs_lnum_t n,
388 const cs_real_t x[]);
389
390/*----------------------------------------------------------------------------*/
391/*
392 * \brief Return the dot product of 2 vectors: x.y using cuBLAS.
393 *
394 * \param[in] stream associated CUDA stream
395 * \param[in] n size of arrays x and y
396 * \param[in] x array of floating-point values (on device)
397 * \param[in] y array of floating-point values (on device)
398 *
399 * \return dot product
400 */
401/*----------------------------------------------------------------------------*/
402
403double
404cs_blas_cublas_dot(cudaStream_t stream,
405 cs_lnum_t n,
406 const cs_real_t x[],
407 const cs_real_t y[]);
408
409#endif /* defined(HAVE_CUBLAS) */
410
411/*----------------------------------------------------------------------------
412 * Compute x <- alpha.x
413 *
414 * This function may be set to use either cuBLAS or a local kernel.
415 *
416 * \param[in] stream associated CUDA stream
417 * \param[in] n number of elements
418 * \param[in] alpha constant value (on device)
419 * \param[in] x vector of elements (on device)
420 *----------------------------------------------------------------------------*/
421
422void
423cs_blas_cuda_scal(cudaStream_t stream,
424 cs_lnum_t n,
425 const cs_real_t *alpha,
426 cs_real_t *x);
427
428#endif /* defined(__CUDACC__) */
429
430/*----------------------------------------------------------------------------*/
431
432#endif /* CS_BLAS_CUDA_H */
cs_base.h
cs_base_cuda.h
cs_blas_cuda_finalize
void cs_blas_cuda_finalize(void)
cs_defs.h
cs_real_t
double cs_real_t
Floating-point value.
Definition: cs_defs.h:332
cs_lnum_t
int cs_lnum_t
local mesh entity id
Definition: cs_defs.h:325
k
@ k
Definition: cs_field_pointer.h:68