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-impl.h>
9 #include <ceed.h>
10 #include <ceed/backend.h>
11 #include <assert.h>
12 #include <math.h>
13 #include <stdbool.h>
14 #include <stdint.h>
15 #include <stdio.h>
16
17 /// @file
18 /// Implementation of public CeedVector interfaces
19
20 /// @cond DOXYGEN_SKIP
21 static struct CeedVector_private ceed_vector_active;
22 static struct CeedVector_private ceed_vector_none;
23 /// @endcond
24
25 /// @addtogroup CeedVectorUser
26 /// @{
27
28 /// Indicate that vector will be provided as an explicit argument to @ref CeedOperatorApply().
29 const CeedVector CEED_VECTOR_ACTIVE = &ceed_vector_active;
30
31 /// Indicate that no vector is applicable (i.e., for @ref CEED_EVAL_WEIGHT).
32 const CeedVector CEED_VECTOR_NONE = &ceed_vector_none;
33
34 /// @}
35
36 /// ----------------------------------------------------------------------------
37 /// CeedVector Internal Functions
38 /// ----------------------------------------------------------------------------
39 /// @addtogroup CeedVectorDeveloper
40 /// @{
41
42 /**
43 @brief View a `CeedVector` passed as a `CeedObject`
44
45 @param[in] vec `CeedVector` to view
46 @param[in] stream Filestream to write to
47
48 @return An error code: 0 - success, otherwise - failure
49
50 @ref Developer
51 **/
CeedVectorView_Object(CeedObject vec,FILE * stream)52 static int CeedVectorView_Object(CeedObject vec, FILE *stream) {
53 CeedCall(CeedVectorView((CeedVector)vec, "%12.8f", stream));
54 return CEED_ERROR_SUCCESS;
55 }
56
57 /**
58 @brief Destroy a `CeedVector` passed as a `CeedObject`
59
60 @param[in,out] vec Address of `CeedVector` to destroy
61
62 @return An error code: 0 - success, otherwise - failure
63
64 @ref Developer
65 **/
CeedVectorDestroy_Object(CeedObject * vec)66 static int CeedVectorDestroy_Object(CeedObject *vec) {
67 CeedCall(CeedVectorDestroy((CeedVector *)vec));
68 return CEED_ERROR_SUCCESS;
69 }
70
71 /// @}
72
73 /// ----------------------------------------------------------------------------
74 /// CeedVector Backend API
75 /// ----------------------------------------------------------------------------
76 /// @addtogroup CeedVectorBackend
77 /// @{
78
79 /**
80 @brief Check for valid data in a `CeedVector`
81
82 @param[in] vec `CeedVector` to check validity
83 @param[out] has_valid_array Variable to store validity
84
85 @return An error code: 0 - success, otherwise - failure
86
87 @ref Backend
88 **/
CeedVectorHasValidArray(CeedVector vec,bool * has_valid_array)89 int CeedVectorHasValidArray(CeedVector vec, bool *has_valid_array) {
90 CeedSize length;
91
92 CeedCheck(vec->HasValidArray, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedVectorHasValidArray");
93 CeedCall(CeedVectorGetLength(vec, &length));
94 if (length == 0) {
95 *has_valid_array = true;
96 return CEED_ERROR_SUCCESS;
97 }
98 CeedCall(vec->HasValidArray(vec, has_valid_array));
99 return CEED_ERROR_SUCCESS;
100 }
101
102 /**
103 @brief Check for borrowed array of a specific @ref CeedMemType in a `CeedVector`
104
105 @param[in] vec `CeedVector` to check
106 @param[in] mem_type Memory type to check
107 @param[out] has_borrowed_array_of_type Variable to store result
108
109 @return An error code: 0 - success, otherwise - failure
110
111 @ref Backend
112 **/
CeedVectorHasBorrowedArrayOfType(CeedVector vec,CeedMemType mem_type,bool * has_borrowed_array_of_type)113 int CeedVectorHasBorrowedArrayOfType(CeedVector vec, CeedMemType mem_type, bool *has_borrowed_array_of_type) {
114 CeedCheck(vec->HasBorrowedArrayOfType, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED,
115 "Backend does not support CeedVectorHasBorrowedArrayOfType");
116 CeedCall(vec->HasBorrowedArrayOfType(vec, mem_type, has_borrowed_array_of_type));
117 return CEED_ERROR_SUCCESS;
118 }
119
120 /**
121 @brief Get the state of a `CeedVector`
122
123 @param[in] vec `CeedVector` to retrieve state
124 @param[out] state Variable to store state
125
126 @return An error code: 0 - success, otherwise - failure
127
128 @ref Backend
129 **/
CeedVectorGetState(CeedVector vec,uint64_t * state)130 int CeedVectorGetState(CeedVector vec, uint64_t *state) {
131 *state = vec->state;
132 return CEED_ERROR_SUCCESS;
133 }
134
135 /**
136 @brief Get the backend data of a `CeedVector`
137
138 @param[in] vec `CeedVector` to retrieve state
139 @param[out] data Variable to store data
140
141 @return An error code: 0 - success, otherwise - failure
142
143 @ref Backend
144 **/
CeedVectorGetData(CeedVector vec,void * data)145 int CeedVectorGetData(CeedVector vec, void *data) {
146 *(void **)data = vec->data;
147 return CEED_ERROR_SUCCESS;
148 }
149
150 /**
151 @brief Set the backend data of a `CeedVector`
152
153 @param[in,out] vec `CeedVector` to retrieve state
154 @param[in] data Data to set
155
156 @return An error code: 0 - success, otherwise - failure
157
158 @ref Backend
159 **/
CeedVectorSetData(CeedVector vec,void * data)160 int CeedVectorSetData(CeedVector vec, void *data) {
161 vec->data = data;
162 return CEED_ERROR_SUCCESS;
163 }
164
165 /**
166 @brief Increment the reference counter for a `CeedVector`
167
168 @param[in,out] vec `CeedVector` to increment the reference counter
169
170 @return An error code: 0 - success, otherwise - failure
171
172 @ref Backend
173 **/
CeedVectorReference(CeedVector vec)174 int CeedVectorReference(CeedVector vec) {
175 CeedCall(CeedObjectReference((CeedObject)vec));
176 return CEED_ERROR_SUCCESS;
177 }
178
179 /// @}
180
181 /// ----------------------------------------------------------------------------
182 /// CeedVector Public API
183 /// ----------------------------------------------------------------------------
184 /// @addtogroup CeedVectorUser
185 /// @{
186
187 /**
188 @brief Create a `CeedVector` of the specified length (does not allocate memory)
189
190 @param[in] ceed `Ceed` object used to create the `CeedVector`
191 @param[in] length Length of vector
192 @param[out] vec Address of the variable where the newly created `CeedVector` will be stored
193
194 @return An error code: 0 - success, otherwise - failure
195
196 @ref User
197 **/
CeedVectorCreate(Ceed ceed,CeedSize length,CeedVector * vec)198 int CeedVectorCreate(Ceed ceed, CeedSize length, CeedVector *vec) {
199 CeedCheck(length >= 0, ceed, CEED_ERROR_UNSUPPORTED, "CeedVector must have length >= 0, received %" CeedSize_FMT, length);
200 if (!ceed->VectorCreate) {
201 Ceed delegate;
202
203 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Vector"));
204 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement VectorCreate");
205 CeedCall(CeedVectorCreate(delegate, length, vec));
206 CeedCall(CeedDestroy(&delegate));
207 return CEED_ERROR_SUCCESS;
208 }
209
210 CeedCall(CeedCalloc(1, vec));
211 CeedCall(CeedObjectCreate(ceed, CeedVectorView_Object, CeedVectorDestroy_Object, &(*vec)->obj));
212 (*vec)->length = length;
213 (*vec)->state = 0;
214 CeedCall(ceed->VectorCreate(length, *vec));
215 return CEED_ERROR_SUCCESS;
216 }
217
218 /**
219 @brief Copy the pointer to a `CeedVector`.
220
221 Both pointers should be destroyed with @ref CeedVectorDestroy().
222
223 Note: If the value of `*vec_copy` passed to this function is non-`NULL`, then it is assumed that `*vec_copy` is a pointer to a `CeedVector`.
224 This `CeedVector` will be destroyed if `*vec_copy` is the only reference to this `CeedVector`.
225
226 @param[in] vec `CeedVector` to copy reference to
227 @param[in,out] vec_copy Variable to store copied reference
228
229 @return An error code: 0 - success, otherwise - failure
230
231 @ref User
232 **/
CeedVectorReferenceCopy(CeedVector vec,CeedVector * vec_copy)233 int CeedVectorReferenceCopy(CeedVector vec, CeedVector *vec_copy) {
234 if (vec != CEED_VECTOR_ACTIVE && vec != CEED_VECTOR_NONE) CeedCall(CeedVectorReference(vec));
235 CeedCall(CeedVectorDestroy(vec_copy));
236 *vec_copy = vec;
237 return CEED_ERROR_SUCCESS;
238 }
239
240 /**
241 @brief Copy a `CeedVector` into a different `CeedVector`.
242
243 @param[in] vec `CeedVector` to copy
244 @param[in,out] vec_copy `CeedVector` to copy array into
245
246 @return An error code: 0 - success, otherwise - failure
247
248 @ref User
249 **/
CeedVectorCopy(CeedVector vec,CeedVector vec_copy)250 int CeedVectorCopy(CeedVector vec, CeedVector vec_copy) {
251 CeedMemType mem_type, mem_type_copy;
252 CeedScalar *array;
253
254 // Get the preferred memory types
255 {
256 Ceed ceed;
257
258 CeedCall(CeedVectorGetCeed(vec, &ceed));
259 CeedCall(CeedGetPreferredMemType(ceed, &mem_type));
260 CeedCall(CeedDestroy(&ceed));
261
262 CeedCall(CeedVectorGetCeed(vec_copy, &ceed));
263 CeedCall(CeedGetPreferredMemType(ceed, &mem_type_copy));
264 CeedCall(CeedDestroy(&ceed));
265 }
266
267 // Check that both have same memory type
268 if (mem_type != mem_type_copy) mem_type = CEED_MEM_HOST;
269
270 // Check compatible lengths
271 {
272 CeedSize length_vec, length_copy;
273
274 CeedCall(CeedVectorGetLength(vec, &length_vec));
275 CeedCall(CeedVectorGetLength(vec_copy, &length_copy));
276 CeedCheck(length_vec == length_copy, CeedVectorReturnCeed(vec), CEED_ERROR_INCOMPATIBLE, "CeedVectors must have the same length to copy");
277 }
278
279 // Copy the values from vec to vec_copy
280 CeedCall(CeedVectorGetArray(vec, mem_type, &array));
281 CeedCall(CeedVectorSetArray(vec_copy, mem_type, CEED_COPY_VALUES, array));
282
283 CeedCall(CeedVectorRestoreArray(vec, &array));
284 return CEED_ERROR_SUCCESS;
285 }
286
287 /**
288 @brief Copy a strided portion of `CeedVector` contents into a different `CeedVector`
289
290 @param[in] vec `CeedVector` to copy
291 @param[in] start First index to copy in the range `[start, stop)`
292 @param[in] stop One past the last element to copy in the range, or `-1` for `length`
293 @param[in] step Stride between indices to copy
294 @param[in,out] vec_copy `CeedVector` to copy values to
295
296 @return An error code: 0 - success, otherwise - failure
297
298 @ref User
299 **/
CeedVectorCopyStrided(CeedVector vec,CeedSize start,CeedSize stop,CeedSize step,CeedVector vec_copy)300 int CeedVectorCopyStrided(CeedVector vec, CeedSize start, CeedSize stop, CeedSize step, CeedVector vec_copy) {
301 CeedSize length;
302 const CeedScalar *array = NULL;
303 CeedScalar *array_copy = NULL;
304
305 // Check length
306 {
307 CeedSize length_vec, length_copy;
308
309 CeedCall(CeedVectorGetLength(vec, &length_vec));
310 CeedCall(CeedVectorGetLength(vec_copy, &length_copy));
311 if (length_vec <= 0 || length_copy <= 0) return CEED_ERROR_SUCCESS;
312 length = length_vec < length_copy ? length_vec : length_copy;
313 }
314 CeedCheck(stop >= -1 && stop <= length, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
315 "Invalid value for stop %" CeedSize_FMT ", must be in the range [-1, length]", stop);
316 CeedCheck(start >= 0 && start <= length && (start <= stop || stop == -1), CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
317 "Invalid value for start %" CeedSize_FMT ", must be in the range [0, stop]", start);
318
319 // Backend version
320 if (vec->CopyStrided && vec_copy->CopyStrided) {
321 CeedCall(vec->CopyStrided(vec, start, stop, step, vec_copy));
322 vec_copy->state += 2;
323 return CEED_ERROR_SUCCESS;
324 }
325
326 // Copy
327 CeedCall(CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &array));
328 CeedCall(CeedVectorGetArray(vec_copy, CEED_MEM_HOST, &array_copy));
329 if (stop == -1) stop = length;
330 for (CeedSize i = start; i < stop; i += step) array_copy[i] = array[i];
331
332 // Cleanup
333 CeedCall(CeedVectorRestoreArrayRead(vec, &array));
334 CeedCall(CeedVectorRestoreArray(vec_copy, &array_copy));
335 return CEED_ERROR_SUCCESS;
336 }
337
338 /**
339 @brief Set the array used by a `CeedVector`, freeing any previously allocated array if applicable.
340
341 The backend may copy values to a different @ref CeedMemType, such as during @ref CeedOperatorApply().
342 See also @ref CeedVectorSyncArray() and @ref CeedVectorTakeArray().
343
344 @param[in,out] vec `CeedVector`
345 @param[in] mem_type Memory type of the array being passed
346 @param[in] copy_mode Copy mode for the array
347 @param[in] array Array to be used, or `NULL` with @ref CEED_COPY_VALUES to have the library allocate
348
349 @return An error code: 0 - success, otherwise - failure
350
351 @ref User
352 **/
CeedVectorSetArray(CeedVector vec,CeedMemType mem_type,CeedCopyMode copy_mode,CeedScalar * array)353 int CeedVectorSetArray(CeedVector vec, CeedMemType mem_type, CeedCopyMode copy_mode, CeedScalar *array) {
354 CeedSize length;
355
356 CeedCheck(vec->SetArray, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED, "Backend does not support VectorSetArray");
357 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
358 "Cannot grant CeedVector array access, the access lock is already in use");
359 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot grant CeedVector array access, a process has read access");
360
361 CeedCall(CeedVectorGetLength(vec, &length));
362 if (length > 0) CeedCall(vec->SetArray(vec, mem_type, copy_mode, array));
363 vec->state += 2;
364 return CEED_ERROR_SUCCESS;
365 }
366
367 /**
368 @brief Set the `CeedVector` to a constant value
369
370 @param[in,out] vec `CeedVector`
371 @param[in] value Value to be used
372
373 @return An error code: 0 - success, otherwise - failure
374
375 @ref User
376 **/
CeedVectorSetValue(CeedVector vec,CeedScalar value)377 int CeedVectorSetValue(CeedVector vec, CeedScalar value) {
378 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
379 "Cannot grant CeedVector array access, the access lock is already in use");
380 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot grant CeedVector array access, a process has read access");
381
382 if (vec->SetValue) {
383 CeedCall(vec->SetValue(vec, value));
384 vec->state += 2;
385 } else {
386 CeedSize length;
387 CeedScalar *array;
388
389 CeedCall(CeedVectorGetArrayWrite(vec, CEED_MEM_HOST, &array));
390 CeedCall(CeedVectorGetLength(vec, &length));
391 for (CeedSize i = 0; i < length; i++) array[i] = value;
392 CeedCall(CeedVectorRestoreArray(vec, &array));
393 }
394 return CEED_ERROR_SUCCESS;
395 }
396
397 /**
398 @brief Set a portion of a `CeedVector` to a constant value.
399
400 Note: The `CeedVector` must already have valid data set via @ref CeedVectorSetArray() or similar.
401
402 @param[in,out] vec `CeedVector`
403 @param[in] start First index to set in range `[start, stop)`
404 @param[in] stop One past the last element to set in the range, or `-1` for `length`
405 @param[in] step Stride between indices to set
406 @param[in] value Value to be used
407
408 @return An error code: 0 - success, otherwise - failure
409
410 @ref User
411 **/
CeedVectorSetValueStrided(CeedVector vec,CeedSize start,CeedSize stop,CeedSize step,CeedScalar value)412 int CeedVectorSetValueStrided(CeedVector vec, CeedSize start, CeedSize stop, CeedSize step, CeedScalar value) {
413 CeedSize length;
414
415 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
416 "Cannot grant CeedVector array access, the access lock is already in use");
417 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot grant CeedVector array access, a process has read access");
418 CeedCall(CeedVectorGetLength(vec, &length));
419 CeedCheck(stop >= -1 && stop <= length, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
420 "Invalid value for stop %" CeedSize_FMT ", must be in the range [-1, length]", stop);
421
422 if (vec->SetValueStrided) {
423 CeedCall(vec->SetValueStrided(vec, start, stop, step, value));
424 vec->state += 2;
425 } else {
426 CeedScalar *array;
427
428 if (length <= 0) return CEED_ERROR_SUCCESS;
429 if (stop == -1) stop = length;
430 CeedCall(CeedVectorGetArray(vec, CEED_MEM_HOST, &array));
431 for (CeedSize i = start; i < stop; i += step) array[i] = value;
432 CeedCall(CeedVectorRestoreArray(vec, &array));
433 }
434 return CEED_ERROR_SUCCESS;
435 }
436
437 /**
438 @brief Sync the `CeedVector` to a specified `mem_type`.
439
440 This function is used to force synchronization of arrays set with @ref CeedVectorSetArray().
441 If the requested `mem_type` is already synchronized, this function results in a no-op.
442
443 @param[in,out] vec `CeedVector`
444 @param[in] mem_type @ref CeedMemType to be synced
445
446 @return An error code: 0 - success, otherwise - failure
447
448 @ref User
449 **/
CeedVectorSyncArray(CeedVector vec,CeedMemType mem_type)450 int CeedVectorSyncArray(CeedVector vec, CeedMemType mem_type) {
451 CeedSize length;
452
453 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot sync CeedVector, the access lock is already in use");
454
455 // Don't sync empty array
456 CeedCall(CeedVectorGetLength(vec, &length));
457 if (length == 0) return CEED_ERROR_SUCCESS;
458
459 if (vec->SyncArray) {
460 CeedCall(vec->SyncArray(vec, mem_type));
461 } else {
462 const CeedScalar *array;
463
464 CeedCall(CeedVectorGetArrayRead(vec, mem_type, &array));
465 CeedCall(CeedVectorRestoreArrayRead(vec, &array));
466 }
467 return CEED_ERROR_SUCCESS;
468 }
469
470 /**
471 @brief Take ownership of the `CeedVector` array set by @ref CeedVectorSetArray() with @ref CEED_USE_POINTER and remove the array from the `CeedVector`.
472
473 The caller is responsible for managing and freeing the array.
474 This function will error if @ref CeedVectorSetArray() was not previously called with @ref CEED_USE_POINTER for the corresponding mem_type.
475
476 @param[in,out] vec `CeedVector`
477 @param[in] mem_type Memory type on which to take the array.
478 If the backend uses a different memory type, this will perform a copy.
479 @param[out] array Array on memory type `mem_type`, or `NULL` if array pointer is not required
480
481 @return An error code: 0 - success, otherwise - failure
482
483 @ref User
484 **/
CeedVectorTakeArray(CeedVector vec,CeedMemType mem_type,CeedScalar ** array)485 int CeedVectorTakeArray(CeedVector vec, CeedMemType mem_type, CeedScalar **array) {
486 CeedSize length;
487 CeedScalar *temp_array = NULL;
488
489 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot take CeedVector array, the access lock is already in use");
490 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot take CeedVector array, a process has read access");
491
492 CeedCall(CeedVectorGetLength(vec, &length));
493 if (length > 0) {
494 bool has_borrowed_array_of_type = true, has_valid_array = true;
495
496 CeedCall(CeedVectorHasBorrowedArrayOfType(vec, mem_type, &has_borrowed_array_of_type));
497 CeedCheck(has_borrowed_array_of_type, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
498 "CeedVector has no borrowed %s array, must set array with CeedVectorSetArray", CeedMemTypes[mem_type]);
499
500 CeedCall(CeedVectorHasValidArray(vec, &has_valid_array));
501 CeedCheck(has_valid_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
502 "CeedVector has no valid data to take, must set data with CeedVectorSetValue or CeedVectorSetArray");
503
504 CeedCall(vec->TakeArray(vec, mem_type, &temp_array));
505 }
506 if (array) (*array) = temp_array;
507 return CEED_ERROR_SUCCESS;
508 }
509
510 /**
511 @brief Get read/write access to a `CeedVector` via the specified memory type.
512
513 Restore access with @ref CeedVectorRestoreArray().
514
515 @param[in,out] vec `CeedVector` to access
516 @param[in] mem_type Memory type on which to access the array.
517 If the backend uses a different memory type, this will perform a copy.
518 @param[out] array Array on memory type `mem_type`
519
520 @note The @ref CeedVectorGetArray() and @ref CeedVectorRestoreArray() functions provide access to array pointers in the desired memory space.
521 Pairing get/restore allows the `CeedVector` to track access, thus knowing if norms or other operations may need to be recomputed.
522
523 @return An error code: 0 - success, otherwise - failure
524
525 @ref User
526 **/
CeedVectorGetArray(CeedVector vec,CeedMemType mem_type,CeedScalar ** array)527 int CeedVectorGetArray(CeedVector vec, CeedMemType mem_type, CeedScalar **array) {
528 CeedSize length;
529
530 CeedCheck(vec->GetArray, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED, "Backend does not support GetArray");
531 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
532 "Cannot grant CeedVector array access, the access lock is already in use");
533 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot grant CeedVector array access, a process has read access");
534
535 CeedCall(CeedVectorGetLength(vec, &length));
536 if (length > 0) {
537 bool has_valid_array = true;
538
539 CeedCall(CeedVectorHasValidArray(vec, &has_valid_array));
540 CeedCheck(has_valid_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
541 "CeedVector has no valid data to read, must set data with CeedVectorSetValue or CeedVectorSetArray");
542
543 CeedCall(vec->GetArray(vec, mem_type, array));
544 } else {
545 *array = NULL;
546 }
547 vec->state++;
548 return CEED_ERROR_SUCCESS;
549 }
550
551 /**
552 @brief Get read-only access to a `CeedVector` via the specified memory type.
553
554 Restore access with @ref CeedVectorRestoreArrayRead().
555
556 @param[in] vec `CeedVector` to access
557 @param[in] mem_type Memory type on which to access the array.
558 If the backend uses a different memory type, this will perform a copy (possibly cached).
559 @param[out] array Array on memory type `mem_type`
560
561 @return An error code: 0 - success, otherwise - failure
562
563 @ref User
564 **/
CeedVectorGetArrayRead(CeedVector vec,CeedMemType mem_type,const CeedScalar ** array)565 int CeedVectorGetArrayRead(CeedVector vec, CeedMemType mem_type, const CeedScalar **array) {
566 CeedSize length;
567
568 CeedCheck(vec->GetArrayRead, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED, "Backend does not support GetArrayRead");
569 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
570 "Cannot grant CeedVector read-only array access, the access lock is already in use");
571
572 CeedCall(CeedVectorGetLength(vec, &length));
573 if (length > 0) {
574 bool has_valid_array = true;
575
576 CeedCall(CeedVectorHasValidArray(vec, &has_valid_array));
577 CeedCheck(has_valid_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
578 "CeedVector has no valid data to read, must set data with CeedVectorSetValue or CeedVectorSetArray");
579
580 CeedCall(vec->GetArrayRead(vec, mem_type, array));
581 } else {
582 *array = NULL;
583 }
584 vec->num_readers++;
585 return CEED_ERROR_SUCCESS;
586 }
587
588 /**
589 @brief Get write access to a `CeedVector` via the specified memory type.
590
591 Restore access with @ref CeedVectorRestoreArray().
592 All old values should be assumed to be invalid.
593
594 @param[in,out] vec `CeedVector` to access
595 @param[in] mem_type Memory type on which to access the array.
596 @param[out] array Array on memory type `mem_type`
597
598 @return An error code: 0 - success, otherwise - failure
599
600 @ref User
601 **/
CeedVectorGetArrayWrite(CeedVector vec,CeedMemType mem_type,CeedScalar ** array)602 int CeedVectorGetArrayWrite(CeedVector vec, CeedMemType mem_type, CeedScalar **array) {
603 CeedSize length;
604
605 CeedCheck(vec->GetArrayWrite, CeedVectorReturnCeed(vec), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedVectorGetArrayWrite");
606 CeedCheck(vec->state % 2 == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
607 "Cannot grant CeedVector array access, the access lock is already in use");
608 CeedCheck(vec->num_readers == 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot grant CeedVector array access, a process has read access");
609
610 CeedCall(CeedVectorGetLength(vec, &length));
611 if (length > 0) {
612 CeedCall(vec->GetArrayWrite(vec, mem_type, array));
613 } else {
614 *array = NULL;
615 }
616 vec->state++;
617 return CEED_ERROR_SUCCESS;
618 }
619
620 /**
621 @brief Restore an array obtained using @ref CeedVectorGetArray() or @ref CeedVectorGetArrayWrite()
622
623 @param[in,out] vec `CeedVector` to restore
624 @param[in,out] array Array of vector data
625
626 @return An error code: 0 - success, otherwise - failure
627
628 @ref User
629 **/
CeedVectorRestoreArray(CeedVector vec,CeedScalar ** array)630 int CeedVectorRestoreArray(CeedVector vec, CeedScalar **array) {
631 CeedSize length;
632
633 CeedCheck(vec->state % 2 == 1, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS, "Cannot restore CeedVector array access, access was not granted");
634 CeedCall(CeedVectorGetLength(vec, &length));
635 if (length > 0 && vec->RestoreArray) CeedCall(vec->RestoreArray(vec));
636 *array = NULL;
637 vec->state++;
638 return CEED_ERROR_SUCCESS;
639 }
640
641 /**
642 @brief Restore an array obtained using @ref CeedVectorGetArrayRead()
643
644 @param[in] vec `CeedVector` to restore
645 @param[in,out] array Array of vector data
646
647 @return An error code: 0 - success, otherwise - failure
648
649 @ref User
650 **/
CeedVectorRestoreArrayRead(CeedVector vec,const CeedScalar ** array)651 int CeedVectorRestoreArrayRead(CeedVector vec, const CeedScalar **array) {
652 CeedSize length;
653
654 CeedCheck(vec->num_readers > 0, CeedVectorReturnCeed(vec), CEED_ERROR_ACCESS,
655 "Cannot restore CeedVector array read access, access was not granted");
656 vec->num_readers--;
657 CeedCall(CeedVectorGetLength(vec, &length));
658 if (length > 0 && vec->num_readers == 0 && vec->RestoreArrayRead) CeedCall(vec->RestoreArrayRead(vec));
659 *array = NULL;
660 return CEED_ERROR_SUCCESS;
661 }
662
663 /**
664 @brief Get the norm of a `CeedVector`.
665
666 Note: This operation is local to the `CeedVector`.
667 This function will likely not provide the desired results for the norm of the libCEED portion of a parallel vector or a `CeedVector` with duplicated or hanging nodes.
668
669 @param[in] vec `CeedVector` to retrieve maximum value
670 @param[in] norm_type Norm type @ref CEED_NORM_1, @ref CEED_NORM_2, or @ref CEED_NORM_MAX
671 @param[out] norm Variable to store norm value
672
673 @return An error code: 0 - success, otherwise - failure
674
675 @ref User
676 **/
CeedVectorNorm(CeedVector vec,CeedNormType norm_type,CeedScalar * norm)677 int CeedVectorNorm(CeedVector vec, CeedNormType norm_type, CeedScalar *norm) {
678 bool has_valid_array = true;
679 CeedSize length;
680
681 CeedCall(CeedVectorHasValidArray(vec, &has_valid_array));
682 CeedCheck(has_valid_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
683 "CeedVector has no valid data to compute norm, must set data with CeedVectorSetValue or CeedVectorSetArray");
684
685 CeedCall(CeedVectorGetLength(vec, &length));
686 if (length == 0) {
687 *norm = 0;
688 return CEED_ERROR_SUCCESS;
689 }
690
691 // Backend impl for GPU, if added
692 if (vec->Norm) {
693 CeedCall(vec->Norm(vec, norm_type, norm));
694 return CEED_ERROR_SUCCESS;
695 }
696
697 const CeedScalar *array;
698 CeedCall(CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &array));
699 assert(array);
700
701 *norm = 0.;
702 switch (norm_type) {
703 case CEED_NORM_1:
704 for (CeedSize i = 0; i < length; i++) {
705 *norm += fabs(array[i]);
706 }
707 break;
708 case CEED_NORM_2:
709 for (CeedSize i = 0; i < length; i++) {
710 *norm += fabs(array[i]) * fabs(array[i]);
711 }
712 break;
713 case CEED_NORM_MAX:
714 for (CeedSize i = 0; i < length; i++) {
715 const CeedScalar abs_v_i = fabs(array[i]);
716 *norm = *norm > abs_v_i ? *norm : abs_v_i;
717 }
718 }
719 if (norm_type == CEED_NORM_2) *norm = sqrt(*norm);
720
721 CeedCall(CeedVectorRestoreArrayRead(vec, &array));
722 return CEED_ERROR_SUCCESS;
723 }
724
725 /**
726 @brief Compute `x = alpha x`
727
728 @param[in,out] x `CeedVector` for scaling
729 @param[in] alpha scaling factor
730
731 @return An error code: 0 - success, otherwise - failure
732
733 @ref User
734 **/
CeedVectorScale(CeedVector x,CeedScalar alpha)735 int CeedVectorScale(CeedVector x, CeedScalar alpha) {
736 bool has_valid_array = true;
737 CeedSize length;
738 CeedScalar *x_array = NULL;
739
740 CeedCall(CeedVectorHasValidArray(x, &has_valid_array));
741 CeedCheck(has_valid_array, CeedVectorReturnCeed(x), CEED_ERROR_BACKEND,
742 "CeedVector has no valid data to scale, must set data with CeedVectorSetValue or CeedVectorSetArray");
743
744 // Return early for empty vector
745 CeedCall(CeedVectorGetLength(x, &length));
746 if (length == 0) return CEED_ERROR_SUCCESS;
747
748 // Backend implementation
749 if (x->Scale) return x->Scale(x, alpha);
750
751 // Default implementation
752 CeedCall(CeedVectorGetArray(x, CEED_MEM_HOST, &x_array));
753 assert(x_array);
754 for (CeedSize i = 0; i < length; i++) x_array[i] *= alpha;
755 CeedCall(CeedVectorRestoreArray(x, &x_array));
756 return CEED_ERROR_SUCCESS;
757 }
758
759 /**
760 @brief Compute `y = alpha x + y`
761
762 @param[in,out] y target `CeedVector` for sum
763 @param[in] alpha scaling factor
764 @param[in] x second `CeedVector`, must be different than ``y`
765
766 @return An error code: 0 - success, otherwise - failure
767
768 @ref User
769 **/
CeedVectorAXPY(CeedVector y,CeedScalar alpha,CeedVector x)770 int CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x) {
771 bool has_valid_array_x = true, has_valid_array_y = true;
772 CeedSize length_x, length_y;
773 CeedScalar *y_array = NULL;
774 CeedScalar const *x_array = NULL;
775
776 CeedCall(CeedVectorGetLength(y, &length_y));
777 CeedCall(CeedVectorGetLength(x, &length_x));
778 CeedCheck(length_x == length_y, CeedVectorReturnCeed(y), CEED_ERROR_UNSUPPORTED,
779 "Cannot add vector of different lengths."
780 " x length: %" CeedSize_FMT " y length: %" CeedSize_FMT,
781 length_x, length_y);
782 CeedCheck(x != y, CeedVectorReturnCeed(y), CEED_ERROR_UNSUPPORTED, "Cannot use same vector for x and y in CeedVectorAXPY");
783
784 CeedCall(CeedVectorHasValidArray(x, &has_valid_array_x));
785 CeedCheck(has_valid_array_x, CeedVectorReturnCeed(y), CEED_ERROR_BACKEND,
786 "CeedVector x has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
787 CeedCall(CeedVectorHasValidArray(y, &has_valid_array_y));
788 CeedCheck(has_valid_array_y, CeedVectorReturnCeed(y), CEED_ERROR_BACKEND,
789 "CeedVector y has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
790
791 {
792 Ceed ceed_x, ceed_y, ceed_parent_x, ceed_parent_y;
793
794 CeedCall(CeedVectorGetCeed(y, &ceed_y));
795 CeedCall(CeedVectorGetCeed(x, &ceed_x));
796 CeedCall(CeedGetParent(ceed_x, &ceed_parent_x));
797 CeedCall(CeedGetParent(ceed_y, &ceed_parent_y));
798 CeedCall(CeedDestroy(&ceed_x));
799 CeedCall(CeedDestroy(&ceed_y));
800 CeedCheck(ceed_parent_x == ceed_parent_y, CeedVectorReturnCeed(y), CEED_ERROR_INCOMPATIBLE,
801 "Vectors x and y must be created by the same Ceed context");
802 CeedCall(CeedDestroy(&ceed_parent_x));
803 CeedCall(CeedDestroy(&ceed_parent_y));
804 }
805
806 // Return early for empty vectors
807 if (length_y == 0) return CEED_ERROR_SUCCESS;
808
809 // Backend implementation
810 if (y->AXPY) {
811 CeedCall(y->AXPY(y, alpha, x));
812 return CEED_ERROR_SUCCESS;
813 }
814
815 // Default implementation
816 CeedCall(CeedVectorGetArray(y, CEED_MEM_HOST, &y_array));
817 CeedCall(CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array));
818
819 assert(x_array);
820 assert(y_array);
821
822 for (CeedSize i = 0; i < length_y; i++) y_array[i] += alpha * x_array[i];
823
824 CeedCall(CeedVectorRestoreArray(y, &y_array));
825 CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
826 return CEED_ERROR_SUCCESS;
827 }
828
829 /**
830 @brief Compute `y = alpha x + beta y`
831
832 @param[in,out] y target `CeedVector` for sum
833 @param[in] alpha first scaling factor
834 @param[in] beta second scaling factor
835 @param[in] x second `CeedVector`, must be different than `y`
836
837 @return An error code: 0 - success, otherwise - failure
838
839 @ref User
840 **/
CeedVectorAXPBY(CeedVector y,CeedScalar alpha,CeedScalar beta,CeedVector x)841 int CeedVectorAXPBY(CeedVector y, CeedScalar alpha, CeedScalar beta, CeedVector x) {
842 bool has_valid_array_x = true, has_valid_array_y = true;
843 CeedSize length_x, length_y;
844 CeedScalar *y_array = NULL;
845 CeedScalar const *x_array = NULL;
846
847 CeedCall(CeedVectorGetLength(y, &length_y));
848 CeedCall(CeedVectorGetLength(x, &length_x));
849 CeedCheck(length_x == length_y, CeedVectorReturnCeed(y), CEED_ERROR_UNSUPPORTED,
850 "Cannot add vector of different lengths."
851 " x length: %" CeedSize_FMT " y length: %" CeedSize_FMT,
852 length_x, length_y);
853 CeedCheck(x != y, CeedVectorReturnCeed(y), CEED_ERROR_UNSUPPORTED, "Cannot use same vector for x and y in CeedVectorAXPBY");
854
855 CeedCall(CeedVectorHasValidArray(x, &has_valid_array_x));
856 CeedCheck(has_valid_array_x, CeedVectorReturnCeed(y), CEED_ERROR_BACKEND,
857 "CeedVector x has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
858 CeedCall(CeedVectorHasValidArray(y, &has_valid_array_y));
859 CeedCheck(has_valid_array_y, CeedVectorReturnCeed(y), CEED_ERROR_BACKEND,
860 "CeedVector y has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
861
862 {
863 Ceed ceed_x, ceed_y, ceed_parent_x, ceed_parent_y;
864
865 CeedCall(CeedVectorGetCeed(y, &ceed_y));
866 CeedCall(CeedVectorGetCeed(x, &ceed_x));
867 CeedCall(CeedGetParent(ceed_x, &ceed_parent_x));
868 CeedCall(CeedGetParent(ceed_y, &ceed_parent_y));
869 CeedCall(CeedDestroy(&ceed_x));
870 CeedCall(CeedDestroy(&ceed_y));
871 CeedCheck(ceed_parent_x == ceed_parent_y, CeedVectorReturnCeed(y), CEED_ERROR_INCOMPATIBLE,
872 "Vectors x and y must be created by the same Ceed context");
873 CeedCall(CeedDestroy(&ceed_parent_x));
874 CeedCall(CeedDestroy(&ceed_parent_y));
875 }
876
877 // Return early for empty vectors
878 if (length_y == 0) return CEED_ERROR_SUCCESS;
879
880 // Backend implementation
881 if (y->AXPBY) {
882 CeedCall(y->AXPBY(y, alpha, beta, x));
883 return CEED_ERROR_SUCCESS;
884 }
885
886 // Default implementation
887 CeedCall(CeedVectorGetArray(y, CEED_MEM_HOST, &y_array));
888 CeedCall(CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array));
889
890 assert(x_array);
891 assert(y_array);
892
893 for (CeedSize i = 0; i < length_y; i++) y_array[i] = alpha * x_array[i] + beta * y_array[i];
894
895 CeedCall(CeedVectorRestoreArray(y, &y_array));
896 CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
897 return CEED_ERROR_SUCCESS;
898 }
899
900 /**
901 @brief Compute the pointwise multiplication \f$w = x .* y\f$.
902
903 Any subset of `x`, `y`, and `w` may be the same `CeedVector`.
904
905 @param[out] w target `CeedVector` for the product
906 @param[in] x first `CeedVector` for product
907 @param[in] y second `CeedVector` for the product
908
909 @return An error code: 0 - success, otherwise - failure
910
911 @ref User
912 **/
CeedVectorPointwiseMult(CeedVector w,CeedVector x,CeedVector y)913 int CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y) {
914 bool has_valid_array_x = true, has_valid_array_y = true;
915 CeedScalar *w_array = NULL;
916 CeedScalar const *x_array = NULL, *y_array = NULL;
917 CeedSize length_w, length_x, length_y;
918
919 CeedCall(CeedVectorGetLength(w, &length_w));
920 CeedCall(CeedVectorGetLength(x, &length_x));
921 CeedCall(CeedVectorGetLength(y, &length_y));
922 CeedCheck(length_x >= length_w && length_y >= length_w, CeedVectorReturnCeed(w), CEED_ERROR_UNSUPPORTED,
923 "Cannot pointwise multiply vectors of incompatible lengths."
924 " w length: %" CeedSize_FMT " x length: %" CeedSize_FMT " y length: %" CeedSize_FMT,
925 length_w, length_x, length_y);
926
927 {
928 Ceed ceed_w, ceed_x, ceed_y, ceed_parent_w, ceed_parent_x, ceed_parent_y;
929
930 CeedCall(CeedVectorGetCeed(w, &ceed_w));
931 CeedCall(CeedVectorGetCeed(x, &ceed_x));
932 CeedCall(CeedVectorGetCeed(y, &ceed_y));
933 CeedCall(CeedGetParent(ceed_w, &ceed_parent_w));
934 CeedCall(CeedGetParent(ceed_x, &ceed_parent_x));
935 CeedCall(CeedGetParent(ceed_y, &ceed_parent_y));
936 CeedCall(CeedDestroy(&ceed_w));
937 CeedCall(CeedDestroy(&ceed_x));
938 CeedCall(CeedDestroy(&ceed_y));
939 CeedCheck(ceed_parent_w == ceed_parent_x && ceed_parent_w == ceed_parent_y, CeedVectorReturnCeed(w), CEED_ERROR_INCOMPATIBLE,
940 "Vectors w, x, and y must be created by the same Ceed context");
941 CeedCall(CeedDestroy(&ceed_parent_w));
942 CeedCall(CeedDestroy(&ceed_parent_x));
943 CeedCall(CeedDestroy(&ceed_parent_y));
944 }
945
946 CeedCall(CeedVectorHasValidArray(x, &has_valid_array_x));
947 CeedCheck(has_valid_array_x, CeedVectorReturnCeed(w), CEED_ERROR_BACKEND,
948 "CeedVector x has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
949 CeedCall(CeedVectorHasValidArray(y, &has_valid_array_y));
950 CeedCheck(has_valid_array_y, CeedVectorReturnCeed(w), CEED_ERROR_BACKEND,
951 "CeedVector y has no valid data, must set data with CeedVectorSetValue or CeedVectorSetArray");
952
953 // Return early for empty vectors
954 if (length_w == 0) return CEED_ERROR_SUCCESS;
955
956 // Backend implementation
957 if (w->PointwiseMult) {
958 CeedCall(w->PointwiseMult(w, x, y));
959 return CEED_ERROR_SUCCESS;
960 }
961
962 // Default implementation
963 if (x == w || y == w) {
964 CeedCall(CeedVectorGetArray(w, CEED_MEM_HOST, &w_array));
965 } else {
966 CeedCall(CeedVectorGetArrayWrite(w, CEED_MEM_HOST, &w_array));
967 }
968 if (x != w) {
969 CeedCall(CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array));
970 } else {
971 x_array = w_array;
972 }
973 if (y != w && y != x) {
974 CeedCall(CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array));
975 } else if (y == x) {
976 y_array = x_array;
977 } else if (y == w) {
978 y_array = w_array;
979 }
980
981 assert(w_array);
982 assert(x_array);
983 assert(y_array);
984
985 for (CeedSize i = 0; i < length_w; i++) w_array[i] = x_array[i] * y_array[i];
986
987 if (y != w && y != x) CeedCall(CeedVectorRestoreArrayRead(y, &y_array));
988 if (x != w) CeedCall(CeedVectorRestoreArrayRead(x, &x_array));
989 CeedCall(CeedVectorRestoreArray(w, &w_array));
990 return CEED_ERROR_SUCCESS;
991 }
992
993 /**
994 @brief Take the reciprocal of a `CeedVector`.
995
996 @param[in,out] vec `CeedVector` to take reciprocal
997
998 @return An error code: 0 - success, otherwise - failure
999
1000 @ref User
1001 **/
CeedVectorReciprocal(CeedVector vec)1002 int CeedVectorReciprocal(CeedVector vec) {
1003 bool has_valid_array = true;
1004 CeedSize length;
1005 CeedScalar *array;
1006
1007 CeedCall(CeedVectorHasValidArray(vec, &has_valid_array));
1008 CeedCheck(has_valid_array, CeedVectorReturnCeed(vec), CEED_ERROR_BACKEND,
1009 "CeedVector has no valid data to compute reciprocal, must set data with CeedVectorSetValue or CeedVectorSetArray");
1010
1011 // Check if vector data set
1012 CeedCheck(vec->state > 0, CeedVectorReturnCeed(vec), CEED_ERROR_INCOMPLETE, "CeedVector must have data set to take reciprocal");
1013
1014 // Return early for empty vector
1015 CeedCall(CeedVectorGetLength(vec, &length));
1016 if (length == 0) return CEED_ERROR_SUCCESS;
1017
1018 // Backend impl for GPU, if added
1019 if (vec->Reciprocal) {
1020 CeedCall(vec->Reciprocal(vec));
1021 return CEED_ERROR_SUCCESS;
1022 }
1023
1024 CeedCall(CeedVectorGetArray(vec, CEED_MEM_HOST, &array));
1025 for (CeedSize i = 0; i < length; i++) {
1026 if (fabs(array[i]) > CEED_EPSILON) array[i] = 1. / array[i];
1027 }
1028
1029 CeedCall(CeedVectorRestoreArray(vec, &array));
1030 return CEED_ERROR_SUCCESS;
1031 }
1032
1033 /**
1034 @brief Set the number of tabs to indent for @ref CeedVectorView() output
1035
1036 @param[in] vec `CeedVector` to set the number of view tabs
1037 @param[in] num_tabs Number of view tabs to set
1038
1039 @return Error code: 0 - success, otherwise - failure
1040
1041 @ref User
1042 **/
CeedVectorSetNumViewTabs(CeedVector vec,CeedInt num_tabs)1043 int CeedVectorSetNumViewTabs(CeedVector vec, CeedInt num_tabs) {
1044 CeedCall(CeedObjectSetNumViewTabs((CeedObject)vec, num_tabs));
1045 return CEED_ERROR_SUCCESS;
1046 }
1047
1048 /**
1049 @brief Get the number of tabs to indent for @ref CeedVectorView() output
1050
1051 @param[in] vec `CeedVector` to get the number of view tabs
1052 @param[out] num_tabs Number of view tabs
1053
1054 @return Error code: 0 - success, otherwise - failure
1055
1056 @ref User
1057 **/
CeedVectorGetNumViewTabs(CeedVector vec,CeedInt * num_tabs)1058 int CeedVectorGetNumViewTabs(CeedVector vec, CeedInt *num_tabs) {
1059 CeedCall(CeedObjectGetNumViewTabs((CeedObject)vec, num_tabs));
1060 return CEED_ERROR_SUCCESS;
1061 }
1062
1063 /**
1064 @brief View a `CeedVector`
1065
1066 Note: It is safe to use any unsigned values for `start` or `stop` and any nonzero integer for `step`.
1067 Any portion of the provided range that is outside the range of valid indices for the `CeedVector` will be ignored.
1068
1069 @param[in] vec `CeedVector` to view
1070 @param[in] start Index of first `CeedVector` entry to view in the range `[start, stop)`
1071 @param[in] stop One past the last element to view in the range, or `-1` for `length`
1072 @param[in] step Step between `CeedVector` entries to view
1073 @param[in] fp_fmt Printing format
1074 @param[in] stream Filestream to write to
1075
1076 @return An error code: 0 - success, otherwise - failure
1077
1078 @ref User
1079 **/
CeedVectorViewRange(CeedVector vec,CeedSize start,CeedSize stop,CeedInt step,const char * fp_fmt,FILE * stream)1080 int CeedVectorViewRange(CeedVector vec, CeedSize start, CeedSize stop, CeedInt step, const char *fp_fmt, FILE *stream) {
1081 char fmt[1024];
1082 char *tabs = NULL;
1083 CeedSize length;
1084 const CeedScalar *x;
1085
1086 CeedCheck(step != 0, CeedVectorReturnCeed(vec), CEED_ERROR_MINOR, "View range 'step' must be nonzero");
1087
1088 {
1089 CeedInt num_tabs = 0;
1090
1091 CeedCall(CeedVectorGetNumViewTabs(vec, &num_tabs));
1092 CeedCall(CeedCalloc(CEED_TAB_WIDTH * num_tabs + 1, &tabs));
1093 for (CeedInt i = 0; i < CEED_TAB_WIDTH * num_tabs; i++) tabs[i] = ' ';
1094 }
1095
1096 CeedCall(CeedVectorGetLength(vec, &length));
1097 fprintf(stream, "%sCeedVector length %" CeedSize_FMT "\n", tabs, length);
1098 if (start != 0 || stop != length || step != 1) {
1099 fprintf(stream, "%s start: %" CeedSize_FMT "\n%s stop: %" CeedSize_FMT "\n%s step: %" CeedInt_FMT "\n", tabs, start, tabs, stop, tabs, step);
1100 }
1101 if (start > length) start = length;
1102 if (stop == -1 || stop > length) stop = length;
1103
1104 snprintf(fmt, sizeof fmt, "%s %s\n", tabs, fp_fmt ? fp_fmt : "%g");
1105 CeedCall(CeedVectorGetArrayRead(vec, CEED_MEM_HOST, &x));
1106 for (CeedSize i = start; step > 0 ? (i < stop) : (i > stop); i += step) fprintf(stream, fmt, x[i]);
1107 CeedCall(CeedVectorRestoreArrayRead(vec, &x));
1108 if (stop != length) fprintf(stream, "%s ...\n", tabs);
1109 CeedCall(CeedFree(&tabs));
1110 return CEED_ERROR_SUCCESS;
1111 }
1112
1113 /**
1114 @brief View a `CeedVector`
1115
1116 @param[in] vec `CeedVector` to view
1117 @param[in] fp_fmt Printing format
1118 @param[in] stream Filestream to write to
1119
1120 @return An error code: 0 - success, otherwise - failure
1121
1122 @ref User
1123 **/
CeedVectorView(CeedVector vec,const char * fp_fmt,FILE * stream)1124 int CeedVectorView(CeedVector vec, const char *fp_fmt, FILE *stream) {
1125 CeedSize length;
1126
1127 CeedCall(CeedVectorGetLength(vec, &length));
1128 CeedCall(CeedVectorViewRange(vec, 0, length, 1, fp_fmt, stream));
1129 return CEED_ERROR_SUCCESS;
1130 }
1131
1132 /**
1133 @brief Get the `Ceed` associated with a `CeedVector`
1134
1135 @param[in] vec `CeedVector` to retrieve state
1136 @param[out] ceed Variable to store `Ceed`
1137
1138 @return An error code: 0 - success, otherwise - failure
1139
1140 @ref Advanced
1141 **/
CeedVectorGetCeed(CeedVector vec,Ceed * ceed)1142 int CeedVectorGetCeed(CeedVector vec, Ceed *ceed) {
1143 CeedCall(CeedObjectGetCeed((CeedObject)vec, ceed));
1144 return CEED_ERROR_SUCCESS;
1145 }
1146
1147 /**
1148 @brief Return the `Ceed` associated with a `CeedVector`
1149
1150 @param[in] vec `CeedVector` to retrieve state
1151
1152 @return `Ceed` associated with the `vec`
1153
1154 @ref Advanced
1155 **/
CeedVectorReturnCeed(CeedVector vec)1156 Ceed CeedVectorReturnCeed(CeedVector vec) { return CeedObjectReturnCeed((CeedObject)vec); }
1157
1158 /**
1159 @brief Get the length of a `CeedVector`
1160
1161 @param[in] vec `CeedVector` to retrieve length
1162 @param[out] length Variable to store length
1163
1164 @return An error code: 0 - success, otherwise - failure
1165
1166 @ref User
1167 **/
CeedVectorGetLength(CeedVector vec,CeedSize * length)1168 int CeedVectorGetLength(CeedVector vec, CeedSize *length) {
1169 *length = vec->length;
1170 return CEED_ERROR_SUCCESS;
1171 }
1172
1173 /**
1174 @brief Destroy a `CeedVector`
1175
1176 @param[in,out] vec `CeedVector` to destroy
1177
1178 @return An error code: 0 - success, otherwise - failure
1179
1180 @ref User
1181 **/
CeedVectorDestroy(CeedVector * vec)1182 int CeedVectorDestroy(CeedVector *vec) {
1183 if (!*vec || *vec == CEED_VECTOR_ACTIVE || *vec == CEED_VECTOR_NONE || CeedObjectDereference((CeedObject)*vec) > 0) {
1184 *vec = NULL;
1185 return CEED_ERROR_SUCCESS;
1186 }
1187 CeedCheck((*vec)->state % 2 == 0, CeedVectorReturnCeed(*vec), CEED_ERROR_ACCESS, "Cannot destroy CeedVector, the writable access lock is in use");
1188 CeedCheck((*vec)->num_readers == 0, CeedVectorReturnCeed(*vec), CEED_ERROR_ACCESS, "Cannot destroy CeedVector, a process has read access");
1189
1190 if ((*vec)->Destroy) CeedCall((*vec)->Destroy(*vec));
1191 CeedCall(CeedObjectDestroy_Private(&(*vec)->obj));
1192 CeedCall(CeedFree(vec));
1193 return CEED_ERROR_SUCCESS;
1194 }
1195
1196 /// @}
1197