1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <cuda_runtime.h>
11 #include <math.h>
12 #include <stdbool.h>
13 #include <string.h>
14
15 #include "../cuda/ceed-cuda-common.h"
16 #include "ceed-cuda-ref.h"
17
18 //------------------------------------------------------------------------------
19 // Check if host/device sync is needed
20 //------------------------------------------------------------------------------
CeedVectorNeedSync_Cuda(const CeedVector vec,CeedMemType mem_type,bool * need_sync)21 static inline int CeedVectorNeedSync_Cuda(const CeedVector vec, CeedMemType mem_type, bool *need_sync) {
22 bool has_valid_array = false;
23 CeedVector_Cuda *impl;
24
25 CeedCallBackend(CeedVectorGetData(vec, &impl));
26 CeedCallBackend(CeedVectorHasValidArray(vec, &has_valid_array));
27 switch (mem_type) {
28 case CEED_MEM_HOST:
29 *need_sync = has_valid_array && !impl->h_array;
30 break;
31 case CEED_MEM_DEVICE:
32 *need_sync = has_valid_array && !impl->d_array;
33 break;
34 }
35 return CEED_ERROR_SUCCESS;
36 }
37
38 //------------------------------------------------------------------------------
39 // Sync host to device
40 //------------------------------------------------------------------------------
CeedVectorSyncH2D_Cuda(const CeedVector vec)41 static inline int CeedVectorSyncH2D_Cuda(const CeedVector vec) {
42 CeedSize length;
43 size_t bytes;
44 CeedVector_Cuda *impl;
45
46 CeedCallBackend(CeedVectorGetData(vec, &impl));
47
48 CeedCheck(impl->h_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND, "No valid host data to sync to device");
49
50 CeedCallBackend(CeedVectorGetLength(vec, &length));
51 bytes = length * sizeof(CeedScalar);
52 if (impl->d_array_borrowed) {
53 impl->d_array = impl->d_array_borrowed;
54 } else if (impl->d_array_owned) {
55 impl->d_array = impl->d_array_owned;
56 } else {
57 CeedCallCuda(CeedVectorReturnCeed(vec), cudaMalloc((void **)&impl->d_array_owned, bytes));
58 impl->d_array = impl->d_array_owned;
59 }
60 CeedCallCuda(CeedVectorReturnCeed(vec), cudaMemcpy(impl->d_array, impl->h_array, bytes, cudaMemcpyHostToDevice));
61 return CEED_ERROR_SUCCESS;
62 }
63
64 //------------------------------------------------------------------------------
65 // Sync device to host
66 //------------------------------------------------------------------------------
CeedVectorSyncD2H_Cuda(const CeedVector vec)67 static inline int CeedVectorSyncD2H_Cuda(const CeedVector vec) {
68 CeedSize length;
69 CeedVector_Cuda *impl;
70
71 CeedCallBackend(CeedVectorGetData(vec, &impl));
72
73 CeedCheck(impl->d_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND, "No valid device data to sync to host");
74
75 if (impl->h_array_borrowed) {
76 impl->h_array = impl->h_array_borrowed;
77 } else if (impl->h_array_owned) {
78 impl->h_array = impl->h_array_owned;
79 } else {
80 CeedSize length;
81
82 CeedCallBackend(CeedVectorGetLength(vec, &length));
83 CeedCallBackend(CeedCalloc(length, &impl->h_array_owned));
84 impl->h_array = impl->h_array_owned;
85 }
86
87 CeedCallBackend(CeedVectorGetLength(vec, &length));
88 size_t bytes = length * sizeof(CeedScalar);
89
90 CeedCallCuda(CeedVectorReturnCeed(vec), cudaMemcpy(impl->h_array, impl->d_array, bytes, cudaMemcpyDeviceToHost));
91 return CEED_ERROR_SUCCESS;
92 }
93
94 //------------------------------------------------------------------------------
95 // Sync arrays
96 //------------------------------------------------------------------------------
CeedVectorSyncArray_Cuda(const CeedVector vec,CeedMemType mem_type)97 static int CeedVectorSyncArray_Cuda(const CeedVector vec, CeedMemType mem_type) {
98 bool need_sync = false;
99
100 // Check whether device/host sync is needed
101 CeedCallBackend(CeedVectorNeedSync_Cuda(vec, mem_type, &need_sync));
102 if (!need_sync) return CEED_ERROR_SUCCESS;
103
104 switch (mem_type) {
105 case CEED_MEM_HOST:
106 return CeedVectorSyncD2H_Cuda(vec);
107 case CEED_MEM_DEVICE:
108 return CeedVectorSyncH2D_Cuda(vec);
109 }
110 return CEED_ERROR_UNSUPPORTED;
111 }
112
113 //------------------------------------------------------------------------------
114 // Set all pointers as invalid
115 //------------------------------------------------------------------------------
CeedVectorSetAllInvalid_Cuda(const CeedVector vec)116 static inline int CeedVectorSetAllInvalid_Cuda(const CeedVector vec) {
117 CeedVector_Cuda *impl;
118
119 CeedCallBackend(CeedVectorGetData(vec, &impl));
120 impl->h_array = NULL;
121 impl->d_array = NULL;
122 return CEED_ERROR_SUCCESS;
123 }
124
125 //------------------------------------------------------------------------------
126 // Check if CeedVector has any valid pointer
127 //------------------------------------------------------------------------------
CeedVectorHasValidArray_Cuda(const CeedVector vec,bool * has_valid_array)128 static inline int CeedVectorHasValidArray_Cuda(const CeedVector vec, bool *has_valid_array) {
129 CeedVector_Cuda *impl;
130
131 CeedCallBackend(CeedVectorGetData(vec, &impl));
132 *has_valid_array = impl->h_array || impl->d_array;
133 return CEED_ERROR_SUCCESS;
134 }
135
136 //------------------------------------------------------------------------------
137 // Check if has array of given type
138 //------------------------------------------------------------------------------
CeedVectorHasArrayOfType_Cuda(const CeedVector vec,CeedMemType mem_type,bool * has_array_of_type)139 static inline int CeedVectorHasArrayOfType_Cuda(const CeedVector vec, CeedMemType mem_type, bool *has_array_of_type) {
140 CeedVector_Cuda *impl;
141
142 CeedCallBackend(CeedVectorGetData(vec, &impl));
143 switch (mem_type) {
144 case CEED_MEM_HOST:
145 *has_array_of_type = impl->h_array_borrowed || impl->h_array_owned;
146 break;
147 case CEED_MEM_DEVICE:
148 *has_array_of_type = impl->d_array_borrowed || impl->d_array_owned;
149 break;
150 }
151 return CEED_ERROR_SUCCESS;
152 }
153
154 //------------------------------------------------------------------------------
155 // Check if has borrowed array of given type
156 //------------------------------------------------------------------------------
CeedVectorHasBorrowedArrayOfType_Cuda(const CeedVector vec,CeedMemType mem_type,bool * has_borrowed_array_of_type)157 static inline int CeedVectorHasBorrowedArrayOfType_Cuda(const CeedVector vec, CeedMemType mem_type, bool *has_borrowed_array_of_type) {
158 CeedVector_Cuda *impl;
159
160 CeedCallBackend(CeedVectorGetData(vec, &impl));
161 switch (mem_type) {
162 case CEED_MEM_HOST:
163 *has_borrowed_array_of_type = impl->h_array_borrowed;
164 break;
165 case CEED_MEM_DEVICE:
166 *has_borrowed_array_of_type = impl->d_array_borrowed;
167 break;
168 }
169 return CEED_ERROR_SUCCESS;
170 }
171
172 //------------------------------------------------------------------------------
173 // Set array from host
174 //------------------------------------------------------------------------------
CeedVectorSetArrayHost_Cuda(const CeedVector vec,const CeedCopyMode copy_mode,CeedScalar * array)175 static int CeedVectorSetArrayHost_Cuda(const CeedVector vec, const CeedCopyMode copy_mode, CeedScalar *array) {
176 CeedSize length;
177 CeedVector_Cuda *impl;
178
179 CeedCallBackend(CeedVectorGetData(vec, &impl));
180 CeedCallBackend(CeedVectorGetLength(vec, &length));
181
182 CeedCallBackend(CeedSetHostCeedScalarArray(array, copy_mode, length, (const CeedScalar **)&impl->h_array_owned,
183 (const CeedScalar **)&impl->h_array_borrowed, (const CeedScalar **)&impl->h_array));
184 return CEED_ERROR_SUCCESS;
185 }
186
187 //------------------------------------------------------------------------------
188 // Set array from device
189 //------------------------------------------------------------------------------
CeedVectorSetArrayDevice_Cuda(const CeedVector vec,const CeedCopyMode copy_mode,CeedScalar * array)190 static int CeedVectorSetArrayDevice_Cuda(const CeedVector vec, const CeedCopyMode copy_mode, CeedScalar *array) {
191 CeedSize length;
192 Ceed ceed;
193 CeedVector_Cuda *impl;
194
195 CeedCallBackend(CeedVectorGetCeed(vec, &ceed));
196 CeedCallBackend(CeedVectorGetData(vec, &impl));
197 CeedCallBackend(CeedVectorGetLength(vec, &length));
198
199 CeedCallBackend(CeedSetDeviceCeedScalarArray_Cuda(ceed, array, copy_mode, length, (const CeedScalar **)&impl->d_array_owned,
200 (const CeedScalar **)&impl->d_array_borrowed, (const CeedScalar **)&impl->d_array));
201 CeedCallBackend(CeedDestroy(&ceed));
202 return CEED_ERROR_SUCCESS;
203 }
204
205 //------------------------------------------------------------------------------
206 // Set the array used by a vector,
207 // freeing any previously allocated array if applicable
208 //------------------------------------------------------------------------------
CeedVectorSetArray_Cuda(const CeedVector vec,const CeedMemType mem_type,const CeedCopyMode copy_mode,CeedScalar * array)209 static int CeedVectorSetArray_Cuda(const CeedVector vec, const CeedMemType mem_type, const CeedCopyMode copy_mode, CeedScalar *array) {
210 CeedVector_Cuda *impl;
211
212 CeedCallBackend(CeedVectorGetData(vec, &impl));
213 CeedCallBackend(CeedVectorSetAllInvalid_Cuda(vec));
214 switch (mem_type) {
215 case CEED_MEM_HOST:
216 return CeedVectorSetArrayHost_Cuda(vec, copy_mode, array);
217 case CEED_MEM_DEVICE:
218 return CeedVectorSetArrayDevice_Cuda(vec, copy_mode, array);
219 }
220 return CEED_ERROR_UNSUPPORTED;
221 }
222
223 //------------------------------------------------------------------------------
224 // Copy host array to value strided
225 //------------------------------------------------------------------------------
CeedHostCopyStrided_Cuda(CeedScalar * h_array,CeedSize start,CeedSize stop,CeedSize step,CeedScalar * h_copy_array)226 static int CeedHostCopyStrided_Cuda(CeedScalar *h_array, CeedSize start, CeedSize stop, CeedSize step, CeedScalar *h_copy_array) {
227 for (CeedSize i = start; i < stop; i += step) h_copy_array[i] = h_array[i];
228 return CEED_ERROR_SUCCESS;
229 }
230
231 //------------------------------------------------------------------------------
232 // Copy device array to value strided (impl in .cu file)
233 //------------------------------------------------------------------------------
234 int CeedDeviceCopyStrided_Cuda(CeedScalar *d_array, CeedSize start, CeedSize stop, CeedSize step, CeedScalar *d_copy_array);
235
236 //------------------------------------------------------------------------------
237 // Copy a vector to a value strided
238 //------------------------------------------------------------------------------
CeedVectorCopyStrided_Cuda(CeedVector vec,CeedSize start,CeedSize stop,CeedSize step,CeedVector vec_copy)239 static int CeedVectorCopyStrided_Cuda(CeedVector vec, CeedSize start, CeedSize stop, CeedSize step, CeedVector vec_copy) {
240 CeedSize length;
241 CeedVector_Cuda *impl;
242
243 CeedCallBackend(CeedVectorGetData(vec, &impl));
244 {
245 CeedSize length_vec, length_copy;
246
247 CeedCallBackend(CeedVectorGetLength(vec, &length_vec));
248 CeedCallBackend(CeedVectorGetLength(vec_copy, &length_copy));
249 length = length_vec < length_copy ? length_vec : length_copy;
250 }
251 if (stop == -1) stop = length;
252 // Set value for synced device/host array
253 if (impl->d_array) {
254 CeedScalar *copy_array;
255
256 CeedCallBackend(CeedVectorGetArray(vec_copy, CEED_MEM_DEVICE, ©_array));
257 #if (CUDA_VERSION >= 12000)
258 cublasHandle_t handle;
259 Ceed ceed;
260
261 CeedCallBackend(CeedVectorGetCeed(vec, &ceed));
262 CeedCallBackend(CeedGetCublasHandle_Cuda(ceed, &handle));
263 #if defined(CEED_SCALAR_IS_FP32)
264 CeedCallCublas(ceed, cublasScopy_64(handle, (int64_t)(stop - start), impl->d_array + start, (int64_t)step, copy_array + start, (int64_t)step));
265 #else /* CEED_SCALAR */
266 CeedCallCublas(ceed, cublasDcopy_64(handle, (int64_t)(stop - start), impl->d_array + start, (int64_t)step, copy_array + start, (int64_t)step));
267 #endif /* CEED_SCALAR */
268 CeedCallBackend(CeedDestroy(&ceed));
269 #else /* CUDA_VERSION */
270 CeedCallBackend(CeedDeviceCopyStrided_Cuda(impl->d_array, start, stop, step, copy_array));
271 #endif /* CUDA_VERSION */
272 CeedCallBackend(CeedVectorRestoreArray(vec_copy, ©_array));
273 impl->h_array = NULL;
274 } else if (impl->h_array) {
275 CeedScalar *copy_array;
276
277 CeedCallBackend(CeedVectorGetArray(vec_copy, CEED_MEM_HOST, ©_array));
278 CeedCallBackend(CeedHostCopyStrided_Cuda(impl->h_array, start, stop, step, copy_array));
279 CeedCallBackend(CeedVectorRestoreArray(vec_copy, ©_array));
280 impl->d_array = NULL;
281 } else {
282 return CeedError(CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND, "CeedVector must have valid data set");
283 }
284 return CEED_ERROR_SUCCESS;
285 }
286
287 //------------------------------------------------------------------------------
288 // Set host array to value
289 //------------------------------------------------------------------------------
CeedHostSetValue_Cuda(CeedScalar * h_array,CeedSize length,CeedScalar val)290 static int CeedHostSetValue_Cuda(CeedScalar *h_array, CeedSize length, CeedScalar val) {
291 for (CeedSize i = 0; i < length; i++) h_array[i] = val;
292 return CEED_ERROR_SUCCESS;
293 }
294
295 //------------------------------------------------------------------------------
296 // Set device array to value (impl in .cu file)
297 //------------------------------------------------------------------------------
298 int CeedDeviceSetValue_Cuda(CeedScalar *d_array, CeedSize length, CeedScalar val);
299
300 //------------------------------------------------------------------------------
301 // Set a vector to a value
302 //------------------------------------------------------------------------------
CeedVectorSetValue_Cuda(CeedVector vec,CeedScalar val)303 static int CeedVectorSetValue_Cuda(CeedVector vec, CeedScalar val) {
304 CeedSize length;
305 CeedVector_Cuda *impl;
306
307 CeedCallBackend(CeedVectorGetData(vec, &impl));
308 CeedCallBackend(CeedVectorGetLength(vec, &length));
309 // Set value for synced device/host array
310 if (!impl->d_array && !impl->h_array) {
311 if (impl->d_array_borrowed) {
312 impl->d_array = impl->d_array_borrowed;
313 } else if (impl->h_array_borrowed) {
314 impl->h_array = impl->h_array_borrowed;
315 } else if (impl->d_array_owned) {
316 impl->d_array = impl->d_array_owned;
317 } else if (impl->h_array_owned) {
318 impl->h_array = impl->h_array_owned;
319 } else {
320 CeedCallBackend(CeedVectorSetArray(vec, CEED_MEM_DEVICE, CEED_COPY_VALUES, NULL));
321 }
322 }
323 if (impl->d_array) {
324 if (val == 0) {
325 CeedCallCuda(CeedVectorReturnCeed(vec), cudaMemset(impl->d_array, 0, length * sizeof(CeedScalar)));
326 } else {
327 CeedCallBackend(CeedDeviceSetValue_Cuda(impl->d_array, length, val));
328 }
329 impl->h_array = NULL;
330 } else if (impl->h_array) {
331 CeedCallBackend(CeedHostSetValue_Cuda(impl->h_array, length, val));
332 impl->d_array = NULL;
333 }
334 return CEED_ERROR_SUCCESS;
335 }
336
337 //------------------------------------------------------------------------------
338 // Set host array to value strided
339 //------------------------------------------------------------------------------
CeedHostSetValueStrided_Cuda(CeedScalar * h_array,CeedSize start,CeedSize stop,CeedSize step,CeedScalar val)340 static int CeedHostSetValueStrided_Cuda(CeedScalar *h_array, CeedSize start, CeedSize stop, CeedSize step, CeedScalar val) {
341 for (CeedSize i = start; i < stop; i += step) h_array[i] = val;
342 return CEED_ERROR_SUCCESS;
343 }
344
345 //------------------------------------------------------------------------------
346 // Set device array to value strided (impl in .cu file)
347 //------------------------------------------------------------------------------
348 int CeedDeviceSetValueStrided_Cuda(CeedScalar *d_array, CeedSize start, CeedSize stop, CeedSize step, CeedScalar val);
349
350 //------------------------------------------------------------------------------
351 // Set a vector to a value strided
352 //------------------------------------------------------------------------------
CeedVectorSetValueStrided_Cuda(CeedVector vec,CeedSize start,CeedSize stop,CeedSize step,CeedScalar val)353 static int CeedVectorSetValueStrided_Cuda(CeedVector vec, CeedSize start, CeedSize stop, CeedSize step, CeedScalar val) {
354 CeedSize length;
355 CeedVector_Cuda *impl;
356
357 CeedCallBackend(CeedVectorGetData(vec, &impl));
358 CeedCallBackend(CeedVectorGetLength(vec, &length));
359 // Set value for synced device/host array
360 if (stop == -1) stop = length;
361 if (impl->d_array) {
362 CeedCallBackend(CeedDeviceSetValueStrided_Cuda(impl->d_array, start, stop, step, val));
363 impl->h_array = NULL;
364 } else if (impl->h_array) {
365 CeedCallBackend(CeedHostSetValueStrided_Cuda(impl->h_array, start, stop, step, val));
366 impl->d_array = NULL;
367 } else {
368 return CeedError(CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND, "CeedVector must have valid data set");
369 }
370 return CEED_ERROR_SUCCESS;
371 }
372
373 //------------------------------------------------------------------------------
374 // Vector Take Array
375 //------------------------------------------------------------------------------
CeedVectorTakeArray_Cuda(CeedVector vec,CeedMemType mem_type,CeedScalar ** array)376 static int CeedVectorTakeArray_Cuda(CeedVector vec, CeedMemType mem_type, CeedScalar **array) {
377 CeedVector_Cuda *impl;
378
379 CeedCallBackend(CeedVectorGetData(vec, &impl));
380 // Sync array to requested mem_type
381 CeedCallBackend(CeedVectorSyncArray(vec, mem_type));
382 // Update pointer
383 switch (mem_type) {
384 case CEED_MEM_HOST:
385 (*array) = impl->h_array_borrowed;
386 impl->h_array_borrowed = NULL;
387 impl->h_array = NULL;
388 break;
389 case CEED_MEM_DEVICE:
390 (*array) = impl->d_array_borrowed;
391 impl->d_array_borrowed = NULL;
392 impl->d_array = NULL;
393 break;
394 }
395 return CEED_ERROR_SUCCESS;
396 }
397
398 //------------------------------------------------------------------------------
399 // Core logic for array syncronization for GetArray.
400 // If a different memory type is most up to date, this will perform a copy
401 //------------------------------------------------------------------------------
CeedVectorGetArrayCore_Cuda(const CeedVector vec,const CeedMemType mem_type,CeedScalar ** array)402 static int CeedVectorGetArrayCore_Cuda(const CeedVector vec, const CeedMemType mem_type, CeedScalar **array) {
403 CeedVector_Cuda *impl;
404
405 CeedCallBackend(CeedVectorGetData(vec, &impl));
406 // Sync array to requested mem_type
407 CeedCallBackend(CeedVectorSyncArray(vec, mem_type));
408 // Update pointer
409 switch (mem_type) {
410 case CEED_MEM_HOST:
411 *array = impl->h_array;
412 break;
413 case CEED_MEM_DEVICE:
414 *array = impl->d_array;
415 break;
416 }
417 return CEED_ERROR_SUCCESS;
418 }
419
420 //------------------------------------------------------------------------------
421 // Get read-only access to a vector via the specified mem_type
422 //------------------------------------------------------------------------------
CeedVectorGetArrayRead_Cuda(const CeedVector vec,const CeedMemType mem_type,const CeedScalar ** array)423 static int CeedVectorGetArrayRead_Cuda(const CeedVector vec, const CeedMemType mem_type, const CeedScalar **array) {
424 return CeedVectorGetArrayCore_Cuda(vec, mem_type, (CeedScalar **)array);
425 }
426
427 //------------------------------------------------------------------------------
428 // Get read/write access to a vector via the specified mem_type
429 //------------------------------------------------------------------------------
CeedVectorGetArray_Cuda(const CeedVector vec,const CeedMemType mem_type,CeedScalar ** array)430 static int CeedVectorGetArray_Cuda(const CeedVector vec, const CeedMemType mem_type, CeedScalar **array) {
431 CeedVector_Cuda *impl;
432
433 CeedCallBackend(CeedVectorGetData(vec, &impl));
434 CeedCallBackend(CeedVectorGetArrayCore_Cuda(vec, mem_type, array));
435 CeedCallBackend(CeedVectorSetAllInvalid_Cuda(vec));
436 switch (mem_type) {
437 case CEED_MEM_HOST:
438 impl->h_array = *array;
439 break;
440 case CEED_MEM_DEVICE:
441 impl->d_array = *array;
442 break;
443 }
444 return CEED_ERROR_SUCCESS;
445 }
446
447 //------------------------------------------------------------------------------
448 // Get write access to a vector via the specified mem_type
449 //------------------------------------------------------------------------------
CeedVectorGetArrayWrite_Cuda(const CeedVector vec,const CeedMemType mem_type,CeedScalar ** array)450 static int CeedVectorGetArrayWrite_Cuda(const CeedVector vec, const CeedMemType mem_type, CeedScalar **array) {
451 bool has_array_of_type = true;
452 CeedVector_Cuda *impl;
453
454 CeedCallBackend(CeedVectorGetData(vec, &impl));
455 CeedCallBackend(CeedVectorHasArrayOfType_Cuda(vec, mem_type, &has_array_of_type));
456 if (!has_array_of_type) {
457 // Allocate if array is not yet allocated
458 CeedCallBackend(CeedVectorSetArray(vec, mem_type, CEED_COPY_VALUES, NULL));
459 } else {
460 // Select dirty array
461 switch (mem_type) {
462 case CEED_MEM_HOST:
463 if (impl->h_array_borrowed) impl->h_array = impl->h_array_borrowed;
464 else impl->h_array = impl->h_array_owned;
465 break;
466 case CEED_MEM_DEVICE:
467 if (impl->d_array_borrowed) impl->d_array = impl->d_array_borrowed;
468 else impl->d_array = impl->d_array_owned;
469 }
470 }
471 return CeedVectorGetArray_Cuda(vec, mem_type, array);
472 }
473
474 //------------------------------------------------------------------------------
475 // Get the norm of a CeedVector
476 //------------------------------------------------------------------------------
CeedVectorNorm_Cuda(CeedVector vec,CeedNormType type,CeedScalar * norm)477 static int CeedVectorNorm_Cuda(CeedVector vec, CeedNormType type, CeedScalar *norm) {
478 Ceed ceed;
479 CeedSize length;
480 #if (CUDA_VERSION < 12000)
481 CeedSize num_calls;
482 #endif /* CUDA_VERSION */
483 const CeedScalar *d_array;
484 CeedVector_Cuda *impl;
485 cublasHandle_t handle;
486
487 CeedCallBackend(CeedVectorGetCeed(vec, &ceed));
488 CeedCallBackend(CeedVectorGetData(vec, &impl));
489 CeedCallBackend(CeedVectorGetLength(vec, &length));
490 CeedCallBackend(CeedGetCublasHandle_Cuda(ceed, &handle));
491
492 #if (CUDA_VERSION < 12000)
493 // With CUDA 12, we can use the 64-bit integer interface. Prior to that,
494 // we need to check if the vector is too long to handle with int32,
495 // and if so, divide it into subsections for repeated cuBLAS calls.
496 num_calls = length / INT_MAX;
497 if (length % INT_MAX > 0) num_calls += 1;
498 #endif /* CUDA_VERSION */
499
500 // Compute norm
501 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &d_array));
502 switch (type) {
503 case CEED_NORM_1: {
504 *norm = 0.0;
505 #if defined(CEED_SCALAR_IS_FP32)
506 #if (CUDA_VERSION >= 12000) // We have CUDA 12, and can use 64-bit integers
507 CeedCallCublas(ceed, cublasSasum_64(handle, (int64_t)length, (float *)d_array, 1, (float *)norm));
508 #else /* CUDA_VERSION */
509 float sub_norm = 0.0;
510 float *d_array_start;
511
512 for (CeedInt i = 0; i < num_calls; i++) {
513 d_array_start = (float *)d_array + (CeedSize)(i)*INT_MAX;
514 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
515 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
516
517 CeedCallCublas(ceed, cublasSasum(handle, (CeedInt)sub_length, (float *)d_array_start, 1, &sub_norm));
518 *norm += sub_norm;
519 }
520 #endif /* CUDA_VERSION */
521 #else /* CEED_SCALAR */
522 #if (CUDA_VERSION >= 12000)
523 CeedCallCublas(ceed, cublasDasum_64(handle, (int64_t)length, (double *)d_array, 1, (double *)norm));
524 #else /* CUDA_VERSION */
525 double sub_norm = 0.0;
526 double *d_array_start;
527
528 for (CeedInt i = 0; i < num_calls; i++) {
529 d_array_start = (double *)d_array + (CeedSize)(i)*INT_MAX;
530 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
531 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
532
533 CeedCallCublas(ceed, cublasDasum(handle, (CeedInt)sub_length, (double *)d_array_start, 1, &sub_norm));
534 *norm += sub_norm;
535 }
536 #endif /* CUDA_VERSION */
537 #endif /* CEED_SCALAR */
538 break;
539 }
540 case CEED_NORM_2: {
541 #if defined(CEED_SCALAR_IS_FP32)
542 #if (CUDA_VERSION >= 12000)
543 CeedCallCublas(ceed, cublasSnrm2_64(handle, (int64_t)length, (float *)d_array, 1, (float *)norm));
544 #else /* CUDA_VERSION */
545 float sub_norm = 0.0, norm_sum = 0.0;
546 float *d_array_start;
547
548 for (CeedInt i = 0; i < num_calls; i++) {
549 d_array_start = (float *)d_array + (CeedSize)(i)*INT_MAX;
550 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
551 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
552
553 CeedCallCublas(ceed, cublasSnrm2(handle, (CeedInt)sub_length, (float *)d_array_start, 1, &sub_norm));
554 norm_sum += sub_norm * sub_norm;
555 }
556 *norm = sqrt(norm_sum);
557 #endif /* CUDA_VERSION */
558 #else /* CEED_SCALAR */
559 #if (CUDA_VERSION >= 12000)
560 CeedCallCublas(ceed, cublasDnrm2_64(handle, (int64_t)length, (double *)d_array, 1, (double *)norm));
561 #else /* CUDA_VERSION */
562 double sub_norm = 0.0, norm_sum = 0.0;
563 double *d_array_start;
564
565 for (CeedInt i = 0; i < num_calls; i++) {
566 d_array_start = (double *)d_array + (CeedSize)(i)*INT_MAX;
567 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
568 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
569
570 CeedCallCublas(ceed, cublasDnrm2(handle, (CeedInt)sub_length, (double *)d_array_start, 1, &sub_norm));
571 norm_sum += sub_norm * sub_norm;
572 }
573 *norm = sqrt(norm_sum);
574 #endif /* CUDA_VERSION */
575 #endif /* CEED_SCALAR */
576 break;
577 }
578 case CEED_NORM_MAX: {
579 #if defined(CEED_SCALAR_IS_FP32)
580 #if (CUDA_VERSION >= 12000)
581 int64_t index;
582 CeedScalar norm_no_abs;
583
584 CeedCallCublas(ceed, cublasIsamax_64(handle, (int64_t)length, (float *)d_array, 1, &index));
585 CeedCallCuda(ceed, cudaMemcpy(&norm_no_abs, impl->d_array + index - 1, sizeof(CeedScalar), cudaMemcpyDeviceToHost));
586 *norm = fabs(norm_no_abs);
587 #else /* CUDA_VERSION */
588 CeedInt index;
589 float sub_max = 0.0, current_max = 0.0;
590 float *d_array_start;
591
592 for (CeedInt i = 0; i < num_calls; i++) {
593 d_array_start = (float *)d_array + (CeedSize)(i)*INT_MAX;
594 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
595 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
596
597 CeedCallCublas(ceed, cublasIsamax(handle, (CeedInt)sub_length, (float *)d_array_start, 1, &index));
598 CeedCallCuda(ceed, cudaMemcpy(&sub_max, d_array_start + index - 1, sizeof(CeedScalar), cudaMemcpyDeviceToHost));
599 if (fabs(sub_max) > current_max) current_max = fabs(sub_max);
600 }
601 *norm = current_max;
602 #endif /* CUDA_VERSION */
603 #else /* CEED_SCALAR */
604 #if (CUDA_VERSION >= 12000)
605 int64_t index;
606 CeedScalar norm_no_abs;
607
608 CeedCallCublas(ceed, cublasIdamax_64(handle, (int64_t)length, (double *)d_array, 1, &index));
609 CeedCallCuda(ceed, cudaMemcpy(&norm_no_abs, impl->d_array + index - 1, sizeof(CeedScalar), cudaMemcpyDeviceToHost));
610 *norm = fabs(norm_no_abs);
611 #else /* CUDA_VERSION */
612 CeedInt index;
613 double sub_max = 0.0, current_max = 0.0;
614 double *d_array_start;
615
616 for (CeedInt i = 0; i < num_calls; i++) {
617 d_array_start = (double *)d_array + (CeedSize)(i)*INT_MAX;
618 CeedSize remaining_length = length - (CeedSize)(i)*INT_MAX;
619 CeedInt sub_length = (i == num_calls - 1) ? (CeedInt)(remaining_length) : INT_MAX;
620
621 CeedCallCublas(ceed, cublasIdamax(handle, (CeedInt)sub_length, (double *)d_array_start, 1, &index));
622 CeedCallCuda(ceed, cudaMemcpy(&sub_max, d_array_start + index - 1, sizeof(CeedScalar), cudaMemcpyDeviceToHost));
623 if (fabs(sub_max) > current_max) current_max = fabs(sub_max);
624 }
625 *norm = current_max;
626 #endif /* CUDA_VERSION */
627 #endif /* CEED_SCALAR */
628 break;
629 }
630 }
631 CeedCallBackend(CeedVectorRestoreArrayRead(vec, &d_array));
632 CeedCallBackend(CeedDestroy(&ceed));
633 return CEED_ERROR_SUCCESS;
634 }
635
636 //------------------------------------------------------------------------------
637 // Take reciprocal of a vector on host
638 //------------------------------------------------------------------------------
CeedHostReciprocal_Cuda(CeedScalar * h_array,CeedSize length)639 static int CeedHostReciprocal_Cuda(CeedScalar *h_array, CeedSize length) {
640 for (CeedSize i = 0; i < length; i++) {
641 if (fabs(h_array[i]) > CEED_EPSILON) h_array[i] = 1. / h_array[i];
642 }
643 return CEED_ERROR_SUCCESS;
644 }
645
646 //------------------------------------------------------------------------------
647 // Take reciprocal of a vector on device (impl in .cu file)
648 //------------------------------------------------------------------------------
649 int CeedDeviceReciprocal_Cuda(CeedScalar *d_array, CeedSize length);
650
651 //------------------------------------------------------------------------------
652 // Take reciprocal of a vector
653 //------------------------------------------------------------------------------
CeedVectorReciprocal_Cuda(CeedVector vec)654 static int CeedVectorReciprocal_Cuda(CeedVector vec) {
655 CeedSize length;
656 CeedVector_Cuda *impl;
657
658 CeedCallBackend(CeedVectorGetData(vec, &impl));
659 CeedCallBackend(CeedVectorGetLength(vec, &length));
660 // Set value for synced device/host array
661 if (impl->d_array) CeedCallBackend(CeedDeviceReciprocal_Cuda(impl->d_array, length));
662 if (impl->h_array) CeedCallBackend(CeedHostReciprocal_Cuda(impl->h_array, length));
663 return CEED_ERROR_SUCCESS;
664 }
665
666 //------------------------------------------------------------------------------
667 // Compute x = alpha x on the host
668 //------------------------------------------------------------------------------
CeedHostScale_Cuda(CeedScalar * x_array,CeedScalar alpha,CeedSize length)669 static int CeedHostScale_Cuda(CeedScalar *x_array, CeedScalar alpha, CeedSize length) {
670 for (CeedSize i = 0; i < length; i++) x_array[i] *= alpha;
671 return CEED_ERROR_SUCCESS;
672 }
673
674 //------------------------------------------------------------------------------
675 // Compute x = alpha x on device (impl in .cu file)
676 //------------------------------------------------------------------------------
677 int CeedDeviceScale_Cuda(CeedScalar *x_array, CeedScalar alpha, CeedSize length);
678
679 //------------------------------------------------------------------------------
680 // Compute x = alpha x
681 //------------------------------------------------------------------------------
CeedVectorScale_Cuda(CeedVector x,CeedScalar alpha)682 static int CeedVectorScale_Cuda(CeedVector x, CeedScalar alpha) {
683 CeedSize length;
684 CeedVector_Cuda *impl;
685
686 CeedCallBackend(CeedVectorGetData(x, &impl));
687 CeedCallBackend(CeedVectorGetLength(x, &length));
688 // Set value for synced device/host array
689 if (impl->d_array) {
690 #if (CUDA_VERSION >= 12000)
691 cublasHandle_t handle;
692
693 CeedCallBackend(CeedGetCublasHandle_Cuda(CeedVectorReturnCeed(x), &handle));
694 #if defined(CEED_SCALAR_IS_FP32)
695 CeedCallCublas(CeedVectorReturnCeed(x), cublasSscal_64(handle, (int64_t)length, &alpha, impl->d_array, 1));
696 #else /* CEED_SCALAR */
697 CeedCallCublas(CeedVectorReturnCeed(x), cublasDscal_64(handle, (int64_t)length, &alpha, impl->d_array, 1));
698 #endif /* CEED_SCALAR */
699 #else /* CUDA_VERSION */
700 CeedCallBackend(CeedDeviceScale_Cuda(impl->d_array, alpha, length));
701 #endif /* CUDA_VERSION */
702 impl->h_array = NULL;
703 } else if (impl->h_array) {
704 CeedCallBackend(CeedHostScale_Cuda(impl->h_array, alpha, length));
705 impl->d_array = NULL;
706 }
707 return CEED_ERROR_SUCCESS;
708 }
709
710 //------------------------------------------------------------------------------
711 // Compute y = alpha x + y on the host
712 //------------------------------------------------------------------------------
CeedHostAXPY_Cuda(CeedScalar * y_array,CeedScalar alpha,CeedScalar * x_array,CeedSize length)713 static int CeedHostAXPY_Cuda(CeedScalar *y_array, CeedScalar alpha, CeedScalar *x_array, CeedSize length) {
714 for (CeedSize i = 0; i < length; i++) y_array[i] += alpha * x_array[i];
715 return CEED_ERROR_SUCCESS;
716 }
717
718 //------------------------------------------------------------------------------
719 // Compute y = alpha x + y on device (impl in .cu file)
720 //------------------------------------------------------------------------------
721 int CeedDeviceAXPY_Cuda(CeedScalar *y_array, CeedScalar alpha, CeedScalar *x_array, CeedSize length);
722
723 //------------------------------------------------------------------------------
724 // Compute y = alpha x + y
725 //------------------------------------------------------------------------------
CeedVectorAXPY_Cuda(CeedVector y,CeedScalar alpha,CeedVector x)726 static int CeedVectorAXPY_Cuda(CeedVector y, CeedScalar alpha, CeedVector x) {
727 CeedSize length;
728 CeedVector_Cuda *y_impl, *x_impl;
729
730 CeedCallBackend(CeedVectorGetData(y, &y_impl));
731 CeedCallBackend(CeedVectorGetData(x, &x_impl));
732 CeedCallBackend(CeedVectorGetLength(y, &length));
733 // Set value for synced device/host array
734 if (y_impl->d_array) {
735 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_DEVICE));
736 #if (CUDA_VERSION >= 12000)
737 cublasHandle_t handle;
738
739 CeedCallBackend(CeedGetCublasHandle_Cuda(CeedVectorReturnCeed(y), &handle));
740 #if defined(CEED_SCALAR_IS_FP32)
741 CeedCallCublas(CeedVectorReturnCeed(y), cublasSaxpy_64(handle, (int64_t)length, &alpha, x_impl->d_array, 1, y_impl->d_array, 1));
742 #else /* CEED_SCALAR */
743 CeedCallCublas(CeedVectorReturnCeed(y), cublasDaxpy_64(handle, (int64_t)length, &alpha, x_impl->d_array, 1, y_impl->d_array, 1));
744 #endif /* CEED_SCALAR */
745 #else /* CUDA_VERSION */
746 CeedCallBackend(CeedDeviceAXPY_Cuda(y_impl->d_array, alpha, x_impl->d_array, length));
747 #endif /* CUDA_VERSION */
748 y_impl->h_array = NULL;
749 } else if (y_impl->h_array) {
750 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_HOST));
751 CeedCallBackend(CeedHostAXPY_Cuda(y_impl->h_array, alpha, x_impl->h_array, length));
752 y_impl->d_array = NULL;
753 }
754 return CEED_ERROR_SUCCESS;
755 }
756
757 //------------------------------------------------------------------------------
758 // Compute y = alpha x + beta y on the host
759 //------------------------------------------------------------------------------
CeedHostAXPBY_Cuda(CeedScalar * y_array,CeedScalar alpha,CeedScalar beta,CeedScalar * x_array,CeedSize length)760 static int CeedHostAXPBY_Cuda(CeedScalar *y_array, CeedScalar alpha, CeedScalar beta, CeedScalar *x_array, CeedSize length) {
761 for (CeedSize i = 0; i < length; i++) y_array[i] = alpha * x_array[i] + beta * y_array[i];
762 return CEED_ERROR_SUCCESS;
763 }
764
765 //------------------------------------------------------------------------------
766 // Compute y = alpha x + beta y on device (impl in .cu file)
767 //------------------------------------------------------------------------------
768 int CeedDeviceAXPBY_Cuda(CeedScalar *y_array, CeedScalar alpha, CeedScalar beta, CeedScalar *x_array, CeedSize length);
769
770 //------------------------------------------------------------------------------
771 // Compute y = alpha x + beta y
772 //------------------------------------------------------------------------------
CeedVectorAXPBY_Cuda(CeedVector y,CeedScalar alpha,CeedScalar beta,CeedVector x)773 static int CeedVectorAXPBY_Cuda(CeedVector y, CeedScalar alpha, CeedScalar beta, CeedVector x) {
774 CeedSize length;
775 CeedVector_Cuda *y_impl, *x_impl;
776
777 CeedCallBackend(CeedVectorGetData(y, &y_impl));
778 CeedCallBackend(CeedVectorGetData(x, &x_impl));
779 CeedCallBackend(CeedVectorGetLength(y, &length));
780 // Set value for synced device/host array
781 if (y_impl->d_array) {
782 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_DEVICE));
783 CeedCallBackend(CeedDeviceAXPBY_Cuda(y_impl->d_array, alpha, beta, x_impl->d_array, length));
784 }
785 if (y_impl->h_array) {
786 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_HOST));
787 CeedCallBackend(CeedHostAXPBY_Cuda(y_impl->h_array, alpha, beta, x_impl->h_array, length));
788 }
789 return CEED_ERROR_SUCCESS;
790 }
791
792 //------------------------------------------------------------------------------
793 // Compute the pointwise multiplication w = x .* y on the host
794 //------------------------------------------------------------------------------
CeedHostPointwiseMult_Cuda(CeedScalar * w_array,CeedScalar * x_array,CeedScalar * y_array,CeedSize length)795 static int CeedHostPointwiseMult_Cuda(CeedScalar *w_array, CeedScalar *x_array, CeedScalar *y_array, CeedSize length) {
796 for (CeedSize i = 0; i < length; i++) w_array[i] = x_array[i] * y_array[i];
797 return CEED_ERROR_SUCCESS;
798 }
799
800 //------------------------------------------------------------------------------
801 // Compute the pointwise multiplication w = x .* y on device (impl in .cu file)
802 //------------------------------------------------------------------------------
803 int CeedDevicePointwiseMult_Cuda(CeedScalar *w_array, CeedScalar *x_array, CeedScalar *y_array, CeedSize length);
804
805 //------------------------------------------------------------------------------
806 // Compute the pointwise multiplication w = x .* y
807 //------------------------------------------------------------------------------
CeedVectorPointwiseMult_Cuda(CeedVector w,CeedVector x,CeedVector y)808 static int CeedVectorPointwiseMult_Cuda(CeedVector w, CeedVector x, CeedVector y) {
809 CeedSize length;
810 CeedVector_Cuda *w_impl, *x_impl, *y_impl;
811
812 CeedCallBackend(CeedVectorGetData(w, &w_impl));
813 CeedCallBackend(CeedVectorGetData(x, &x_impl));
814 CeedCallBackend(CeedVectorGetData(y, &y_impl));
815 CeedCallBackend(CeedVectorGetLength(w, &length));
816 // Set value for synced device/host array
817 if (!w_impl->d_array && !w_impl->h_array) {
818 CeedCallBackend(CeedVectorSetValue(w, 0.0));
819 }
820 if (w_impl->d_array) {
821 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_DEVICE));
822 CeedCallBackend(CeedVectorSyncArray(y, CEED_MEM_DEVICE));
823 CeedCallBackend(CeedDevicePointwiseMult_Cuda(w_impl->d_array, x_impl->d_array, y_impl->d_array, length));
824 }
825 if (w_impl->h_array) {
826 CeedCallBackend(CeedVectorSyncArray(x, CEED_MEM_HOST));
827 CeedCallBackend(CeedVectorSyncArray(y, CEED_MEM_HOST));
828 CeedCallBackend(CeedHostPointwiseMult_Cuda(w_impl->h_array, x_impl->h_array, y_impl->h_array, length));
829 }
830 return CEED_ERROR_SUCCESS;
831 }
832
833 //------------------------------------------------------------------------------
834 // Destroy the vector
835 //------------------------------------------------------------------------------
CeedVectorDestroy_Cuda(const CeedVector vec)836 static int CeedVectorDestroy_Cuda(const CeedVector vec) {
837 CeedVector_Cuda *impl;
838
839 CeedCallBackend(CeedVectorGetData(vec, &impl));
840 CeedCallCuda(CeedVectorReturnCeed(vec), cudaFree(impl->d_array_owned));
841 CeedCallBackend(CeedFree(&impl->h_array_owned));
842 CeedCallBackend(CeedFree(&impl));
843 return CEED_ERROR_SUCCESS;
844 }
845
846 //------------------------------------------------------------------------------
847 // Create a vector of the specified length (does not allocate memory)
848 //------------------------------------------------------------------------------
CeedVectorCreate_Cuda(CeedSize n,CeedVector vec)849 int CeedVectorCreate_Cuda(CeedSize n, CeedVector vec) {
850 CeedVector_Cuda *impl;
851 Ceed ceed;
852
853 CeedCallBackend(CeedVectorGetCeed(vec, &ceed));
854 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "HasValidArray", CeedVectorHasValidArray_Cuda));
855 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "HasBorrowedArrayOfType", CeedVectorHasBorrowedArrayOfType_Cuda));
856 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "SetArray", CeedVectorSetArray_Cuda));
857 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "TakeArray", CeedVectorTakeArray_Cuda));
858 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "CopyStrided", CeedVectorCopyStrided_Cuda));
859 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "SetValue", CeedVectorSetValue_Cuda));
860 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "SetValueStrided", CeedVectorSetValueStrided_Cuda));
861 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "SyncArray", CeedVectorSyncArray_Cuda));
862 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "GetArray", CeedVectorGetArray_Cuda));
863 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "GetArrayRead", CeedVectorGetArrayRead_Cuda));
864 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "GetArrayWrite", CeedVectorGetArrayWrite_Cuda));
865 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Norm", CeedVectorNorm_Cuda));
866 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Reciprocal", CeedVectorReciprocal_Cuda));
867 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Scale", CeedVectorScale_Cuda));
868 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "AXPY", CeedVectorAXPY_Cuda));
869 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "AXPBY", CeedVectorAXPBY_Cuda));
870 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "PointwiseMult", CeedVectorPointwiseMult_Cuda));
871 CeedCallBackend(CeedSetBackendFunction(ceed, "Vector", vec, "Destroy", CeedVectorDestroy_Cuda));
872 CeedCallBackend(CeedDestroy(&ceed));
873 CeedCallBackend(CeedCalloc(1, &impl));
874 CeedCallBackend(CeedVectorSetData(vec, impl));
875 return CEED_ERROR_SUCCESS;
876 }
877
878 //------------------------------------------------------------------------------
879