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