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