SIMD CG transform and related improvements

.github/workflows/nalgebra-ci-build.yml

@@ -64,7 +64,7 @@ jobs:
           override: true
       - uses: actions/checkout@v4
       - name: check
-        run: cargo check --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser,rayon;
+        run: cargo check --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,wide-tests,rkyv-safe-deser,rayon;
   test-nalgebra:
     runs-on: ubuntu-latest
     #    env:
@@ -81,7 +81,7 @@ jobs:
           override: true
       - uses: actions/checkout@v4
       - name: test
-        run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,rkyv-safe-deser,rayon;
+        run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests,wide-tests,rkyv-safe-deser,rayon;
   test-nalgebra-glm:
     runs-on: ubuntu-latest
     steps:

Cargo.toml

@@ -98,6 +98,7 @@ arbitrary = ["quickcheck"]
 proptest-support = ["proptest"]
 slow-tests = []
 rkyv-safe-deser = ["rkyv-serialize", "rkyv/validation"]
+wide-tests = ["simba/wide", "simba/rand"]
 
 [dependencies]
 nalgebra-macros = { version = "0.3.0", path = "nalgebra-macros", optional = true }
@@ -168,14 +169,15 @@ required-features = ["compare"]
 name = "nalgebra_bench"
 harness = false
 path = "benches/lib.rs"
-required-features = ["rand"]
+required-features = ["rand", "wide-tests"]
 
 #[profile.bench]
 #opt-level = 0
 #lto = false
 
 [profile.bench]
 lto = true
+codegen-units = 1
 
 [package.metadata.docs.rs]
 # Enable all the features when building the docs on docs.rs

benches/core/cg.rs

@@ -0,0 +1,85 @@
+use na::{Matrix3, Matrix4, Orthographic3, Point2, Point3, Vector2, Vector3};
+use rand::{Rng, SeedableRng};
+use rand_isaac::IsaacRng;
+use simba::simd::WideF32x4;
+
+#[path = "../common/macros.rs"]
+mod macros;
+
+bench_binop_ref!(
+    mat3_transform_vector2,
+    Matrix3<f32>,
+    Vector2<f32>,
+    transform_vector
+);
+bench_binop_ref!(
+    mat4_transform_vector3,
+    Matrix4<f32>,
+    Vector3<f32>,
+    transform_vector
+);
+bench_binop_ref!(
+    mat3_transform_point2,
+    Matrix3<f32>,
+    Point2<f32>,
+    transform_point
+);
+bench_binop_ref!(
+    mat4_transform_point3,
+    Matrix4<f32>,
+    Point3<f32>,
+    transform_point
+);
+
+bench_binop_ref!(
+    mat3_transform_vector2_x4wide,
+    Matrix3<WideF32x4>,
+    Vector2<WideF32x4>,
+    simd_transform_vector
+);
+bench_binop_ref!(
+    mat4_transform_vector3_x4wide,
+    Matrix4<WideF32x4>,
+    Vector3<WideF32x4>,
+    simd_transform_vector
+);
+bench_binop_ref!(
+    mat3_transform_point2_x4wide,
+    Matrix3<WideF32x4>,
+    Point2<WideF32x4>,
+    simd_transform_point
+);
+bench_binop_ref!(
+    mat4_transform_point3_x4wide,
+    Matrix4<WideF32x4>,
+    Point3<WideF32x4>,
+    simd_transform_point
+);
+
+fn mat4_transform_vector3_no_division(bench: &mut criterion::Criterion) {
+    let mut rng = IsaacRng::seed_from_u64(0);
+    let orthographic = Orthographic3::from_fov(
+        rng.random_range(0.5..2.0),
+        rng.random_range(30.0..90.0),
+        rng.random_range(0.5..1.5),
+        rng.random_range(2.0..1000.0),
+    )
+    .to_homogeneous();
+    let vector = rng.random();
+    bench.bench_function("mat4_transform_vector3_no_division", move |bh| {
+        bh.iter(|| orthographic.transform_vector(&vector))
+    });
+}
+
+criterion_group!(
+    cg,
+    mat3_transform_vector2,
+    mat4_transform_vector3,
+    mat3_transform_point2,
+    mat4_transform_point3,
+    mat3_transform_vector2_x4wide,
+    mat4_transform_vector3_x4wide,
+    mat3_transform_point2_x4wide,
+    mat4_transform_point3_x4wide,
+    mat4_transform_vector3_no_division,
+);

benches/core/mod.rs

@@ -1,5 +1,7 @@
+pub use self::cg::cg;
 pub use self::matrix::matrix;
 pub use self::vector::vector;
 
+mod cg;
 mod matrix;
 mod vector;

benches/lib.rs

@@ -21,6 +21,7 @@ fn reproducible_dmatrix(nrows: usize, ncols: usize) -> DMatrix<f64> {
 }
 
 criterion_main!(
+    core::cg,
     core::matrix,
     core::vector,
     geometry::quaternion,

src/base/cg.rs

@@ -20,6 +20,7 @@ use crate::geometry::{
 };
 
 use simba::scalar::{ClosedAddAssign, ClosedMulAssign, RealField};
+use simba::simd::SimdRealField;
 
 /// # Translation and scaling in any dimension
 impl<T, D: DimName> OMatrix<T, D, D>
@@ -402,6 +403,9 @@ where
         + Allocator<DimNameDiff<D, U1>, DimNameDiff<D, U1>>,
 {
     /// Transforms the given vector, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting vector is divided by the last component of the homogeneous
+    /// coordinates if it is not zero or returned unchanged otherwise.
     #[inline]
     pub fn transform_vector(
         &self,
@@ -423,8 +427,40 @@ where
     }
 }
 
+/// # Transformation of vectors and points
+impl<T: SimdRealField, D: DimNameSub<U1>, S: Storage<T, D, D>> SquareMatrix<T, D, S>
+where
+    DefaultAllocator: Allocator<D, D>
+        + Allocator<DimNameDiff<D, U1>>
+        + Allocator<DimNameDiff<D, U1>, DimNameDiff<D, U1>>,
+{
+    /// Transforms the given vector, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting vector is divided by the last component of the homogeneous
+    /// coordinates if it is not zero or returned unchanged otherwise.
+    #[inline]
+    pub fn simd_transform_vector(
+        &self,
+        v: &OVector<T, DimNameDiff<D, U1>>,
+    ) -> OVector<T, DimNameDiff<D, U1>> {
+        let transform = self.generic_view(
+            (0, 0),
+            (DimNameDiff::<D, U1>::name(), DimNameDiff::<D, U1>::name()),
+        );
+        let normalizer =
+            self.generic_view((D::DIM - 1, 0), (Const::<1>, DimNameDiff::<D, U1>::name()));
+        let n = normalizer.tr_dot(v);
+
+        let n = n.clone().select(n.simd_ne(T::zero()), T::one());
+        transform * (v / n)
+    }
+}
+
 impl<T: RealField, S: Storage<T, Const<3>, Const<3>>> SquareMatrix<T, Const<3>, S> {
     /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting point is divided by the `z` component if it is not zero
+    /// or returned unchanged otherwise.
     #[inline]
     pub fn transform_point(&self, pt: &Point<T, 2>) -> Point<T, 2> {
         let transform = self.fixed_view::<2, 2>(0, 0);
@@ -440,8 +476,28 @@ impl<T: RealField, S: Storage<T, Const<3>, Const<3>>> SquareMatrix<T, Const<3>,
     }
 }
 
+impl<T: SimdRealField, S: Storage<T, Const<3>, Const<3>>> SquareMatrix<T, Const<3>, S> {
+    /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting point is divided by the `z` component if it is not zero
+    /// or returned unchanged otherwise.
+    #[inline]
+    pub fn simd_transform_point(&self, pt: &Point<T, 2>) -> Point<T, 2> {
+        let transform = self.fixed_view::<2, 2>(0, 0);
+        let translation = self.fixed_view::<2, 1>(0, 2);
+        let normalizer = self.fixed_view::<1, 2>(2, 0);
+        let n = normalizer.tr_dot(&pt.coords) + unsafe { self.get_unchecked((2, 2)).clone() };
+
+        let n = n.clone().select(n.simd_ne(T::zero()), T::one());
+        (transform * pt + translation) / n
+    }
+}
+
 impl<T: RealField, S: Storage<T, Const<4>, Const<4>>> SquareMatrix<T, Const<4>, S> {
     /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting vector is divided by the `w` component if it is not zero
+    /// or returned unchanged otherwise.
     #[inline]
     pub fn transform_point(&self, pt: &Point<T, 3>) -> Point<T, 3> {
         let transform = self.fixed_view::<3, 3>(0, 0);
@@ -456,3 +512,20 @@ impl<T: RealField, S: Storage<T, Const<4>, Const<4>>> SquareMatrix<T, Const<4>,
         }
     }
 }
+
+impl<T: SimdRealField, S: Storage<T, Const<4>, Const<4>>> SquareMatrix<T, Const<4>, S> {
+    /// Transforms the given point, assuming the matrix `self` uses homogeneous coordinates.
+    ///
+    /// Each component of the resulting vector is divided by the `w` component if it is not zero
+    /// or returned unchanged otherwise.
+    #[inline]
+    pub fn simd_transform_point(&self, pt: &Point<T, 3>) -> Point<T, 3> {
+        let transform = self.fixed_view::<3, 3>(0, 0);
+        let translation = self.fixed_view::<3, 1>(0, 3);
+        let normalizer = self.fixed_view::<1, 3>(3, 0);
+        let n = normalizer.tr_dot(&pt.coords) + unsafe { self.get_unchecked((3, 3)).clone() };
+
+        let n = n.clone().select(n.simd_ne(T::zero()), T::one());
+        (transform * pt + translation) / n
+    }
+}

src/geometry/perspective.rs

@@ -233,6 +233,8 @@ impl<T: RealField> Perspective3<T> {
 
     // TODO: when we get specialization, specialize the Mul impl instead.
     /// Projects a point. Faster than matrix multiplication.
+    ///
+    /// Each component of the resulting point will be NaN or Inf if `p.z` is zero.
     #[inline]
     #[must_use]
     pub fn project_point(&self, p: &Point3<T>) -> Point3<T> {
@@ -261,6 +263,13 @@ impl<T: RealField> Perspective3<T> {
 
     // TODO: when we get specialization, specialize the Mul impl instead.
     /// Projects a vector. Faster than matrix multiplication.
+    ///
+    /// `x` and `y` components of the result is a projection as a 2D vector.
+    ///
+    /// Components `x` and `y` of the resulting vector will be NaN or Inf if `p.z` is zero.
+    ///
+    /// `z` component of the resulting vector always equals to `-self.m33` and does
+    /// not depend on the components of `p`.
     #[inline]
     #[must_use]
     pub fn project_vector<SB>(&self, p: &Vector<T, U3, SB>) -> Vector3<T>
@@ -271,7 +280,7 @@ impl<T: RealField> Perspective3<T> {
         Vector3::new(
             self.matrix[(0, 0)].clone() * p[0].clone() * inverse_denom.clone(),
             self.matrix[(1, 1)].clone() * p[1].clone() * inverse_denom,
-            self.matrix[(2, 2)].clone(),
+            -self.matrix[(2, 2)].clone(),
         )
     }