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