1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <ceed-impl.h> 20 #include <math.h> 21 #include <stdint.h> 22 #include <stdio.h> 23 24 /// @file 25 /// Implementation of public CeedVector interfaces 26 27 /// @cond DOXYGEN_SKIP 28 static struct CeedVector_private ceed_vector_active; 29 static struct CeedVector_private ceed_vector_none; 30 /// @endcond 31 32 /// @addtogroup CeedVectorUser 33 /// @{ 34 35 /// Indicate that vector will be provided as an explicit argument to 36 /// CeedOperatorApply(). 37 const CeedVector CEED_VECTOR_ACTIVE = &ceed_vector_active; 38 39 /// Indicate that no vector is applicable (i.e., for @ref CEED_EVAL_WEIGHT). 40 const CeedVector CEED_VECTOR_NONE = &ceed_vector_none; 41 42 /// @} 43 44 /// ---------------------------------------------------------------------------- 45 /// CeedVector Backend API 46 /// ---------------------------------------------------------------------------- 47 /// @addtogroup CeedVectorBackend 48 /// @{ 49 50 /** 51 @brief Get the state of a CeedVector 52 53 @param vec CeedVector to retrieve state 54 @param[out] state Variable to store state 55 56 @return An error code: 0 - success, otherwise - failure 57 58 @ref Backend 59 **/ 60 int CeedVectorGetState(CeedVector vec, uint64_t *state) { 61 *state = vec->state; 62 return CEED_ERROR_SUCCESS; 63 } 64 65 /** 66 @brief Add a reference to a CeedVector 67 68 @param[out] vec CeedVector to increment reference counter 69 70 @return An error code: 0 - success, otherwise - failure 71 72 @ref Backend 73 **/ 74 int CeedVectorAddReference(CeedVector vec) { 75 vec->ref_count++; 76 return CEED_ERROR_SUCCESS; 77 } 78 79 /** 80 @brief Get the backend data of a CeedVector 81 82 @param vec CeedVector to retrieve state 83 @param[out] data Variable to store data 84 85 @return An error code: 0 - success, otherwise - failure 86 87 @ref Backend 88 **/ 89 int CeedVectorGetData(CeedVector vec, void *data) { 90 *(void **)data = vec->data; 91 return CEED_ERROR_SUCCESS; 92 } 93 94 /** 95 @brief Set the backend data of a CeedVector 96 97 @param[out] vec CeedVector to retrieve state 98 @param data Data to set 99 100 @return An error code: 0 - success, otherwise - failure 101 102 @ref Backend 103 **/ 104 int CeedVectorSetData(CeedVector vec, void *data) { 105 vec->data = data; 106 return CEED_ERROR_SUCCESS; 107 } 108 109 /** 110 @brief Increment the reference counter for a CeedVector 111 112 @param vec CeedVector to increment the reference counter 113 114 @return An error code: 0 - success, otherwise - failure 115 116 @ref Backend 117 **/ 118 int CeedVectorReference(CeedVector vec) { 119 vec->ref_count++; 120 return CEED_ERROR_SUCCESS; 121 } 122 123 /// @} 124 125 /// ---------------------------------------------------------------------------- 126 /// CeedVector Public API 127 /// ---------------------------------------------------------------------------- 128 /// @addtogroup CeedVectorUser 129 /// @{ 130 131 /** 132 @brief Create a CeedVector of the specified length (does not allocate memory) 133 134 @param ceed Ceed object where the CeedVector will be created 135 @param length Length of vector 136 @param[out] vec Address of the variable where the newly created 137 CeedVector will be stored 138 139 @return An error code: 0 - success, otherwise - failure 140 141 @ref User 142 **/ 143 int CeedVectorCreate(Ceed ceed, CeedInt length, CeedVector *vec) { 144 int ierr; 145 146 if (!ceed->VectorCreate) { 147 Ceed delegate; 148 ierr = CeedGetObjectDelegate(ceed, &delegate, "Vector"); CeedChk(ierr); 149 150 if (!delegate) 151 // LCOV_EXCL_START 152 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 153 "Backend does not support VectorCreate"); 154 // LCOV_EXCL_STOP 155 156 ierr = CeedVectorCreate(delegate, length, vec); CeedChk(ierr); 157 return CEED_ERROR_SUCCESS; 158 } 159 160 ierr = CeedCalloc(1, vec); CeedChk(ierr); 161 (*vec)->ceed = ceed; 162 ierr = CeedReference(ceed); CeedChk(ierr); 163 (*vec)->ref_count = 1; 164 (*vec)->length = length; 165 (*vec)->state = 0; 166 ierr = ceed->VectorCreate(length, *vec); CeedChk(ierr); 167 return CEED_ERROR_SUCCESS; 168 } 169 170 /** 171 @brief Copy the pointer to a CeedVector. Both pointers should 172 be destroyed with `CeedVectorDestroy()`; 173 Note: If `*vec_copy` is non-NULL, then it is assumed that 174 `*vec_copy` is a pointer to a CeedVector. This 175 CeedVector will be destroyed if `*vec_copy` is the only 176 reference to this CeedVector. 177 178 @param vec CeedVector to copy reference to 179 @param[out] vec_copy Variable to store copied reference 180 181 @return An error code: 0 - success, otherwise - failure 182 183 @ref User 184 **/ 185 int CeedVectorReferenceCopy(CeedVector vec, CeedVector *vec_copy) { 186 int ierr; 187 188 ierr = CeedVectorReference(vec); CeedChk(ierr); 189 ierr = CeedVectorDestroy(vec_copy); CeedChk(ierr); 190 *vec_copy = vec; 191 return CEED_ERROR_SUCCESS; 192 } 193 194 /** 195 @brief Set the array used by a CeedVector, freeing any previously allocated 196 array if applicable. The backend may copy values to a different 197 memtype, such as during @ref CeedOperatorApply(). 198 See also @ref CeedVectorSyncArray() and @ref CeedVectorTakeArray(). 199 200 @param vec CeedVector 201 @param mem_type Memory type of the array being passed 202 @param copy_mode Copy mode for the array 203 @param array Array to be used, or NULL with @ref CEED_COPY_VALUES to have the 204 library allocate 205 206 @return An error code: 0 - success, otherwise - failure 207 208 @ref User 209 **/ 210 int CeedVectorSetArray(CeedVector vec, CeedMemType mem_type, 211 CeedCopyMode copy_mode, 212 CeedScalar *array) { 213 int ierr; 214 215 if (!vec->SetArray) 216 // LCOV_EXCL_START 217 return CeedError(vec->ceed, CEED_ERROR_UNSUPPORTED, 218 "Backend does not support VectorSetArray"); 219 // LCOV_EXCL_STOP 220 221 if (vec->state % 2 == 1) 222 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 223 "Cannot grant CeedVector array access, the " 224 "access lock is already in use"); 225 226 if (vec->num_readers > 0) 227 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 228 "Cannot grant CeedVector array access, a " 229 "process has read access"); 230 231 ierr = vec->SetArray(vec, mem_type, copy_mode, array); CeedChk(ierr); 232 vec->state += 2; 233 return CEED_ERROR_SUCCESS; 234 } 235 236 /** 237 @brief Set the CeedVector to a constant value 238 239 @param vec CeedVector 240 @param[in] value Value to be used 241 242 @return An error code: 0 - success, otherwise - failure 243 244 @ref User 245 **/ 246 int CeedVectorSetValue(CeedVector vec, CeedScalar value) { 247 int ierr; 248 249 if (vec->state % 2 == 1) 250 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 251 "Cannot grant CeedVector array access, the " 252 "access lock is already in use"); 253 254 if (vec->SetValue) { 255 ierr = vec->SetValue(vec, value); CeedChk(ierr); 256 } else { 257 CeedScalar *array; 258 ierr = CeedVectorGetArray(vec, CEED_MEM_HOST, &array); CeedChk(ierr); 259 for (int i=0; i<vec->length; i++) array[i] = value; 260 ierr = CeedVectorRestoreArray(vec, &array); CeedChk(ierr); 261 } 262 vec->state += 2; 263 return CEED_ERROR_SUCCESS; 264 } 265 266 /** 267 @brief Sync the CeedVector to a specified memtype. This function is used to 268 force synchronization of arrays set with @ref CeedVectorSetArray(). 269 If the requested memtype is already synchronized, this function 270 results in a no-op. 271 272 @param vec CeedVector 273 @param mem_type Memtype to be synced 274 275 @return An error code: 0 - success, otherwise - failure 276 277 @ref User 278 **/ 279 int CeedVectorSyncArray(CeedVector vec, CeedMemType mem_type) { 280 int ierr; 281 282 if (vec->state % 2 == 1) 283 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 284 "Cannot sync CeedVector, the access lock is " 285 "already in use"); 286 287 if (vec->SyncArray) { 288 ierr = vec->SyncArray(vec, mem_type); CeedChk(ierr); 289 } else { 290 const CeedScalar *array; 291 ierr = CeedVectorGetArrayRead(vec, mem_type, &array); CeedChk(ierr); 292 ierr = CeedVectorRestoreArrayRead(vec, &array); CeedChk(ierr); 293 } 294 return CEED_ERROR_SUCCESS; 295 } 296 297 /** 298 @brief Take ownership of the CeedVector array and remove the array from the 299 CeedVector. The caller is responsible for managing and freeing 300 the array. 301 302 @param vec CeedVector 303 @param mem_type Memory type on which to take the array. If the backend 304 uses a different memory type, this will perform a copy. 305 @param[out] array Array on memory type mem_type, or NULL if array pointer is 306 not required 307 308 @return An error code: 0 - success, otherwise - failure 309 310 @ref User 311 **/ 312 int CeedVectorTakeArray(CeedVector vec, CeedMemType mem_type, 313 CeedScalar **array) { 314 int ierr; 315 316 if (vec->state % 2 == 1) 317 // LCOV_EXCL_START 318 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 319 "Cannot take CeedVector array, the access " 320 "lock is already in use"); 321 // LCOV_EXCL_STOP 322 if (vec->num_readers > 0) 323 // LCOV_EXCL_START 324 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 325 "Cannot take CeedVector array, a process " 326 "has read access"); 327 // LCOV_EXCL_STOP 328 329 CeedScalar *temp_array = NULL; 330 ierr = vec->TakeArray(vec, mem_type, &temp_array); CeedChk(ierr); 331 if (array) (*array) = temp_array; 332 return CEED_ERROR_SUCCESS; 333 } 334 335 /** 336 @brief Get read/write access to a CeedVector via the specified memory type. 337 Restore access with @ref CeedVectorRestoreArray(). 338 339 @param vec CeedVector to access 340 @param mem_type Memory type on which to access the array. If the backend 341 uses a different memory type, this will perform a copy. 342 @param[out] array Array on memory type mem_type 343 344 @note The CeedVectorGetArray* and CeedVectorRestoreArray* functions provide 345 access to array pointers in the desired memory space. Pairing get/restore 346 allows the Vector to track access, thus knowing if norms or other 347 operations may need to be recomputed. 348 349 @return An error code: 0 - success, otherwise - failure 350 351 @ref User 352 **/ 353 int CeedVectorGetArray(CeedVector vec, CeedMemType mem_type, 354 CeedScalar **array) { 355 int ierr; 356 357 if (!vec->GetArray) 358 // LCOV_EXCL_START 359 return CeedError(vec->ceed, CEED_ERROR_UNSUPPORTED, 360 "Backend does not support GetArray"); 361 // LCOV_EXCL_STOP 362 363 if (vec->state % 2 == 1) 364 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 365 "Cannot grant CeedVector array access, the " 366 "access lock is already in use"); 367 368 if (vec->num_readers > 0) 369 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 370 "Cannot grant CeedVector array access, a " 371 "process has read access"); 372 373 ierr = vec->GetArray(vec, mem_type, array); CeedChk(ierr); 374 vec->state += 1; 375 return CEED_ERROR_SUCCESS; 376 } 377 378 /** 379 @brief Get read-only access to a CeedVector via the specified memory type. 380 Restore access with @ref CeedVectorRestoreArrayRead(). 381 382 @param vec CeedVector to access 383 @param mem_type Memory type on which to access the array. If the backend 384 uses a different memory type, this will perform a copy 385 (possibly cached). 386 @param[out] array Array on memory type mem_type 387 388 @return An error code: 0 - success, otherwise - failure 389 390 @ref User 391 **/ 392 int CeedVectorGetArrayRead(CeedVector vec, CeedMemType mem_type, 393 const CeedScalar **array) { 394 int ierr; 395 396 if (!vec->GetArrayRead) 397 // LCOV_EXCL_START 398 return CeedError(vec->ceed, CEED_ERROR_UNSUPPORTED, 399 "Backend does not support GetArrayRead"); 400 // LCOV_EXCL_STOP 401 402 if (vec->state % 2 == 1) 403 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 404 "Cannot grant CeedVector read-only array " 405 "access, the access lock is already in use"); 406 407 ierr = vec->GetArrayRead(vec, mem_type, array); CeedChk(ierr); 408 vec->num_readers++; 409 return CEED_ERROR_SUCCESS; 410 } 411 412 /** 413 @brief Restore an array obtained using @ref CeedVectorGetArray() 414 415 @param vec CeedVector to restore 416 @param array Array of vector data 417 418 @return An error code: 0 - success, otherwise - failure 419 420 @ref User 421 **/ 422 int CeedVectorRestoreArray(CeedVector vec, CeedScalar **array) { 423 int ierr; 424 425 if (!vec->RestoreArray) 426 // LCOV_EXCL_START 427 return CeedError(vec->ceed, CEED_ERROR_UNSUPPORTED, 428 "Backend does not support RestoreArray"); 429 // LCOV_EXCL_STOP 430 431 if (vec->state % 2 != 1) 432 return CeedError(vec->ceed, CEED_ERROR_ACCESS, 433 "Cannot restore CeedVector array access, " 434 "access was not granted"); 435 436 ierr = vec->RestoreArray(vec); CeedChk(ierr); 437 *array = NULL; 438 vec->state += 1; 439 return CEED_ERROR_SUCCESS; 440 } 441 442 /** 443 @brief Restore an array obtained using @ref CeedVectorGetArrayRead() 444 445 @param vec CeedVector to restore 446 @param array Array of vector data 447 448 @return An error code: 0 - success, otherwise - failure 449 450 @ref User 451 **/ 452 int CeedVectorRestoreArrayRead(CeedVector vec, const CeedScalar **array) { 453 int ierr; 454 455 if (!vec->RestoreArrayRead) 456 // LCOV_EXCL_START 457 return CeedError(vec->ceed, CEED_ERROR_UNSUPPORTED, 458 "Backend does not support RestoreArrayRead"); 459 // LCOV_EXCL_STOP 460 461 ierr = vec->RestoreArrayRead(vec); CeedChk(ierr); 462 *array = NULL; 463 vec->num_readers--; 464 return CEED_ERROR_SUCCESS; 465 } 466 467 /** 468 @brief Get the norm of a CeedVector. 469 470 Note: This operation is local to the CeedVector. This function will likely 471 not provide the desired results for the norm of the libCEED portion 472 of a parallel vector or a CeedVector with duplicated or hanging nodes. 473 474 @param vec CeedVector to retrieve maximum value 475 @param norm_type Norm type @ref CEED_NORM_1, @ref CEED_NORM_2, or @ref CEED_NORM_MAX 476 @param[out] norm Variable to store norm value 477 478 @return An error code: 0 - success, otherwise - failure 479 480 @ref User 481 **/ 482 int CeedVectorNorm(CeedVector vec, CeedNormType norm_type, CeedScalar *norm) { 483 int ierr; 484 485 // Backend impl for GPU, if added 486 if (vec->Norm) { 487 ierr = vec->Norm(vec, norm_type, norm); CeedChk(ierr); 488 return CEED_ERROR_SUCCESS; 489 } 490 491 const CeedScalar *array; 492 ierr = CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &array); CeedChk(ierr); 493 494 *norm = 0.; 495 switch (norm_type) { 496 case CEED_NORM_1: 497 for (int i=0; i<vec->length; i++) { 498 *norm += fabs(array[i]); 499 } 500 break; 501 case CEED_NORM_2: 502 for (int i=0; i<vec->length; i++) { 503 *norm += fabs(array[i])*fabs(array[i]); 504 } 505 break; 506 case CEED_NORM_MAX: 507 for (int i=0; i<vec->length; i++) { 508 const CeedScalar abs_v_i = fabs(array[i]); 509 *norm = *norm > abs_v_i ? *norm : abs_v_i; 510 } 511 } 512 if (norm_type == CEED_NORM_2) 513 *norm = sqrt(*norm); 514 515 ierr = CeedVectorRestoreArrayRead(vec, &array); CeedChk(ierr); 516 return CEED_ERROR_SUCCESS; 517 } 518 519 /** 520 @brief Compute x = alpha x 521 522 @param[in,out] x vector for scaling 523 @param[in] alpha scaling factor 524 525 @return An error code: 0 - success, otherwise - failure 526 527 @ref User 528 **/ 529 int CeedVectorScale(CeedVector x, CeedScalar alpha) { 530 int ierr; 531 CeedScalar *x_array; 532 CeedInt n_x; 533 534 ierr = CeedVectorGetLength(x, &n_x); CeedChk(ierr); 535 536 // Backend implementation 537 if (x->Scale) 538 return x->Scale(x, alpha); 539 540 // Default implementation 541 ierr = CeedVectorGetArray(x, CEED_MEM_HOST, &x_array); CeedChk(ierr); 542 for (CeedInt i=0; i<n_x; i++) 543 x_array[i] *= alpha; 544 ierr = CeedVectorRestoreArray(x, &x_array); CeedChk(ierr); 545 546 return CEED_ERROR_SUCCESS; 547 } 548 549 /** 550 @brief Compute y = alpha x + y 551 552 @param[in,out] y target vector for sum 553 @param[in] alpha scaling factor 554 @param[in] x second vector, must be different than y 555 556 @return An error code: 0 - success, otherwise - failure 557 558 @ref User 559 **/ 560 int CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x) { 561 int ierr; 562 CeedScalar *y_array; 563 CeedScalar const *x_array; 564 CeedInt n_x, n_y; 565 566 ierr = CeedVectorGetLength(y, &n_y); CeedChk(ierr); 567 ierr = CeedVectorGetLength(x, &n_x); CeedChk(ierr); 568 if (n_x != n_y) 569 // LCOV_EXCL_START 570 return CeedError(y->ceed, CEED_ERROR_UNSUPPORTED, 571 "Cannot add vector of different lengths"); 572 // LCOV_EXCL_STOP 573 if (x == y) 574 // LCOV_EXCL_START 575 return CeedError(y->ceed, CEED_ERROR_UNSUPPORTED, 576 "Cannot use same vector for x and y in CeedVectorAXPY"); 577 // LCOV_EXCL_STOP 578 579 Ceed ceed_parent_x, ceed_parent_y; 580 ierr = CeedGetParent(x->ceed, &ceed_parent_x); CeedChk(ierr); 581 ierr = CeedGetParent(y->ceed, &ceed_parent_y); CeedChk(ierr); 582 if (ceed_parent_x != ceed_parent_y) 583 // LCOV_EXCL_START 584 return CeedError(y->ceed, CEED_ERROR_INCOMPATIBLE, 585 "Vectors x and y must be created by the same Ceed context"); 586 // LCOV_EXCL_STOP 587 588 // Backend implementation 589 if (y->AXPY) { 590 ierr = y->AXPY(y, alpha, x); CeedChk(ierr); 591 return CEED_ERROR_SUCCESS; 592 } 593 594 // Default implementation 595 ierr = CeedVectorGetArray(y, CEED_MEM_HOST, &y_array); CeedChk(ierr); 596 ierr = CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); CeedChk(ierr); 597 598 for (CeedInt i=0; i<n_y; i++) 599 y_array[i] += alpha * x_array[i]; 600 601 ierr = CeedVectorRestoreArray(y, &y_array); CeedChk(ierr); 602 ierr = CeedVectorRestoreArrayRead(x, &x_array); CeedChk(ierr); 603 604 return CEED_ERROR_SUCCESS; 605 } 606 607 /** 608 @brief Compute the pointwise multiplication w = x .* y. Any 609 subset of x, y, and w may be the same vector. 610 611 @param[out] w target vector for the product 612 @param[in] x first vector for product 613 @param[in] y second vector for the product 614 615 @return An error code: 0 - success, otherwise - failure 616 617 @ ref User 618 **/ 619 int CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y) { 620 int ierr; 621 CeedScalar *w_array; 622 CeedScalar const *x_array, *y_array; 623 CeedInt n_w, n_x, n_y; 624 625 ierr = CeedVectorGetLength(w, &n_w); CeedChk(ierr); 626 ierr = CeedVectorGetLength(x, &n_x); CeedChk(ierr); 627 ierr = CeedVectorGetLength(y, &n_y); CeedChk(ierr); 628 if (n_w != n_x || n_w != n_y) 629 // LCOV_EXCL_START 630 return CeedError(w->ceed, CEED_ERROR_UNSUPPORTED, 631 "Cannot multiply vectors of different lengths"); 632 // LCOV_EXCL_STOP 633 634 Ceed ceed_parent_w, ceed_parent_x, ceed_parent_y; 635 ierr = CeedGetParent(w->ceed, &ceed_parent_w); CeedChk(ierr); 636 ierr = CeedGetParent(x->ceed, &ceed_parent_x); CeedChk(ierr); 637 ierr = CeedGetParent(y->ceed, &ceed_parent_y); CeedChk(ierr); 638 if ((ceed_parent_w != ceed_parent_y) || 639 (ceed_parent_w != ceed_parent_y)) 640 // LCOV_EXCL_START 641 return CeedError(w->ceed, CEED_ERROR_INCOMPATIBLE, 642 "Vectors w, x, and y must be created by the same Ceed context"); 643 // LCOV_EXCL_STOP 644 645 // Backend implementation 646 if (w->PointwiseMult) { 647 ierr = w->PointwiseMult(w, x, y); CeedChk(ierr); 648 return CEED_ERROR_SUCCESS; 649 } 650 651 // Default implementation 652 ierr = CeedVectorGetArray(w, CEED_MEM_HOST, &w_array); CeedChk(ierr); 653 if (x != w) { 654 ierr = CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); CeedChk(ierr); 655 } else { 656 x_array = w_array; 657 } 658 if (y != w && y != x) { 659 ierr = CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array); CeedChk(ierr); 660 } else if (y != x) { 661 y_array = w_array; 662 } else { 663 y_array = x_array; 664 } 665 666 for (CeedInt i=0; i<n_w; i++) 667 w_array[i] = x_array[i] * y_array[i]; 668 669 if (y != w && y != x) { 670 ierr = CeedVectorRestoreArrayRead(y, &y_array); CeedChk(ierr); 671 } 672 if (x != w) { 673 ierr = CeedVectorRestoreArrayRead(x, &x_array); CeedChk(ierr); 674 } 675 ierr = CeedVectorRestoreArray(w, &w_array); CeedChk(ierr); 676 return CEED_ERROR_SUCCESS; 677 } 678 679 /** 680 @brief Take the reciprocal of a CeedVector. 681 682 @param vec CeedVector to take reciprocal 683 684 @return An error code: 0 - success, otherwise - failure 685 686 @ref User 687 **/ 688 int CeedVectorReciprocal(CeedVector vec) { 689 int ierr; 690 691 // Check if vector data set 692 if (!vec->state) 693 // LCOV_EXCL_START 694 return CeedError(vec->ceed, CEED_ERROR_INCOMPLETE, 695 "CeedVector must have data set to take reciprocal"); 696 // LCOV_EXCL_STOP 697 698 // Backend impl for GPU, if added 699 if (vec->Reciprocal) { 700 ierr = vec->Reciprocal(vec); CeedChk(ierr); 701 return CEED_ERROR_SUCCESS; 702 } 703 704 CeedInt len; 705 ierr = CeedVectorGetLength(vec, &len); CeedChk(ierr); 706 CeedScalar *array; 707 ierr = CeedVectorGetArray(vec, CEED_MEM_HOST, &array); CeedChk(ierr); 708 for (CeedInt i=0; i<len; i++) 709 if (fabs(array[i]) > CEED_EPSILON) 710 array[i] = 1./array[i]; 711 712 ierr = CeedVectorRestoreArray(vec, &array); CeedChk(ierr); 713 return CEED_ERROR_SUCCESS; 714 } 715 716 /** 717 @brief View a CeedVector 718 719 @param[in] vec CeedVector to view 720 @param[in] fp_fmt Printing format 721 @param[in] stream Filestream to write to 722 723 @return An error code: 0 - success, otherwise - failure 724 725 @ref User 726 **/ 727 int CeedVectorView(CeedVector vec, const char *fp_fmt, FILE *stream) { 728 const CeedScalar *x; 729 730 int ierr = CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &x); CeedChk(ierr); 731 732 char fmt[1024]; 733 fprintf(stream, "CeedVector length %ld\n", (long)vec->length); 734 snprintf(fmt, sizeof fmt, " %s\n", fp_fmt ? fp_fmt : "%g"); 735 for (CeedInt i=0; i<vec->length; i++) 736 fprintf(stream, fmt, x[i]); 737 738 ierr = CeedVectorRestoreArrayRead(vec, &x); CeedChk(ierr); 739 return CEED_ERROR_SUCCESS; 740 } 741 742 /** 743 @brief Get the Ceed associated with a CeedVector 744 745 @param vec CeedVector to retrieve state 746 @param[out] ceed Variable to store ceed 747 748 @return An error code: 0 - success, otherwise - failure 749 750 @ref Advanced 751 **/ 752 int CeedVectorGetCeed(CeedVector vec, Ceed *ceed) { 753 *ceed = vec->ceed; 754 return CEED_ERROR_SUCCESS; 755 } 756 757 /** 758 @brief Get the length of a CeedVector 759 760 @param vec CeedVector to retrieve length 761 @param[out] length Variable to store length 762 763 @return An error code: 0 - success, otherwise - failure 764 765 @ref User 766 **/ 767 int CeedVectorGetLength(CeedVector vec, CeedInt *length) { 768 *length = vec->length; 769 return CEED_ERROR_SUCCESS; 770 } 771 772 /** 773 @brief Destroy a CeedVector 774 775 @param vec CeedVector to destroy 776 777 @return An error code: 0 - success, otherwise - failure 778 779 @ref User 780 **/ 781 int CeedVectorDestroy(CeedVector *vec) { 782 int ierr; 783 784 if (!*vec || --(*vec)->ref_count > 0) return CEED_ERROR_SUCCESS; 785 786 if (((*vec)->state % 2) == 1) 787 return CeedError((*vec)->ceed, CEED_ERROR_ACCESS, 788 "Cannot destroy CeedVector, the writable access " 789 "lock is in use"); 790 791 if ((*vec)->num_readers > 0) 792 // LCOV_EXCL_START 793 return CeedError((*vec)->ceed, CEED_ERROR_ACCESS, 794 "Cannot destroy CeedVector, a process has " 795 "read access"); 796 // LCOV_EXCL_STOP 797 798 if ((*vec)->Destroy) { 799 ierr = (*vec)->Destroy(*vec); CeedChk(ierr); 800 } 801 802 ierr = CeedDestroy(&(*vec)->ceed); CeedChk(ierr); 803 ierr = CeedFree(vec); CeedChk(ierr); 804 return CEED_ERROR_SUCCESS; 805 } 806 807 /// @} 808