xref: /libCEED/rust/libceed/src/vector.rs (revision edab612303ce5c583510db4ec9520fcda426d3c1)
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 //! A Ceed Vector constitutes the main data structure and serves as input/output
18 //! for Ceed Operators.
19 
20 use std::{
21     ops::{Deref, DerefMut},
22     os::raw::c_char,
23 };
24 
25 use crate::prelude::*;
26 
27 // -----------------------------------------------------------------------------
28 // CeedVector option
29 // -----------------------------------------------------------------------------
30 #[derive(Debug, Clone, Copy)]
31 pub enum VectorOpt<'a> {
32     Some(&'a Vector<'a>),
33     Active,
34     None,
35 }
36 /// Construct a VectorOpt reference from a Vector reference
37 impl<'a> From<&'a Vector<'_>> for VectorOpt<'a> {
38     fn from(vec: &'a Vector) -> Self {
39         debug_assert!(vec.ptr != unsafe { bind_ceed::CEED_VECTOR_NONE });
40         debug_assert!(vec.ptr != unsafe { bind_ceed::CEED_VECTOR_ACTIVE });
41         Self::Some(vec)
42     }
43 }
44 impl<'a> VectorOpt<'a> {
45     /// Transform a Rust libCEED VectorOpt into C libCEED CeedVector
46     pub(crate) fn to_raw(self) -> bind_ceed::CeedVector {
47         match self {
48             Self::Some(vec) => vec.ptr,
49             Self::None => unsafe { bind_ceed::CEED_VECTOR_NONE },
50             Self::Active => unsafe { bind_ceed::CEED_VECTOR_ACTIVE },
51         }
52     }
53 }
54 
55 // -----------------------------------------------------------------------------
56 // CeedVector context wrapper
57 // -----------------------------------------------------------------------------
58 #[derive(Debug)]
59 pub struct Vector<'a> {
60     ceed: &'a crate::Ceed,
61     pub(crate) ptr: bind_ceed::CeedVector,
62 }
63 impl From<&'_ Vector<'_>> for bind_ceed::CeedVector {
64     fn from(vec: &Vector) -> Self {
65         vec.ptr
66     }
67 }
68 
69 // -----------------------------------------------------------------------------
70 // Destructor
71 // -----------------------------------------------------------------------------
72 impl<'a> Drop for Vector<'a> {
73     fn drop(&mut self) {
74         let not_none_and_active = self.ptr != unsafe { bind_ceed::CEED_VECTOR_NONE }
75             && self.ptr != unsafe { bind_ceed::CEED_VECTOR_ACTIVE };
76 
77         if not_none_and_active {
78             unsafe { bind_ceed::CeedVectorDestroy(&mut self.ptr) };
79         }
80     }
81 }
82 
83 // -----------------------------------------------------------------------------
84 // Display
85 // -----------------------------------------------------------------------------
86 impl<'a> fmt::Display for Vector<'a> {
87     /// View a Vector
88     ///
89     /// ```
90     /// # use libceed::prelude::*;
91     /// # let ceed = libceed::Ceed::default_init();
92     /// let vec = libceed::vector::Vector::from_slice(&ceed, &[1., 2., 3.]).unwrap();
93     /// assert_eq!(
94     ///     vec.to_string(),
95     ///     "CeedVector length 3
96     ///     1.00000000
97     ///     2.00000000
98     ///     3.00000000
99     /// "
100     /// )
101     /// ```
102     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
103         let mut ptr = std::ptr::null_mut();
104         let mut sizeloc = crate::MAX_BUFFER_LENGTH;
105         let format = CString::new("%12.8f").expect("CString::new failed");
106         let format_c: *const c_char = format.into_raw();
107         let cstring = unsafe {
108             let file = bind_ceed::open_memstream(&mut ptr, &mut sizeloc);
109             bind_ceed::CeedVectorView(self.ptr, format_c, file);
110             bind_ceed::fclose(file);
111             CString::from_raw(ptr)
112         };
113         cstring.to_string_lossy().fmt(f)
114     }
115 }
116 
117 // -----------------------------------------------------------------------------
118 // Implementations
119 // -----------------------------------------------------------------------------
120 impl<'a> Vector<'a> {
121     // Constructors
122     pub fn create(ceed: &'a crate::Ceed, n: usize) -> crate::Result<Self> {
123         let n = i32::try_from(n).unwrap();
124         let mut ptr = std::ptr::null_mut();
125         let ierr = unsafe { bind_ceed::CeedVectorCreate(ceed.ptr, n, &mut ptr) };
126         ceed.check_error(ierr)?;
127         Ok(Self { ceed, ptr })
128     }
129 
130     pub(crate) fn from_raw(
131         ceed: &'a crate::Ceed,
132         ptr: bind_ceed::CeedVector,
133     ) -> crate::Result<Self> {
134         Ok(Self { ceed, ptr })
135     }
136 
137     /// Create a Vector from a slice
138     ///
139     /// # arguments
140     ///
141     /// * `slice` - values to initialize vector with
142     ///
143     /// ```
144     /// # use libceed::prelude::*;
145     /// # let ceed = libceed::Ceed::default_init();
146     /// let vec = vector::Vector::from_slice(&ceed, &[1., 2., 3.]).unwrap();
147     /// assert_eq!(vec.length(), 3, "Incorrect length from slice");
148     /// ```
149     pub fn from_slice(ceed: &'a crate::Ceed, v: &[f64]) -> crate::Result<Self> {
150         let mut x = Self::create(ceed, v.len())?;
151         x.set_slice(v)?;
152         Ok(x)
153     }
154 
155     /// Create a Vector from a mutable array reference
156     ///
157     /// # arguments
158     ///
159     /// * `slice` - values to initialize vector with
160     ///
161     /// ```
162     /// # use libceed::prelude::*;
163     /// # let ceed = libceed::Ceed::default_init();
164     /// let mut rust_vec = vec![1., 2., 3.];
165     /// let vec = libceed::vector::Vector::from_array(&ceed, &mut rust_vec).unwrap();
166     ///
167     /// assert_eq!(vec.length(), 3, "Incorrect length from slice");
168     /// ```
169     pub fn from_array(ceed: &'a crate::Ceed, v: &mut [f64]) -> crate::Result<Self> {
170         let x = Self::create(ceed, v.len())?;
171         let (host, user_pointer) = (
172             crate::MemType::Host as bind_ceed::CeedMemType,
173             crate::CopyMode::UsePointer as bind_ceed::CeedCopyMode,
174         );
175         let v = v.as_ptr() as *mut f64;
176         let ierr = unsafe { bind_ceed::CeedVectorSetArray(x.ptr, host, user_pointer, v) };
177         ceed.check_error(ierr)?;
178         Ok(x)
179     }
180 
181     /// Returns the length of a CeedVector
182     ///
183     /// ```
184     /// # use libceed::prelude::*;
185     /// # let ceed = libceed::Ceed::default_init();
186     /// let vec = ceed.vector(10).unwrap();
187     ///
188     /// let n = vec.length();
189     /// assert_eq!(n, 10, "Incorrect length");
190     /// ```
191     pub fn length(&self) -> usize {
192         let mut n = 0;
193         unsafe { bind_ceed::CeedVectorGetLength(self.ptr, &mut n) };
194         usize::try_from(n).unwrap()
195     }
196 
197     /// Returns the length of a CeedVector
198     ///
199     /// ```
200     /// # use libceed::prelude::*;
201     /// # let ceed = libceed::Ceed::default_init();
202     /// let vec = ceed.vector(10).unwrap();
203     /// assert_eq!(vec.len(), 10, "Incorrect length");
204     /// ```
205     pub fn len(&self) -> usize {
206         self.length()
207     }
208 
209     /// Set the CeedVector to a constant value
210     ///
211     /// # arguments
212     ///
213     /// * `val` - Value to be used
214     ///
215     /// ```
216     /// # use libceed::prelude::*;
217     /// # let ceed = libceed::Ceed::default_init();
218     /// let len = 10;
219     /// let mut vec = ceed.vector(len).unwrap();
220     ///
221     /// let val = 42.0;
222     /// vec.set_value(val).unwrap();
223     ///
224     /// vec.view().iter().for_each(|v| {
225     ///     assert_eq!(*v, val, "Value not set correctly");
226     /// });
227     /// ```
228     pub fn set_value(&mut self, value: f64) -> crate::Result<i32> {
229         let ierr = unsafe { bind_ceed::CeedVectorSetValue(self.ptr, value) };
230         self.ceed.check_error(ierr)
231     }
232 
233     /// Set values from a slice of the same length
234     ///
235     /// # arguments
236     ///
237     /// * `slice` - values to into self; length must match
238     ///
239     /// ```
240     /// # use libceed::prelude::*;
241     /// # let ceed = libceed::Ceed::default_init();
242     /// let mut vec = ceed.vector(4).unwrap();
243     /// vec.set_slice(&[10., 11., 12., 13.]).unwrap();
244     ///
245     /// vec.view().iter().enumerate().for_each(|(i, v)| {
246     ///     assert_eq!(*v, 10. + i as f64, "Slice not set correctly");
247     /// });
248     /// ```
249     pub fn set_slice(&mut self, slice: &[f64]) -> crate::Result<i32> {
250         assert_eq!(self.length(), slice.len());
251         let (host, copy_mode) = (
252             crate::MemType::Host as bind_ceed::CeedMemType,
253             crate::CopyMode::CopyValues as bind_ceed::CeedCopyMode,
254         );
255         let ierr = unsafe {
256             bind_ceed::CeedVectorSetArray(self.ptr, host, copy_mode, slice.as_ptr() as *mut f64)
257         };
258         self.ceed.check_error(ierr)
259     }
260 
261     /// Sync the CeedVector to a specified memtype
262     ///
263     /// # arguments
264     ///
265     /// * `mtype` - Memtype to be synced
266     ///
267     /// ```
268     /// # use libceed::prelude::*;
269     /// # let ceed = libceed::Ceed::default_init();
270     /// let len = 10;
271     /// let mut vec = ceed.vector(len).unwrap();
272     ///
273     /// let val = 42.0;
274     /// vec.set_value(val);
275     /// vec.sync(MemType::Host).unwrap();
276     ///
277     /// vec.view().iter().for_each(|v| {
278     ///     assert_eq!(*v, val, "Value not set correctly");
279     /// });
280     /// ```
281     pub fn sync(&self, mtype: crate::MemType) -> crate::Result<i32> {
282         let ierr =
283             unsafe { bind_ceed::CeedVectorSyncArray(self.ptr, mtype as bind_ceed::CeedMemType) };
284         self.ceed.check_error(ierr)
285     }
286 
287     /// Create an immutable view
288     ///
289     /// ```
290     /// # use libceed::prelude::*;
291     /// # let ceed = libceed::Ceed::default_init();
292     /// let vec = ceed.vector_from_slice(&[10., 11., 12., 13.]).unwrap();
293     ///
294     /// let v = vec.view();
295     /// assert_eq!(v[0..2], [10., 11.]);
296     ///
297     /// // It is valid to have multiple immutable views
298     /// let w = vec.view();
299     /// assert_eq!(v[1..], w[1..]);
300     /// ```
301     pub fn view(&self) -> VectorView {
302         VectorView::new(self)
303     }
304 
305     /// Create an mutable view
306     ///
307     /// ```
308     /// # use libceed::prelude::*;
309     /// # let ceed = libceed::Ceed::default_init();
310     /// let mut vec = ceed.vector_from_slice(&[10., 11., 12., 13.]).unwrap();
311     ///
312     /// {
313     ///     let mut v = vec.view_mut();
314     ///     v[2] = 9.;
315     /// }
316     ///
317     /// let w = vec.view();
318     /// assert_eq!(w[2], 9., "View did not mutate data");
319     /// ```
320     pub fn view_mut(&mut self) -> VectorViewMut {
321         VectorViewMut::new(self)
322     }
323 
324     /// Return the norm of a CeedVector
325     ///
326     /// # arguments
327     ///
328     /// * `ntype` - Norm type One, Two, or Max
329     ///
330     /// ```
331     /// # use libceed::prelude::*;
332     /// # let ceed = libceed::Ceed::default_init();
333     /// let vec = ceed.vector_from_slice(&[1., 2., 3., 4.]).unwrap();
334     ///
335     /// let max_norm = vec.norm(NormType::Max).unwrap();
336     /// assert_eq!(max_norm, 4.0, "Incorrect Max norm");
337     ///
338     /// let l1_norm = vec.norm(NormType::One).unwrap();
339     /// assert_eq!(l1_norm, 10., "Incorrect L1 norm");
340     ///
341     /// let l2_norm = vec.norm(NormType::Two).unwrap();
342     /// assert!((l2_norm - 5.477) < 1e-3, "Incorrect L2 norm");
343     /// ```
344     pub fn norm(&self, ntype: crate::NormType) -> crate::Result<f64> {
345         let mut res: f64 = 0.0;
346         let ierr = unsafe {
347             bind_ceed::CeedVectorNorm(self.ptr, ntype as bind_ceed::CeedNormType, &mut res)
348         };
349         self.ceed.check_error(ierr)?;
350         Ok(res)
351     }
352 
353     /// Compute x = alpha x for a CeedVector
354     ///
355     /// # arguments
356     ///
357     /// * `alpha` - scaling factor
358     ///
359     /// ```
360     /// # use libceed::prelude::*;
361     /// # let ceed = libceed::Ceed::default_init();
362     /// let mut vec = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
363     ///
364     /// vec = vec.scale(-1.0).unwrap();
365     /// vec.view().iter().enumerate().for_each(|(i, &v)| {
366     ///     assert_eq!(v, -(i as f64), "Value not set correctly");
367     /// });
368     /// ```
369     #[allow(unused_mut)]
370     pub fn scale(mut self, alpha: f64) -> crate::Result<Self> {
371         let ierr = unsafe { bind_ceed::CeedVectorScale(self.ptr, alpha) };
372         self.ceed.check_error(ierr)?;
373         Ok(self)
374     }
375 
376     /// Compute y = alpha x + y for a pair of CeedVectors
377     ///
378     /// # arguments
379     ///
380     /// * `alpha` - scaling factor
381     /// * `x`     - second vector, must be different than self
382     ///
383     /// ```
384     /// # use libceed::prelude::*;
385     /// # let ceed = libceed::Ceed::default_init();
386     /// let x = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
387     /// let mut y = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
388     ///
389     /// y = y.axpy(-0.5, &x).unwrap();
390     /// y.view().iter().enumerate().for_each(|(i, &v)| {
391     ///     assert_eq!(v, (i as f64) / 2.0, "Value not set correctly");
392     /// });
393     /// ```
394     #[allow(unused_mut)]
395     pub fn axpy(mut self, alpha: f64, x: &crate::Vector) -> crate::Result<Self> {
396         let ierr = unsafe { bind_ceed::CeedVectorAXPY(self.ptr, alpha, x.ptr) };
397         self.ceed.check_error(ierr)?;
398         Ok(self)
399     }
400 
401     /// Compute the pointwise multiplication w = x .* y for CeedVectors
402     ///
403     /// # arguments
404     ///
405     /// * `x` - first vector for product
406     /// * `y` - second vector for product
407     ///
408     /// ```
409     /// # use libceed::prelude::*;
410     /// # let ceed = libceed::Ceed::default_init();
411     /// let mut w = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
412     /// let x = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
413     /// let y = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
414     ///
415     /// w = w.pointwise_mult(&x, &y).unwrap();
416     /// w.view().iter().enumerate().for_each(|(i, &v)| {
417     ///     assert_eq!(v, (i as f64).powf(2.0), "Value not set correctly");
418     /// });
419     /// ```
420     #[allow(unused_mut)]
421     pub fn pointwise_mult(mut self, x: &crate::Vector, y: &crate::Vector) -> crate::Result<Self> {
422         let ierr = unsafe { bind_ceed::CeedVectorPointwiseMult(self.ptr, x.ptr, y.ptr) };
423         self.ceed.check_error(ierr)?;
424         Ok(self)
425     }
426 
427     /// Compute the pointwise multiplication w = w .* x for CeedVectors
428     ///
429     /// # arguments
430     ///
431     /// * `x` - second vector for product
432     ///
433     /// ```
434     /// # use libceed::prelude::*;
435     /// # let ceed = libceed::Ceed::default_init();
436     /// let mut w = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
437     /// let x = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
438     ///
439     /// w = w.pointwise_scale(&x).unwrap();
440     /// w.view().iter().enumerate().for_each(|(i, &v)| {
441     ///     assert_eq!(v, (i as f64).powf(2.0), "Value not set correctly");
442     /// });
443     /// ```
444     #[allow(unused_mut)]
445     pub fn pointwise_scale(mut self, x: &crate::Vector) -> crate::Result<Self> {
446         let ierr = unsafe { bind_ceed::CeedVectorPointwiseMult(self.ptr, self.ptr, x.ptr) };
447         self.ceed.check_error(ierr)?;
448         Ok(self)
449     }
450 
451     /// Compute the pointwise multiplication w = w .* w for a CeedVector
452     ///
453     /// ```
454     /// # use libceed::prelude::*;
455     /// # let ceed = libceed::Ceed::default_init();
456     /// let mut w = ceed.vector_from_slice(&[0., 1., 2., 3., 4.]).unwrap();
457     ///
458     /// w = w.pointwise_square().unwrap();
459     /// w.view().iter().enumerate().for_each(|(i, &v)| {
460     ///     assert_eq!(v, (i as f64).powf(2.0), "Value not set correctly");
461     /// });
462     /// ```
463     #[allow(unused_mut)]
464     pub fn pointwise_square(mut self) -> crate::Result<Self> {
465         let ierr = unsafe { bind_ceed::CeedVectorPointwiseMult(self.ptr, self.ptr, self.ptr) };
466         self.ceed.check_error(ierr)?;
467         Ok(self)
468     }
469 }
470 
471 // -----------------------------------------------------------------------------
472 // Vector Viewer
473 // -----------------------------------------------------------------------------
474 /// A (host) view of a Vector with Deref to slice.  We can't make
475 /// Vector itself Deref to slice because we can't handle the drop to
476 /// call bind_ceed::CeedVectorRestoreArrayRead().
477 #[derive(Debug)]
478 pub struct VectorView<'a> {
479     vec: &'a Vector<'a>,
480     array: *const f64,
481 }
482 
483 impl<'a> VectorView<'a> {
484     /// Construct a VectorView from a Vector reference
485     fn new(vec: &'a Vector) -> Self {
486         let mut array = std::ptr::null();
487         unsafe {
488             bind_ceed::CeedVectorGetArrayRead(
489                 vec.ptr,
490                 crate::MemType::Host as bind_ceed::CeedMemType,
491                 &mut array,
492             );
493         }
494         Self {
495             vec: vec,
496             array: array,
497         }
498     }
499 }
500 
501 // Destructor
502 impl<'a> Drop for VectorView<'a> {
503     fn drop(&mut self) {
504         unsafe {
505             bind_ceed::CeedVectorRestoreArrayRead(self.vec.ptr, &mut self.array);
506         }
507     }
508 }
509 
510 // Data access
511 impl<'a> Deref for VectorView<'a> {
512     type Target = [f64];
513     fn deref(&self) -> &[f64] {
514         unsafe { std::slice::from_raw_parts(self.array, self.vec.len()) }
515     }
516 }
517 
518 // Viewing
519 impl<'a> fmt::Display for VectorView<'a> {
520     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
521         write!(f, "VectorView({:?})", self.deref())
522     }
523 }
524 
525 // -----------------------------------------------------------------------------
526 // Vector Viewer Mutable
527 // -----------------------------------------------------------------------------
528 /// A mutable (host) view of a Vector with Deref to slice.
529 #[derive(Debug)]
530 pub struct VectorViewMut<'a> {
531     vec: &'a Vector<'a>,
532     array: *mut f64,
533 }
534 
535 impl<'a> VectorViewMut<'a> {
536     /// Construct a VectorViewMut from a Vector reference
537     fn new(vec: &'a mut Vector) -> Self {
538         let mut ptr = std::ptr::null_mut();
539         unsafe {
540             bind_ceed::CeedVectorGetArray(
541                 vec.ptr,
542                 crate::MemType::Host as bind_ceed::CeedMemType,
543                 &mut ptr,
544             );
545         }
546         Self {
547             vec: vec,
548             array: ptr,
549         }
550     }
551 }
552 
553 // Destructor
554 impl<'a> Drop for VectorViewMut<'a> {
555     fn drop(&mut self) {
556         unsafe {
557             bind_ceed::CeedVectorRestoreArray(self.vec.ptr, &mut self.array);
558         }
559     }
560 }
561 
562 // Data access
563 impl<'a> Deref for VectorViewMut<'a> {
564     type Target = [f64];
565     fn deref(&self) -> &[f64] {
566         unsafe { std::slice::from_raw_parts(self.array, self.vec.len()) }
567     }
568 }
569 
570 // Mutable data access
571 impl<'a> DerefMut for VectorViewMut<'a> {
572     fn deref_mut(&mut self) -> &mut [f64] {
573         unsafe { std::slice::from_raw_parts_mut(self.array, self.vec.len()) }
574     }
575 }
576 
577 // Viewing
578 impl<'a> fmt::Display for VectorViewMut<'a> {
579     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
580         write!(f, "VectorViewMut({:?})", self.deref())
581     }
582 }
583 
584 // -----------------------------------------------------------------------------
585