xref: /libCEED/interface/ceed-vector.c (revision 43e1b16f0c7c5dbfae04543c0675e255f265ff67)
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