Merge pull request #1104 from dimforge/dev

Release v0.31.0

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

@@ -38,8 +38,12 @@ jobs:
       run: cargo build --features serde-serialize
     - name: Build nalgebra-lapack
       run: cd nalgebra-lapack; cargo build;
-    - name: Build nalgebra-sparse
+    - name: Build nalgebra-sparse --no-default-features
+      run: cd nalgebra-sparse; cargo build --no-default-features;
+    - name: Build nalgebra-sparse (default features)
       run: cd nalgebra-sparse; cargo build;
+    - name: Build nalgebra-sparse --all-features
+      run: cd nalgebra-sparse; cargo build --all-features;
   # Run this on it’s own job because it alone takes a lot of time.
   # So it’s best to let it run in parallel to the other jobs.
   build-nalgebra-all-features:
@@ -57,13 +61,13 @@ jobs:
     steps:
       - uses: actions/checkout@v2
       - name: test
-        run: cargo test --features arbitrary,rand,serde-serialize,abomonation-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests;
+        run: cargo test --features arbitrary,rand,serde-serialize,sparse,debug,io,compare,libm,proptest-support,slow-tests;
   test-nalgebra-glm:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v2
       - name: test nalgebra-glm
-        run: cargo test -p nalgebra-glm --features arbitrary,serde-serialize,abomonation-serialize;
+        run: cargo test -p nalgebra-glm --features arbitrary,serde-serialize;
   test-nalgebra-sparse:
     runs-on: ubuntu-latest
     steps:
@@ -71,10 +75,10 @@ jobs:
       - name: test nalgebra-sparse
         # Manifest-path is necessary because cargo otherwise won't correctly forward features
         # We increase number of proptest cases to hopefully catch more potential bugs
-        run: PROPTEST_CASES=10000 cargo test --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io
+        run: PROPTEST_CASES=10000 cargo test --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io,serde-serialize
       - name: test nalgebra-sparse (slow tests)
         # Unfortunately, the "slow-tests" take so much time that we need to run them with --release
-        run: PROPTEST_CASES=10000 cargo test --release --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io,slow-tests slow
+        run: PROPTEST_CASES=10000 cargo test --release --manifest-path=nalgebra-sparse/Cargo.toml --features compare,proptest-support,io,serde-serialize,slow-tests slow
   test-nalgebra-macros:
     runs-on: ubuntu-latest
     steps:
@@ -120,6 +124,8 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: Jimver/cuda-toolkit@v0.2.4
+        with:
+          cuda: '11.2.2'
       - name: Install nightly-2021-12-04
         uses: actions-rs/toolchain@v1
         with:
@@ -128,4 +134,6 @@ jobs:
       - uses: actions/checkout@v2
       - run: rustup target add nvptx64-nvidia-cuda
       - run: cargo build --no-default-features --features cuda
-      - run: cargo build --no-default-features --features cuda --target=nvptx64-nvidia-cuda
+      - run: cargo build --no-default-features --features cuda --target=nvptx64-nvidia-cuda
+        env: 
+          CUDA_ARCH: "350"

CHANGELOG.md

@@ -4,12 +4,42 @@ documented here.
 
 This project adheres to [Semantic Versioning](https://semver.org/).
 
+## [0.31.0] (30 Apr. 2022)
+
+### Breaking changes
+- Switch to `cust` 0.3 (for CUDA support).
+- Switch to `rkyv` 0.7
+- Remove support for serialization based on `abomonation`.
+- Remove support for conversions between `nalgebra` types and `glam` 0.13.
+
+### Modified
+- The aliases for `Const` types have been simplified to help `rust-analyzer`.
+
+### Added
+- Add `TryFrom` conversion between `UnitVector2/3/4` and `glam`’s `Vec2/3/4`.
+- `nalgebra-sparse`: added support for serialization of sparse matrices with `serde`.
+- `nalgebra-sparse`: add a CSC matrix constructor from unsorted (but valid) data.
+- `nalgebra-lapack`: add generalized eigenvalues/eigenvectors calculation + QZ decomposition.
+
+### Fixed
+- Improve stability of SVD.
+- Fix slerp for `UnitComplex`.
+
+## [0.30.1] (09 Jan. 2022)
+
+### Added
+- Add conversion from/to types of `glam` 0.19 and 0.20.
+
 ## [0.30.0] (02 Jan. 2022)
 
 ### Breaking changes
 - The `Dim` trait is now marked as unsafe.
 - The `Matrix::pow` and `Matrix::pow_mut` methods only allow positive integer exponents now. To compute negative
   exponents, the user is free to invert the matrix before calling `pow` with the exponent’s absolute value.
+- Remove the `Bounded` requirement from `RealField`. Replace it by methods returning `Option<Self>` so that they can
+  still be implemented by unbounded types (by returning `None`).
+- The `ComplexField` trait derives  from `FromPrimitive` again. We can actually keep this because all its methods
+  return `Option<Self>`, meaning that it could be implemented by any type.
 
 ### Modified
 - Use more concise debug impls for matrices and geometric transformation types.

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name    = "nalgebra"
-version = "0.30.0"
+version = "0.31.0"
 authors = [ "Sébastien Crozet <developer@crozet.re>" ]
 
 description = "General-purpose linear algebra library with transformations and statically-sized or dynamically-sized matrices."
@@ -32,17 +32,18 @@ compare = [ "matrixcompare-core" ]
 libm    = [ "simba/libm" ]
 libm-force = [ "simba/libm_force" ]
 macros = [ "nalgebra-macros" ]
-cuda   = [ "cust", "simba/cuda" ]
+cuda   = [ "cust_core", "simba/cuda" ]
 
 # Conversion
 convert-mint = [ "mint" ]
 convert-bytemuck = [ "bytemuck" ]
-convert-glam013 = [ "glam013" ]
 convert-glam014 = [ "glam014" ]
 convert-glam015 = [ "glam015" ]
 convert-glam016 = [ "glam016" ]
 convert-glam017 = [ "glam017" ]
 convert-glam018 = [ "glam018" ]
+convert-glam019 = [ "glam019" ]
+convert-glam020 = [ "glam020" ]
 
 # Serialization
 ## To use serde in a #[no-std] environment, enable the
@@ -51,9 +52,8 @@ convert-glam018 = [ "glam018" ]
 ## `serde-serialize`.
 serde-serialize-no-std = [ "serde", "num-complex/serde" ]
 serde-serialize        = [ "serde-serialize-no-std", "serde/std" ]
-abomonation-serialize  = [ "abomonation" ]
 rkyv-serialize-no-std  = [ "rkyv" ]
-rkyv-serialize         = [ "rkyv-serialize-no-std", "rkyv/std" ]
+rkyv-serialize         = [ "rkyv-serialize-no-std", "rkyv/std", "bytecheck" ]
 
 # Randomness
 ## To use rand in a #[no-std] environment, enable the
@@ -79,8 +79,8 @@ alga           = { version = "0.9", default-features = false, optional = true }
 rand_distr     = { version = "0.4", default-features = false, optional = true }
 matrixmultiply = { version = "0.3", optional = true }
 serde          = { version = "1.0", default-features = false, features = [ "derive" ], optional = true }
-abomonation    = { version = "0.7", optional = true }
+rkyv           = { version = "~0.7.1", optional = true }
-rkyv           = { version = "~0.6.4", default-features = false, features = ["const_generics"], optional = true }
+bytecheck      = { version = "~0.6.1", optional = true }
 mint           = { version = "0.5", optional = true }
 quickcheck     = { version = "1", optional = true }
 pest           = { version = "2", optional = true }
@@ -88,15 +88,14 @@ pest_derive    = { version = "2", optional = true }
 bytemuck       = { version = "1.5", optional = true }
 matrixcompare-core = { version = "0.1", optional = true }
 proptest       = { version = "1", optional = true, default-features = false, features = ["std"] }
-glam013        = { package = "glam", version = "0.13", optional = true }
 glam014        = { package = "glam", version = "0.14", optional = true }
 glam015        = { package = "glam", version = "0.15", optional = true }
 glam016        = { package = "glam", version = "0.16", optional = true }
 glam017        = { package = "glam", version = "0.17", optional = true }
 glam018        = { package = "glam", version = "0.18", optional = true }
-
+glam019        = { package = "glam", version = "0.19", optional = true }
-[target.'cfg(not(target_os = "cuda"))'.dependencies]
+glam020        = { package = "glam", version = "0.20", optional = true }
-cust           = { version = "0.2", optional = true }
+cust_core      = { version = "0.1", optional = true }
 
 
 [dev-dependencies]

README.md

@@ -42,6 +42,9 @@ And our gold sponsors:
 
 <p>
   <a href="https://fragcolor.com">
-    <img src="https://dimforge.com/img/fragcolor_logo1_color_black.svg" width="151px">
+    <img src="https://dimforge.com/img/fragcolor_logo2_color_black.svg" width="300px">
+  </a>
+  <a href="https://resolutiongames.com/">
+    <img src="https://dimforge.com/img/logo_resolution_games.png" width="300px" />
   </a>
 </p>

examples/cargo/Cargo.toml

@@ -4,7 +4,7 @@ version = "0.0.0"
 authors = [ "You" ]
 
 [dependencies]
-nalgebra = "0.30.0"
+nalgebra = "0.31.0"
 
 [[bin]]
 name = "example"

nalgebra-glm/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "nalgebra-glm"
-version = "0.16.0"
+version = "0.17.0"
 authors = ["sebcrozet <developer@crozet.re>"]
 
 description = "A computer-graphics oriented API for nalgebra, inspired by the C++ GLM library."
@@ -21,13 +21,11 @@ default               = [ "std" ]
 std                   = [ "nalgebra/std", "simba/std" ]
 arbitrary             = [ "nalgebra/arbitrary" ]
 serde-serialize       = [ "nalgebra/serde-serialize-no-std" ]
-abomonation-serialize = [ "nalgebra/abomonation-serialize" ]
 cuda                  = [ "nalgebra/cuda" ]
 
 # Conversion
 convert-mint = [ "nalgebra/mint" ]
 convert-bytemuck = [ "nalgebra/bytemuck" ]
-convert-glam013 = [ "nalgebra/glam013" ]
 convert-glam014 = [ "nalgebra/glam014" ]
 convert-glam015 = [ "nalgebra/glam015" ]
 convert-glam016 = [ "nalgebra/glam016" ]
@@ -38,4 +36,4 @@ convert-glam018 = [ "nalgebra/glam018" ]
 num-traits = { version = "0.2", default-features = false }
 approx = { version = "0.5", default-features = false }
 simba = { version = "0.7", default-features = false }
-nalgebra = { path =  "..", version = "0.30", default-features = false }
+nalgebra = { path =  "..", version = "0.31", default-features = false }

nalgebra-lapack/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name    = "nalgebra-lapack"
-version = "0.21.0"
+version = "0.22.0"
 authors = [ "Sébastien Crozet <developer@crozet.re>", "Andrew Straw <strawman@astraw.com>" ]
 
 description   = "Matrix decompositions using nalgebra matrices and Lapack bindings."
@@ -29,7 +29,7 @@ accelerate = ["lapack-src/accelerate"]
 intel-mkl  = ["lapack-src/intel-mkl"]
 
 [dependencies]
-nalgebra      = { version = "0.30", path = ".." }
+nalgebra      = { version = "0.31", path = ".." }
 num-traits    = "0.2"
 num-complex   = { version = "0.4", default-features = false }
 simba         = "0.7"
@@ -39,7 +39,7 @@ lapack-src    = { version = "0.8", default-features = false }
 # clippy = "*"
 
 [dev-dependencies]
-nalgebra   = { version = "0.30", features = [ "arbitrary", "rand" ], path = ".." }
+nalgebra   = { version = "0.31", features = [ "arbitrary", "rand" ], path = ".." }
 proptest   = { version = "1", default-features = false, features = ["std"] }
 quickcheck = "1"
 approx     = "0.5"

nalgebra-lapack/src/generalized_eigenvalues.rs

@@ -0,0 +1,350 @@
+#[cfg(feature = "serde-serialize")]
+use serde::{Deserialize, Serialize};
+
+use num::Zero;
+use num_complex::Complex;
+
+use simba::scalar::RealField;
+
+use crate::ComplexHelper;
+use na::allocator::Allocator;
+use na::dimension::{Const, Dim};
+use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar};
+
+use lapack;
+
+/// Generalized eigenvalues and generalized eigenvectors (left and right) of a pair of N*N real square matrices.
+///
+/// Each generalized eigenvalue (lambda) satisfies determinant(A - lambda*B) = 0
+///
+/// The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
+/// of (A,B) satisfies
+///
+/// A * v(j) = lambda(j) * B * v(j).
+///
+/// The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
+/// of (A,B) satisfies
+///
+/// u(j)**H * A  = lambda(j) * u(j)**H * B .
+/// where u(j)**H is the conjugate-transpose of u(j).
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(serialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+         OVector<T, D>: Serialize,
+         OMatrix<T, D, D>: Serialize")
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(deserialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+         OVector<T, D>: Deserialize<'de>,
+         OMatrix<T, D, D>: Deserialize<'de>")
+    )
+)]
+#[derive(Clone, Debug)]
+pub struct GeneralizedEigen<T: Scalar, D: Dim>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+{
+    alphar: OVector<T, D>,
+    alphai: OVector<T, D>,
+    beta: OVector<T, D>,
+    vsl: OMatrix<T, D, D>,
+    vsr: OMatrix<T, D, D>,
+}
+
+impl<T: Scalar + Copy, D: Dim> Copy for GeneralizedEigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+    OMatrix<T, D, D>: Copy,
+    OVector<T, D>: Copy,
+{
+}
+
+impl<T: GeneralizedEigenScalar + RealField + Copy, D: Dim> GeneralizedEigen<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+{
+    /// Attempts to compute the generalized eigenvalues, and left and right associated eigenvectors
+    /// via the raw returns from LAPACK's dggev and sggev routines
+    ///
+    /// Panics if the method did not converge.
+    pub fn new(a: OMatrix<T, D, D>, b: OMatrix<T, D, D>) -> Self {
+        Self::try_new(a, b).expect("Calculation of generalized eigenvalues failed.")
+    }
+
+    /// Attempts to compute the generalized eigenvalues (and eigenvectors) via the raw returns from LAPACK's
+    /// dggev and sggev routines
+    ///
+    /// Returns `None` if the method did not converge.
+    pub fn try_new(mut a: OMatrix<T, D, D>, mut b: OMatrix<T, D, D>) -> Option<Self> {
+        assert!(
+            a.is_square() && b.is_square(),
+            "Unable to compute the generalized eigenvalues of non-square matrices."
+        );
+
+        assert!(
+            a.shape_generic() ==  b.shape_generic(),
+            "Unable to compute the generalized eigenvalues of two square matrices of different dimensions."
+        );
+
+        let (nrows, ncols) = a.shape_generic();
+        let n = nrows.value();
+
+        let mut info = 0;
+
+        let mut alphar = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut alphai = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut beta = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut vsl = Matrix::zeros_generic(nrows, ncols);
+        let mut vsr = Matrix::zeros_generic(nrows, ncols);
+
+        let lwork = T::xggev_work_size(
+            b'V',
+            b'V',
+            n as i32,
+            a.as_mut_slice(),
+            n as i32,
+            b.as_mut_slice(),
+            n as i32,
+            alphar.as_mut_slice(),
+            alphai.as_mut_slice(),
+            beta.as_mut_slice(),
+            vsl.as_mut_slice(),
+            n as i32,
+            vsr.as_mut_slice(),
+            n as i32,
+            &mut info,
+        );
+        lapack_check!(info);
+
+        let mut work = vec![T::zero(); lwork as usize];
+
+        T::xggev(
+            b'V',
+            b'V',
+            n as i32,
+            a.as_mut_slice(),
+            n as i32,
+            b.as_mut_slice(),
+            n as i32,
+            alphar.as_mut_slice(),
+            alphai.as_mut_slice(),
+            beta.as_mut_slice(),
+            vsl.as_mut_slice(),
+            n as i32,
+            vsr.as_mut_slice(),
+            n as i32,
+            &mut work,
+            lwork,
+            &mut info,
+        );
+        lapack_check!(info);
+
+        Some(GeneralizedEigen {
+            alphar,
+            alphai,
+            beta,
+            vsl,
+            vsr,
+        })
+    }
+
+    /// Calculates the generalized eigenvectors (left and right) associated with the generalized eigenvalues
+    ///
+    /// Outputs two matrices.
+    /// The first output matrix contains the left eigenvectors of the generalized eigenvalues
+    /// as columns.
+    /// The second matrix contains the right eigenvectors of the generalized eigenvalues
+    /// as columns.
+    pub fn eigenvectors(&self) -> (OMatrix<Complex<T>, D, D>, OMatrix<Complex<T>, D, D>)
+    where
+        DefaultAllocator:
+            Allocator<Complex<T>, D, D> + Allocator<Complex<T>, D> + Allocator<(Complex<T>, T), D>,
+    {
+        /*
+         How the eigenvectors are built up:
+
+         Since the input entries are all real, the generalized eigenvalues if complex come in pairs
+         as a consequence of the [complex conjugate root thorem](https://en.wikipedia.org/wiki/Complex_conjugate_root_theorem)
+         The Lapack routine output reflects this by expecting the user to unpack the real and complex eigenvalues associated
+         eigenvectors from the real matrix output via the following procedure
+
+         (Note: VL stands for the lapack real matrix output containing the left eigenvectors as columns,
+                VR stands for the lapack real matrix output containing the right eigenvectors as columns)
+
+         If the j-th and (j+1)-th eigenvalues form a complex conjugate pair,
+         then
+
+          u(j) = VL(:,j)+i*VL(:,j+1)
+          u(j+1) = VL(:,j)-i*VL(:,j+1)
+
+         and
+
+          u(j) = VR(:,j)+i*VR(:,j+1)
+          v(j+1) = VR(:,j)-i*VR(:,j+1).
+        */
+
+        let n = self.vsl.shape().0;
+
+        let mut l = self.vsl.map(|x| Complex::new(x, T::RealField::zero()));
+
+        let mut r = self.vsr.map(|x| Complex::new(x, T::RealField::zero()));
+
+        let eigenvalues = self.raw_eigenvalues();
+
+        let mut c = 0;
+
+        while c < n {
+            if eigenvalues[c].0.im.abs() != T::RealField::zero() && c + 1 < n {
+                // taking care of the left eigenvector matrix
+                l.column_mut(c).zip_apply(&self.vsl.column(c + 1), |r, i| {
+                    *r = Complex::new(r.re.clone(), i.clone());
+                });
+                l.column_mut(c + 1).zip_apply(&self.vsl.column(c), |i, r| {
+                    *i = Complex::new(r.clone(), -i.re.clone());
+                });
+
+                // taking care of the right eigenvector matrix
+                r.column_mut(c).zip_apply(&self.vsr.column(c + 1), |r, i| {
+                    *r = Complex::new(r.re.clone(), i.clone());
+                });
+                r.column_mut(c + 1).zip_apply(&self.vsr.column(c), |i, r| {
+                    *i = Complex::new(r.clone(), -i.re.clone());
+                });
+
+                c += 2;
+            } else {
+                c += 1;
+            }
+        }
+
+        (l, r)
+    }
+
+    /// Outputs the unprocessed (almost) version of  generalized eigenvalues ((alphar, alphai), beta)
+    /// straight from LAPACK
+    #[must_use]
+    pub fn raw_eigenvalues(&self) -> OVector<(Complex<T>, T), D>
+    where
+        DefaultAllocator: Allocator<(Complex<T>, T), D>,
+    {
+        let mut out = Matrix::from_element_generic(
+            self.vsl.shape_generic().0,
+            Const::<1>,
+            (Complex::zero(), T::RealField::zero()),
+        );
+
+        for i in 0..out.len() {
+            out[i] = (Complex::new(self.alphar[i], self.alphai[i]), self.beta[i])
+        }
+
+        out
+    }
+}
+
+/*
+ *
+ * Lapack functions dispatch.
+ *
+ */
+/// Trait implemented by scalars for which Lapack implements the RealField GeneralizedEigen decomposition.
+pub trait GeneralizedEigenScalar: Scalar {
+    #[allow(missing_docs)]
+    fn xggev(
+        jobvsl: u8,
+        jobvsr: u8,
+        n: i32,
+        a: &mut [Self],
+        lda: i32,
+        b: &mut [Self],
+        ldb: i32,
+        alphar: &mut [Self],
+        alphai: &mut [Self],
+        beta: &mut [Self],
+        vsl: &mut [Self],
+        ldvsl: i32,
+        vsr: &mut [Self],
+        ldvsr: i32,
+        work: &mut [Self],
+        lwork: i32,
+        info: &mut i32,
+    );
+
+    #[allow(missing_docs)]
+    fn xggev_work_size(
+        jobvsl: u8,
+        jobvsr: u8,
+        n: i32,
+        a: &mut [Self],
+        lda: i32,
+        b: &mut [Self],
+        ldb: i32,
+        alphar: &mut [Self],
+        alphai: &mut [Self],
+        beta: &mut [Self],
+        vsl: &mut [Self],
+        ldvsl: i32,
+        vsr: &mut [Self],
+        ldvsr: i32,
+        info: &mut i32,
+    ) -> i32;
+}
+
+macro_rules! generalized_eigen_scalar_impl (
+    ($N: ty, $xggev: path) => (
+        impl GeneralizedEigenScalar for $N {
+            #[inline]
+            fn xggev(jobvsl:  u8,
+                     jobvsr:  u8,
+                     n:      i32,
+                     a:      &mut [$N],
+                     lda:    i32,
+                     b:      &mut [$N],
+                     ldb:    i32,
+                     alphar: &mut [$N],
+                     alphai: &mut [$N],
+                     beta  : &mut [$N],
+                     vsl:    &mut [$N],
+                     ldvsl:  i32,
+                     vsr:    &mut [$N],
+                     ldvsr:  i32,
+                     work:   &mut [$N],
+                     lwork:  i32,
+                     info:   &mut i32) {
+                unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, info); }
+            }
+
+
+            #[inline]
+            fn xggev_work_size(jobvsl:  u8,
+                               jobvsr:  u8,
+                               n:      i32,
+                               a:      &mut [$N],
+                               lda:    i32,
+                               b:      &mut [$N],
+                               ldb:    i32,
+                               alphar: &mut [$N],
+                               alphai: &mut [$N],
+                               beta  : &mut [$N],
+                               vsl:    &mut [$N],
+                               ldvsl:  i32,
+                               vsr:    &mut [$N],
+                               ldvsr:  i32,
+                               info:   &mut i32)
+                               -> i32 {
+                let mut work = [ Zero::zero() ];
+                let lwork    = -1 as i32;
+
+                unsafe { $xggev(jobvsl, jobvsr, n, a, lda, b, ldb, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, info); }
+                ComplexHelper::real_part(work[0]) as i32
+            }
+        }
+    )
+);
+
+generalized_eigen_scalar_impl!(f32, lapack::sggev);
+generalized_eigen_scalar_impl!(f64, lapack::dggev);

nalgebra-lapack/src/lib.rs

@@ -83,9 +83,11 @@ mod lapack_check;
 
 mod cholesky;
 mod eigen;
+mod generalized_eigenvalues;
 mod hessenberg;
 mod lu;
 mod qr;
+mod qz;
 mod schur;
 mod svd;
 mod symmetric_eigen;
@@ -94,9 +96,11 @@ use num_complex::Complex;
 
 pub use self::cholesky::{Cholesky, CholeskyScalar};
 pub use self::eigen::Eigen;
+pub use self::generalized_eigenvalues::GeneralizedEigen;
 pub use self::hessenberg::Hessenberg;
 pub use self::lu::{LUScalar, LU};
 pub use self::qr::QR;
+pub use self::qz::QZ;
 pub use self::schur::Schur;
 pub use self::svd::SVD;
 pub use self::symmetric_eigen::SymmetricEigen;

nalgebra-lapack/src/qz.rs

@@ -0,0 +1,321 @@
+#[cfg(feature = "serde-serialize")]
+use serde::{Deserialize, Serialize};
+
+use num::Zero;
+use num_complex::Complex;
+
+use simba::scalar::RealField;
+
+use crate::ComplexHelper;
+use na::allocator::Allocator;
+use na::dimension::{Const, Dim};
+use na::{DefaultAllocator, Matrix, OMatrix, OVector, Scalar};
+
+use lapack;
+
+/// QZ decomposition of a pair of N*N square matrices.
+///
+/// Retrieves the left and right matrices of Schur Vectors (VSL and VSR)
+/// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the
+/// decomposed input matrix `a`  equals `VSL * S * VSL.transpose()` and
+/// decomposed input matrix `b`  equals `VSL * T * VSL.transpose()`.
+#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(serialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+         OVector<T, D>: Serialize,
+         OMatrix<T, D, D>: Serialize")
+    )
+)]
+#[cfg_attr(
+    feature = "serde-serialize",
+    serde(
+        bound(deserialize = "DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+         OVector<T, D>: Deserialize<'de>,
+         OMatrix<T, D, D>: Deserialize<'de>")
+    )
+)]
+#[derive(Clone, Debug)]
+pub struct QZ<T: Scalar, D: Dim>
+where
+    DefaultAllocator: Allocator<T, D> + Allocator<T, D, D>,
+{
+    alphar: OVector<T, D>,
+    alphai: OVector<T, D>,
+    beta: OVector<T, D>,
+    vsl: OMatrix<T, D, D>,
+    s: OMatrix<T, D, D>,
+    vsr: OMatrix<T, D, D>,
+    t: OMatrix<T, D, D>,
+}
+
+impl<T: Scalar + Copy, D: Dim> Copy for QZ<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+    OMatrix<T, D, D>: Copy,
+    OVector<T, D>: Copy,
+{
+}
+
+impl<T: QZScalar + RealField, D: Dim> QZ<T, D>
+where
+    DefaultAllocator: Allocator<T, D, D> + Allocator<T, D>,
+{
+    /// Attempts to compute the QZ decomposition of input real square matrices `a` and `b`.
+    ///
+    /// i.e retrieves the left and right matrices of Schur Vectors (VSL and VSR)
+    /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the
+    /// decomposed matrix `a`  equals `VSL * S * VSL.transpose()` and
+    /// decomposed matrix `b`  equals `VSL * T * VSL.transpose()`.
+    ///
+    /// Panics if the method did not converge.
+    pub fn new(a: OMatrix<T, D, D>, b: OMatrix<T, D, D>) -> Self {
+        Self::try_new(a, b).expect("QZ decomposition: convergence failed.")
+    }
+
+    /// Computes the decomposition of input matrices `a` and `b` into a pair of matrices of Schur vectors
+    /// , a quasi-upper triangular matrix and an upper-triangular matrix .
+    ///
+    /// Returns `None` if the method did not converge.
+    pub fn try_new(mut a: OMatrix<T, D, D>, mut b: OMatrix<T, D, D>) -> Option<Self> {
+        assert!(
+            a.is_square() && b.is_square(),
+            "Unable to compute the qz decomposition of non-square matrices."
+        );
+
+        assert!(
+            a.shape_generic() ==  b.shape_generic(),
+            "Unable to compute the qz decomposition of two square matrices of different dimensions."
+        );
+
+        let (nrows, ncols) = a.shape_generic();
+        let n = nrows.value();
+
+        let mut info = 0;
+
+        let mut alphar = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut alphai = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut beta = Matrix::zeros_generic(nrows, Const::<1>);
+        let mut vsl = Matrix::zeros_generic(nrows, ncols);
+        let mut vsr = Matrix::zeros_generic(nrows, ncols);
+        // Placeholders:
+        let mut bwork = [0i32];
+        let mut unused = 0;
+
+        let lwork = T::xgges_work_size(
+            b'V',
+            b'V',
+            b'N',
+            n as i32,
+            a.as_mut_slice(),
+            n as i32,
+            b.as_mut_slice(),
+            n as i32,
+            &mut unused,
+            alphar.as_mut_slice(),
+            alphai.as_mut_slice(),
+            beta.as_mut_slice(),
+            vsl.as_mut_slice(),
+            n as i32,
+            vsr.as_mut_slice(),
+            n as i32,
+            &mut bwork,
+            &mut info,
+        );
+        lapack_check!(info);
+
+        let mut work = vec![T::zero(); lwork as usize];
+
+        T::xgges(
+            b'V',
+            b'V',
+            b'N',
+            n as i32,
+            a.as_mut_slice(),
+            n as i32,
+            b.as_mut_slice(),
+            n as i32,
+            &mut unused,
+            alphar.as_mut_slice(),
+            alphai.as_mut_slice(),
+            beta.as_mut_slice(),
+            vsl.as_mut_slice(),
+            n as i32,
+            vsr.as_mut_slice(),
+            n as i32,
+            &mut work,
+            lwork,
+            &mut bwork,
+            &mut info,
+        );
+        lapack_check!(info);
+
+        Some(QZ {
+            alphar,
+            alphai,
+            beta,
+            vsl,
+            s: a,
+            vsr,
+            t: b,
+        })
+    }
+
+    /// Retrieves the left and right matrices of Schur Vectors (VSL and VSR)
+    /// the upper-quasitriangular matrix `S` and upper triangular matrix `T` such that the
+    /// decomposed input matrix `a`  equals `VSL * S * VSL.transpose()` and
+    /// decomposed input matrix `b`  equals `VSL * T * VSL.transpose()`.
+    pub fn unpack(
+        self,
+    ) -> (
+        OMatrix<T, D, D>,
+        OMatrix<T, D, D>,
+        OMatrix<T, D, D>,
+        OMatrix<T, D, D>,
+    ) {
+        (self.vsl, self.s, self.t, self.vsr)
+    }
+
+    /// outputs the unprocessed (almost) version of  generalized eigenvalues ((alphar, alpai), beta)
+    /// straight from LAPACK
+    #[must_use]
+    pub fn raw_eigenvalues(&self) -> OVector<(Complex<T>, T), D>
+    where
+        DefaultAllocator: Allocator<(Complex<T>, T), D>,
+    {
+        let mut out = Matrix::from_element_generic(
+            self.vsl.shape_generic().0,
+            Const::<1>,
+            (Complex::zero(), T::RealField::zero()),
+        );
+
+        for i in 0..out.len() {
+            out[i] = (
+                Complex::new(self.alphar[i].clone(), self.alphai[i].clone()),
+                self.beta[i].clone(),
+            )
+        }
+
+        out
+    }
+}
+
+/*
+ *
+ * Lapack functions dispatch.
+ *
+ */
+/// Trait implemented by scalars for which Lapack implements the RealField QZ decomposition.
+pub trait QZScalar: Scalar {
+    #[allow(missing_docs)]
+    fn xgges(
+        jobvsl: u8,
+        jobvsr: u8,
+        sort: u8,
+        // select: ???
+        n: i32,
+        a: &mut [Self],
+        lda: i32,
+        b: &mut [Self],
+        ldb: i32,
+        sdim: &mut i32,
+        alphar: &mut [Self],
+        alphai: &mut [Self],
+        beta: &mut [Self],
+        vsl: &mut [Self],
+        ldvsl: i32,
+        vsr: &mut [Self],
+        ldvsr: i32,
+        work: &mut [Self],
+        lwork: i32,
+        bwork: &mut [i32],
+        info: &mut i32,
+    );
+
+    #[allow(missing_docs)]
+    fn xgges_work_size(
+        jobvsl: u8,
+        jobvsr: u8,
+        sort: u8,
+        // select: ???
+        n: i32,
+        a: &mut [Self],
+        lda: i32,
+        b: &mut [Self],
+        ldb: i32,
+        sdim: &mut i32,
+        alphar: &mut [Self],
+        alphai: &mut [Self],
+        beta: &mut [Self],
+        vsl: &mut [Self],
+        ldvsl: i32,
+        vsr: &mut [Self],
+        ldvsr: i32,
+        bwork: &mut [i32],
+        info: &mut i32,
+    ) -> i32;
+}
+
+macro_rules! qz_scalar_impl (
+    ($N: ty, $xgges: path) => (
+        impl QZScalar for $N {
+            #[inline]
+            fn xgges(jobvsl:  u8,
+                     jobvsr:  u8,
+                     sort:   u8,
+                     // select: ???
+                     n:      i32,
+                     a:      &mut [$N],
+                     lda:    i32,
+                     b:      &mut [$N],
+                     ldb:    i32,
+                     sdim:   &mut i32,
+                     alphar: &mut [$N],
+                     alphai: &mut [$N],
+                     beta  : &mut [$N],
+                     vsl:    &mut [$N],
+                     ldvsl:  i32,
+                     vsr:    &mut [$N],
+                     ldvsr:  i32,
+                     work:   &mut [$N],
+                     lwork:  i32,
+                     bwork:  &mut [i32],
+                     info:   &mut i32) {
+                unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, bwork, info); }
+            }
+
+
+            #[inline]
+            fn xgges_work_size(jobvsl:  u8,
+                               jobvsr:  u8,
+                               sort:   u8,
+                               // select: ???
+                               n:      i32,
+                               a:      &mut [$N],
+                               lda:    i32,
+                               b:      &mut [$N],
+                               ldb:    i32,
+                               sdim:   &mut i32,
+                               alphar: &mut [$N],
+                               alphai: &mut [$N],
+                               beta  : &mut [$N],
+                               vsl:    &mut [$N],
+                               ldvsl:  i32,
+                               vsr:    &mut [$N],
+                               ldvsr:  i32,
+                               bwork:  &mut [i32],
+                               info:   &mut i32)
+                               -> i32 {
+                let mut work = [ Zero::zero() ];
+                let lwork    = -1 as i32;
+
+                unsafe { $xgges(jobvsl, jobvsr, sort, None, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, &mut work, lwork, bwork, info); }
+                ComplexHelper::real_part(work[0]) as i32
+            }
+        }
+    )
+);
+
+qz_scalar_impl!(f32, lapack::sgges);
+qz_scalar_impl!(f64, lapack::dgges);

nalgebra-lapack/tests/linalg/generalized_eigenvalues.rs

@@ -0,0 +1,72 @@
+use na::dimension::Const;
+use na::{DMatrix, OMatrix};
+use nl::GeneralizedEigen;
+use num_complex::Complex;
+use simba::scalar::ComplexField;
+
+use crate::proptest::*;
+use proptest::{prop_assert, prop_compose, proptest};
+
+prop_compose! {
+    fn f64_dynamic_dim_squares()
+        (n in PROPTEST_MATRIX_DIM)
+        (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){
+    (a,b)
+}}
+
+proptest! {
+    #[test]
+    fn ge((a,b) in f64_dynamic_dim_squares()){
+
+        let a_c = a.clone().map(|x| Complex::new(x, 0.0));
+        let b_c = b.clone().map(|x| Complex::new(x, 0.0));
+        let n = a.shape_generic().0;
+
+        let ge = GeneralizedEigen::new(a.clone(), b.clone());
+        let (vsl,vsr) = ge.clone().eigenvectors();
+
+
+        for (i,(alpha,beta)) in ge.raw_eigenvalues().iter().enumerate() {
+            let l_a = a_c.clone() * Complex::new(*beta, 0.0);
+            let l_b = b_c.clone() * *alpha;
+
+            prop_assert!(
+                relative_eq!(
+                    ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
+                    OMatrix::zeros_generic(n, Const::<1>),
+                    epsilon = 1.0e-5));
+
+            prop_assert!(
+                relative_eq!(
+                    (vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
+                    OMatrix::zeros_generic(Const::<1>, n),
+                    epsilon = 1.0e-5))
+        };
+    }
+
+    #[test]
+    fn ge_static(a in matrix4(), b in matrix4()) {
+
+        let ge = GeneralizedEigen::new(a.clone(), b.clone());
+        let a_c =a.clone().map(|x| Complex::new(x, 0.0));
+        let b_c = b.clone().map(|x| Complex::new(x, 0.0));
+        let (vsl,vsr) = ge.eigenvectors();
+        let eigenvalues = ge.raw_eigenvalues();
+
+        for (i,(alpha,beta)) in eigenvalues.iter().enumerate() {
+            let l_a = a_c.clone() * Complex::new(*beta, 0.0);
+            let l_b = b_c.clone() * *alpha;
+
+            prop_assert!(
+                relative_eq!(
+                    ((&l_a - &l_b)*vsr.column(i)).map(|x| x.modulus()),
+                    OMatrix::zeros_generic(Const::<4>, Const::<1>),
+                    epsilon = 1.0e-5));
+            prop_assert!(
+                relative_eq!((vsl.column(i).adjoint()*(&l_a - &l_b)).map(|x| x.modulus()),
+                             OMatrix::zeros_generic(Const::<1>, Const::<4>),
+                             epsilon = 1.0e-5))
+        }
+    }
+
+}

nalgebra-lapack/tests/linalg/mod.rs

@@ -1,6 +1,8 @@
 mod cholesky;
+mod generalized_eigenvalues;
 mod lu;
 mod qr;
+mod qz;
 mod real_eigensystem;
 mod schur;
 mod svd;

nalgebra-lapack/tests/linalg/qz.rs

@@ -0,0 +1,34 @@
+use na::DMatrix;
+use nl::QZ;
+
+use crate::proptest::*;
+use proptest::{prop_assert, prop_compose, proptest};
+
+prop_compose! {
+    fn f64_dynamic_dim_squares()
+        (n in PROPTEST_MATRIX_DIM)
+        (a in matrix(PROPTEST_F64,n,n), b in matrix(PROPTEST_F64,n,n)) -> (DMatrix<f64>, DMatrix<f64>){
+    (a,b)
+}}
+
+proptest! {
+    #[test]
+    fn qz((a,b) in f64_dynamic_dim_squares()) {
+
+        let qz = QZ::new(a.clone(), b.clone());
+        let (vsl,s,t,vsr) = qz.clone().unpack();
+
+        prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7));
+
+    }
+
+    #[test]
+    fn qz_static(a in matrix4(), b in matrix4()) {
+        let qz = QZ::new(a.clone(), b.clone());
+        let (vsl,s,t,vsr) = qz.unpack();
+
+        prop_assert!(relative_eq!(&vsl * s * vsr.transpose(), a, epsilon = 1.0e-7));
+        prop_assert!(relative_eq!(vsl * t * vsr.transpose(), b, epsilon = 1.0e-7));
+    }
+}

nalgebra-macros/Cargo.toml

@@ -21,5 +21,5 @@ quote = "1.0"
 proc-macro2 = "1.0"
 
 [dev-dependencies]
-nalgebra = { version = "0.30.0", path = ".." }
+nalgebra = { version = "0.31.0", path = ".." }
 trybuild = "1.0.42"

nalgebra-macros/src/lib.rs

@@ -111,7 +111,7 @@ impl Parse for Matrix {
 
 /// Construct a fixed-size matrix directly from data.
 ///
-/// **Note: Requires the `macro` feature to be enabled (enabled by default)**.
+/// **Note: Requires the `macros` feature to be enabled (enabled by default)**.
 ///
 /// This macro facilitates easy construction of matrices when the entries of the matrix are known
 /// (either as constants or expressions). This macro produces an instance of `SMatrix`. This means
@@ -125,7 +125,6 @@ impl Parse for Matrix {
 /// (`;`) designates that a new row begins.
 ///
 /// # Examples
-///
 /// ```
 /// use nalgebra::matrix;
 ///
@@ -164,12 +163,13 @@ pub fn matrix(stream: TokenStream) -> TokenStream {
 
 /// Construct a dynamic matrix directly from data.
 ///
-/// **Note: Requires the `macro` feature to be enabled (enabled by default)**.
+/// **Note: Requires the `macros` feature to be enabled (enabled by default)**.
 ///
 /// The syntax is exactly the same as for [`matrix!`], but instead of producing instances of
 /// `SMatrix`, it produces instances of `DMatrix`. At the moment it is not usable
 /// in `const fn` contexts.
 ///
+/// # Example
 /// ```
 /// use nalgebra::dmatrix;
 ///
@@ -233,7 +233,7 @@ impl Parse for Vector {
 
 /// Construct a fixed-size column vector directly from data.
 ///
-/// **Note: Requires the `macro` feature to be enabled (enabled by default)**.
+/// **Note: Requires the `macros` feature to be enabled (enabled by default)**.
 ///
 /// Similarly to [`matrix!`], this macro facilitates easy construction of fixed-size vectors.
 /// However, whereas the [`matrix!`] macro expects each row to be separated by a semi-colon,
@@ -243,8 +243,7 @@ impl Parse for Vector {
 /// `vector!` is intended to be the most readable and performant way of constructing small,
 /// fixed-size vectors, and it is usable in `const fn` contexts.
 ///
-/// ## Examples
+/// # Example
-///
 /// ```
 /// use nalgebra::vector;
 ///
@@ -265,12 +264,13 @@ pub fn vector(stream: TokenStream) -> TokenStream {
 
 /// Construct a dynamic column vector directly from data.
 ///
-/// **Note: Requires the `macro` feature to be enabled (enabled by default)**.
+/// **Note: Requires the `macros` feature to be enabled (enabled by default)**.
 ///
 /// The syntax is exactly the same as for [`vector!`], but instead of producing instances of
 /// `SVector`, it produces instances of `DVector`. At the moment it is not usable
 /// in `const fn` contexts.
 ///
+/// # Example
 /// ```
 /// use nalgebra::dvector;
 ///
@@ -294,15 +294,14 @@ pub fn dvector(stream: TokenStream) -> TokenStream {
 
 /// Construct a fixed-size point directly from data.
 ///
-/// **Note: Requires the `macro` feature to be enabled (enabled by default)**.
+/// **Note: Requires the `macros` feature to be enabled (enabled by default)**.
 ///
 /// Similarly to [`vector!`], this macro facilitates easy construction of points.
 ///
 /// `point!` is intended to be the most readable and performant way of constructing small,
 /// points, and it is usable in `const fn` contexts.
 ///
-/// ## Examples
+/// # Example
-///
 /// ```
 /// use nalgebra::point;
 ///

nalgebra-sparse/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "nalgebra-sparse"
-version = "0.6.0"
+version = "0.7.0"
 authors = [ "Andreas Longva", "Sébastien Crozet <developer@crozet.re>" ]
 edition = "2018"
 description = "Sparse matrix computation based on nalgebra."
@@ -15,6 +15,7 @@ license = "Apache-2.0"
 [features]
 proptest-support = ["proptest", "nalgebra/proptest-support"]
 compare = [ "matrixcompare-core" ]
+serde-serialize = [ "serde/std" ]
 
 # Enable matrix market I/O
 io      = [ "pest", "pest_derive" ]
@@ -23,18 +24,20 @@ io      = [ "pest", "pest_derive" ]
 slow-tests = []
 
 [dependencies]
-nalgebra = { version="0.30", path = "../" }
+nalgebra = { version="0.31", path = "../" }
 num-traits = { version = "0.2", default-features = false }
 proptest = { version = "1.0", optional = true }
 matrixcompare-core = { version = "0.1.0", optional = true }
 pest           = { version = "2", optional = true }
 pest_derive    = { version = "2", optional = true }
+serde = { version = "1.0", default-features = false, features = [ "derive" ], optional = true }
 
 [dev-dependencies]
 itertools = "0.10"
 matrixcompare = { version = "0.3.0", features = [ "proptest-support" ] }
-nalgebra = { version="0.30", path = "../", features = ["compare"] }
+nalgebra = { version="0.31", path = "../", features = ["compare"] }
+serde_json = "1.0"
 
 [package.metadata.docs.rs]
 # Enable certain features when building docs for docs.rs
-features = [ "proptest-support", "compare" ]
+features = [ "proptest-support", "compare" ]

nalgebra-sparse/src/convert/serial.rs

@@ -14,6 +14,7 @@ use crate::coo::CooMatrix;
 use crate::cs;
 use crate::csc::CscMatrix;
 use crate::csr::CsrMatrix;
+use crate::utils::{apply_permutation, compute_sort_permutation};
 
 /// Converts a dense matrix to [`CooMatrix`].
 pub fn convert_dense_coo<T, R, C, S>(dense: &Matrix<T, R, C, S>) -> CooMatrix<T>
@@ -376,29 +377,12 @@ fn sort_lane<T: Clone>(
     assert_eq!(values.len(), workspace.len());
 
     let permutation = workspace;
-    // Set permutation to identity
+    compute_sort_permutation(permutation, minor_idx);
-    for (i, p) in permutation.iter_mut().enumerate() {
-        *p = i;
-    }
-
-    // Compute permutation needed to bring minor indices into sorted order
-    // Note: Using sort_unstable here avoids internal allocations, which is crucial since
-    // each lane might have a small number of elements
-    permutation.sort_unstable_by_key(|idx| minor_idx[*idx]);
 
     apply_permutation(minor_idx_result, minor_idx, permutation);
     apply_permutation(values_result, values, permutation);
 }
 
-// TODO: Move this into `utils` or something?
-fn apply_permutation<T: Clone>(out_slice: &mut [T], in_slice: &[T], permutation: &[usize]) {
-    assert_eq!(out_slice.len(), in_slice.len());
-    assert_eq!(out_slice.len(), permutation.len());
-    for (out_element, old_pos) in out_slice.iter_mut().zip(permutation) {
-        *out_element = in_slice[*old_pos].clone();
-    }
-}
-
 /// Given *sorted* indices and corresponding scalar values, combines duplicates with the given
 /// associative combiner and calls the provided produce methods with combined indices and values.
 fn combine_duplicates<T: Clone>(

nalgebra-sparse/src/coo.rs

@@ -1,5 +1,8 @@
 //! An implementation of the COO sparse matrix format.
 
+#[cfg(feature = "serde-serialize")]
+mod coo_serde;
+
 use crate::SparseFormatError;
 
 /// A COO representation of a sparse matrix.

nalgebra-sparse/src/coo/coo_serde.rs

@@ -0,0 +1,65 @@
+use crate::coo::CooMatrix;
+use serde::{de, Deserialize, Deserializer, Serialize, Serializer};
+
+/// This is an intermediate type for (de)serializing `CooMatrix`.
+///
+/// Deserialization requires using a `try_from_*` function for validation. We could have used
+/// the `remote = "Self"` trick (https://github.com/serde-rs/serde/issues/1220) which allows
+/// to directly serialize/deserialize the original fields and combine it with validation.
+/// However, this would lead to nested serialization of the `CsMatrix` and `SparsityPattern`
+/// types. Instead, we decided that we want a more human-readable serialization format using
+/// field names like `row_indices` and `col_indices`. The easiest way to achieve this is to
+/// introduce an intermediate type. It also allows the serialization format to stay constant
+/// even if the internal layout in `nalgebra` changes.
+///
+/// We want to avoid unnecessary copies when serializing (i.e. cloning slices into owned
+/// storage). Therefore, we use generic arguments to allow using slices during serialization and
+/// owned storage (i.e. `Vec`) during deserialization. Without a major update of serde, slices
+/// and `Vec`s should always (de)serialize identically.
+#[derive(Serialize, Deserialize)]
+struct CooMatrixSerializationData<Indices, Values> {
+    nrows: usize,
+    ncols: usize,
+    row_indices: Indices,
+    col_indices: Indices,
+    values: Values,
+}
+
+impl<T> Serialize for CooMatrix<T>
+where
+    T: Serialize + Clone,
+{
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        CooMatrixSerializationData::<&[usize], &[T]> {
+            nrows: self.nrows(),
+            ncols: self.ncols(),
+            row_indices: self.row_indices(),
+            col_indices: self.col_indices(),
+            values: self.values(),
+        }
+        .serialize(serializer)
+    }
+}
+
+impl<'de, T> Deserialize<'de> for CooMatrix<T>
+where
+    T: Deserialize<'de> + Clone,
+{
+    fn deserialize<D>(deserializer: D) -> Result<CooMatrix<T>, D::Error>
+    where
+        D: Deserializer<'de>,
+    {
+        let de = CooMatrixSerializationData::<Vec<usize>, Vec<T>>::deserialize(deserializer)?;
+        CooMatrix::try_from_triplets(
+            de.nrows,
+            de.ncols,
+            de.row_indices,
+            de.col_indices,
+            de.values,
+        )
+        .map_err(|e| de::Error::custom(e))
+    }
+}

nalgebra-sparse/src/cs.rs

@@ -6,7 +6,8 @@ use num_traits::One;
 use nalgebra::Scalar;
 
 use crate::pattern::SparsityPattern;
-use crate::{SparseEntry, SparseEntryMut};
+use crate::utils::{apply_permutation, compute_sort_permutation};
+use crate::{SparseEntry, SparseEntryMut, SparseFormatError, SparseFormatErrorKind};
 
 /// An abstract compressed matrix.
 ///
@@ -543,3 +544,151 @@ pub fn convert_counts_to_offsets(counts: &mut [usize]) {
         offset += count;
     }
 }
+
+/// Validates cs data, optionally sorts minor indices and values
+pub(crate) fn validate_and_optionally_sort_cs_data<T>(
+    major_dim: usize,
+    minor_dim: usize,
+    major_offsets: &[usize],
+    minor_indices: &mut [usize],
+    values: Option<&mut [T]>,
+    sort: bool,
+) -> Result<(), SparseFormatError>
+where
+    T: Scalar,
+{
+    let mut values_option = values;
+
+    if let Some(values) = values_option.as_mut() {
+        if minor_indices.len() != values.len() {
+            return Err(SparseFormatError::from_kind_and_msg(
+                SparseFormatErrorKind::InvalidStructure,
+                "Number of values and minor indices must be the same.",
+            ));
+        }
+    } else if sort {
+        unreachable!("Internal error: Sorting currently not supported if no values are present.");
+    }
+    if major_offsets.len() == 0 {
+        return Err(SparseFormatError::from_kind_and_msg(
+            SparseFormatErrorKind::InvalidStructure,
+            "Number of offsets should be greater than 0.",
+        ));
+    }
+    if major_offsets.len() != major_dim + 1 {
+        return Err(SparseFormatError::from_kind_and_msg(
+            SparseFormatErrorKind::InvalidStructure,
+            "Length of offset array is not equal to (major_dim + 1).",
+        ));
+    }
+
+    // Check that the first and last offsets conform to the specification
+    {
+        let first_offset_ok = *major_offsets.first().unwrap() == 0;
+        let last_offset_ok = *major_offsets.last().unwrap() == minor_indices.len();
+        if !first_offset_ok || !last_offset_ok {
+            return Err(SparseFormatError::from_kind_and_msg(
+                SparseFormatErrorKind::InvalidStructure,
+                "First or last offset is incompatible with format.",
+            ));
+        }
+    }
+
+    // Set up required buffers up front
+    let mut minor_idx_buffer: Vec<usize> = Vec::new();
+    let mut values_buffer: Vec<T> = Vec::new();
+    let mut minor_index_permutation: Vec<usize> = Vec::new();
+
+    // Test that each lane has strictly monotonically increasing minor indices, i.e.
+    // minor indices within a lane are sorted, unique. Sort minor indices within a lane if needed.
+    // In addition, each minor index must be in bounds with respect to the minor dimension.
+    {
+        for lane_idx in 0..major_dim {
+            let range_start = major_offsets[lane_idx];
+            let range_end = major_offsets[lane_idx + 1];
+
+            // Test that major offsets are monotonically increasing
+            if range_start > range_end {
+                return Err(SparseFormatError::from_kind_and_msg(
+                    SparseFormatErrorKind::InvalidStructure,
+                    "Offsets are not monotonically increasing.",
+                ));
+            }
+
+            let minor_idx_in_lane = minor_indices.get(range_start..range_end).ok_or(
+                SparseFormatError::from_kind_and_msg(
+                    SparseFormatErrorKind::IndexOutOfBounds,
+                    "A major offset is out of bounds.",
+                ),
+            )?;
+
+            // We test for in-bounds, uniqueness and monotonicity at the same time
+            // to ensure that we only visit each minor index once
+            let mut prev = None;
+            let mut monotonic = true;
+
+            for &minor_idx in minor_idx_in_lane {
+                if minor_idx >= minor_dim {
+                    return Err(SparseFormatError::from_kind_and_msg(
+                        SparseFormatErrorKind::IndexOutOfBounds,
+                        "A minor index is out of bounds.",
+                    ));
+                }
+
+                if let Some(prev) = prev {
+                    if prev >= minor_idx {
+                        if !sort {
+                            return Err(SparseFormatError::from_kind_and_msg(
+                                SparseFormatErrorKind::InvalidStructure,
+                                "Minor indices are not strictly monotonically increasing in each lane.",
+                            ));
+                        }
+                        monotonic = false;
+                    }
+                }
+                prev = Some(minor_idx);
+            }
+
+            // sort if indices are not monotonic and sorting is expected
+            if !monotonic && sort {
+                let range_size = range_end - range_start;
+                minor_index_permutation.resize(range_size, 0);
+                compute_sort_permutation(&mut minor_index_permutation, &minor_idx_in_lane);
+                minor_idx_buffer.clear();
+                minor_idx_buffer.extend_from_slice(&minor_idx_in_lane);
+                apply_permutation(
+                    &mut minor_indices[range_start..range_end],
+                    &minor_idx_buffer,
+                    &minor_index_permutation,
+                );
+
+                // check duplicates
+                prev = None;
+                for &minor_idx in &minor_indices[range_start..range_end] {
+                    if let Some(prev) = prev {
+                        if prev == minor_idx {
+                            return Err(SparseFormatError::from_kind_and_msg(
+                                SparseFormatErrorKind::DuplicateEntry,
+                                "Input data contains duplicate entries.",
+                            ));
+                        }
+                    }
+                    prev = Some(minor_idx);
+                }
+
+                // sort values if they exist
+                if let Some(values) = values_option.as_mut() {
+                    values_buffer.clear();
+                    values_buffer.extend_from_slice(&values[range_start..range_end]);
+                    apply_permutation(
+                        &mut values[range_start..range_end],
+                        &values_buffer,
+                        &minor_index_permutation,
+                    );
+                }
+            }
+        }
+    }
+
+    Ok(())
+}

nalgebra-sparse/src/csc.rs

@@ -3,6 +3,10 @@
 //! This is the module-level documentation. See [`CscMatrix`] for the main documentation of the
 //! CSC implementation.
 
+#[cfg(feature = "serde-serialize")]
+mod csc_serde;
+
+use crate::cs;
 use crate::cs::{CsLane, CsLaneIter, CsLaneIterMut, CsLaneMut, CsMatrix};
 use crate::csr::CsrMatrix;
 use crate::pattern::{SparsityPattern, SparsityPatternFormatError, SparsityPatternIter};
@@ -170,6 +174,50 @@ impl<T> CscMatrix<T> {
         Self::try_from_pattern_and_values(pattern, values)
     }
 
+    /// Try to construct a CSC matrix from raw CSC data with unsorted row indices.
+    ///
+    /// It is assumed that each column contains unique row indices that are in
+    /// bounds with respect to the number of rows in the matrix. If this is not the case,
+    /// an error is returned to indicate the failure.
+    ///
+    /// An error is returned if the data given does not conform to the CSC storage format
+    /// with the exception of having unsorted row indices and values.
+    /// See the documentation for [CscMatrix](struct.CscMatrix.html) for more information.
+    pub fn try_from_unsorted_csc_data(
+        num_rows: usize,
+        num_cols: usize,
+        col_offsets: Vec<usize>,
+        mut row_indices: Vec<usize>,
+        mut values: Vec<T>,
+    ) -> Result<Self, SparseFormatError>
+    where
+        T: Scalar,
+    {
+        let result = cs::validate_and_optionally_sort_cs_data(
+            num_cols,
+            num_rows,
+            &col_offsets,
+            &mut row_indices,
+            Some(&mut values),
+            true,
+        );
+
+        match result {
+            Ok(()) => {
+                let pattern = unsafe {
+                    SparsityPattern::from_offset_and_indices_unchecked(
+                        num_cols,
+                        num_rows,
+                        col_offsets,
+                        row_indices,
+                    )
+                };
+                Self::try_from_pattern_and_values(pattern, values)
+            }
+            Err(err) => Err(err),
+        }
+    }
+
     /// Try to construct a CSC matrix from a sparsity pattern and associated non-zero values.
     ///
     /// Returns an error if the number of values does not match the number of minor indices

nalgebra-sparse/src/csc/csc_serde.rs

@@ -0,0 +1,65 @@
+use crate::CscMatrix;
+use serde::{de, Deserialize, Deserializer, Serialize, Serializer};
+
+/// This is an intermediate type for (de)serializing `CscMatrix`.
+///
+/// Deserialization requires using a `try_from_*` function for validation. We could have used
+/// the `remote = "Self"` trick (https://github.com/serde-rs/serde/issues/1220) which allows
+/// to directly serialize/deserialize the original fields and combine it with validation.
+/// However, this would lead to nested serialization of the `CsMatrix` and `SparsityPattern`
+/// types. Instead, we decided that we want a more human-readable serialization format using
+/// field names like `col_offsets` and `row_indices`. The easiest way to achieve this is to
+/// introduce an intermediate type. It also allows the serialization format to stay constant
+/// even if the internal layout in `nalgebra` changes.
+///
+/// We want to avoid unnecessary copies when serializing (i.e. cloning slices into owned
+/// storage). Therefore, we use generic arguments to allow using slices during serialization and
+/// owned storage (i.e. `Vec`) during deserialization. Without a major update of serde, slices
+/// and `Vec`s should always (de)serialize identically.
+#[derive(Serialize, Deserialize)]
+struct CscMatrixSerializationData<Indices, Values> {
+    nrows: usize,
+    ncols: usize,
+    col_offsets: Indices,
+    row_indices: Indices,
+    values: Values,
+}
+
+impl<T> Serialize for CscMatrix<T>
+where
+    T: Serialize + Clone,
+{
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        CscMatrixSerializationData::<&[usize], &[T]> {
+            nrows: self.nrows(),
+            ncols: self.ncols(),
+            col_offsets: self.col_offsets(),
+            row_indices: self.row_indices(),
+            values: self.values(),
+        }
+        .serialize(serializer)
+    }
+}
+
+impl<'de, T> Deserialize<'de> for CscMatrix<T>
+where
+    T: Deserialize<'de> + Clone,
+{
+    fn deserialize<D>(deserializer: D) -> Result<CscMatrix<T>, D::Error>
+    where
+        D: Deserializer<'de>,
+    {
+        let de = CscMatrixSerializationData::<Vec<usize>, Vec<T>>::deserialize(deserializer)?;
+        CscMatrix::try_from_csc_data(
+            de.nrows,
+            de.ncols,
+            de.col_offsets,
+            de.row_indices,
+            de.values,
+        )
+        .map_err(|e| de::Error::custom(e))
+    }
+}

nalgebra-sparse/src/csr.rs

@@ -2,6 +2,11 @@
 //!
 //! This is the module-level documentation. See [`CsrMatrix`] for the main documentation of the
 //! CSC implementation.
+
+#[cfg(feature = "serde-serialize")]
+mod csr_serde;
+
+use crate::cs;
 use crate::cs::{CsLane, CsLaneIter, CsLaneIterMut, CsLaneMut, CsMatrix};
 use crate::csc::CscMatrix;
 use crate::pattern::{SparsityPattern, SparsityPatternFormatError, SparsityPatternIter};
@@ -10,7 +15,6 @@ use crate::{SparseEntry, SparseEntryMut, SparseFormatError, SparseFormatErrorKin
 use nalgebra::Scalar;
 use num_traits::One;
 
-use std::iter::FromIterator;
 use std::slice::{Iter, IterMut};
 
 /// A CSR representation of a sparse matrix.
@@ -184,62 +188,35 @@ impl<T> CsrMatrix<T> {
         num_rows: usize,
         num_cols: usize,
         row_offsets: Vec<usize>,
-        col_indices: Vec<usize>,
+        mut col_indices: Vec<usize>,
-        values: Vec<T>,
+        mut values: Vec<T>,
     ) -> Result<Self, SparseFormatError>
     where
         T: Scalar,
     {
-        use SparsityPatternFormatError::*;
+        let result = cs::validate_and_optionally_sort_cs_data(
-        let count = col_indices.len();
-        let mut p: Vec<usize> = (0..count).collect();
-
-        if col_indices.len() != values.len() {
-            return Err(SparseFormatError::from_kind_and_msg(
-                SparseFormatErrorKind::InvalidStructure,
-                "Number of values and column indices must be the same",
-            ));
-        }
-
-        if row_offsets.len() == 0 {
-            return Err(SparseFormatError::from_kind_and_msg(
-                SparseFormatErrorKind::InvalidStructure,
-                "Number of offsets should be greater than 0",
-            ));
-        }
-
-        for (index, &offset) in row_offsets[0..row_offsets.len() - 1].iter().enumerate() {
-            let next_offset = row_offsets[index + 1];
-            if next_offset > count {
-                return Err(SparseFormatError::from_kind_and_msg(
-                    SparseFormatErrorKind::InvalidStructure,
-                    "No row offset should be greater than the number of column indices",
-                ));
-            }
-            if offset > next_offset {
-                return Err(NonmonotonicOffsets).map_err(pattern_format_error_to_csr_error);
-            }
-            p[offset..next_offset].sort_by(|a, b| {
-                let x = &col_indices[*a];
-                let y = &col_indices[*b];
-                x.partial_cmp(y).unwrap()
-            });
-        }
-
-        // permute indices
-        let sorted_col_indices: Vec<usize> =
-            Vec::from_iter((p.iter().map(|i| &col_indices[*i])).cloned());
-
-        // permute values
-        let sorted_values: Vec<T> = Vec::from_iter((p.iter().map(|i| &values[*i])).cloned());
-
-        return Self::try_from_csr_data(
             num_rows,
             num_cols,
-            row_offsets,
+            &row_offsets,
-            sorted_col_indices,
+            &mut col_indices,
-            sorted_values,
+            Some(&mut values),
+            true,
         );
+
+        match result {
+            Ok(()) => {
+                let pattern = unsafe {
+                    SparsityPattern::from_offset_and_indices_unchecked(
+                        num_rows,
+                        num_cols,
+                        row_offsets,
+                        col_indices,
+                    )
+                };
+                Self::try_from_pattern_and_values(pattern, values)
+            }
+            Err(err) => Err(err),
+        }
     }
 
     /// Try to construct a CSR matrix from a sparsity pattern and associated non-zero values.

nalgebra-sparse/src/csr/csr_serde.rs

@@ -0,0 +1,65 @@
+use crate::CsrMatrix;
+use serde::{de, Deserialize, Deserializer, Serialize, Serializer};
+
+/// This is an intermediate type for (de)serializing `CsrMatrix`.
+///
+/// Deserialization requires using a `try_from_*` function for validation. We could have used
+/// the `remote = "Self"` trick (https://github.com/serde-rs/serde/issues/1220) which allows
+/// to directly serialize/deserialize the original fields and combine it with validation.
+/// However, this would lead to nested serialization of the `CsMatrix` and `SparsityPattern`
+/// types. Instead, we decided that we want a more human-readable serialization format using
+/// field names like `row_offsets` and `cal_indices`. The easiest way to achieve this is to
+/// introduce an intermediate type. It also allows the serialization format to stay constant
+/// even if the internal layout in `nalgebra` changes.
+///
+/// We want to avoid unnecessary copies when serializing (i.e. cloning slices into owned
+/// storage). Therefore, we use generic arguments to allow using slices during serialization and
+/// owned storage (i.e. `Vec`) during deserialization. Without a major update of serde, slices
+/// and `Vec`s should always (de)serialize identically.
+#[derive(Serialize, Deserialize)]
+struct CsrMatrixSerializationData<Indices, Values> {
+    nrows: usize,
+    ncols: usize,
+    row_offsets: Indices,
+    col_indices: Indices,
+    values: Values,
+}
+
+impl<T> Serialize for CsrMatrix<T>
+where
+    T: Serialize + Clone,
+{
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        CsrMatrixSerializationData::<&[usize], &[T]> {
+            nrows: self.nrows(),
+            ncols: self.ncols(),
+            row_offsets: self.row_offsets(),
+            col_indices: self.col_indices(),
+            values: self.values(),
+        }
+        .serialize(serializer)
+    }
+}
+
+impl<'de, T> Deserialize<'de> for CsrMatrix<T>
+where
+    T: Deserialize<'de> + Clone,
+{
+    fn deserialize<D>(deserializer: D) -> Result<CsrMatrix<T>, D::Error>
+    where
+        D: Deserializer<'de>,
+    {
+        let de = CsrMatrixSerializationData::<Vec<usize>, Vec<T>>::deserialize(deserializer)?;
+        CsrMatrix::try_from_csr_data(
+            de.nrows,
+            de.ncols,
+            de.row_offsets,
+            de.col_indices,
+            de.values,
+        )
+        .map_err(|e| de::Error::custom(e))
+    }
+}

nalgebra-sparse/src/lib.rs

@@ -160,6 +160,7 @@ pub mod ops;
 pub mod pattern;
 
 pub(crate) mod cs;
+pub(crate) mod utils;
 
 #[cfg(feature = "proptest-support")]
 pub mod proptest;

nalgebra-sparse/src/pattern.rs

@@ -1,4 +1,8 @@
 //! Sparsity patterns for CSR and CSC matrices.
+
+#[cfg(feature = "serde-serialize")]
+mod pattern_serde;
+
 use crate::cs::transpose_cs;
 use crate::SparseFormatError;
 use std::error::Error;
@@ -184,6 +188,35 @@ impl SparsityPattern {
         })
     }
 
+    /// Try to construct a sparsity pattern from the given dimensions, major offsets
+    /// and minor indices.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the number of major offsets is not exactly one greater than the major dimension
+    /// or if major offsets do not start with 0 and end with the number of minor indices.
+    pub unsafe fn from_offset_and_indices_unchecked(
+        major_dim: usize,
+        minor_dim: usize,
+        major_offsets: Vec<usize>,
+        minor_indices: Vec<usize>,
+    ) -> Self {
+        assert_eq!(major_offsets.len(), major_dim + 1);
+
+        // Check that the first and last offsets conform to the specification
+        {
+            let first_offset_ok = *major_offsets.first().unwrap() == 0;
+            let last_offset_ok = *major_offsets.last().unwrap() == minor_indices.len();
+            assert!(first_offset_ok && last_offset_ok);
+        }
+
+        Self {
+            major_offsets,
+            minor_indices,
+            minor_dim,
+        }
+    }
+
     /// An iterator over the explicitly stored "non-zero" entries (i, j).
     ///
     /// The iteration happens in a lane-major fashion, meaning that the lane index i

nalgebra-sparse/src/pattern/pattern_serde.rs

@@ -0,0 +1,56 @@
+use crate::pattern::SparsityPattern;
+use serde::{de, Deserialize, Deserializer, Serialize, Serializer};
+
+/// This is an intermediate type for (de)serializing `SparsityPattern`.
+///
+/// Deserialization requires using a `try_from_*` function for validation. We could have used
+/// the `remote = "Self"` trick (https://github.com/serde-rs/serde/issues/1220) which allows
+/// to directly serialize/deserialize the original fields and combine it with validation.
+/// However, this would lead to nested serialization of the `CsMatrix` and `SparsityPattern`
+/// types. Instead, we decided that we want a more human-readable serialization format using
+/// field names like `major_offsets` and `minor_indices`. The easiest way to achieve this is to
+/// introduce an intermediate type. It also allows the serialization format to stay constant
+/// even when the internal layout in `nalgebra` changes.
+///
+/// We want to avoid unnecessary copies when serializing (i.e. cloning slices into owned
+/// storage). Therefore, we use generic arguments to allow using slices during serialization and
+/// owned storage (i.e. `Vec`) during deserialization. Without a major update of serde, slices
+/// and `Vec`s should always (de)serialize identically.
+#[derive(Serialize, Deserialize)]
+struct SparsityPatternSerializationData<Indices> {
+    major_dim: usize,
+    minor_dim: usize,
+    major_offsets: Indices,
+    minor_indices: Indices,
+}
+
+impl Serialize for SparsityPattern {
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        SparsityPatternSerializationData::<&[usize]> {
+            major_dim: self.major_dim(),
+            minor_dim: self.minor_dim(),
+            major_offsets: self.major_offsets(),
+            minor_indices: self.minor_indices(),
+        }
+        .serialize(serializer)
+    }
+}
+
+impl<'de> Deserialize<'de> for SparsityPattern {
+    fn deserialize<D>(deserializer: D) -> Result<SparsityPattern, D::Error>
+    where
+        D: Deserializer<'de>,
+    {
+        let de = SparsityPatternSerializationData::<Vec<usize>>::deserialize(deserializer)?;
+        SparsityPattern::try_from_offsets_and_indices(
+            de.major_dim,
+            de.minor_dim,
+            de.major_offsets,
+            de.minor_indices,
+        )
+        .map_err(|e| de::Error::custom(e))
+    }
+}

nalgebra-sparse/src/utils.rs

@@ -0,0 +1,26 @@
+//! Helper functions for sparse matrix computations
+
+/// permutes entries of in_slice according to permutation slice and puts them to out_slice
+#[inline]
+pub fn apply_permutation<T: Clone>(out_slice: &mut [T], in_slice: &[T], permutation: &[usize]) {
+    assert_eq!(out_slice.len(), in_slice.len());
+    assert_eq!(out_slice.len(), permutation.len());
+    for (out_element, old_pos) in out_slice.iter_mut().zip(permutation) {
+        *out_element = in_slice[*old_pos].clone();
+    }
+}
+
+/// computes permutation by using provided indices as keys
+#[inline]
+pub fn compute_sort_permutation(permutation: &mut [usize], indices: &[usize]) {
+    assert_eq!(permutation.len(), indices.len());
+    // Set permutation to identity
+    for (i, p) in permutation.iter_mut().enumerate() {
+        *p = i;
+    }
+
+    // Compute permutation needed to bring minor indices into sorted order
+    // Note: Using sort_unstable here avoids internal allocations, which is crucial since
+    // each lane might have a small number of elements
+    permutation.sort_unstable_by_key(|idx| indices[*idx]);
+}

nalgebra-sparse/tests/serde.rs

@@ -0,0 +1,206 @@
+#![cfg(feature = "serde-serialize")]
+//! Serialization tests
+#[cfg(any(not(feature = "proptest-support"), not(feature = "compare")))]
+compile_error!("Tests must be run with features `proptest-support` and `compare`");
+
+#[macro_use]
+pub mod common;
+
+use nalgebra_sparse::coo::CooMatrix;
+use nalgebra_sparse::csc::CscMatrix;
+use nalgebra_sparse::csr::CsrMatrix;
+use nalgebra_sparse::pattern::SparsityPattern;
+
+use proptest::prelude::*;
+use serde::{Deserialize, Serialize};
+
+use crate::common::{csc_strategy, csr_strategy};
+
+fn json_roundtrip<T: Serialize + for<'a> Deserialize<'a>>(csr: &T) -> T {
+    let serialized = serde_json::to_string(csr).unwrap();
+    let deserialized: T = serde_json::from_str(&serialized).unwrap();
+    deserialized
+}
+
+#[test]
+fn pattern_roundtrip() {
+    {
+        // A pattern with zero explicitly stored entries
+        let pattern =
+            SparsityPattern::try_from_offsets_and_indices(3, 2, vec![0, 0, 0, 0], Vec::new())
+                .unwrap();
+
+        assert_eq!(json_roundtrip(&pattern), pattern);
+    }
+
+    {
+        // Arbitrary pattern
+        let offsets = vec![0, 2, 2, 5];
+        let indices = vec![0, 5, 1, 2, 3];
+        let pattern =
+            SparsityPattern::try_from_offsets_and_indices(3, 6, offsets.clone(), indices.clone())
+                .unwrap();
+
+        assert_eq!(json_roundtrip(&pattern), pattern);
+    }
+}
+
+#[test]
+#[rustfmt::skip]
+fn pattern_deserialize_invalid() {
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0,2,2,5],"minor_indices":[0,5,1,2,3]}"#).is_ok());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":0,"minor_dim":0,"major_offsets":[],"minor_indices":[]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 3, 5],"minor_indices":[0, 1, 2, 3, 5]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[1, 2, 2, 5],"minor_indices":[0, 5, 1, 2, 3]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 2, 2, 4],"minor_indices":[0, 5, 1, 2, 3]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 2, 2],"minor_indices":[0, 5, 1, 2, 3]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 3, 2, 5],"minor_indices":[0, 1, 2, 3, 4]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 2, 2, 5],"minor_indices":[0, 2, 3, 1, 4]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 2, 2, 5],"minor_indices":[0, 6, 1, 2, 3]}"#).is_err());
+    assert!(serde_json::from_str::<SparsityPattern>(r#"{"major_dim":3,"minor_dim":6,"major_offsets":[0, 2, 2, 5],"minor_indices":[0, 5, 2, 2, 3]}"#).is_err());
+}
+
+#[test]
+fn coo_roundtrip() {
+    {
+        // A COO matrix without entries
+        let matrix =
+            CooMatrix::<i32>::try_from_triplets(3, 2, Vec::new(), Vec::new(), Vec::new()).unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+
+    {
+        // Arbitrary COO matrix, no duplicates
+        let i = vec![0, 1, 0, 0, 2];
+        let j = vec![0, 2, 1, 3, 3];
+        let v = vec![2, 3, 7, 3, 1];
+        let matrix =
+            CooMatrix::<i32>::try_from_triplets(3, 5, i.clone(), j.clone(), v.clone()).unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+}
+
+#[test]
+fn coo_deserialize_invalid() {
+    // Valid matrix: {"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,1]}
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,1]}"#).is_ok());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":0,"ncols":0,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":-3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,4,5]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,8,3],"values":[2,3,7,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0],"col_indices":[0,2,1,8,3],"values":[2,3,7,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,10,0,0,2],"col_indices":[0,2,1,3,3],"values":[2,3,7,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CooMatrix<i32>>(r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2],"col_indices":[0,2,1,30,3],"values":[2,3,7,3,4]}"#).is_err());
+}
+
+#[test]
+fn coo_deserialize_duplicates() {
+    assert_eq!(
+        serde_json::from_str::<CooMatrix<i32>>(
+            r#"{"nrows":3,"ncols":5,"row_indices":[0,1,0,0,2,0,1],"col_indices":[0,2,1,3,3,0,2],"values":[2,3,7,3,1,5,6]}"#
+        ).unwrap(),
+        CooMatrix::<i32>::try_from_triplets(
+            3,
+            5,
+            vec![0, 1, 0, 0, 2, 0, 1],
+            vec![0, 2, 1, 3, 3, 0, 2],
+            vec![2, 3, 7, 3, 1, 5, 6]
+        )
+        .unwrap()
+    );
+}
+
+#[test]
+fn csc_roundtrip() {
+    {
+        // A CSC matrix with zero explicitly stored entries
+        let offsets = vec![0, 0, 0, 0];
+        let indices = vec![];
+        let values = Vec::<i32>::new();
+        let matrix = CscMatrix::try_from_csc_data(2, 3, offsets, indices, values).unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+
+    {
+        // An arbitrary CSC matrix
+        let offsets = vec![0, 2, 2, 5];
+        let indices = vec![0, 5, 1, 2, 3];
+        let values = vec![0, 1, 2, 3, 4];
+        let matrix =
+            CscMatrix::try_from_csc_data(6, 3, offsets.clone(), indices.clone(), values.clone())
+                .unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+}
+
+#[test]
+fn csc_deserialize_invalid() {
+    // Valid matrix: {"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_ok());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":0,"ncols":0,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":-6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4,5]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,8,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3,1,1],"values":[0,1,2,3,4]}"#).is_err());
+    // The following actually panics ('range end index 10 out of range for slice of length 5', nalgebra-sparse\src\pattern.rs:156:38)
+    //assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,10,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CscMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+}
+
+#[test]
+fn csr_roundtrip() {
+    {
+        // A CSR matrix with zero explicitly stored entries
+        let offsets = vec![0, 0, 0, 0];
+        let indices = vec![];
+        let values = Vec::<i32>::new();
+        let matrix = CsrMatrix::try_from_csr_data(3, 2, offsets, indices, values).unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+
+    {
+        // An arbitrary CSR matrix
+        let offsets = vec![0, 2, 2, 5];
+        let indices = vec![0, 5, 1, 2, 3];
+        let values = vec![0, 1, 2, 3, 4];
+        let matrix =
+            CsrMatrix::try_from_csr_data(3, 6, offsets.clone(), indices.clone(), values.clone())
+                .unwrap();
+
+        assert_eq!(json_roundtrip(&matrix), matrix);
+    }
+}
+
+#[test]
+fn csr_deserialize_invalid() {
+    // Valid matrix: {"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_ok());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":0,"ncols":0,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":-3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4,5]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,8,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,2,2,5],"col_indices":[0,5,1,2,3,1,1],"values":[0,1,2,3,4]}"#).is_err());
+    // The following actually panics ('range end index 10 out of range for slice of length 5', nalgebra-sparse\src\pattern.rs:156:38)
+    //assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":3,"ncols":6,"row_offsets":[0,10,2,5],"col_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+    assert!(serde_json::from_str::<CsrMatrix<i32>>(r#"{"nrows":6,"ncols":3,"col_offsets":[0,2,2,5],"row_indices":[0,5,1,2,3],"values":[0,1,2,3,4]}"#).is_err());
+}
+
+proptest! {
+    #[test]
+    fn csc_roundtrip_proptest(csc in csc_strategy()) {
+        prop_assert_eq!(json_roundtrip(&csc), csc);
+    }
+
+    #[test]
+    fn csr_roundtrip_proptest(csr in csr_strategy()) {
+        prop_assert_eq!(json_roundtrip(&csr), csr);
+    }
+}

nalgebra-sparse/tests/unit_tests/csc.rs

@@ -5,6 +5,8 @@ use nalgebra_sparse::{SparseEntry, SparseEntryMut, SparseFormatErrorKind};
 use proptest::prelude::*;
 use proptest::sample::subsequence;
 
+use super::test_data_examples::{InvalidCsDataExamples, ValidCsDataExamples};
+
 use crate::assert_panics;
 use crate::common::csc_strategy;
 
@@ -171,11 +173,26 @@ fn csc_matrix_valid_data() {
     }
 }
 
+#[test]
+fn csc_matrix_valid_data_unsorted_column_indices() {
+    let valid_data: ValidCsDataExamples = ValidCsDataExamples::new();
+
+    let (offsets, indices, values) = valid_data.valid_unsorted_cs_data;
+    let csc = CscMatrix::try_from_unsorted_csc_data(5, 4, offsets, indices, values).unwrap();
+
+    let (offsets2, indices2, values2) = valid_data.valid_cs_data;
+    let expected_csc = CscMatrix::try_from_csc_data(5, 4, offsets2, indices2, values2).unwrap();
+
+    assert_eq!(csc, expected_csc);
+}
+
 #[test]
 fn csc_matrix_try_from_invalid_csc_data() {
+    let invalid_data: InvalidCsDataExamples = InvalidCsDataExamples::new();
     {
         // Empty offset array (invalid length)
-        let matrix = CscMatrix::try_from_csc_data(0, 0, Vec::new(), Vec::new(), Vec::<u32>::new());
+        let (offsets, indices, values) = invalid_data.empty_offset_array;
+        let matrix = CscMatrix::try_from_csc_data(0, 0, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
             &SparseFormatErrorKind::InvalidStructure
@@ -184,10 +201,8 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Offset array invalid length for arbitrary data
-        let offsets = vec![0, 3, 5];
+        let (offsets, indices, values) =
-        let indices = vec![0, 1, 2, 3, 5];
+            invalid_data.offset_array_invalid_length_for_arbitrary_data;
-        let values = vec![0, 1, 2, 3, 4];
-
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -197,9 +212,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Invalid first entry in offsets array
-        let offsets = vec![1, 2, 2, 5];
+        let (offsets, indices, values) = invalid_data.invalid_first_entry_in_offsets_array;
-        let indices = vec![0, 5, 1, 2, 3];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -209,9 +222,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Invalid last entry in offsets array
-        let offsets = vec![0, 2, 2, 4];
+        let (offsets, indices, values) = invalid_data.invalid_last_entry_in_offsets_array;
-        let indices = vec![0, 5, 1, 2, 3];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -221,9 +232,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Invalid length of offsets array
-        let offsets = vec![0, 2, 2];
+        let (offsets, indices, values) = invalid_data.invalid_length_of_offsets_array;
-        let indices = vec![0, 5, 1, 2, 3];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -233,9 +242,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Nonmonotonic offsets
-        let offsets = vec![0, 3, 2, 5];
+        let (offsets, indices, values) = invalid_data.nonmonotonic_offsets;
-        let indices = vec![0, 1, 2, 3, 4];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -257,9 +264,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Minor index out of bounds
-        let offsets = vec![0, 2, 2, 5];
+        let (offsets, indices, values) = invalid_data.minor_index_out_of_bounds;
-        let indices = vec![0, 6, 1, 2, 3];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -269,9 +274,7 @@ fn csc_matrix_try_from_invalid_csc_data() {
 
     {
         // Duplicate entry
-        let offsets = vec![0, 2, 2, 5];
+        let (offsets, indices, values) = invalid_data.duplicate_entry;
-        let indices = vec![0, 5, 2, 2, 3];
-        let values = vec![0, 1, 2, 3, 4];
         let matrix = CscMatrix::try_from_csc_data(6, 3, offsets, indices, values);
         assert_eq!(
             matrix.unwrap_err().kind(),
@@ -280,6 +283,121 @@ fn csc_matrix_try_from_invalid_csc_data() {
     }
 }
 
+#[test]
+fn csc_matrix_try_from_unsorted_invalid_csc_data() {
+    let invalid_data: InvalidCsDataExamples = InvalidCsDataExamples::new();
+    {
+        // Empty offset array (invalid length)
+        let (offsets, indices, values) = invalid_data.empty_offset_array;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(0, 0, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Offset array invalid length for arbitrary data
+        let (offsets, indices, values) =
+            invalid_data.offset_array_invalid_length_for_arbitrary_data;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Invalid first entry in offsets array
+        let (offsets, indices, values) = invalid_data.invalid_first_entry_in_offsets_array;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Invalid last entry in offsets array
+        let (offsets, indices, values) = invalid_data.invalid_last_entry_in_offsets_array;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Invalid length of offsets array
+        let (offsets, indices, values) = invalid_data.invalid_length_of_offsets_array;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Nonmonotonic offsets
+        let (offsets, indices, values) = invalid_data.nonmonotonic_offsets;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+
+    {
+        // Major offset out of bounds
+        let (offsets, indices, values) = invalid_data.major_offset_out_of_bounds;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::IndexOutOfBounds
+        );
+    }
+
+    {
+        // Minor index out of bounds
+        let (offsets, indices, values) = invalid_data.minor_index_out_of_bounds;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::IndexOutOfBounds
+        );
+    }
+
+    {
+        // Duplicate entry
+        let (offsets, indices, values) = invalid_data.duplicate_entry;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::DuplicateEntry
+        );
+    }
+
+    {
+        // Duplicate entry in unsorted lane
+        let (offsets, indices, values) = invalid_data.duplicate_entry_unsorted;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::DuplicateEntry
+        );
+    }
+
+    {
+        // Wrong values length
+        let (offsets, indices, values) = invalid_data.wrong_values_length;
+        let matrix = CscMatrix::try_from_unsorted_csc_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
+}
+
 #[test]
 fn csc_disassemble_avoids_clone_when_owned() {
     // Test that disassemble avoids cloning the sparsity pattern when it holds the sole reference

nalgebra-sparse/tests/unit_tests/csr.rs

@@ -5,7 +5,7 @@ use nalgebra_sparse::{SparseEntry, SparseEntryMut, SparseFormatErrorKind};
 use proptest::prelude::*;
 use proptest::sample::subsequence;
 
-use super::test_data_examples::InvalidCsrDataExamples;
+use super::test_data_examples::{InvalidCsDataExamples, ValidCsDataExamples};
 
 use crate::assert_panics;
 use crate::common::csr_strategy;
@@ -175,30 +175,20 @@ fn csr_matrix_valid_data() {
 
 #[test]
 fn csr_matrix_valid_data_unsorted_column_indices() {
-    let csr = CsrMatrix::try_from_unsorted_csr_data(
+    let valid_data: ValidCsDataExamples = ValidCsDataExamples::new();
-        4,
-        5,
-        vec![0, 3, 5, 8, 11],
-        vec![4, 1, 3, 3, 1, 2, 3, 0, 3, 4, 1],
-        vec![5, 1, 4, 7, 4, 2, 3, 1, 8, 9, 6],
-    )
-    .unwrap();
 
-    let expected_csr = CsrMatrix::try_from_csr_data(
+    let (offsets, indices, values) = valid_data.valid_unsorted_cs_data;
-        4,
+    let csr = CsrMatrix::try_from_unsorted_csr_data(4, 5, offsets, indices, values).unwrap();
-        5,
+
-        vec![0, 3, 5, 8, 11],
+    let (offsets2, indices2, values2) = valid_data.valid_cs_data;
-        vec![1, 3, 4, 1, 3, 0, 2, 3, 1, 3, 4],
+    let expected_csr = CsrMatrix::try_from_csr_data(4, 5, offsets2, indices2, values2).unwrap();
-        vec![1, 4, 5, 4, 7, 1, 2, 3, 6, 8, 9],
-    )
-    .unwrap();
 
     assert_eq!(csr, expected_csr);
 }
 
 #[test]
 fn csr_matrix_try_from_invalid_csr_data() {
-    let invalid_data: InvalidCsrDataExamples = InvalidCsrDataExamples::new();
+    let invalid_data: InvalidCsDataExamples = InvalidCsDataExamples::new();
     {
         // Empty offset array (invalid length)
         let (offsets, indices, values) = invalid_data.empty_offset_array;
@@ -293,7 +283,7 @@ fn csr_matrix_try_from_invalid_csr_data() {
 
 #[test]
 fn csr_matrix_try_from_unsorted_invalid_csr_data() {
-    let invalid_data: InvalidCsrDataExamples = InvalidCsrDataExamples::new();
+    let invalid_data: InvalidCsDataExamples = InvalidCsDataExamples::new();
     {
         // Empty offset array (invalid length)
         let (offsets, indices, values) = invalid_data.empty_offset_array;
@@ -355,6 +345,16 @@ fn csr_matrix_try_from_unsorted_invalid_csr_data() {
         );
     }
 
+    {
+        // Major offset out of bounds
+        let (offsets, indices, values) = invalid_data.major_offset_out_of_bounds;
+        let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::IndexOutOfBounds
+        );
+    }
+
     {
         // Minor index out of bounds
         let (offsets, indices, values) = invalid_data.minor_index_out_of_bounds;
@@ -374,6 +374,26 @@ fn csr_matrix_try_from_unsorted_invalid_csr_data() {
             &SparseFormatErrorKind::DuplicateEntry
         );
     }
+
+    {
+        // Duplicate entry in unsorted lane
+        let (offsets, indices, values) = invalid_data.duplicate_entry_unsorted;
+        let matrix = CsrMatrix::try_from_unsorted_csr_data(3, 6, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::DuplicateEntry
+        );
+    }
+
+    {
+        // Wrong values length
+        let (offsets, indices, values) = invalid_data.wrong_values_length;
+        let matrix = CsrMatrix::try_from_unsorted_csr_data(6, 3, offsets, indices, values);
+        assert_eq!(
+            matrix.unwrap_err().kind(),
+            &SparseFormatErrorKind::InvalidStructure
+        );
+    }
 }
 
 #[test]

nalgebra-sparse/tests/unit_tests/test_data_examples.rs

@@ -1,5 +1,31 @@
-/// Examples of *invalid* raw CSR data `(offsets, indices, values)`.
+/// Examples of *valid* raw CS data `(offsets, indices, values)`.
-pub struct InvalidCsrDataExamples {
+pub struct ValidCsDataExamples {
+    pub valid_cs_data: (Vec<usize>, Vec<usize>, Vec<i32>),
+    pub valid_unsorted_cs_data: (Vec<usize>, Vec<usize>, Vec<i32>),
+}
+
+impl ValidCsDataExamples {
+    pub fn new() -> Self {
+        let valid_cs_data = (
+            vec![0, 3, 5, 8, 11],
+            vec![1, 3, 4, 1, 3, 0, 2, 3, 1, 3, 4],
+            vec![1, 4, 5, 4, 7, 1, 2, 3, 6, 8, 9],
+        );
+        let valid_unsorted_cs_data = (
+            vec![0, 3, 5, 8, 11],
+            vec![4, 1, 3, 3, 1, 2, 3, 0, 3, 4, 1],
+            vec![5, 1, 4, 7, 4, 2, 3, 1, 8, 9, 6],
+        );
+
+        return Self {
+            valid_cs_data,
+            valid_unsorted_cs_data,
+        };
+    }
+}
+
+/// Examples of *invalid* raw CS data `(offsets, indices, values)`.
+pub struct InvalidCsDataExamples {
     pub empty_offset_array: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub offset_array_invalid_length_for_arbitrary_data: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub invalid_first_entry_in_offsets_array: (Vec<usize>, Vec<usize>, Vec<i32>),
@@ -7,11 +33,14 @@ pub struct InvalidCsrDataExamples {
     pub invalid_length_of_offsets_array: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub nonmonotonic_offsets: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub nonmonotonic_minor_indices: (Vec<usize>, Vec<usize>, Vec<i32>),
+    pub major_offset_out_of_bounds: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub minor_index_out_of_bounds: (Vec<usize>, Vec<usize>, Vec<i32>),
     pub duplicate_entry: (Vec<usize>, Vec<usize>, Vec<i32>),
+    pub duplicate_entry_unsorted: (Vec<usize>, Vec<usize>, Vec<i32>),
+    pub wrong_values_length: (Vec<usize>, Vec<usize>, Vec<i32>),
 }
 
-impl InvalidCsrDataExamples {
+impl InvalidCsDataExamples {
     pub fn new() -> Self {
         let empty_offset_array = (Vec::<usize>::new(), Vec::<usize>::new(), Vec::<i32>::new());
         let offset_array_invalid_length_for_arbitrary_data =
@@ -25,9 +54,13 @@ impl InvalidCsrDataExamples {
         let nonmonotonic_offsets = (vec![0, 3, 2, 5], vec![0, 1, 2, 3, 4], vec![0, 1, 2, 3, 4]);
         let nonmonotonic_minor_indices =
             (vec![0, 2, 2, 5], vec![0, 2, 3, 1, 4], vec![0, 1, 2, 3, 4]);
+        let major_offset_out_of_bounds =
+            (vec![0, 7, 2, 5], vec![0, 2, 3, 1, 4], vec![0, 1, 2, 3, 4]);
         let minor_index_out_of_bounds =
             (vec![0, 2, 2, 5], vec![0, 6, 1, 2, 3], vec![0, 1, 2, 3, 4]);
-        let duplicate_entry = (vec![0, 2, 2, 5], vec![0, 5, 2, 2, 3], vec![0, 1, 2, 3, 4]);
+        let duplicate_entry = (vec![0, 1, 2, 5], vec![1, 3, 2, 3, 3], vec![0, 1, 2, 3, 4]);
+        let duplicate_entry_unsorted = (vec![0, 1, 4, 5], vec![1, 3, 2, 3, 3], vec![0, 1, 2, 3, 4]);
+        let wrong_values_length = (vec![0, 1, 2, 5], vec![1, 3, 2, 3, 0], vec![5, 4]);
 
         return Self {
             empty_offset_array,
@@ -37,8 +70,11 @@ impl InvalidCsrDataExamples {
             invalid_length_of_offsets_array,
             nonmonotonic_minor_indices,
             nonmonotonic_offsets,
+            major_offset_out_of_bounds,
             minor_index_out_of_bounds,
             duplicate_entry,
+            duplicate_entry_unsorted,
+            wrong_values_length,
         };
     }
 }

src/base/array_storage.rs

@@ -1,7 +1,5 @@
 use std::fmt::{self, Debug, Formatter};
 // use std::hash::{Hash, Hasher};
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 use std::ops::Mul;
 
 #[cfg(feature = "serde-serialize-no-std")]
@@ -13,9 +11,6 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer};
 #[cfg(feature = "serde-serialize-no-std")]
 use std::marker::PhantomData;
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use crate::base::allocator::Allocator;
 use crate::base::default_allocator::DefaultAllocator;
 use crate::base::dimension::{Const, ToTypenum};
@@ -32,10 +27,12 @@ use std::mem;
 /// A array-based statically sized matrix data storage.
 #[repr(transparent)]
 #[derive(Copy, Clone, PartialEq, Eq, Hash)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct ArrayStorage<T, const R: usize, const C: usize>(pub [[T; R]; C]);
 
 impl<T, const R: usize, const C: usize> ArrayStorage<T, R, C> {
@@ -281,71 +278,3 @@ unsafe impl<T: Scalar + Copy + bytemuck::Pod, const R: usize, const C: usize> by
     for ArrayStorage<T, R, C>
 {
 }
-
-#[cfg(feature = "abomonation-serialize")]
-impl<T, const R: usize, const C: usize> Abomonation for ArrayStorage<T, R, C>
-where
-    T: Scalar + Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        for element in self.as_slice() {
-            element.entomb(writer)?;
-        }
-
-        Ok(())
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, mut bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        for element in self.as_mut_slice() {
-            let temp = bytes;
-            bytes = element.exhume(temp)?
-        }
-        Some(bytes)
-    }
-
-    fn extent(&self) -> usize {
-        self.as_slice().iter().fold(0, |acc, e| acc + e.extent())
-    }
-}
-
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::ArrayStorage;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive, const R: usize, const C: usize> Archive for ArrayStorage<T, R, C> {
-        type Archived = ArrayStorage<T::Archived, R, C>;
-        type Resolver = <[[T; R]; C] as Archive>::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.0.resolve(
-                pos + offset_of!(Self::Archived, 0),
-                resolver,
-                project_struct!(out: Self::Archived => 0),
-            );
-        }
-    }
-
-    impl<T: Serialize<S>, S: Fallible + ?Sized, const R: usize, const C: usize> Serialize<S>
-        for ArrayStorage<T, R, C>
-    {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            self.0.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, D: Fallible + ?Sized, const R: usize, const C: usize>
-        Deserialize<ArrayStorage<T, R, C>, D> for ArrayStorage<T::Archived, R, C>
-    where
-        T::Archived: Deserialize<T, D>,
-    {
-        fn deserialize(&self, deserializer: &mut D) -> Result<ArrayStorage<T, R, C>, D::Error> {
-            Ok(ArrayStorage(self.0.deserialize(deserializer)?))
-        }
-    }
-}

src/base/blas.rs

@@ -175,8 +175,7 @@ where
     /// Note that this is **not** the matrix multiplication as in, e.g., numpy. For matrix
     /// multiplication, use one of: `.gemm`, `.mul_to`, `.mul`, the `*` operator.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Vector3, Matrix2x3};
     /// let vec1 = Vector3::new(1.0, 2.0, 3.0);
@@ -207,8 +206,7 @@ where
     /// Note that this is **not** the matrix multiplication as in, e.g., numpy. For matrix
     /// multiplication, use one of: `.gemm`, `.mul_to`, `.mul`, the `*` operator.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Vector2, Complex};
     /// let vec1 = Vector2::new(Complex::new(1.0, 2.0), Complex::new(3.0, 4.0));
@@ -232,8 +230,7 @@ where
 
     /// The dot product between the transpose of `self` and `rhs`.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Vector3, RowVector3, Matrix2x3, Matrix3x2};
     /// let vec1 = Vector3::new(1.0, 2.0, 3.0);
@@ -285,8 +282,7 @@ where
     ///
     /// If `b` is zero, `self` is never read from.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Vector3;
     /// let mut vec1 = Vector3::new(1.0, 2.0, 3.0);
@@ -308,8 +304,7 @@ where
     ///
     /// If `b` is zero, `self` is never read from.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Vector3;
     /// let mut vec1 = Vector3::new(1.0, 2.0, 3.0);
@@ -333,8 +328,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2};
     /// let mut vec1 = Vector2::new(1.0, 2.0);
@@ -425,8 +419,7 @@ where
     /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part
     /// (including the diagonal) is actually read.
     ///
-    /// # Examples:
+    /// # Examples
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2};
     /// let mat = Matrix2::new(1.0, 2.0,
@@ -468,8 +461,7 @@ where
     /// If `beta` is zero, `self` is never read. If `self` is read, only its lower-triangular part
     /// (including the diagonal) is actually read.
     ///
-    /// # Examples:
+    /// # Examples
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2, Complex};
     /// let mat = Matrix2::new(Complex::new(1.0, 0.0), Complex::new(2.0, -0.1),
@@ -552,8 +544,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2};
     /// let mat = Matrix2::new(1.0, 3.0,
@@ -587,8 +578,7 @@ where
     /// For real matrices, this is the same as `.gemv_tr`.
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2, Complex};
     /// let mat = Matrix2::new(Complex::new(1.0, 2.0), Complex::new(3.0, 4.0),
@@ -656,8 +646,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2x3, Vector2, Vector3};
     /// let mut mat = Matrix2x3::repeat(4.0);
@@ -688,8 +677,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix2x3, Vector2, Vector3, Complex};
@@ -722,8 +710,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix2x3, Matrix3x4, Matrix2x4};
@@ -763,8 +750,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix3x2, Matrix3x4, Matrix2x4};
@@ -821,8 +807,7 @@ where
     ///
     /// If `beta` is zero, `self` is never read.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix3x2, Matrix3x4, Matrix2x4, Complex};
@@ -921,8 +906,7 @@ where
     /// If `beta` is zero, `self` is never read. The result is symmetric. Only the lower-triangular
     /// (including the diagonal) part of `self` is read/written.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2};
     /// let mut mat = Matrix2::identity();
@@ -934,6 +918,7 @@ where
     /// mat.ger_symm(10.0, &vec1, &vec2, 5.0);
     /// assert_eq!(mat.lower_triangle(), expected.lower_triangle());
     /// assert_eq!(mat.m12, 99999.99999); // This was untouched.
+    /// ```
     #[inline]
     #[deprecated(note = "This is renamed `syger` to match the original BLAS terminology.")]
     pub fn ger_symm<D2: Dim, D3: Dim, SB, SC>(
@@ -958,8 +943,7 @@ where
     /// If `beta` is zero, `self` is never read. The result is symmetric. Only the lower-triangular
     /// (including the diagonal) part of `self` is read/written.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2};
     /// let mut mat = Matrix2::identity();
@@ -971,6 +955,7 @@ where
     /// mat.syger(10.0, &vec1, &vec2, 5.0);
     /// assert_eq!(mat.lower_triangle(), expected.lower_triangle());
     /// assert_eq!(mat.m12, 99999.99999); // This was untouched.
+    /// ```
     #[inline]
     pub fn syger<D2: Dim, D3: Dim, SB, SC>(
         &mut self,
@@ -993,8 +978,7 @@ where
     /// If `beta` is zero, `self` is never read. The result is symmetric. Only the lower-triangular
     /// (including the diagonal) part of `self` is read/written.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Matrix2, Vector2, Complex};
     /// let mut mat = Matrix2::identity();
@@ -1006,6 +990,7 @@ where
     /// mat.hegerc(Complex::new(10.0, 20.0), &vec1, &vec2, Complex::new(5.0, 15.0));
     /// assert_eq!(mat.lower_triangle(), expected.lower_triangle());
     /// assert_eq!(mat.m12, Complex::new(99999.99999, 88888.88888)); // This was untouched.
+    /// ```
     #[inline]
     pub fn hegerc<D2: Dim, D3: Dim, SB, SC>(
         &mut self,
@@ -1031,8 +1016,7 @@ where
     ///
     /// This uses the provided workspace `work` to avoid allocations for intermediate results.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{DMatrix, DVector};
@@ -1053,6 +1037,7 @@ where
     ///
     /// mat.quadform_tr_with_workspace(&mut workspace, 10.0, &lhs, &mid, 5.0);
     /// assert_relative_eq!(mat, expected);
+    /// ```
     pub fn quadform_tr_with_workspace<D2, S2, R3, C3, S3, D4, S4>(
         &mut self,
         work: &mut Vector<T, D2, S2>,
@@ -1085,8 +1070,7 @@ where
     /// If `D1` is a type-level integer, then the allocation is performed on the stack.
     /// Use `.quadform_tr_with_workspace(...)` instead to avoid allocations.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix2, Matrix3, Matrix2x3, Vector2};
@@ -1100,6 +1084,7 @@ where
     ///
     /// mat.quadform_tr(10.0, &lhs, &mid, 5.0);
     /// assert_relative_eq!(mat, expected);
+    /// ```
     pub fn quadform_tr<R3, C3, S3, D4, S4>(
         &mut self,
         alpha: T,
@@ -1124,6 +1109,7 @@ where
     ///
     /// This uses the provided workspace `work` to avoid allocations for intermediate results.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{DMatrix, DVector};
@@ -1145,6 +1131,7 @@ where
     ///
     /// mat.quadform_with_workspace(&mut workspace, 10.0, &mid, &rhs, 5.0);
     /// assert_relative_eq!(mat, expected);
+    /// ```
     pub fn quadform_with_workspace<D2, S2, D3, S3, R4, C4, S4>(
         &mut self,
         work: &mut Vector<T, D2, S2>,
@@ -1180,6 +1167,7 @@ where
     /// If `D2` is a type-level integer, then the allocation is performed on the stack.
     /// Use `.quadform_with_workspace(...)` instead to avoid allocations.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix2, Matrix3x2, Matrix3};
@@ -1194,6 +1182,7 @@ where
     ///
     /// mat.quadform(10.0, &mid, &rhs, 5.0);
     /// assert_relative_eq!(mat, expected);
+    /// ```
     pub fn quadform<D2, S2, R3, C3, S3>(
         &mut self,
         alpha: T,

src/base/dimension.rs

@@ -13,10 +13,12 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
 /// Dim of dynamically-sized algebraic entities.
 #[derive(Clone, Copy, Eq, PartialEq, Debug)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Dynamic {
     value: usize,
 }
@@ -202,9 +204,11 @@ dim_ops!(
 
 #[derive(Debug, Copy, Clone, PartialEq, Eq, Hash)]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Const<const R: usize>;
 
 /// Trait implemented exclusively by type-level integers.
@@ -239,37 +243,6 @@ impl<'de, const D: usize> Deserialize<'de> for Const<D> {
     }
 }
 
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Const;
-    use rkyv::{Archive, Deserialize, Fallible, Serialize};
-
-    impl<const R: usize> Archive for Const<R> {
-        type Archived = Self;
-        type Resolver = ();
-
-        fn resolve(
-            &self,
-            _: usize,
-            _: Self::Resolver,
-            _: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-        }
-    }
-
-    impl<S: Fallible + ?Sized, const R: usize> Serialize<S> for Const<R> {
-        fn serialize(&self, _: &mut S) -> Result<Self::Resolver, S::Error> {
-            Ok(())
-        }
-    }
-
-    impl<D: Fallible + ?Sized, const R: usize> Deserialize<Self, D> for Const<R> {
-        fn deserialize(&self, _: &mut D) -> Result<Self, D::Error> {
-            Ok(Const)
-        }
-    }
-}
-
 pub trait ToConst {
     type Const: DimName;
 }
@@ -309,24 +282,24 @@ impl<const T: usize> DimName for Const<T> {
 
 pub type U1 = Const<1>;
 
-impl ToTypenum for Const<{ typenum::U1::USIZE }> {
+impl ToTypenum for Const<1> {
     type Typenum = typenum::U1;
 }
 
 impl ToConst for typenum::U1 {
-    type Const = Const<{ typenum::U1::USIZE }>;
+    type Const = Const<1>;
 }
 
 macro_rules! from_to_typenum (
-    ($($D: ident),* $(,)*) => {$(
+    ($($D: ident, $VAL: expr);* $(;)*) => {$(
-        pub type $D = Const<{ typenum::$D::USIZE }>;
+        pub type $D = Const<$VAL>;
 
-        impl ToTypenum for Const<{ typenum::$D::USIZE }> {
+        impl ToTypenum for Const<$VAL> {
             type Typenum = typenum::$D;
         }
 
         impl ToConst for typenum::$D {
-            type Const = Const<{ typenum::$D::USIZE }>;
+            type Const = Const<$VAL>;
         }
 
         impl IsNotStaticOne for $D { }
@@ -334,12 +307,12 @@ macro_rules! from_to_typenum (
 );
 
 from_to_typenum!(
-    U0, /*U1,*/ U2, U3, U4, U5, U6, U7, U8, U9, U10, U11, U12, U13, U14, U15, U16, U17, U18,
+    U0, 0; /*U1,1;*/ U2, 2; U3, 3; U4, 4; U5, 5; U6, 6; U7, 7; U8, 8; U9, 9; U10, 10; U11, 11; U12, 12; U13, 13; U14, 14; U15, 15; U16, 16; U17, 17; U18, 18;
-    U19, U20, U21, U22, U23, U24, U25, U26, U27, U28, U29, U30, U31, U32, U33, U34, U35, U36, U37,
+    U19, 19; U20, 20; U21, 21; U22, 22; U23, 23; U24, 24; U25, 25; U26, 26; U27, 27; U28, 28; U29, 29; U30, 30; U31, 31; U32, 32; U33, 33; U34, 34; U35, 35; U36, 36; U37, 37;
-    U38, U39, U40, U41, U42, U43, U44, U45, U46, U47, U48, U49, U50, U51, U52, U53, U54, U55, U56,
+    U38, 38; U39, 39; U40, 40; U41, 41; U42, 42; U43, 43; U44, 44; U45, 45; U46, 46; U47, 47; U48, 48; U49, 49; U50, 50; U51, 51; U52, 52; U53, 53; U54, 54; U55, 55; U56, 56;
-    U57, U58, U59, U60, U61, U62, U63, U64, U65, U66, U67, U68, U69, U70, U71, U72, U73, U74, U75,
+    U57, 57; U58, 58; U59, 59; U60, 60; U61, 61; U62, 62; U63, 63; U64, 64; U65, 65; U66, 66; U67, 67; U68, 68; U69, 69; U70, 70; U71, 71; U72, 72; U73, 73; U74, 74; U75, 75;
-    U76, U77, U78, U79, U80, U81, U82, U83, U84, U85, U86, U87, U88, U89, U90, U91, U92, U93, U94,
+    U76, 76; U77, 77; U78, 78; U79, 79; U80, 80; U81, 81; U82, 82; U83, 83; U84, 84; U85, 85; U86, 86; U87, 87; U88, 88; U89, 89; U90, 90; U91, 91; U92, 92; U93, 93; U94, 94;
-    U95, U96, U97, U98, U99, U100, U101, U102, U103, U104, U105, U106, U107, U108, U109, U110,
+    U95, 95; U96, 96; U97, 97; U98, 98; U99, 99; U100, 100; U101, 101; U102, 102; U103, 103; U104, 104; U105, 105; U106, 106; U107, 107; U108, 108; U109, 109; U110, 110;
-    U111, U112, U113, U114, U115, U116, U117, U118, U119, U120, U121, U122, U123, U124, U125, U126,
+    U111, 111; U112, 112; U113, 113; U114, 114; U115, 115; U116, 116; U117, 117; U118, 118; U119, 119; U120, 120; U121, 121; U122, 122; U123, 123; U124, 124; U125, 125; U126, 126;
-    U127
+    U127, 127
 );

src/base/matrix.rs

@@ -1,6 +1,4 @@
 use num::{One, Zero};
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use std::any::TypeId;
@@ -13,9 +11,6 @@ use std::mem;
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, Field, SupersetOf};
 use simba::simd::SimdPartialOrd;
 
@@ -155,10 +150,12 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
 /// some concrete types for `T` and a compatible data storage type `S`).
 #[repr(C)]
 #[derive(Clone, Copy)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Matrix<T, R, C, S> {
     /// The data storage that contains all the matrix components. Disappointed?
     ///
@@ -254,21 +251,6 @@ where
     }
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T: Scalar, R: Dim, C: Dim, S: Abomonation> Abomonation for Matrix<T, R, C, S> {
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.data.entomb(writer)
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.data.exhume(bytes)
-    }
-
-    fn extent(&self) -> usize {
-        self.data.extent()
-    }
-}
-
 #[cfg(feature = "compare")]
 impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> matrixcompare_core::Matrix<T>
     for Matrix<T, R, C, S>
@@ -311,53 +293,6 @@ where
 {
 }
 
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Matrix;
-    use core::marker::PhantomData;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive, R: Archive, C: Archive, S: Archive> Archive for Matrix<T, R, C, S> {
-        type Archived = Matrix<T::Archived, R::Archived, C::Archived, S::Archived>;
-        type Resolver = S::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.data.resolve(
-                pos + offset_of!(Self::Archived, data),
-                resolver,
-                project_struct!(out: Self::Archived => data),
-            );
-        }
-    }
-
-    impl<T: Archive, R: Archive, C: Archive, S: Serialize<_S>, _S: Fallible + ?Sized> Serialize<_S>
-        for Matrix<T, R, C, S>
-    {
-        fn serialize(&self, serializer: &mut _S) -> Result<Self::Resolver, _S::Error> {
-            self.data.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, R: Archive, C: Archive, S: Archive, D: Fallible + ?Sized>
-        Deserialize<Matrix<T, R, C, S>, D>
-        for Matrix<T::Archived, R::Archived, C::Archived, S::Archived>
-    where
-        S::Archived: Deserialize<S, D>,
-    {
-        fn deserialize(&self, deserializer: &mut D) -> Result<Matrix<T, R, C, S>, D::Error> {
-            Ok(Matrix {
-                data: self.data.deserialize(deserializer)?,
-                _phantoms: PhantomData,
-            })
-        }
-    }
-}
-
 impl<T, R, C, S> Matrix<T, R, C, S> {
     /// Creates a new matrix with the given data without statically checking that the matrix
     /// dimension matches the storage dimension.
@@ -434,8 +369,6 @@ where
 {
     /// Assumes a matrix's entries to be initialized. This operation should be near zero-cost.
     ///
-    /// For the similar method that operates on matrix slices, see [`slice_assume_init`].
-    ///
     /// # Safety
     /// The user must make sure that every single entry of the buffer has been initialized,
     /// or Undefined Behavior will immediately occur.
@@ -456,12 +389,12 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
 
     /// The shape of this matrix returned as the tuple (number of rows, number of columns).
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Matrix3x4;
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.shape(), (3, 4));
+    /// ```
     #[inline]
     #[must_use]
     pub fn shape(&self) -> (usize, usize) {
@@ -478,12 +411,12 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
 
     /// The number of rows of this matrix.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Matrix3x4;
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.nrows(), 3);
+    /// ```
     #[inline]
     #[must_use]
     pub fn nrows(&self) -> usize {
@@ -492,12 +425,12 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
 
     /// The number of columns of this matrix.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Matrix3x4;
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.ncols(), 4);
+    /// ```
     #[inline]
     #[must_use]
     pub fn ncols(&self) -> usize {
@@ -506,14 +439,14 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
 
     /// The strides (row stride, column stride) of this matrix.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::DMatrix;
     /// let mat = DMatrix::<f32>::zeros(10, 10);
     /// let slice = mat.slice_with_steps((0, 0), (5, 3), (1, 2));
     /// // The column strides is the number of steps (here 2) multiplied by the corresponding dimension.
     /// assert_eq!(mat.strides(), (1, 10));
+    /// ```
     #[inline]
     #[must_use]
     pub fn strides(&self) -> (usize, usize) {
@@ -1108,8 +1041,7 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
 impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
     /// Iterates through this matrix coordinates in column-major order.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::Matrix2x3;
     /// let mat = Matrix2x3::new(11, 12, 13,
@@ -1122,6 +1054,7 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(*it.next().unwrap(), 13);
     /// assert_eq!(*it.next().unwrap(), 23);
     /// assert!(it.next().is_none());
+    /// ```
     #[inline]
     pub fn iter(&self) -> MatrixIter<'_, T, R, C, S> {
         MatrixIter::new(&self.data)
@@ -1144,6 +1077,7 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
     }
 
     /// Iterate through the columns of this matrix.
+    ///
     /// # Example
     /// ```
     /// # use nalgebra::Matrix2x3;

src/base/unit.rs

@@ -1,14 +1,9 @@
 use std::fmt;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 use std::ops::Deref;
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use crate::allocator::Allocator;
 use crate::base::DefaultAllocator;
 use crate::storage::RawStorage;
@@ -26,10 +21,12 @@ use crate::{Dim, Matrix, OMatrix, RealField, Scalar, SimdComplexField, SimdRealF
 /// in their documentation, read their dedicated pages directly.
 #[repr(transparent)]
 #[derive(Clone, Hash, Copy)]
-// #[cfg_attr(
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-//     all(not(target_os = "cuda"), feature = "cuda"),
+#[cfg_attr(
-//     derive(cust::DeviceCopy)
+    feature = "rkyv-serialize-no-std",
-// )]
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
+)]
+// #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Unit<T> {
     pub(crate) value: T,
 }
@@ -66,65 +63,8 @@ impl<'de, T: Deserialize<'de>> Deserialize<'de> for Unit<T> {
     }
 }
 
-#[cfg(feature = "abomonation-serialize")]
+#[cfg(feature = "cuda")]
-impl<T: Abomonation> Abomonation for Unit<T> {
+unsafe impl<T: cust_core::DeviceCopy, R, C, S> cust_core::DeviceCopy for Unit<Matrix<T, R, C, S>>
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.value.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.value.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.value.exhume(bytes)
-    }
-}
-
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Unit;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive> Archive for Unit<T> {
-        type Archived = Unit<T::Archived>;
-        type Resolver = T::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut ::core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.value.resolve(
-                pos + offset_of!(Self::Archived, value),
-                resolver,
-                project_struct!(out: Self::Archived => value),
-            );
-        }
-    }
-
-    impl<T: Serialize<S>, S: Fallible + ?Sized> Serialize<S> for Unit<T> {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            self.value.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, D: Fallible + ?Sized> Deserialize<Unit<T>, D> for Unit<T::Archived>
-    where
-        T::Archived: Deserialize<T, D>,
-    {
-        fn deserialize(&self, deserializer: &mut D) -> Result<Unit<T>, D::Error> {
-            Ok(Unit {
-                value: self.value.deserialize(deserializer)?,
-            })
-        }
-    }
-}
-
-#[cfg(all(not(target_os = "cuda"), feature = "cuda"))]
-unsafe impl<T: cust::memory::DeviceCopy, R, C, S> cust::memory::DeviceCopy
-    for Unit<Matrix<T, R, C, S>>
 where
     T: Scalar,
     R: Dim,

src/base/vec_storage.rs

@@ -1,6 +1,3 @@
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
-
 #[cfg(all(feature = "alloc", not(feature = "std")))]
 use alloc::vec::Vec;
 
@@ -18,8 +15,6 @@ use serde::{
 };
 
 use crate::Storage;
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
 use std::mem::MaybeUninit;
 
 /*
@@ -402,21 +397,6 @@ where
     }
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T: Abomonation, R: Dim, C: Dim> Abomonation for VecStorage<T, R, C> {
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.data.entomb(writer)
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.data.exhume(bytes)
-    }
-
-    fn extent(&self) -> usize {
-        self.data.extent()
-    }
-}
-
 impl<T, R: Dim> Extend<T> for VecStorage<T, R, Dynamic> {
     /// Extends the number of columns of the `VecStorage` with elements
     /// from the given iterator.

src/geometry/dual_quaternion.rs

@@ -19,6 +19,7 @@ use simba::scalar::{ClosedNeg, RealField};
 /// `DualQuaternions` are stored as \[..real, ..dual\].
 /// Both of the quaternion components are laid out in `i, j, k, w` order.
 ///
+/// # Example
 /// ```
 /// # use nalgebra::{DualQuaternion, Quaternion};
 ///
@@ -39,10 +40,12 @@ use simba::scalar::{ClosedNeg, RealField};
 ///  See <https://github.com/dimforge/nalgebra/issues/487>
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct DualQuaternion<T> {
     /// The real component of the quaternion
     pub real: Quaternion<T>,
@@ -623,6 +626,7 @@ where
     ///     dq.rotation().euler_angles().0, std::f32::consts::FRAC_PI_2, epsilon = 1.0e-6
     /// );
     /// assert_relative_eq!(dq.translation().vector.y, 3.0, epsilon = 1.0e-6);
+    /// ```
     #[inline]
     #[must_use]
     pub fn sclerp(&self, other: &Self, t: T) -> Self
@@ -713,6 +717,7 @@ where
 
     /// Return the rotation part of this unit dual quaternion.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3};
@@ -733,6 +738,7 @@ where
 
     /// Return the translation part of this unit dual quaternion.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3};
@@ -758,6 +764,7 @@ where
 
     /// Builds an isometry from this unit dual quaternion.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3};
@@ -783,6 +790,7 @@ where
     ///
     /// This is the same as the multiplication `self * pt`.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3, Point3};
@@ -807,6 +815,7 @@ where
     ///
     /// This is the same as the multiplication `self * v`.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3};
@@ -831,6 +840,7 @@ where
     /// This may be cheaper than inverting the unit dual quaternion and
     /// transforming the point.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3, Point3};
@@ -856,6 +866,7 @@ where
     /// This may be cheaper than inverting the unit dual quaternion and
     /// transforming the vector.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3};
@@ -880,6 +891,7 @@ where
     /// cheaper than inverting the unit dual quaternion and transforming the
     /// vector.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Unit, Vector3};
@@ -909,6 +921,7 @@ where
     /// Converts this unit dual quaternion interpreted as an isometry
     /// into its equivalent homogeneous transformation matrix.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Matrix4, UnitDualQuaternion, UnitQuaternion, Vector3};

src/geometry/dual_quaternion_construction.rs

@@ -27,7 +27,6 @@ impl<T: Scalar> DualQuaternion<T> {
     /// The dual quaternion multiplicative identity.
     ///
     /// # Example
-    ///
     /// ```
     /// # use nalgebra::{DualQuaternion, Quaternion};
     ///
@@ -134,6 +133,7 @@ impl<T: SimdRealField> UnitDualQuaternion<T> {
     /// The unit dual quaternion multiplicative identity, which also represents
     /// the identity transformation as an isometry.
     ///
+    /// # Example
     /// ```
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3, Point3};
     /// let ident = UnitDualQuaternion::identity();
@@ -171,6 +171,7 @@ where
     /// Return a dual quaternion representing the translation and orientation
     /// given by the provided rotation quaternion and translation vector.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{UnitDualQuaternion, UnitQuaternion, Vector3, Point3};
@@ -196,6 +197,7 @@ where
     /// Return a unit dual quaternion representing the translation and orientation
     /// given by the provided isometry.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::{Isometry3, UnitDualQuaternion, UnitQuaternion, Vector3, Point3};

src/geometry/isometry.rs

@@ -1,15 +1,10 @@
 use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use std::fmt;
 use std::hash;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Serialize};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::{RealField, SubsetOf};
 use simba::simd::SimdRealField;
 
@@ -55,10 +50,7 @@ use crate::geometry::{AbstractRotation, Point, Translation};
 ///
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
-#[cfg_attr(
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
-    all(not(target_os = "cuda"), feature = "cuda"),
-    derive(cust::DeviceCopy)
-)]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
     feature = "serde-serialize-no-std",
@@ -74,6 +66,11 @@ use crate::geometry::{AbstractRotation, Point, Translation};
                        Owned<T, Const<D>>: Deserialize<'de>,
                        T: Scalar"))
 )]
+#[cfg_attr(
+    feature = "rkyv-serialize-no-std",
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
+)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 pub struct Isometry<T, R, const D: usize> {
     /// The pure rotational part of this isometry.
     pub rotation: R,
@@ -81,89 +78,6 @@ pub struct Isometry<T, R, const D: usize> {
     pub translation: Translation<T, D>,
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T, R, const D: usize> Abomonation for Isometry<T, R, D>
-where
-    T: SimdRealField,
-    R: Abomonation,
-    Translation<T, D>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.rotation.entomb(writer)?;
-        self.translation.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.rotation.extent() + self.translation.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.rotation
-            .exhume(bytes)
-            .and_then(|bytes| self.translation.exhume(bytes))
-    }
-}
-
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Isometry;
-    use crate::{base::Scalar, geometry::Translation};
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Scalar + Archive, R: Archive, const D: usize> Archive for Isometry<T, R, D>
-    where
-        T::Archived: Scalar,
-    {
-        type Archived = Isometry<T::Archived, R::Archived, D>;
-        type Resolver = (R::Resolver, <Translation<T, D> as Archive>::Resolver);
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.rotation.resolve(
-                pos + offset_of!(Self::Archived, rotation),
-                resolver.0,
-                project_struct!(out: Self::Archived => rotation),
-            );
-            self.translation.resolve(
-                pos + offset_of!(Self::Archived, translation),
-                resolver.1,
-                project_struct!(out: Self::Archived => translation),
-            );
-        }
-    }
-
-    impl<T: Scalar + Serialize<S>, R: Serialize<S>, S: Fallible + ?Sized, const D: usize>
-        Serialize<S> for Isometry<T, R, D>
-    where
-        T::Archived: Scalar,
-    {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            Ok((
-                self.rotation.serialize(serializer)?,
-                self.translation.serialize(serializer)?,
-            ))
-        }
-    }
-
-    impl<T: Scalar + Archive, R: Archive, _D: Fallible + ?Sized, const D: usize>
-        Deserialize<Isometry<T, R, D>, _D> for Isometry<T::Archived, R::Archived, D>
-    where
-        T::Archived: Scalar + Deserialize<T, _D>,
-        R::Archived: Scalar + Deserialize<R, _D>,
-    {
-        fn deserialize(&self, deserializer: &mut _D) -> Result<Isometry<T, R, D>, _D::Error> {
-            Ok(Isometry {
-                rotation: self.rotation.deserialize(deserializer)?,
-                translation: self.translation.deserialize(deserializer)?,
-            })
-        }
-    }
-}
-
 impl<T: Scalar + hash::Hash, R: hash::Hash, const D: usize> hash::Hash for Isometry<T, R, D>
 where
     Owned<T, Const<D>>: hash::Hash,

src/geometry/orthographic.rs

@@ -19,10 +19,12 @@ use crate::geometry::{Point3, Projective3};
 
 /// A 3D orthographic projection stored as a homogeneous 4x4 matrix.
 #[repr(C)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Orthographic3<T> {
     matrix: Matrix4<T>,
@@ -319,6 +321,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The left offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -336,6 +339,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The right offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -353,6 +357,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The bottom offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -370,6 +375,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The top offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -387,6 +393,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The near plane offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -404,6 +411,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// The far plane offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -526,6 +534,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the left offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -545,6 +554,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the right offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -564,6 +574,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the bottom offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -583,6 +594,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the top offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -602,6 +614,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the near plane offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -621,6 +634,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the far plane offset of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -640,6 +654,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the view cuboid offsets along the `x` axis.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -665,6 +680,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the view cuboid offsets along the `y` axis.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;
@@ -690,6 +706,7 @@ impl<T: RealField> Orthographic3<T> {
 
     /// Sets the near and far plane offsets of the view cuboid.
     ///
+    /// # Example
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use nalgebra::Orthographic3;

src/geometry/perspective.rs

@@ -20,10 +20,12 @@ use crate::geometry::{Point3, Projective3};
 
 /// A 3D perspective projection stored as a homogeneous 4x4 matrix.
 #[repr(C)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Perspective3<T> {
     matrix: Matrix4<T>,

src/geometry/point.rs

@@ -3,15 +3,10 @@ use num::One;
 use std::cmp::Ordering;
 use std::fmt;
 use std::hash;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::simd::SimdPartialOrd;
 
 use crate::base::allocator::Allocator;
@@ -41,6 +36,11 @@ use std::mem::MaybeUninit;
 /// of said transformations for details.
 #[repr(C)]
 #[derive(Clone)]
+#[cfg_attr(
+    feature = "rkyv-serialize-no-std",
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
+)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 pub struct OPoint<T: Scalar, D: DimName>
 where
     DefaultAllocator: Allocator<T, D>,
@@ -74,12 +74,11 @@ where
 {
 }
 
-#[cfg(all(not(target_os = "cuda"), feature = "cuda"))]
+#[cfg(feature = "cuda")]
-unsafe impl<T: Scalar + cust::memory::DeviceCopy, D: DimName> cust::memory::DeviceCopy
+unsafe impl<T: Scalar + cust_core::DeviceCopy, D: DimName> cust_core::DeviceCopy for OPoint<T, D>
-    for OPoint<T, D>
 where
     DefaultAllocator: Allocator<T, D>,
-    OVector<T, D>: cust::memory::DeviceCopy,
+    OVector<T, D>: cust_core::DeviceCopy,
 {
 }
 
@@ -130,26 +129,6 @@ where
     }
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T, D: DimName> Abomonation for OPoint<T, D>
-where
-    T: Scalar,
-    OVector<T, D>: Abomonation,
-    DefaultAllocator: Allocator<T, D>,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.coords.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.coords.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.coords.exhume(bytes)
-    }
-}
-
 impl<T: Scalar, D: DimName> OPoint<T, D>
 where
     DefaultAllocator: Allocator<T, D>,
@@ -292,6 +271,7 @@ where
     /// assert_eq!(it.next(), Some(2.0));
     /// assert_eq!(it.next(), Some(3.0));
     /// assert_eq!(it.next(), None);
+    /// ```
     #[inline]
     pub fn iter(
         &self,
@@ -318,6 +298,7 @@ where
     /// }
     ///
     /// assert_eq!(p, Point3::new(10.0, 20.0, 30.0));
+    /// ```
     #[inline]
     pub fn iter_mut(
         &mut self,

src/geometry/quaternion.rs

@@ -2,17 +2,12 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use num::Zero;
 use std::fmt;
 use std::hash::{Hash, Hasher};
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use crate::base::storage::Owned;
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::{ClosedNeg, RealField};
 use simba::simd::{SimdBool, SimdOption, SimdRealField};
 
@@ -28,10 +23,12 @@ use crate::geometry::{Point3, Rotation};
 /// that may be used as a rotation.
 #[repr(C)]
 #[derive(Copy, Clone)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Quaternion<T> {
     /// This quaternion as a 4D vector of coordinates in the `[ x, y, z, w ]` storage order.
     pub coords: Vector4<T>,
@@ -77,24 +74,6 @@ where
 {
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T: Scalar> Abomonation for Quaternion<T>
-where
-    Vector4<T>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.coords.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.coords.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.coords.exhume(bytes)
-    }
-}
-
 #[cfg(feature = "serde-serialize-no-std")]
 impl<T: Scalar> Serialize for Quaternion<T>
 where
@@ -123,48 +102,6 @@ where
     }
 }
 
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Quaternion;
-    use crate::base::Vector4;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive> Archive for Quaternion<T> {
-        type Archived = Quaternion<T::Archived>;
-        type Resolver = <Vector4<T> as Archive>::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.coords.resolve(
-                pos + offset_of!(Self::Archived, coords),
-                resolver,
-                project_struct!(out: Self::Archived => coords),
-            );
-        }
-    }
-
-    impl<T: Serialize<S>, S: Fallible + ?Sized> Serialize<S> for Quaternion<T> {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            self.coords.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, D: Fallible + ?Sized> Deserialize<Quaternion<T>, D> for Quaternion<T::Archived>
-    where
-        T::Archived: Deserialize<T, D>,
-    {
-        fn deserialize(&self, deserializer: &mut D) -> Result<Quaternion<T>, D::Error> {
-            Ok(Quaternion {
-                coords: self.coords.deserialize(deserializer)?,
-            })
-        }
-    }
-}
-
 impl<T: SimdRealField> Quaternion<T>
 where
     T::Element: SimdRealField,
@@ -428,6 +365,7 @@ where
     /// let expected = Quaternion::new(-20.0, 0.0, 0.0, 0.0);
     /// let result = a.inner(&b);
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-5);
+    /// ```
     #[inline]
     #[must_use]
     pub fn inner(&self, other: &Self) -> Self {
@@ -1068,8 +1006,8 @@ impl<T: RealField + fmt::Display> fmt::Display for Quaternion<T> {
 /// A unit quaternions. May be used to represent a rotation.
 pub type UnitQuaternion<T> = Unit<Quaternion<T>>;
 
-#[cfg(all(not(target_os = "cuda"), feature = "cuda"))]
+#[cfg(feature = "cuda")]
-unsafe impl<T: cust::memory::DeviceCopy> cust::memory::DeviceCopy for UnitQuaternion<T> {}
+unsafe impl<T: cust_core::DeviceCopy> cust_core::DeviceCopy for UnitQuaternion<T> {}
 
 impl<T: Scalar + ClosedNeg + PartialEq> PartialEq for UnitQuaternion<T> {
     #[inline]
@@ -1253,8 +1191,7 @@ where
     /// Panics if the angle between both quaternion is 180 degrees (in which case the interpolation
     /// is not well-defined). Use `.try_slerp` instead to avoid the panic.
     ///
-    /// # Examples:
+    /// # Example
-    ///
     /// ```
     /// # use nalgebra::geometry::UnitQuaternion;
     ///
@@ -1476,7 +1413,6 @@ where
     /// Builds a rotation matrix from this unit quaternion.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1559,7 +1495,6 @@ where
     /// Converts this unit quaternion into its equivalent homogeneous transformation matrix.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1583,7 +1518,6 @@ where
     /// This is the same as the multiplication `self * pt`.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1604,7 +1538,6 @@ where
     /// This is the same as the multiplication `self * v`.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1625,7 +1558,6 @@ where
     /// point.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1648,7 +1580,6 @@ where
     /// vector.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -1669,7 +1600,6 @@ where
     /// vector.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;

src/geometry/rotation.rs

@@ -2,8 +2,6 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use num::{One, Zero};
 use std::fmt;
 use std::hash;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
@@ -11,9 +9,6 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer};
 #[cfg(feature = "serde-serialize-no-std")]
 use crate::base::storage::Owned;
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::RealField;
 use simba::simd::SimdRealField;
 
@@ -54,10 +49,12 @@ use crate::geometry::Point;
 /// * [Conversion to a matrix <span style="float:right;">`matrix`, `to_homogeneous`…</span>](#conversion-to-a-matrix)
 ///
 #[repr(C)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Rotation<T, const D: usize> {
     matrix: SMatrix<T, D, D>,
@@ -94,25 +91,6 @@ where
 {
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T, const D: usize> Abomonation for Rotation<T, D>
-where
-    T: Scalar,
-    SMatrix<T, D, D>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.matrix.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.matrix.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.matrix.exhume(bytes)
-    }
-}
-
 #[cfg(feature = "serde-serialize-no-std")]
 impl<T: Scalar, const D: usize> Serialize for Rotation<T, D>
 where

src/geometry/scale.rs

@@ -2,15 +2,10 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use num::{One, Zero};
 use std::fmt;
 use std::hash;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use crate::base::allocator::Allocator;
 use crate::base::dimension::{DimNameAdd, DimNameSum, U1};
 use crate::base::storage::Owned;
@@ -22,10 +17,12 @@ use crate::geometry::Point;
 
 /// A scale which supports non-uniform scaling.
 #[repr(C)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Scale<T, const D: usize> {
     /// The scale coordinates, i.e., how much is multiplied to a point's coordinates when it is
@@ -64,25 +61,6 @@ where
 {
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T, const D: usize> Abomonation for Scale<T, D>
-where
-    T: Scalar,
-    SVector<T, D>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.vector.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.vector.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.vector.exhume(bytes)
-    }
-}
-
 #[cfg(feature = "serde-serialize-no-std")]
 impl<T: Scalar, const D: usize> Serialize for Scale<T, D>
 where
@@ -111,49 +89,6 @@ where
     }
 }
 
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Scale;
-    use crate::base::SVector;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive, const D: usize> Archive for Scale<T, D> {
-        type Archived = Scale<T::Archived, D>;
-        type Resolver = <SVector<T, D> as Archive>::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.vector.resolve(
-                pos + offset_of!(Self::Archived, vector),
-                resolver,
-                project_struct!(out: Self::Archived => vector),
-            );
-        }
-    }
-
-    impl<T: Serialize<S>, S: Fallible + ?Sized, const D: usize> Serialize<S> for Scale<T, D> {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            self.vector.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, _D: Fallible + ?Sized, const D: usize> Deserialize<Scale<T, D>, _D>
-        for Scale<T::Archived, D>
-    where
-        T::Archived: Deserialize<T, _D>,
-    {
-        fn deserialize(&self, deserializer: &mut _D) -> Result<Scale<T, D>, _D::Error> {
-            Ok(Scale {
-                vector: self.vector.deserialize(deserializer)?,
-            })
-        }
-    }
-}
-
 impl<T: Scalar, const D: usize> Scale<T, D> {
     /// Inverts `self`.
     ///

src/geometry/similarity.rs

@@ -3,15 +3,9 @@ use num::Zero;
 use std::fmt;
 use std::hash;
 
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
-
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Serialize};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::{RealField, SubsetOf};
 use simba::simd::SimdRealField;
 
@@ -24,10 +18,7 @@ use crate::geometry::{AbstractRotation, Isometry, Point, Translation};
 /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation.
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
-#[cfg_attr(
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
-    all(not(target_os = "cuda"), feature = "cuda"),
-    derive(cust::DeviceCopy)
-)]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
     feature = "serde-serialize-no-std",
@@ -43,30 +34,17 @@ use crate::geometry::{AbstractRotation, Isometry, Point, Translation};
                        DefaultAllocator: Allocator<T, Const<D>>,
                        Owned<T, Const<D>>: Deserialize<'de>"))
 )]
+#[cfg_attr(
+    feature = "rkyv-serialize-no-std",
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
+)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 pub struct Similarity<T, R, const D: usize> {
     /// The part of this similarity that does not include the scaling factor.
     pub isometry: Isometry<T, R, D>,
     scaling: T,
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T: Scalar, R, const D: usize> Abomonation for Similarity<T, R, D>
-where
-    Isometry<T, R, D>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.isometry.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.isometry.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.isometry.exhume(bytes)
-    }
-}
-
 impl<T: Scalar + hash::Hash, R: hash::Hash, const D: usize> hash::Hash for Similarity<T, R, D>
 where
     Owned<T, Const<D>>: hash::Hash,

src/geometry/similarity_construction.rs

@@ -38,7 +38,6 @@ where
     /// Creates a new identity similarity.
     ///
     /// # Example
-    ///
     /// ```
     /// # use nalgebra::{Similarity2, Point2, Similarity3, Point3};
     ///
@@ -95,7 +94,6 @@ where
     /// its axis passing through the point `p`.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -146,7 +144,6 @@ where
     /// Creates a new similarity from a translation, a rotation, and an uniform scaling factor.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -188,7 +185,6 @@ where
     /// Creates a new similarity from a translation and a rotation angle.
     ///
     /// # Example
-    ///
     /// ```
     /// # #[macro_use] extern crate approx;
     /// # use std::f32;
@@ -232,7 +228,6 @@ macro_rules! similarity_construction_impl(
             /// factor.
             ///
             /// # Example
-            ///
             /// ```
             /// # #[macro_use] extern crate approx;
             /// # use std::f32;
@@ -288,7 +283,6 @@ macro_rules! similarity_construction_impl(
             ///   to `eye - at`. Non-collinearity is not checked.
             ///
             /// # Example
-            ///
             /// ```
             /// # #[macro_use] extern crate approx;
             /// # use std::f32;
@@ -316,7 +310,7 @@ macro_rules! similarity_construction_impl(
                 Self::from_isometry(Isometry::<_, $Rot<T>, 3>::face_towards(eye, target, up), scaling)
             }
 
-            /// Deprecated: Use [`SimilarityMatrix3::face_towards`] instead.
+            /// Deprecated: Use [`SimilarityMatrix3::face_towards`](Self::face_towards) instead.
             #[deprecated(note="renamed to `face_towards`")]
             pub fn new_observer_frames(eye:    &Point3<T>,
                                        target: &Point3<T>,
@@ -338,7 +332,6 @@ macro_rules! similarity_construction_impl(
             ///   requirement of this parameter is to not be collinear to `target - eye`.
             ///
             /// # Example
-            ///
             /// ```
             /// # #[macro_use] extern crate approx;
             /// # use std::f32;
@@ -376,7 +369,6 @@ macro_rules! similarity_construction_impl(
             ///   requirement of this parameter is to not be collinear to `target - eye`.
             ///
             /// # Example
-            ///
             /// ```
             /// # #[macro_use] extern crate approx;
             /// # use std::f32;

src/geometry/transform.rs

@@ -60,26 +60,17 @@ where
 
 /// Tag representing the most general (not necessarily inversible) `Transform` type.
 #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-#[cfg_attr(
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
-    all(not(target_os = "cuda"), feature = "cuda"),
-    derive(cust::DeviceCopy)
-)]
 pub enum TGeneral {}
 
 /// Tag representing the most general inversible `Transform` type.
 #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-#[cfg_attr(
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
-    all(not(target_os = "cuda"), feature = "cuda"),
-    derive(cust::DeviceCopy)
-)]
 pub enum TProjective {}
 
 /// Tag representing an affine `Transform`. Its bottom-row is equal to `(0, 0 ... 0, 1)`.
 #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-#[cfg_attr(
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
-    all(not(target_os = "cuda"), feature = "cuda"),
-    derive(cust::DeviceCopy)
-)]
 pub enum TAffine {}
 
 impl TCategory for TGeneral {
@@ -207,13 +198,13 @@ where
 {
 }
 
-#[cfg(all(not(target_os = "cuda"), feature = "cuda"))]
+#[cfg(feature = "cuda")]
-unsafe impl<T: RealField + cust::memory::DeviceCopy, C: TCategory, const D: usize>
+unsafe impl<T: RealField + cust_core::DeviceCopy, C: TCategory, const D: usize>
-    cust::memory::DeviceCopy for Transform<T, C, D>
+    cust_core::DeviceCopy for Transform<T, C, D>
 where
     Const<D>: DimNameAdd<U1>,
     DefaultAllocator: Allocator<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>,
-    Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>: cust::memory::DeviceCopy,
+    Owned<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>: cust_core::DeviceCopy,
 {
 }
 

src/geometry/translation.rs

@@ -2,15 +2,10 @@ use approx::{AbsDiffEq, RelativeEq, UlpsEq};
 use num::{One, Zero};
 use std::fmt;
 use std::hash;
-#[cfg(feature = "abomonation-serialize")]
-use std::io::{Result as IOResult, Write};
 
 #[cfg(feature = "serde-serialize-no-std")]
 use serde::{Deserialize, Deserializer, Serialize, Serializer};
 
-#[cfg(feature = "abomonation-serialize")]
-use abomonation::Abomonation;
-
 use simba::scalar::{ClosedAdd, ClosedNeg, ClosedSub};
 
 use crate::base::allocator::Allocator;
@@ -22,10 +17,12 @@ use crate::geometry::Point;
 
 /// A translation.
 #[repr(C)]
+#[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
 #[cfg_attr(
-    all(not(target_os = "cuda"), feature = "cuda"),
+    feature = "rkyv-serialize-no-std",
-    derive(cust::DeviceCopy)
+    derive(rkyv::Archive, rkyv::Serialize, rkyv::Deserialize)
 )]
+#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Translation<T, const D: usize> {
     /// The translation coordinates, i.e., how much is added to a point's coordinates when it is
@@ -64,25 +61,6 @@ where
 {
 }
 
-#[cfg(feature = "abomonation-serialize")]
-impl<T, const D: usize> Abomonation for Translation<T, D>
-where
-    T: Scalar,
-    SVector<T, D>: Abomonation,
-{
-    unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> {
-        self.vector.entomb(writer)
-    }
-
-    fn extent(&self) -> usize {
-        self.vector.extent()
-    }
-
-    unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> {
-        self.vector.exhume(bytes)
-    }
-}
-
 #[cfg(feature = "serde-serialize-no-std")]
 impl<T: Scalar, const D: usize> Serialize for Translation<T, D>
 where
@@ -111,49 +89,6 @@ where
     }
 }
 
-#[cfg(feature = "rkyv-serialize-no-std")]
-mod rkyv_impl {
-    use super::Translation;
-    use crate::base::SVector;
-    use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize};
-
-    impl<T: Archive, const D: usize> Archive for Translation<T, D> {
-        type Archived = Translation<T::Archived, D>;
-        type Resolver = <SVector<T, D> as Archive>::Resolver;
-
-        fn resolve(
-            &self,
-            pos: usize,
-            resolver: Self::Resolver,
-            out: &mut core::mem::MaybeUninit<Self::Archived>,
-        ) {
-            self.vector.resolve(
-                pos + offset_of!(Self::Archived, vector),
-                resolver,
-                project_struct!(out: Self::Archived => vector),
-            );
-        }
-    }
-
-    impl<T: Serialize<S>, S: Fallible + ?Sized, const D: usize> Serialize<S> for Translation<T, D> {
-        fn serialize(&self, serializer: &mut S) -> Result<Self::Resolver, S::Error> {
-            self.vector.serialize(serializer)
-        }
-    }
-
-    impl<T: Archive, _D: Fallible + ?Sized, const D: usize> Deserialize<Translation<T, D>, _D>
-        for Translation<T::Archived, D>
-    where
-        T::Archived: Deserialize<T, _D>,
-    {
-        fn deserialize(&self, deserializer: &mut _D) -> Result<Translation<T, D>, _D::Error> {
-            Ok(Translation {
-                vector: self.vector.deserialize(deserializer)?,
-            })
-        }
-    }
-}
-
 impl<T: Scalar, const D: usize> Translation<T, D> {
     /// Creates a new translation from the given vector.
     #[inline]
@@ -255,6 +190,7 @@ impl<T: Scalar + ClosedAdd, const D: usize> Translation<T, D> {
     /// let t = Translation3::new(1.0, 2.0, 3.0);
     /// let transformed_point = t.transform_point(&Point3::new(4.0, 5.0, 6.0));
     /// assert_eq!(transformed_point, Point3::new(5.0, 7.0, 9.0));
+    /// ```
     #[inline]
     #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
@@ -271,6 +207,7 @@ impl<T: Scalar + ClosedSub, const D: usize> Translation<T, D> {
     /// let t = Translation3::new(1.0, 2.0, 3.0);
     /// let transformed_point = t.inverse_transform_point(&Point3::new(4.0, 5.0, 6.0));
     /// assert_eq!(transformed_point, Point3::new(3.0, 3.0, 3.0));
+    /// ```
     #[inline]
     #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {

src/geometry/unit_complex.rs

@@ -31,8 +31,8 @@ use std::cmp::{Eq, PartialEq};
 /// * [Conversion to a matrix <span style="float:right;">`to_rotation_matrix`, `to_homogeneous`…</span>](#conversion-to-a-matrix)
 pub type UnitComplex<T> = Unit<Complex<T>>;
 
-#[cfg(all(not(target_os = "cuda"), feature = "cuda"))]
+#[cfg(feature = "cuda")]
-unsafe impl<T: cust::memory::DeviceCopy> cust::memory::DeviceCopy for UnitComplex<T> {}
+unsafe impl<T: cust_core::DeviceCopy> cust_core::DeviceCopy for UnitComplex<T> {}
 
 impl<T: Scalar + PartialEq> PartialEq for UnitComplex<T> {
     #[inline]
@@ -410,7 +410,8 @@ where
     #[inline]
     #[must_use]
     pub fn slerp(&self, other: &Self, t: T) -> Self {
-        Self::new(self.angle() * (T::one() - t.clone()) + other.angle() * t)
+        let delta = other / self;
+        self * Self::new(delta.angle() * t)
     }
 }
 

src/lib.rs

@@ -78,12 +78,13 @@ an optimized set of tools for computer graphics and physics. Those features incl
     unused_mut,
     unused_parens,
     unused_qualifications,
-    unused_results,
     rust_2018_idioms,
     rust_2018_compatibility,
     future_incompatible,
     missing_copy_implementations
 )]
+#![cfg_attr(feature = "rkyv-serialize-no-std", warn(unused_results))] // TODO: deny this once bytecheck stops generating warnings.
+#![cfg_attr(not(feature = "rkyv-serialize-no-std"), deny(unused_results))]
 #![doc(
     html_favicon_url = "https://nalgebra.org/img/favicon.ico",
     html_root_url = "https://docs.rs/nalgebra/0.25.0"
@@ -246,8 +247,8 @@ pub fn min<T: Ord>(a: T, b: T) -> T {
 
 /// The absolute value of `a`.
 ///
-/// Deprecated: Use [`Matrix::abs`] or [`RealField::abs`] instead.
+/// Deprecated: Use [`Matrix::abs`] or [`ComplexField::abs`] instead.
-#[deprecated(note = "use the inherent method `Matrix::abs` or `RealField::abs` instead")]
+#[deprecated(note = "use the inherent method `Matrix::abs` or `ComplexField::abs` instead")]
 #[inline]
 pub fn abs<T: Signed>(a: &T) -> T {
     a.abs()

src/linalg/svd.rs

@@ -98,7 +98,7 @@ where
             matrix,
             compute_u,
             compute_v,
-            T::RealField::default_epsilon(),
+            T::RealField::default_epsilon() * crate::convert(5.0),
             0,
         )
         .unwrap()
@@ -888,13 +888,13 @@ fn compute_2x2_uptrig_svd<T: RealField>(
             v_t = Some(csv.clone());
         }
 
-        if compute_u {
+        let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone();
-            let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone();
+        let su = (m22 * csv.s()) / v1.clone();
-            let su = (m22 * csv.s()) / v1.clone();
+        let (csu, sgn_u) = GivensRotation::new(cu, su);
-            let (csu, sgn_u) = GivensRotation::new(cu, su);
+        v1 *= sgn_u.clone();
+        v2 *= sgn_u;
 
-            v1 *= sgn_u.clone();
+        if compute_u {
-            v2 *= sgn_u;
             u = Some(csu);
         }
     }

src/third_party/glam/common/glam_isometry.rs

@@ -1,5 +1,5 @@
 use super::glam::{DMat3, DMat4, DQuat, DVec2, DVec3, Mat3, Mat4, Quat, Vec2, Vec3};
-use crate::{Isometry2, Isometry3, Matrix3, Matrix4, Translation3, UnitQuaternion, Vector2};
+use crate::{Isometry2, Isometry3, Matrix3, Matrix4, Translation3, UnitQuaternion};
 use std::convert::TryFrom;
 
 impl From<Isometry2<f32>> for Mat3 {
@@ -36,18 +36,18 @@ impl From<Isometry3<f64>> for (DVec3, DQuat) {
     }
 }
 
-impl From<Isometry2<f32>> for (Vec3, Quat) {
+impl From<Isometry2<f32>> for (Vec2, f32) {
-    fn from(iso: Isometry2<f32>) -> (Vec3, Quat) {
+    fn from(iso: Isometry2<f32>) -> (Vec2, f32) {
-        let tra = Vec3::new(iso.translation.x, iso.translation.y, 0.0);
+        let tra = Vec2::new(iso.translation.x, iso.translation.y);
-        let rot = Quat::from_axis_angle(Vec3::Z, iso.rotation.angle());
+        let rot = iso.rotation.angle();
         (tra, rot)
     }
 }
 
-impl From<Isometry2<f64>> for (DVec3, DQuat) {
+impl From<Isometry2<f64>> for (DVec2, f64) {
-    fn from(iso: Isometry2<f64>) -> (DVec3, DQuat) {
+    fn from(iso: Isometry2<f64>) -> (DVec2, f64) {
-        let tra = DVec3::new(iso.translation.x, iso.translation.y, 0.0);
+        let tra = DVec2::new(iso.translation.x, iso.translation.y);
-        let rot = DQuat::from_axis_angle(DVec3::Z, iso.rotation.angle());
+        let rot = iso.rotation.angle();
         (tra, rot)
     }
 }
@@ -64,30 +64,6 @@ impl From<(DVec3, DQuat)> for Isometry3<f64> {
     }
 }
 
-impl From<(Vec3, Quat)> for Isometry2<f32> {
-    fn from((tra, rot): (Vec3, Quat)) -> Self {
-        Isometry2::new([tra.x, tra.y].into(), rot.to_axis_angle().1)
-    }
-}
-
-impl From<(DVec3, DQuat)> for Isometry2<f64> {
-    fn from((tra, rot): (DVec3, DQuat)) -> Self {
-        Isometry2::new([tra.x, tra.y].into(), rot.to_axis_angle().1)
-    }
-}
-
-impl From<(Vec2, Quat)> for Isometry2<f32> {
-    fn from((tra, rot): (Vec2, Quat)) -> Self {
-        Isometry2::new(tra.into(), rot.to_axis_angle().1)
-    }
-}
-
-impl From<(DVec2, DQuat)> for Isometry2<f64> {
-    fn from((tra, rot): (DVec2, DQuat)) -> Self {
-        Isometry2::new(tra.into(), rot.to_axis_angle().1)
-    }
-}
-
 impl From<(Vec2, f32)> for Isometry2<f32> {
     fn from((tra, rot): (Vec2, f32)) -> Self {
         Isometry2::new(tra.into(), rot)
@@ -112,18 +88,6 @@ impl From<DQuat> for Isometry3<f64> {
     }
 }
 
-impl From<Quat> for Isometry2<f32> {
-    fn from(rot: Quat) -> Self {
-        Isometry2::new(Vector2::zeros(), rot.to_axis_angle().1)
-    }
-}
-
-impl From<DQuat> for Isometry2<f64> {
-    fn from(rot: DQuat) -> Self {
-        Isometry2::new(Vector2::zeros(), rot.to_axis_angle().1)
-    }
-}
-
 impl From<Vec3> for Isometry3<f32> {
     fn from(tra: Vec3) -> Self {
         Isometry3::from_parts(tra.into(), UnitQuaternion::identity())
@@ -148,18 +112,6 @@ impl From<DVec2> for Isometry2<f64> {
     }
 }
 
-impl From<Vec3> for Isometry2<f32> {
-    fn from(tra: Vec3) -> Self {
-        Isometry2::new([tra.x, tra.y].into(), 0.0)
-    }
-}
-
-impl From<DVec3> for Isometry2<f64> {
-    fn from(tra: DVec3) -> Self {
-        Isometry2::new([tra.x, tra.y].into(), 0.0)
-    }
-}
-
 impl TryFrom<Mat3> for Isometry2<f32> {
     type Error = ();
 

src/third_party/glam/common/glam_matrix.rs

@@ -3,7 +3,11 @@ use super::glam::{
     Mat4, UVec2, UVec3, UVec4, Vec2, Vec3, Vec3A, Vec4,
 };
 use crate::storage::RawStorage;
-use crate::{Matrix, Matrix2, Matrix3, Matrix4, Vector, Vector2, Vector3, Vector4, U2, U3, U4};
+use crate::{
+    Matrix, Matrix2, Matrix3, Matrix4, Unit, UnitVector2, UnitVector3, UnitVector4, Vector,
+    Vector2, Vector3, Vector4, U2, U3, U4,
+};
+use std::convert::TryFrom;
 
 macro_rules! impl_vec_conversion(
     ($N: ty, $Vec2: ty, $Vec3: ty, $Vec4: ty) => {
@@ -66,6 +70,63 @@ impl_vec_conversion!(i32, IVec2, IVec3, IVec4);
 impl_vec_conversion!(u32, UVec2, UVec3, UVec4);
 impl_vec_conversion!(bool, BVec2, BVec3, BVec4);
 
+const ERR: &'static str = "Normalization failed.";
+
+macro_rules! impl_unit_vec_conversion(
+    ($N: ty, $Vec2: ty, $Vec3: ty, $Vec4: ty) => {
+        impl TryFrom<$Vec2> for UnitVector2<$N> {
+            type Error = &'static str;
+            #[inline]
+            fn try_from(e: $Vec2) -> Result<Self, Self::Error> {
+                Unit::try_new(e.into(), 0.0).ok_or(ERR)
+            }
+        }
+
+        impl From<UnitVector2<$N>> for $Vec2
+        {
+            #[inline]
+            fn from(e: UnitVector2<$N>) -> $Vec2 {
+                e.into_inner().into()
+            }
+        }
+
+        impl TryFrom<$Vec3> for UnitVector3<$N> {
+            type Error = &'static str;
+            #[inline]
+            fn try_from(e: $Vec3) -> Result<Self, Self::Error> {
+                Unit::try_new(e.into(), 0.0).ok_or(ERR)
+            }
+        }
+
+        impl From<UnitVector3<$N>> for $Vec3
+        {
+            #[inline]
+            fn from(e: UnitVector3<$N>) -> $Vec3 {
+                e.into_inner().into()
+            }
+        }
+
+        impl TryFrom<$Vec4> for UnitVector4<$N> {
+            type Error = &'static str;
+            #[inline]
+            fn try_from(e: $Vec4) -> Result<Self, Self::Error> {
+                Unit::try_new(e.into(), 0.0).ok_or(ERR)
+            }
+        }
+
+        impl From<UnitVector4<$N>> for $Vec4
+        {
+            #[inline]
+            fn from(e: UnitVector4<$N>) -> $Vec4 {
+                e.into_inner().into()
+            }
+        }
+    }
+);
+
+impl_unit_vec_conversion!(f32, Vec2, Vec3, Vec4);
+impl_unit_vec_conversion!(f64, DVec2, DVec3, DVec4);
+
 impl From<Vec3A> for Vector3<f32> {
     #[inline]
     fn from(e: Vec3A) -> Vector3<f32> {
@@ -83,6 +144,21 @@ where
     }
 }
 
+impl TryFrom<Vec3A> for UnitVector3<f32> {
+    type Error = &'static str;
+    #[inline]
+    fn try_from(e: Vec3A) -> Result<Self, Self::Error> {
+        Unit::try_new(e.into(), 0.0).ok_or(ERR)
+    }
+}
+
+impl From<UnitVector3<f32>> for Vec3A {
+    #[inline]
+    fn from(e: UnitVector3<f32>) -> Vec3A {
+        e.into_inner().into()
+    }
+}
+
 impl From<Mat2> for Matrix2<f32> {
     #[inline]
     fn from(e: Mat2) -> Matrix2<f32> {

src/third_party/glam/mod.rs

@@ -1,5 +1,3 @@
-#[cfg(feature = "glam013")]
-mod v013;
 #[cfg(feature = "glam014")]
 mod v014;
 #[cfg(feature = "glam015")]
@@ -10,3 +8,7 @@ mod v016;
 mod v017;
 #[cfg(feature = "glam018")]
 mod v018;
+#[cfg(feature = "glam019")]
+mod v019;
+#[cfg(feature = "glam020")]
+mod v020;

src/third_party/glam/v019/mod.rs

@@ -15,4 +15,4 @@ mod glam_translation;
 #[path = "../common/glam_unit_complex.rs"]
 mod glam_unit_complex;
 
-pub(self) use glam013 as glam;
+pub(self) use glam019 as glam;

src/third_party/glam/v020/mod.rs

@@ -0,0 +1,18 @@
+#[path = "../common/glam_isometry.rs"]
+mod glam_isometry;
+#[path = "../common/glam_matrix.rs"]
+mod glam_matrix;
+#[path = "../common/glam_point.rs"]
+mod glam_point;
+#[path = "../common/glam_quaternion.rs"]
+mod glam_quaternion;
+#[path = "../common/glam_rotation.rs"]
+mod glam_rotation;
+#[path = "../common/glam_similarity.rs"]
+mod glam_similarity;
+#[path = "../common/glam_translation.rs"]
+mod glam_translation;
+#[path = "../common/glam_unit_complex.rs"]
+mod glam_unit_complex;
+
+pub(self) use glam020 as glam;

tests/core/abomonation.rs

@@ -1,51 +0,0 @@
-use abomonation::{decode, encode, Abomonation};
-use na::{
-    DMatrix, Isometry3, IsometryMatrix3, Matrix3x4, Point3, Quaternion, Rotation3, Similarity3,
-    SimilarityMatrix3, Translation3,
-};
-use rand::random;
-
-#[test]
-fn abomonate_dmatrix() {
-    assert_encode_and_decode(DMatrix::<f32>::new_random(3, 5));
-}
-
-macro_rules! test_abomonation(
-    ($($test: ident, $ty: ty);* $(;)*) => {$(
-        #[test]
-        fn $test() {
-            assert_encode_and_decode(random::<$ty>());
-        }
-    )*}
-);
-
-test_abomonation! {
-    abomonate_matrix3x4, Matrix3x4<f32>;
-    abomonate_point3, Point3<f32>;
-    abomonate_translation3, Translation3<f64>;
-    abomonate_rotation3, Rotation3<f64>;
-    abomonate_isometry3, Isometry3<f32>;
-    abomonate_isometry_matrix3, IsometryMatrix3<f64>;
-    abomonate_similarity3, Similarity3<f32>;
-    abomonate_similarity_matrix3, SimilarityMatrix3<f32>;
-    abomonate_quaternion, Quaternion<f32>;
-}
-
-fn assert_encode_and_decode<T: Abomonation + PartialEq + Clone>(original_data: T) {
-    // Hold on to a clone for later comparison
-    let data = original_data.clone();
-
-    // Encode
-    let mut bytes = Vec::new();
-    unsafe {
-        let _ = encode(&original_data, &mut bytes);
-    }
-
-    // Drop the original, so that dangling pointers are revealed by the test
-    drop(original_data);
-
-    if let Some((result, rest)) = unsafe { decode::<T>(&mut bytes) } {
-        assert!(result == &data);
-        assert!(rest.len() == 0, "binary data was not decoded completely");
-    }
-}

tests/core/mod.rs

@@ -1,5 +1,3 @@
-#[cfg(feature = "abomonation-serialize")]
-mod abomonation;
 mod blas;
 mod cg;
 mod conversion;

tests/geometry/rotation.rs

@@ -32,7 +32,9 @@ fn quaternion_euler_angles_issue_494() {
 
 #[cfg(feature = "proptest-support")]
 mod proptest_tests {
+    use approx::AbsDiffEq;
     use na::{self, Rotation2, Rotation3, Unit};
+    use na::{UnitComplex, UnitQuaternion};
     use simba::scalar::RealField;
     use std::f64;
 
@@ -229,5 +231,74 @@ mod proptest_tests {
                 prop_assert_eq!(r, Rotation3::identity())
             }
         }
+
+        //
+        //In general, `slerp(a,b,t)` should equal `(b/a)^t * a` even though in practice,
+        //we may not use that formula directly for complex numbers or quaternions
+        //
+
+        #[test]
+        fn slerp_powf_agree_2(a in unit_complex(), b in unit_complex(), t in PROPTEST_F64) {
+            let z1 = a.slerp(&b, t);
+            let z2 = (b/a).powf(t) * a;
+            prop_assert!(relative_eq!(z1,z2,epsilon=1e-10));
+        }
+
+        #[test]
+        fn slerp_powf_agree_3(a in unit_quaternion(), b in unit_quaternion(), t in PROPTEST_F64) {
+            if let Some(z1) = a.try_slerp(&b, t, f64::default_epsilon()) {
+                let z2 = (b/a).powf(t) * a;
+                prop_assert!(relative_eq!(z1,z2,epsilon=1e-10));
+            }
+        }
+
+        //
+        //when not antipodal, slerp should always take the shortest path between two orientations
+        //
+
+        #[test]
+        fn slerp_takes_shortest_path_2(
+            z in unit_complex(), dtheta in -f64::pi()..f64::pi(), t in 0.0..1.0f64
+        ) {
+
+            //ambiguous when at ends of angle range, so we don't really care here
+            if dtheta.abs() != f64::pi() {
+
+                //make two complex numbers separated by an angle between -pi and pi
+                let (z1, z2) = (z, z * UnitComplex::new(dtheta));
+                let z3 = z1.slerp(&z2, t);
+
+                //since the angle is no larger than a half-turn, and t is between 0 and 1,
+                //the shortest path just corresponds to adding the scaled angle
+                let a1 = z3.angle();
+                let a2 = na::wrap(z1.angle() + dtheta*t, -f64::pi(), f64::pi());
+
+                prop_assert!(relative_eq!(a1, a2, epsilon=1e-10));
+            }
+
+        }
+
+        #[test]
+        fn slerp_takes_shortest_path_3(
+            q in unit_quaternion(), dtheta in -f64::pi()..f64::pi(), t in 0.0..1.0f64
+        ) {
+
+            //ambiguous when at ends of angle range, so we don't really care here
+            if let Some(axis) = q.axis() {
+
+                //make two quaternions separated by an angle between -pi and pi
+                let (q1, q2) = (q, q * UnitQuaternion::from_axis_angle(&axis, dtheta));
+                let q3 = q1.slerp(&q2, t);
+
+                //since the angle is no larger than a half-turn, and t is between 0 and 1,
+                //the shortest path just corresponds to adding the scaled angle
+                let q4 = q1 * UnitQuaternion::from_axis_angle(&axis, dtheta*t);
+                prop_assert!(relative_eq!(q3, q4, epsilon=1e-10));
+
+            }
+
+        }
+
+
     }
 }

tests/lib.rs

@@ -9,8 +9,6 @@ compile_error!(
      Example: `cargo test --features debug,compare,rand,macros`"
 );
 
-#[cfg(feature = "abomonation-serialize")]
-extern crate abomonation;
 #[cfg(all(feature = "debug", feature = "compare", feature = "rand"))]
 #[macro_use]
 extern crate approx;

tests/linalg/svd.rs

@@ -460,3 +460,42 @@ fn svd_sorted() {
         epsilon = 1.0e-5
     );
 }
+
+#[test]
+// Exercises bug reported in issue #983 of nalgebra (https://github.com/dimforge/nalgebra/issues/983)
+fn svd_regression_issue_983() {
+    let m = nalgebra::dmatrix![
+        10.74785316637712f64, -5.994983325167452, -6.064492921857296;
+        -4.149751381521569, 20.654504205822462, -4.470436210703133;
+        -22.772715014220207, -1.4554372570788008, 18.108113992170573
+    ]
+    .transpose();
+    let svd1 = m.clone().svd(true, true);
+    let svd2 = m.clone().svd(false, true);
+    let svd3 = m.clone().svd(true, false);
+    let svd4 = m.svd(false, false);
+
+    assert_relative_eq!(svd1.singular_values, svd2.singular_values, epsilon = 1e-9);
+    assert_relative_eq!(svd1.singular_values, svd3.singular_values, epsilon = 1e-9);
+    assert_relative_eq!(svd1.singular_values, svd4.singular_values, epsilon = 1e-9);
+    assert_relative_eq!(
+        svd1.singular_values,
+        nalgebra::dvector![3.16188022e+01, 2.23811978e+01, 0.],
+        epsilon = 1e-6
+    );
+}
+
+#[test]
+// Exercises bug reported in issue #1072 of nalgebra (https://github.com/dimforge/nalgebra/issues/1072)
+fn svd_regression_issue_1072() {
+    let x = nalgebra::dmatrix![-6.206610118536945f64, -3.67612186839874; -1.2755730783423473, 6.047238193479124];
+    let mut x_svd = x.svd(true, true);
+    x_svd.singular_values = nalgebra::dvector![1.0, 0.0];
+    let y = x_svd.recompose().unwrap();
+    let y_svd = y.svd(true, true);
+    assert_relative_eq!(
+        y_svd.singular_values,
+        nalgebra::dvector![1.0, 0.0],
+        epsilon = 1e-9
+    );
+}