diffsol_la/scalar/
mod.rs

1use std::{
2    fmt::Display,
3    ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign},
4};
5
6#[cfg(any(feature = "cuda", feature = "cuda-oxide"))]
7pub mod cuda;
8
9use crate::vector::VectorView;
10
11/// A scalar type suitable for numerical computations in ODE solvers.
12///
13/// This trait aggregates the crate-local numeric requirements shared across diffsol.
14///
15/// # Implementations
16/// diffsol provides implementations for `f64` and `f32`.
17///
18/// # Examples
19/// ```
20/// use diffsol_la::Scalar;
21///
22/// fn compute<T: Scalar>(x: T, y: T) -> T {
23///     x * x + y
24/// }
25/// ```
26pub trait Scalar:
27    num_traits::Signed
28    + num_traits::Pow<Self, Output = Self>
29    + num_traits::Pow<i32, Output = Self>
30    + num_traits::FromPrimitive
31    + num_traits::ToPrimitive
32    + Display
33    + std::fmt::Debug
34    + Copy
35    + PartialEq
36    + PartialOrd
37    + AddAssign
38    + SubAssign
39    + MulAssign
40    + std::ops::DivAssign
41    + Send
42    + Sync
43    + 'static
44{
45    /// Machine epsilon for this scalar type (smallest representable positive value such that 1.0 + EPSILON != 1.0).
46    const EPSILON: Self;
47    /// Positive infinity value for this scalar type.
48    const INFINITY: Self;
49    /// Not-a-Number (NaN) value for this scalar type.
50    const NAN: Self;
51    /// Check if this value is NaN.
52    fn is_nan(self) -> bool;
53    /// Square root.
54    fn sqrt(self) -> Self;
55    /// Cube root.
56    fn cbrt(self) -> Self;
57    /// Exponential.
58    fn exp(self) -> Self;
59    /// Sine.
60    fn sin(self) -> Self;
61    /// Cosine.
62    fn cos(self) -> Self;
63    /// Maximum of two values.
64    fn max(self, other: Self) -> Self;
65
66    /// Addition, potentially reassociated/relaxed for speed under the `algebraic-ops` feature.
67    fn algebraic_add(self, other: Self) -> Self {
68        self + other
69    }
70    /// Subtraction, potentially reassociated/relaxed for speed under the `algebraic-ops` feature.
71    fn algebraic_sub(self, other: Self) -> Self {
72        self - other
73    }
74    /// Multiplication, potentially reassociated/relaxed for speed under the `algebraic-ops` feature.
75    fn algebraic_mul(self, other: Self) -> Self {
76        self * other
77    }
78    /// Division, potentially reassociated/relaxed for speed under the `algebraic-ops` feature.
79    fn algebraic_div(self, other: Self) -> Self {
80        self / other
81    }
82}
83
84/// A [`Scalar`] that also satisfies nalgebra's numeric field requirements.
85pub trait NalgebraScalar:
86    Scalar + nalgebra::Scalar + nalgebra::SimdRealField + nalgebra::ComplexField<RealField = Self>
87{
88}
89
90impl<T> NalgebraScalar for T where
91    T: Scalar
92        + nalgebra::Scalar
93        + nalgebra::SimdRealField
94        + nalgebra::ComplexField<RealField = Self>
95{
96}
97
98/// A [`Scalar`] that also satisfies faer's numeric field requirements.
99pub trait FaerScalar: Scalar + faer::traits::ComplexField + faer::traits::RealField {}
100
101impl<T> FaerScalar for T where T: Scalar + faer::traits::ComplexField + faer::traits::RealField {}
102
103/// The index type used throughout diffsol for indexing vectors and matrices.
104pub type IndexType = usize;
105
106impl Scalar for f64 {
107    const EPSILON: Self = f64::EPSILON;
108    const INFINITY: Self = f64::INFINITY;
109    const NAN: Self = f64::NAN;
110    fn is_nan(self) -> bool {
111        self.is_nan()
112    }
113    fn sqrt(self) -> Self {
114        self.sqrt()
115    }
116    fn cbrt(self) -> Self {
117        self.cbrt()
118    }
119    fn exp(self) -> Self {
120        self.exp()
121    }
122    fn sin(self) -> Self {
123        self.sin()
124    }
125    fn cos(self) -> Self {
126        self.cos()
127    }
128    fn max(self, other: Self) -> Self {
129        self.max(other)
130    }
131    #[cfg(feature = "algebraic-ops")]
132    fn algebraic_add(self, other: Self) -> Self {
133        self.algebraic_add(other)
134    }
135    #[cfg(feature = "algebraic-ops")]
136    fn algebraic_sub(self, other: Self) -> Self {
137        self.algebraic_sub(other)
138    }
139    #[cfg(feature = "algebraic-ops")]
140    fn algebraic_mul(self, other: Self) -> Self {
141        self.algebraic_mul(other)
142    }
143    #[cfg(feature = "algebraic-ops")]
144    fn algebraic_div(self, other: Self) -> Self {
145        self.algebraic_div(other)
146    }
147}
148
149impl Scalar for f32 {
150    const EPSILON: Self = f32::EPSILON;
151    const INFINITY: Self = f32::INFINITY;
152    const NAN: Self = f32::NAN;
153    fn is_nan(self) -> bool {
154        self.is_nan()
155    }
156    fn sqrt(self) -> Self {
157        self.sqrt()
158    }
159    fn cbrt(self) -> Self {
160        self.cbrt()
161    }
162    fn exp(self) -> Self {
163        self.exp()
164    }
165    fn sin(self) -> Self {
166        self.sin()
167    }
168    fn cos(self) -> Self {
169        self.cos()
170    }
171    fn max(self, other: Self) -> Self {
172        self.max(other)
173    }
174    #[cfg(feature = "algebraic-ops")]
175    fn algebraic_add(self, other: Self) -> Self {
176        self.algebraic_add(other)
177    }
178    #[cfg(feature = "algebraic-ops")]
179    fn algebraic_sub(self, other: Self) -> Self {
180        self.algebraic_sub(other)
181    }
182    #[cfg(feature = "algebraic-ops")]
183    fn algebraic_mul(self, other: Self) -> Self {
184        self.algebraic_mul(other)
185    }
186    #[cfg(feature = "algebraic-ops")]
187    fn algebraic_div(self, other: Self) -> Self {
188        self.algebraic_div(other)
189    }
190}
191
192impl<T: FaerScalar> From<faer::Scale<T>> for Scale<T> {
193    fn from(s: faer::Scale<T>) -> Self {
194        Scale(s.0)
195    }
196}
197impl<T: FaerScalar> From<Scale<T>> for faer::Scale<T> {
198    fn from(s: Scale<T>) -> Self {
199        faer::Scale(s.value())
200    }
201}
202impl<T: Scalar> From<T> for Scale<T> {
203    fn from(s: T) -> Self {
204        Scale(s)
205    }
206}
207
208/// A wrapper for scalar values used when scaling vectors and matrices.
209#[derive(Copy, Clone, Debug)]
210pub struct Scale<E: Scalar>(pub E);
211
212impl<E: Scalar> Scale<E> {
213    /// Get the underlying scalar value.
214    #[inline]
215    pub fn value(self) -> E {
216        self.0
217    }
218}
219
220/// Create a `Scale` wrapper from a scalar value.
221///
222/// This is a convenience function equivalent to `Scale(value)`.
223#[inline]
224pub fn scale<E: Scalar>(value: E) -> Scale<E> {
225    Scale(value)
226}
227
228macro_rules! impl_bin_op {
229    ($trait:ident, $method:ident, $operator:tt) => {
230        impl<E: Scalar> $trait<Scale<E>> for Scale<E> {
231            type Output = Scale<E>;
232
233            #[inline]
234            fn $method(self, rhs: Scale<E>) -> Self::Output {
235                Scale(self.0 $operator rhs.0)
236            }
237        }
238    };
239}
240
241macro_rules! impl_assign_bin_op {
242    ($trait:ident, $method:ident, $operator:tt) => {
243        impl<E: Scalar> $trait<Scale<E>> for Scale<E> {
244            #[inline]
245            fn $method(&mut self, rhs: Scale<E>) {
246                self.0 = self.0 $operator rhs.0
247            }
248        }
249    };
250}
251
252impl_bin_op!(Mul, mul, *);
253impl_bin_op!(Add, add, +);
254impl_bin_op!(Sub, sub, -);
255
256impl_assign_bin_op!(MulAssign, mul_assign, *);
257impl_assign_bin_op!(AddAssign, add_assign, +);
258impl_assign_bin_op!(SubAssign, sub_assign, -);
259
260impl<E: Scalar> PartialEq for Scale<E> {
261    #[inline]
262    fn eq(&self, rhs: &Self) -> bool {
263        self.0 == rhs.0
264    }
265}
266
267impl<V: VectorView<'static, T = E>, E: Scalar> Mul<V> for Scale<E> {
268    type Output = V::Owned;
269    #[inline]
270    fn mul(self, rhs: V) -> Self::Output {
271        rhs * scale(self.0)
272    }
273}
274
275#[test]
276fn test_scale() {
277    assert_eq!(scale(2.0) * scale(3.0), scale(6.0));
278}
diffsol_la/scalar/mod.rs

diffsol_la/scalar/
mod.rs
