Merge branch 'dimforge:dev' into dev

2024-05-04 19:35:48 +03:00 · 2024-05-04 19:35:48 +03:00 · 02d66f28d0
parent d65b457410 9bc7f8b0d8
commit 02d66f28d0
47 changed files with 786 additions and 519 deletions
--- a/.github/workflows/nalgebra-ci-build.yml
+++ b/.github/workflows/nalgebra-ci-build.yml
@ -49,8 +49,6 @@ jobs:
  build-nalgebra-all-features:
    runs-on: ubuntu-latest
    steps:
-      # Needed because the --all-features build which enables cuda support.
-      - uses: Jimver/cuda-toolkit@v0.2.8
      - uses: actions/checkout@v2
      - run: cargo build --all-features;
      - run: cargo build -p nalgebra-glm --all-features;
@ -120,23 +118,6 @@ jobs:
        run: xargo build --verbose --no-default-features -p nalgebra-glm --target=x86_64-unknown-linux-gnu;
      - name: build thumbv7em-none-eabihf nalgebra-glm
        run: xargo build --verbose --no-default-features -p nalgebra-glm --target=thumbv7em-none-eabihf;
-  build-cuda:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: Jimver/cuda-toolkit@v0.2.8
-        with:
-          cuda: '11.5.0'
-      - name: Install nightly-2021-12-04
-        uses: actions-rs/toolchain@v1
-        with:
-          toolchain: nightly-2021-12-04
-          override: true
-      - 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
-        env: 
-          CUDA_ARCH: "350"
  docs:
    runs-on: ubuntu-latest
    steps:
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@ -1,40 +1,67 @@
 # Change Log
+
 All notable changes to `nalgebra`, starting with the version 0.6.0 will be
 documented here.

 This project adheres to [Semantic Versioning](https://semver.org/).

+## Unreleased
+
+### Removed
+
+- The `cuda` feature has been removed, as the toolchain it depends on
+  is long abandoned.
+
+## [0.32.5] (28 March 2024)
+
+## Fixed
+
+- Fix numerical issue on SVD with near-identity matrix.
+
+## [0.32.4] (19 Feb 2023)
+
+### Added
+
+- Add the `glam-0.25` feature to enable conversion from/to types from `glam` v0.25.
+
 ## [0.32.3] (09 July 2023)

 ### Modified
+
 - Statically sized matrices are now serialized as tuples to match how serde
  serialized plain arrays.
 - Don’t require `Scalar` for matrix `PartialEq` and `Eq`.

 ### Added
+
 - Allow trailing punctuation in macros `vector!`, `matrix!`, `point!`, etc.
 - Add the methods `Matrix1::as_scalar`, `::as_scalar_mut`, `::to_scalar`, `::into_scalar`.
 - Add `Rotation3::euler_angles_ordered`, a generalized euler angles calculation.
 - Add the `glam-0.24` feature to enable conversion from/to types from `glam` v0.24.
+- Add the `glam-0.25` feature to enable conversion from/to types from `glam` v0.25.
 - Add the `lerp` method to points.
 - Implement `Clone` for `MatrixIter`.

 ### Fixed
+
 - Fixed severe catastrophic cancellation issue in variance calculation.

 ## [0.32.2] (07 March 2023)

 ### Added
+
 - Add the `glam-0.23` to enable conversion from/to type from `glam` v0.23.

 ## [0.32.1] (14 Jan. 2023)

 ### Modified
+
 - Updated `nalgebra-macros` to use the new `Dyn`, avoiding macro-generated deprecation warnings.

 ## [0.32.0] (14 Jan. 2023)

 ### Modified
+
 - Renamed all `MatrixSlice` types to `MatrixView`. In general all occurrences of the world `Slice` or `slice` have been
  replaced by `View` or `view`.
 - Deprecated all the types involving `Slice` in its name, in favor of the word `View`.
@ -42,10 +69,12 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - Renamed `Dynamic` to `Dyn` and make `Dyn` a tuple struct.

 ### Added
+
 - Add `Cholesky::ln_determinant` to compute the natural logarithm of the determinant of a matrix decomposed
  with Cholesky. This can be more numerically stable than computing the determinant itself when very small and/or
  large values are involved.
- Added new methods `Matrix::as_view` and `Matrix::as_view_mut`, which are very useful when working with view-based APIs.
+- Added new methods `Matrix::as_view` and `Matrix::as_view_mut`, which are very useful when working with view-based
+  APIs.
 - Added parallel column iterator using `rayon`: `Matrix::par_column_iter` and `Matrix::par_column_iter_mut`. The `rayon`
  feature must be enabled to access these methods.
 - Implement `ReshapableStorage` for matrix slices (only for unit strides at the moment).
@ -53,76 +82,89 @@ This project adheres to [Semantic Versioning](https://semver.org/).
  `Const::<4>` when we need const dimensions.

 ### Fixed
+
 - Fixed the implementation of `Rotation3::euler_angles` to return the angles in the documented order (roll, pitch, yaw).

 ## [0.31.4] (13 Nov. 2022)

 ### Added
- Add a `convert-glam022` feature to enable conversion between `nalgebra` and `glam v0.22`.

+- Add a `convert-glam022` feature to enable conversion between `nalgebra` and `glam v0.22`.

 ## [0.31.3] (30 Oct. 2022)

 ### Added
+
 - Add `Matrix::try_cast` to attempt casting the inner scalar types when that cast may fail.

 ### Fixed
+
 - Fixed the usage of `CheckBytes` with `rkyv`.

 ## [0.31.2] (09 Oct. 2022)

 ### Modified
+
 - Use `#[inline]` on the `Dim` implementation for `Const` to improve opt-level 1 performance.
 - Make the `Point::new` constructions const-fn.

 ### Added
- Add `UnitVector::cast` to change the underlying scalar type.

+- Add `UnitVector::cast` to change the underlying scalar type.

 ## [0.31.1] (31 July 2022)

 ### Modified
+
 - Improve performances of multiplication of two sparse matrices.

 ### Added
+
 - Add `Matrix::from_row_iterator` to build a matrix from an iterator yielding components in row-major order.
 - Add support for conversion from/to types of `glam` 0.21.
 - `nalgebra-sparse`: add support for the matrix-market export of sparse matrices.
 - `nalgebra-lapack`: add a `GE` for solving the generalized eigenvalues problem.

 ### Fixed
+
 - Fix `Rotation3::from_matrix` and `UnitQuaternion::from_matrix` when the input matrix is already a valid
  rotation matrix.

 ## [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.
@ -132,6 +174,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
  return `Option<Self>`, meaning that it could be implemented by any type.

 ### Modified
+
 - Use more concise debug impls for matrices and geometric transformation types.
 - The singular values computed by the SVD are now sorted in increasing order by default. Use `SVD::new_unordered`
  instead to reproduce the older behavior without the sorting overhead.
@ -141,6 +184,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - The `Matrix::pow` and `Matrix::pow_mut` methods will now also work with integer matrices.

 ### Added
+
 - Added the conversion trait `From<Vec<T>>` and method `from_vec_storage` for `RowDVector`.
 - Added implementation of `From` and `Into` for converting between `nalgebra` types and types from
  `glam 0.18`. These can be enabled by enabling the `convert-glam018` cargo features.
@ -161,6 +205,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - `nalgebra-sparse`: added reading MatrixMarked data files to a sparse `CooMatrix`.

 ### Fixed
+
 - Fixed a potential unsoundness with `matrix.get(i)` and `matrix.get_mut(i)` where `i`  is an `usize`, and `matrix`
  is a matrix slice with non-default strides.
 - Fixed potential unsoundness with `vector.perp` where `vector` isn’t actually a 2D vector as expected.
@ -169,9 +214,10 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - Fixed the `no-std` build of `nalgebra-glm`.
 - Fix the `pow` and `pow_mut` functions (the result was incorrect for some exponent values).

-
 ## [0.29.0]
+
 ### Breaking changes
+
 - We updated to the version 0.6 of `simba`. This means that the trait bounds `T: na::RealField`, `na::ComplexField`,
  `na::SimdRealField`, `na:SimdComplexField` no imply that `T: Copy` (they only imply that `T: Clone`). This may affect
  generic code.
@ -183,6 +229,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
  way.

 ### Modified
+
 - `Orthographic3::from_matrix_unchecked` is now `const fn`.
 - `Perspective3::from_matrix_unchecked` is now `const fn`.
 - `Rotation::from_matrix_unchecked` is now `const fn`.
@ -190,6 +237,7 @@ This project adheres to [Semantic Versioning](https://semver.org/).
  `Copy` are now much safer to work with thanks to the refactoring of the `Allocator` system.

 ### Added
+
 - The conversion traits form the `bytemuck` crates are now implemented for the geometric types too.
 - Added operator overloading for `Transform * UnitComplex`, `UnitComplex * Transform`, `Transform ×= UnitComplex`,
  `Transform ÷= UnitComplex`.
@ -200,16 +248,21 @@ This project adheres to [Semantic Versioning](https://semver.org/).
  cargo features.

 ## [0.28.0]
+
 ### Added
+
 - Implement `Hash` for `Transform`.
 - Implement `Borrow` and `BorrowMut` for contiguous slices.

 ### Modified
- The `OPoint<T, D>` type has been added. It takes the dimension number as a type-level integer (e.g. `Const<3>`) instead
+
+- The `OPoint<T, D>` type has been added. It takes the dimension number as a type-level integer (e.g. `Const<3>`)
+  instead
  of a const-generic. The type `Point<T, const D: usize>` is now an alias for `OPoint`. This changes doesn't affect any
  of the existing code using `Point`. However, it will allow the use `OPoint` in a generic context where the dimension
  cannot be easily expressed as a const-generic (because of the current limitation of const-generics in Rust).
- Several clippy warnings were fixed. This results in some method signature changes (e.g. taking `self` instead of `&self`)
+- Several clippy warnings were fixed. This results in some method signature changes (e.g. taking `self` instead
+  of `&self`)
  but this should not have any practical infulances on existing codebase.
 - The `Point::new` constructors are no longer const-fn. This is due to some limitations in const-fn
  not allowing custom trait-bounds. Use the `point!` macro instead to build points in const environments.
@ -217,30 +270,38 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 - Methods returning `Result<(), ()>` now return `bool` instead.

 ### Fixed
+
 - Fixed a potential unsoundess issue when converting a mutable slice to a `&mut[T]`.

 ## [0.27.1]
+
 ### Fixed
+
 - Fixed a bug in the conversion from `glam::Vec2` or `glam::DVec2` to `Isometry2`.

 ## [0.27.0]
+
 This removes the `convert-glam` and `convert-glam-unchecked` optional features.
 Instead, this adds the `convert-glam013`, `convert-glam014`, and `convert-glam015` optional features for
 conversions targeting the versions 0.13, 0.14, and 0.15 of `glam`.

 ### Added
+
 - Add macros `matrix!`, `dmatrix!`, `vector!`, `dvector!`, `point!` for constructing matrices/vectors/points in a
-  more convenient way. See [#886](https://github.com/dimforge/nalgebra/pull/886) and [#899](https://github.com/dimforge/nalgebra/pull/899).
+  more convenient way. See [#886](https://github.com/dimforge/nalgebra/pull/886)
+  and [#899](https://github.com/dimforge/nalgebra/pull/899).
 - Add `CooMatrix::reserve` to `nalgebra-sparse`.
 - Add basic support for serialization using `rkyv`. Can be enabled with the features `rkyv-serialize` or
  `rkyv-serialize-no-std`.

-
 ### Fixed
+
 - Fixed a potential unsoundness issue after deserializing an invalid `DVector` using `serde`.

 ## [0.26.2]
+
 ### Added
+
 - Conversion from an array `[T; D]` to an isometry `Isometry<T, _, D>` (as a translation).
 - Conversion from a static vector `SVector<T; D>` to an isometry `Isometry<T, _, D>` (as a translation).
 - Conversion from a point `Point<T; D>` to an isometry `Isometry<T, _, D>` (as a translation).
@ -250,14 +311,17 @@ conversions targeting the versions 0.13, 0.14, and 0.15 of `glam`.
 - Conversion of a glam type `Vec2/3/4` from/to a `Translation2/3/4`.

 ## [0.26.1]
+
 Fix a regression introduced in 0.26.0 preventing `DVector` from being serialized with `serde`.

 ## [0.26.0]
+
 This release integrates `min-const-generics` to nalgebra. See
 [our blog post](https://www.dimforge.com/blog/2021/04/12/integrating-const-generics-to-nalgebra)
 for details about this release.

 ### Added
+
 - Add type aliases for unit vector, e.g., `UnitVector3`.
 - Add a `pow` and `pow_mut` function to square matrices.
 - Add `Cholesky::determinant` to compute the determinant of a matrix decomposed
@ -265,8 +329,8 @@ for details about this release.
 - Add the `serde-serialize-no-std` feature to enable serialization of static matrices/vectors
  with serde, but without requiring `std`.

-
 ### Modified
+
 - The `serde` crate isn't enabled by default now. Enable the `serde-serialize` or the
  `serde-serialize-no-std` features instead.
 - The `Const<const D: usize>` type has been introduced to represent dimensions known
@ -288,16 +352,22 @@ for details about this release.
  constructor is also a `const fn` now.

 ## [0.25.4]
+
 ### Fixed
+
 - Fix a compilation error when only the `serde-serialize` feature is enabled.

 ## [0.25.3]
+
 ### Added
+
 - The `Vector::simd_cap_magnitude` method to cap the magnitude of the vector with
  SIMD components.

 ## [0.25.2]
+
 ### Added
+
 - A `convert-glam` cargo feature to enable implementations of `From` traits to convert
  between `glam` types and `nalgebra` types.
 - A `convert-glam-unchecked` cargo feature to enable some extra `glam`/`nalgebra` conversions that may
@ -308,56 +378,70 @@ for details about this release.
  type of the components of a given entity. Example: `vector.cast::<f32>()`.

 ## [0.25.1]
+
 This release replaces the version 0.25.0 which has been yanked. The 0.25.0 version
 added significant complication to build `nalgebra` targeting a `#[no-std]` platform
 not supported by `rand`.

 The `rand` dependency is now optional (and disabled by default). You may enable it with:
+
 - The `rand-no-std` cargo feature when targeting a `#[no-std]` environment.
 - The `rand` cargo feature when targeting a `std` environment.

 ## [0.25.0] - Yanked
+
 This updates all the dependencies of nalgebra to their latest version, including:
+
 - rand 0.8
 - proptest 1.0
 - simba 0.4

 ### New crate: nalgebra-sparse
+
 Alongside this release of `nalgebra`, we are releasing `nalgebra-sparse`: a crate dedicated to sparse matrix
 computation with `nalgebra`. The `sparse` module of `nalgebra`itself still exists for backward compatibility,
 but it will be deprecated soon in favor of the `nalgebra-sparse` crate.

 ### Added
+
 * Add `UnitDualQuaternion`, a dual-quaternion with unit magnitude which can be used as an isometry transformation.
 * Add `UDU::new()` and `matrix.udu()` to compute the UDU factorization of a matrix.
 * Add `ColPivQR::new()` and `matrix.col_piv_qr()` to compute the QR decomposition with column pivoting of a matrix.
-* Add `from_basis_unchecked` to all the rotation types. This builds a rotation from a set of basis vectors (representing the columns of the corresponding rotation matrix).
+* Add `from_basis_unchecked` to all the rotation types. This builds a rotation from a set of basis vectors (representing
+  the columns of the corresponding rotation matrix).
 * Add `Matrix::cap_magnitude` to cap the magnitude of a vector.
-* Add `UnitQuaternion::append_axisangle_linearized` to approximately append a rotation represented as an axis-angle to a rotation represented as an unit quaternion.
+* Add `UnitQuaternion::append_axisangle_linearized` to approximately append a rotation represented as an axis-angle to a
+  rotation represented as an unit quaternion.
 * Mark the iterators on matrix components as `DoubleEndedIter`.
 * Re-export `simba::simd::SimdValue` at the root of the `nalgebra` crate.

 ## [0.24.0]

 ### Added
+
 * The `DualQuaternion` type. It is still work-in-progress, but the basics are here:
  creation from its real and dual part, multiplication of two dual quaternions,
  and normalization.

 ### Removed
+
 * There is no blanket `impl<T> PartialEq for Unit<T>` anymore. Instead, it is
  implemented specifically for `UnitComplex`, `UnitQuaternion` and `Unit<Vector>`.

 ## [0.23.2]
+
 In this release, we improved the documentation of some of the geometric types
 by applying changes similar to what we did in the version 0.23.1 for matrices.

 ### Added
+
 * The `Isometry::inv_mul` method which is a more efficient way of doing
  `isometry1.inverse() * isometry2`.

 ## [0.23.1]
+
 In this release we improved the documentation of the matrix and vector types by:
+
 - Grouping `impl` bocks logically, adding a title comment to these impl blocks.
 - Reference these impl blocks docs at the top of the documentation page for `Matrix`.
 - Reduce the depth of type aliasing. Now all vector and matrix types are aliases of `Matrix`
@ -366,13 +450,17 @@ In this release we improved the documentation of the matrix and vector types by:
 ## [0.23.0]

 ### Added
- * The `.inverse_transform_unit_vector(v)` was added to `Rotation2/3`, `Isometry2/3`, `UnitQuaternion`, and `UnitComplex`.
+
+* The `.inverse_transform_unit_vector(v)` was added to `Rotation2/3`, `Isometry2/3`, `UnitQuaternion`,
+  and `UnitComplex`.
  It applies the corresponding rotation to a unit vector `Unit<Vector2/3>`.
- * The `Point.map(f)` and `Point.apply(f)` to apply a function to each component of the point, similarly to `Vector.map(f)`
+* The `Point.map(f)` and `Point.apply(f)` to apply a function to each component of the point, similarly
+  to `Vector.map(f)`
  and `Vector.apply(f)`.
 * The `Quaternion::from([N; 4])` conversion to build a quaternion from an array of four elements.
 * The `Isometry::from(Translation)` conversion to build an isometry from a translation.
- * The `Vector::ith_axis(i)` which build a unit vector, e.g., `Unit<Vector3<f32>>` with its i-th component set to 1.0, and the
+* The `Vector::ith_axis(i)` which build a unit vector, e.g., `Unit<Vector3<f32>>` with its i-th component set to 1.0,
+  and the
  others set to zero.
 * The `Isometry.lerp_slerp` and `Isometry.try_lerp_slerp` methods to interpolate between two isometries using linear
  interpolation for the translational part, and spherical interpolation for the rotational part.
@ -380,12 +468,14 @@ In this release we improved the documentation of the matrix and vector types by:
  spherical interpolation.

 ## [0.22.0]
+
 In this release, we are using the new version 0.2 of simba. One major change of that version is that the
 use of `libm` is now opt-in when building targeting `no-std` environment. If you are using floating-point
 operations with nalgebra in a `no-std` environment, you will need to enable the new `libm` feature
 of nalgebra for your code to compile again.

 ### Added
+
 * The `libm` feature that enables `libm` when building for `no-std` environment.
 * The `libm-force` feature that enables `libm` even when building for a not `no-std` environment.
 * `Cholesky::new_unchecked` which build a Cholesky decomposition without checking that its input is
@ -396,6 +486,7 @@ of nalgebra for your code to compile again.
 * The `Vector::ith(i, x)` that builds a vector filled with zeros except for the `i`-th component set to `x`.

 ## [0.21.0]
+
 In this release, we are no longer relying on traits from the __alga__ crate for our generic code.
 Instead, we use traits from the new [simba](https://crates.io/crates/simba) crate which are both
 simpler, and allow for significant optimizations like AoSoA SIMD.
@ -404,75 +495,98 @@ Refer to the [monthly dimforge blogpost](https://www.dimforge.org/blog/2020/04/0
 for details about this switch and its benefits.

 ### Added
+
 * It is now possible to use SIMD types like `simba::f32x4` as scalar types for nalgebra's matrices and
  geometric types.
+
 ### Modified
+
 * Use of traits like `alga::general::{RealField, ComplexField}` have now been replaced by
  `simba::scalar::{RealField, ComplexField}`.
 * The implementation of traits from the __alga__ crate (and well as the dependency to _alga__) are now
  omitted unless the `alga` cargo feature is activated.
+
 ### Removed
+
 * The `Neg` unary operator is no longer implemented for `UnitComplex` and `UnitQuaternion`. This caused
  hard-to-track errors when we mistakenly write, e.g., `-q * v` instead of `-(q * v)`.
 * The `na::convert_unchecked` is no longer marked as unsafe.

 ## [0.20.0]
+
 ### Added
+
 * `cholesky.rank_one_update(...)` which performs a rank-one update on the cholesky decomposition of a matrix.
 * `From<&Matrix>` is now implemented for matrix slices.
 * `.try_set_magnitude(...)` which sets the magnitude of a vector, while keeping its direction.
 * Implementations of `From` and `Into` for the conversion between matrix slices and standard (`&[N]` `&mut [N]`) slices.

 ### Modified
+
 * We started some major changes in order to allow non-Copy types to be used as scalar types inside of matrices/vectors.

 ## [0.19.0]
+
 ### Added
+
 * `.remove_rows_at` and `remove_columns_at` which removes a set of rows or columns (specified by indices) from a matrix.
 * Several formatting traits have been implemented for all matrices/vectors: `LowerExp`, `UpperExp`, `Octal`, `LowerHex`,
  `UpperHex`, `Binary`, `Pointer`.
 * `UnitQuaternion::quaternions_mean(...)` which computes the mean rotation of a set of unit quaternions. This implements
-    the algorithm from _Oshman, Yaakov, and Avishy Carmi, "Attitude estimation from vector observations using a genetic-algorithm-embedded quaternion particle filter."
+  the algorithm from _Oshman, Yaakov, and Avishy Carmi, "Attitude estimation from vector observations using a
+  genetic-algorithm-embedded quaternion particle filter."

 ### Modified
+
 * It is now possible to get the `min/max` element of unsigned integer matrices.

 ### Added to nalgebra-glm
+
 * Some infinite and reversed perspectives: `::infinite_perspective_rh_no`, `::infinite_perspective_rh_zo`,
  `::reversed_perspective_rh_zo`, and `::reversed_infinite_perspective_rh_zo`.

 ## [0.18.0]
+
 This release adds full complex number support to nalgebra. This includes all common vector/matrix operations as well
 as matrix decomposition. This excludes geometric type (like `Isometry`, `Rotation`, `Translation`, etc.) from the
 `geometry` module.

 ### Added
+
 #### Quaternion and geometric operations
-  * Add trigonometric functions for quaternions: `.cos, .sin, .tan, .acos, .asin, .atan, .cosh, .sinh, .tanh, .acosh, .asinh, .atanh`.
+
+* Add trigonometric functions for
+  quaternions: `.cos, .sin, .tan, .acos, .asin, .atan, .cosh, .sinh, .tanh, .acosh, .asinh, .atanh`.
 * Add geometric algebra operations for quaternions: `.inner, .outer, .project, .rejection`
 * Add `.left_div, .right_div` for quaternions.
 * Add `.renormalize` to `Unit<...>` and `Rotation3` to correct potential drift due to repeated operations.
  Those drifts could cause them not to be pure rotations anymore.

 #### Convolution
+
 * `.convolve_full(kernel)` returns the convolution of `self` by `kernel`.
-  * `.convolve_valid(kernel)` returns the convolution of `self` by `kernel` after removal of all the elements relying on zero-padding.
+* `.convolve_valid(kernel)` returns the convolution of `self` by `kernel` after removal of all the elements relying on
+  zero-padding.
 * `.convolve_same(kernel)` returns the convolution of `self` by `kernel` with a result of the same size as `self`.

 #### Complex number support
+
 * Add the `::from_matrix` constructor too all rotation types to extract a rotation from a raw matrix.
 * Add the `::from_matrix_eps` constructor too all rotation types to extract a rotation from a raw matrix. This takes
  more argument than `::from_matrix` to control the convergence of the underlying optimization algorithm.
 * Add `.camax()` which returns the matrix component with the greatest L1-norm.
 * Add `.camin()` which returns the matrix component with the smallest L1-norm.
 * Add `.ad_mul(b)` for matrix-multiplication of `self.adjoint() * b`.
-  * Add `.ad_mul_to(b)` which is the same as `.ad_mul` but with a provided matrix to be filled with the result of the multiplication.
+* Add `.ad_mul_to(b)` which is the same as `.ad_mul` but with a provided matrix to be filled with the result of the
+  multiplication.
 * Add BLAS operations involving complex conjugation (following similar names as the original BLAS spec):
    * `.dotc(rhs)` equal to  `self.adjoint() * rhs`.
    * `.gerc(alpha, x, y, beta)` equivalent to `self = alpha * x * y.adjoint() + beta * self`
    * `.hegerc` which is like `gerc` but for Hermitian matrices.
-      * `.syger` which is the new name of `.ger_symm` which is equivalent to `self = alpha * x * y.transpose() + beta * self`.
-      * `.sygemv` which is the new name of `.gemv_symm` which is equivalent to `self = alpha * a * x + beta * self` with `a` symmetric.
+    * `.syger` which is the new name of `.ger_symm` which is equivalent
+      to `self = alpha * x * y.transpose() + beta * self`.
+    * `.sygemv` which is the new name of `.gemv_symm` which is equivalent to `self = alpha * a * x + beta * self`
+      with `a` symmetric.
    * `.hegemv(alpha, a, x, beta)` which is like `.sygemv` but with `a` Hermitian.
    * `.gemv_ad(alpha, a, x, beta)` which is equivalent to `self = alpha * a.adjoint() * x + beta * self`.
    * `.gemm_ad(alpha, a, b, beta)` which is equivalent to `self = alpha * a.adjoint() * b + beta * self`.
@ -481,14 +595,17 @@ as matrix decomposition. This excludes geometric type (like `Isometry`, `Rotatio
 Note that all the other BLAS operation will continue to work for all fields, including floats and complex numbers.

 ### Renamed
-  * `RealSchur` has been renamed `Schur` because it can now work with complex matrices.

+* `RealSchur` has been renamed `Schur` because it can now work with complex matrices.

 ## [0.17.0]

 ### Added
-  * Add swizzling up to dimension 3 for vectors. For example, you can do `v.zxy()` as an equivalent to `Vector3::new(v.z, v.x, v.y)`.
-  * Add swizzling up to dimension 3 for points. For example, you can do `p.zxy()` as an equivalent to `Point3::new(p.z, p.x, p.y)`.
+
+* Add swizzling up to dimension 3 for vectors. For example, you can do `v.zxy()` as an equivalent
+  to `Vector3::new(v.z, v.x, v.y)`.
+* Add swizzling up to dimension 3 for points. For example, you can do `p.zxy()` as an equivalent
+  to `Point3::new(p.z, p.x, p.y)`.
 * Add `.copy_from_slice` to copy matrix components from a slice in column-major order.
 * Add `.dot` to quaternions.
 * Add `.zip_zip_map` for iterating on three matrices simultaneously, and applying a closure to them.
@ -499,7 +616,8 @@ Note that all the other BLAS operation will continue to work for all fields, inc
 * Add `From/Into` impls to allow the conversion of any transformation type to a matrix.
 * Add `Into` impls to convert a matrix slice into an owned matrix.
 * Add `Point*::from_slice` to create a point from a slice.
-  * Add `.map_with_location` to matrices to apply a map which passes the component indices to the user-defined closure alongside
+* Add `.map_with_location` to matrices to apply a map which passes the component indices to the user-defined closure
+  alongside
  the component itself.
 * Add impl `From<Vector>` for `Point`.
 * Add impl `From<Vector4>` for `Quaternion`.
@ -509,7 +627,8 @@ Note that all the other BLAS operation will continue to work for all fields, inc
 * Add `.to_homogeneous` to square matrices (and with dimensions higher than 1x1). This will increase their number of row
  and columns by 1. The new column and row are filled with 0, except for the diagonal element which is set to 1.
 * Implement `Extend<Vec>` for matrices with a dynamic storage. The provided `Vec` is assumed to represent a column-major
-    matrix with the same number of rows as the one being extended. This will effectively append new columns on the right of
+  matrix with the same number of rows as the one being extended. This will effectively append new columns on the right
+  of
  the matrix being extended.
 * Implement `Extend<Vec>` for vectors with a dynamic storage. This will concatenate the vector with the given `Vec`.
 * Implement `Extend<Matrix<...>>` for matrices with dynamic storage. This will concatenate the columns of both matrices.
@ -518,10 +637,12 @@ Note that all the other BLAS operation will continue to work for all fields, inc
 * Add a `.len()` method to retrieve the size of a `MatrixVec`.

 ### Modified
+
 * The orthographic projection no longer require that `bottom < top`, that `left < right`, and that `znear < zfar`. The
  only restriction now ith that they must not be equal (in which case the projection would be singular).
 * The `Point::from_coordinates` methods is deprecated. Use `Point::from` instead.
-  * The `.transform_point` and `.transform_vector` methods are now inherent methods for matrices so that the user does not have to
+* The `.transform_point` and `.transform_vector` methods are now inherent methods for matrices so that the user does not
+  have to
  explicitly import the `Transform` trait from the alga crate.
 * Renamed the matrix storage types: `MatrixArray` -> `ArrayStorage` and `MatrixVec` -> `VecStorage`.
 * Renamed `.unwrap()` to `.into_inner()` for geometric types that wrap another type.
@ -529,32 +650,46 @@ Note that all the other BLAS operation will continue to work for all fields, inc
 * Deprecate several functions at the root of the crate (replaced by methods).

 ### Removed
+
    * Remove the `Deref` impl for `MatrixVec` as it could cause hard-to-understand compilation errors.

 ### nalgebra-glm
+
 * Add several alternative projection computations, e.g., `ortho_lh`, `ortho_lh_no`, `perspective_lh`, etc.
-  * Add features matching those of nalgebra, in particular: `serde-serialize`, `abmonation-serialize`, std` (enabled by default).
+* Add features matching those of nalgebra, in particular:`serde-serialize`, `abmonation-serialize`, std` (enabled by
+  default).

 ## [0.16.0]
+
 All dependencies have been updated to their latest versions.

 ## Modified
-  * Adjust `UnitQuaternion`s, `Rotation3`s, and `Rotation2`s generated from the `Standard` distribution to be uniformly distributed.
+
+* Adjust `UnitQuaternion`s, `Rotation3`s, and `Rotation2`s generated from the `Standard` distribution to be uniformly
+  distributed.
+
 ### Added
+
 * Add a feature `stdweb` to activate the dependency feature `rand/stdweb`.
 * Add blas-like methods `.imin()` and `.imax()` that return the index of the minimum and maximum entry of a vector.
 * Add construction of a `Point` from an array by implementing the `From` trait.
 * Add support for generating uniformly distributed random unit column vectors using the `Standard` distribution.

 ## [0.15.0]
+
 The most notable change of this release is the support for using part of the library without the rust standard
-library (i.e. it supports `#![no_std]`). See the corresponding [documentation](https://nalgebra.org/wasm_and_embedded_programming/).
+library (i.e. it supports `#![no_std]`). See the
+corresponding [documentation](https://nalgebra.org/wasm_and_embedded_programming/).
+
 ### Modified
+
 * Rename the `core` module to `base` to avoid conflicts with the `core` crate implicitly imported when
  `#![no_std]` is enabled.
 * Constructors of the `MatrixSlice*` types have been renamed from `new_*` to `from_slice_*`. This was
  necessary to avoid the `incoherent_fundamental_impls` lint that is going to become a hard error.
+
 ### Added
+
 * Add `UnitQuaternion` constructor `::new_eps(...)` and `::from_scaled_axis_eps(...)` that return the
  identity if the magnitude of the input axisangle is smaller than the epsilon provided.
 * Add methods `.rotation_between_axis(...)` and `.scaled_rotation_between_axis(...)` to `UnitComplex`
@ -564,9 +699,13 @@ library (i.e. it supports `#![no_std]`). See the corresponding [documentation](h
 * Add functions to construct a random matrix with a user-defined distribution: `::from_distribution(...)`.

 ## [0.14.0]
+
 ### Modified
+
 * Allow the `Isometry * Unit<Vector>` multiplication.
+
 ### Added
+
 * Add blas-like operations: `.quadform(...)` and `.quadform_tr(...)` to compute respectively
  the quadratic forms `self = alpha * A.transpose() * B * A + beta * self` and
  `alpha * A * B * A.transpose() + beta * self`. Here, `A, B` are matrices with
@ -596,8 +735,6 @@ library (i.e. it supports `#![no_std]`). See the corresponding [documentation](h
  `UnitQuaternion::rotation_between_axis(...)` that take Unit vectors instead of
  Vector as arguments.

-
-
 ## [0.13.0]

 The **nalgebra-lapack** crate has been updated. This now includes a broad range
@ -607,6 +744,7 @@ This adds support for serialization using the
 [abomonation](https://crates.io/crates/abomonation) crate.

 ### Breaking semantic change
+
 * The implementation of slicing with steps now matches the documentation.
  Before, step identified the number to add to pass from one column/row index
  to the next one. This made 0 step invalid. Now (and on the documentation so
@ -615,6 +753,7 @@ This adds support for serialization using the
  example, a step of, say, 3 on previous versions should now bet set to 2.

 ### Modified
+
 * The trait `Axpy` has been replaced by a method `.axpy`.
 * The alias `MatrixNM` is now deprecated. Use `MatrixMN` instead (we
  reordered M and N to be in alphabetical order).
@ -624,6 +763,7 @@ This adds support for serialization using the
  `.component_div_assign(...)` instead.

 ### Added
+
 * `alga::general::Real` is now re-exported by nalgebra.
  elements.)
 * `::zeros(...)` that creates a matrix filled with zeroes.
@ -659,17 +799,20 @@ Pure Rust implementation of some Blas operations:
 * `.ger_symm(...)` is the same as `.ger` except that `self` is assumed symmetric.

 New slicing methods:
+
 * `.rows_range(...)` that retrieves a reference to a range of rows.
 * `.rows_range_mut(...)` that retrieves a mutable reference to a range of rows.
 * `.columns_range(...)` that retrieves a reference to a range of columns.
 * `.columns_range_mut(...)` that retrieves a mutable reference to a range of columns.

 Matrix decompositions implemented in pure Rust:
+
 * Cholesky, SVD, LU, QR, Hessenberg, Schur, Symmetric eigendecompositions,
  Bidiagonal, Symmetric tridiagonal
 * Computation of householder reflectors and givens rotations.

 Matrix edition:
+
 * `.upper_triangle()` extracts the upper triangle of a matrix, including the diagonal.
 * `.lower_triangle()` extracts the lower triangle of a matrix, including the diagonal.
 * `.fill(...)` fills the matrix with a single value.
@ -686,37 +829,45 @@ Matrix edition:
 * `.swap_columns(...)` swaps two columns.

 Column removal:
+
 * `.remove_column(...)` removes one column.
 * `.remove_fixed_columns<D>(...)` removes `D` columns.
 * `.remove_columns(...)` removes a number of columns known at run-time.

 Row removal:
+
 * `.remove_row(...)` removes one row.
 * `.remove_fixed_rows<D>(...)` removes `D` rows.
 * `.remove_rows(...)` removes a number of rows known at run-time.

 Column insertion:
+
 * `.insert_column(...)` adds one column at the given position.
 * `.insert_fixed_columns<D>(...)` adds `D` columns at the given position.
 * `.insert_columns(...)` adds at the given position a number of columns known at run-time.

 Row insertion:
+
 * `.insert_row(...)` adds one row at the given position.
 * `.insert_fixed_rows<D>(...)` adds `D` rows at the given position.
 * `.insert_rows(...)` adds at the given position a number of rows known at run-time.

 ## [0.12.0]
+
 The main change of this release is the update of the dependency serde to 1.0.

 ### Added
+
 * `.trace()` that computes the trace of a matrix (the sum of its diagonal
  elements.)

 ## [0.11.0]
+
 The [website](https://nalgebra.org) has been fully rewritten and gives a good
 overview of all the added/modified features.

 This version is a major rewrite of the library. Major changes are:
+
 * Algebraic traits are now defined by the [alga](https://crates.io/crates/alga) crate.
  All other mathematical traits, except `Axpy` have been removed from
  **nalgebra**.
@ -734,9 +885,11 @@ This version is a major rewrite of the library. Major changes are:
 * Matrix **slices** are now implemented.

 ### Added
+
 Lots of features including rectangular matrices, slices, and Serde
 serialization. Refer to the brand new [website](https://nalgebra.org) for more
 details. The following free-functions have been added as well:
+
 * `::id()` that returns the universal [identity element](https://nalgebra.org/performance_tricks/#the-id-type)
  of type `Id`.
 * `::inf_sup()` that returns both the infimum and supremum of a value at the
@ -746,6 +899,7 @@ details. The following free-functions have been added as well:
  the interval width to it.

 ### Modified
+
 * `::cast`            -> `::convert`
 * `point.as_vector()` -> `point.coords`
 * `na::origin`        -> `P::origin()`
@ -756,6 +910,7 @@ details. The following free-functions have been added as well:
 * `::angle_between`   -> `::angle`

 Componentwise multiplication and division has been replaced by methods:
+
 * multiplication -> `.componentwise_mul`, `.componentwise_mul_mut`.
 * division -> `.componentwise_div`, `.componentwise_div_mut`.

@ -765,8 +920,8 @@ only:
 `.eig`), `::hessenberg`, `::qr`, `::to_homogeneous`, `::to_rotation_matrix`,
 `::transpose`, `::shape`.

-
 The following free-functions are now replaced by static methods only:
+
 * `::householder_matrix` under the name `::new_householder_generic`
 * `::identity`
 * `::new_identity` under the name `::identity`
@ -775,6 +930,7 @@ The following free-functions are now replaced by static methods only:

 The following free-function are now replaced methods accessible through traits
 only:
+
 * `::transform` -> methods `.transform_point` and `.transform_vector` of the `alga::linear::Transformation` trait.
 * `::inverse_transform` -> methods `.inverse_transform_point` and
  `.inverse_transform_vector` of the `alga::linear::ProjectiveTransformation`
@ -791,9 +947,8 @@ only:
  `alga::linear::Rotation` trait.
 * `::is_zero` -> method with the same name from `num::Zero`.

-
-
 ### Removed
+
 * The free functions `::prepend_rotation`, `::append_rotation`,
  `::append_rotation_wrt_center`, `::append_rotation_wrt_point`,
  `::append_transformation`, and `::append_translation ` have been removed.
@ -817,11 +972,14 @@ only:
  `rotation.angle()` and `rotation.axis()`.

 ## [0.10.0]
+
 ### Added
+
 Binary operations are now allowed between references as well. For example
 `Vector3<f32> + &Vector3<f32>` is now possible.

 ### Modified
+
 Removed unused parameters to methods from the `ApproxEq` trait. Those were
 required before rust 1.0 to help type inference. They are not needed any more
 since it now allowed to write for a type `T` that implements `ApproxEq`:
@ -829,7 +987,9 @@ since it now allowed to write for a type `T` that implements `ApproxEq`:
 `ApproxEq::approx_epsilon(None::<T>)`.

 ## [0.9.0]
+
 ### Modified
+
 * Renamed:
    - `::from_col_vector` -> `::from_column_vector`
    - `::from_col_iter` -> `::from_column_iter`
@ -851,11 +1011,13 @@ Other similar trait changes are to be expected in the future, e.g., for the

 Methods marked `unsafe` for reasons unrelated to memory safety are no
 longer unsafe. Instead, their name end with `_unchecked`. In particular:
+
 * `Rotation3::new_with_matrix` -> `Rotation3::from_matrix_unchecked`
 * `PerspectiveMatrix3::new_with_matrix` -> `PerspectiveMatrix3::from_matrix_unchecked`
 * `OrthographicMatrix3::new_with_matrix` -> `OrthographicMatrix3::from_matrix_unchecked`

 ### Added
+
 - A `Unit<T>` type that wraps normalized values. In particular,
  `UnitQuaternion<N>` is now an alias for `Unit<Quaternion<N>>`.
 - `.ln()`, `.exp()` and `.powf(..)` for quaternions and unit quaternions.
@ -872,7 +1034,9 @@ crate for vectors, rotations and points. To enable them, activate the
 `abstract_algebra` feature.

 ## [0.8.0]
+
 ### Modified
+
 * Almost everything (types, methods, and traits) now use fulls names instead
  of abbreviations (e.g. `Vec3` becomes `Vector3`). Most changes are obvious.
  Note however that:
@ -886,19 +1050,24 @@ crate for vectors, rotations and points. To enable them, activate the
  e.g., `vec.rs` becomes `vector.rs`.

 ## [0.7.0]
+
 ### Added
+
 * Added implementation of assignment operators (+=, -=, etc.) for
  everything.
+
 ### Modified
+
 * Points and vectors are now linked to each other with associated types
  (on the PointAsVector trait).

-
 ## [0.6.0]
+
 **Announcement:** a users forum has been created for `nalgebra`, `ncollide`, and `nphysics`. See
 you [there](https://users.nphysics.org)!

 ### Added
+
 * Added a dependency to [generic-array](https://crates.io/crates/generic-array). Feature-gated:
  requires `features="generic_sizes"`.
 * Added statically sized vectors with user-defined sizes: `VectorN`. Feature-gated: requires
@ -907,10 +1076,13 @@ you [there](https://users.nphysics.org)!
  translation): `Similarity2`, `Similarity3`.

 ### Removed
+
 * Removed zero-sized elements `Vector0`, `Point0`.
-  * Removed 4-dimensional transformations `Rotation4` and `Isometry4` (which had an implementation too incomplete to be useful).
+* Removed 4-dimensional transformations `Rotation4` and `Isometry4` (which had an implementation too incomplete to be
+  useful).

 ### Modified
+
 * Vectors are now multipliable with isometries. This will result into a pure rotation (this is how
  vectors differ from point semantically: they design directions, so they are not translatable).
 * `{Isometry3, Rotation3}::look_at` reimplemented and renamed to `::look_at_rh` and `::look_at_lh` to agree
--- a/Cargo.toml
+++ b/Cargo.toml
@ -23,7 +23,7 @@ path = "src/lib.rs"

 [features]
 default = [ "std", "macros" ]
-std     = [ "matrixmultiply", "simba/std" ]
+std     = [ "matrixmultiply", "num-traits/std", "num-complex/std", "num-rational/std", "approx/std", "simba/std" ]
 sparse  = [ ]
 debug   = [ "approx/num-complex", "rand" ]
 alloc   = [ ]
@ -32,7 +32,6 @@ compare = [ "matrixcompare-core" ]
 libm    = [ "simba/libm" ]
 libm-force = [ "simba/libm_force" ]
 macros = [ "nalgebra-macros" ]
-cuda   = [ "cust_core", "simba/cuda" ]


 # Conversion
@ -49,6 +48,7 @@ convert-glam021 = [ "glam021" ]
 convert-glam022 = [ "glam022" ]
 convert-glam023 = [ "glam023" ]
 convert-glam024 = [ "glam024" ]
+convert-glam025 = [ "glam025" ]

 # Serialization
 ## To use serde in a #[no-std] environment, enable the
@ -104,7 +104,7 @@ glam021        = { package = "glam", version = "0.21", optional = true }
 glam022        = { package = "glam", version = "0.22", optional = true }
 glam023        = { package = "glam", version = "0.23", optional = true }
 glam024        = { package = "glam", version = "0.24", optional = true }
-cust_core      = { version = "0.1", optional = true }
+glam025        = { package = "glam", version = "0.25", optional = true }
 rayon          = { version = "1.6", optional = true }

 [dev-dependencies]
--- a/examples/dimensional_genericity.rs
+++ b/examples/dimensional_genericity.rs
@ -28,7 +28,7 @@ where
 }

 /// Reflects a 3D vector wrt. the 3D plane with normal `plane_normal`.
-/// /!\ This is an exact replicate of `reflect_wrt_hyperplane2, but for 3D.
+/// /!\ This is an exact replicate of `reflect_wrt_hyperplane2`, but for 3D.
 fn reflect_wrt_hyperplane3<T>(plane_normal: &Unit<Vector3<T>>, vector: &Vector3<T>) -> Vector3<T>
 where
    T: RealField,
--- a/nalgebra-glm/Cargo.toml
+++ b/nalgebra-glm/Cargo.toml
@ -21,7 +21,6 @@ default               = [ "std" ]
 std                   = [ "nalgebra/std", "simba/std" ]
 arbitrary             = [ "nalgebra/arbitrary" ]
 serde-serialize       = [ "nalgebra/serde-serialize-no-std" ]
-cuda                  = [ "nalgebra/cuda" ]

 # Conversion
 convert-mint = [ "nalgebra/mint" ]
--- a/nalgebra-glm/src/ext/quaternion_common.rs
+++ b/nalgebra-glm/src/ext/quaternion_common.rs
@ -1,4 +1,4 @@
-use na::{self, Unit};
+use na::Unit;

 use crate::aliases::Qua;
 use crate::RealNumber;
--- a/nalgebra-glm/src/gtx/vector_query.rs
+++ b/nalgebra-glm/src/gtx/vector_query.rs
@ -41,7 +41,10 @@ pub fn is_comp_null<T: Number, const D: usize>(v: &TVec<T, D>, epsilon: T) -> TV

 /// Returns `true` if `v` has a magnitude of 1 (up to an epsilon).
 pub fn is_normalized<T: RealNumber, const D: usize>(v: &TVec<T, D>, epsilon: T) -> bool {
-    abs_diff_eq!(v.norm_squared(), T::one(), epsilon = epsilon * epsilon)
+    // sqrt(1 + epsilon_{norm²} = 1 + epsilon_{norm}
+    // ==> epsilon_{norm²} = epsilon_{norm}² + 2*epsilon_{norm}
+    // For small epsilon, epsilon² is basically zero, so use 2*epsilon.
+    abs_diff_eq!(v.norm_squared(), T::one(), epsilon = epsilon + epsilon)
 }

 /// Returns `true` if `v` is zero (up to an epsilon).
--- a/nalgebra-glm/src/traits.rs
+++ b/nalgebra-glm/src/traits.rs
@ -1,7 +1,6 @@
 use approx::AbsDiffEq;
 use num::{Bounded, Signed};

-use core::cmp::PartialOrd;
 use na::Scalar;
 use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, RealField};

--- a/nalgebra-glm/src/vector_relational.rs
+++ b/nalgebra-glm/src/vector_relational.rs
@ -141,7 +141,7 @@ pub fn less_than<T: Number, const D: usize>(x: &TVec<T, D>, y: &TVec<T, D>) -> T
    x.zip_map(y, |x, y| x < y)
 }

-/// Component-wise `>=` comparison.
+/// Component-wise `<=` comparison.
 ///
 /// # Examples:
 ///
--- a/nalgebra-sparse/src/cs.rs
+++ b/nalgebra-sparse/src/cs.rs
@ -1,4 +1,3 @@
-use std::mem::replace;
 use std::ops::Range;

 use num_traits::One;
@ -369,7 +368,7 @@ where
        if let Some(minor_indices) = lane {
            let count = minor_indices.len();

-            let remaining = replace(&mut self.remaining_values, &mut []);
+            let remaining = std::mem::take(&mut self.remaining_values);
            let (values_in_lane, remaining) = remaining.split_at_mut(count);
            self.remaining_values = remaining;
            self.current_lane_idx += 1;
@ -578,7 +577,7 @@ where
    } else if sort {
        unreachable!("Internal error: Sorting currently not supported if no values are present.");
    }
-    if major_offsets.len() == 0 {
+    if major_offsets.is_empty() {
        return Err(SparseFormatError::from_kind_and_msg(
            SparseFormatErrorKind::InvalidStructure,
            "Number of offsets should be greater than 0.",
@ -624,12 +623,12 @@ where
                ));
            }

-            let minor_idx_in_lane = minor_indices.get(range_start..range_end).ok_or(
+            let minor_idx_in_lane = minor_indices.get(range_start..range_end).ok_or_else(|| {
                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
--- a/nalgebra-sparse/src/csc.rs
+++ b/nalgebra-sparse/src/csc.rs
@ -625,6 +625,15 @@ pub struct CscTripletIter<'a, T> {
    values_iter: Iter<'a, T>,
 }

+impl<'a, T> Clone for CscTripletIter<'a, T> {
+    fn clone(&self) -> Self {
+        CscTripletIter {
+            pattern_iter: self.pattern_iter.clone(),
+            values_iter: self.values_iter.clone(),
+        }
+    }
+}
+
 impl<'a, T: Clone> CscTripletIter<'a, T> {
    /// Adapts the triplet iterator to return owned values.
    ///
--- a/nalgebra-sparse/src/csr.rs
+++ b/nalgebra-sparse/src/csr.rs
@ -626,6 +626,15 @@ pub struct CsrTripletIter<'a, T> {
    values_iter: Iter<'a, T>,
 }

+impl<'a, T> Clone for CsrTripletIter<'a, T> {
+    fn clone(&self) -> Self {
+        CsrTripletIter {
+            pattern_iter: self.pattern_iter.clone(),
+            values_iter: self.values_iter.clone(),
+        }
+    }
+}
+
 impl<'a, T: Clone> CsrTripletIter<'a, T> {
    /// Adapts the triplet iterator to return owned values.
    ///
--- a/nalgebra-sparse/src/ops/impl_std_ops.rs
+++ b/nalgebra-sparse/src/ops/impl_std_ops.rs
@ -59,7 +59,7 @@ macro_rules! impl_sp_plus_minus {
            let mut result = $matrix_type::try_from_pattern_and_values(pattern, values)
                .unwrap();
            $spadd_fn(T::zero(), &mut result, T::one(), Op::NoOp(&a)).unwrap();
-            $spadd_fn(T::one(), &mut result, $factor * T::one(), Op::NoOp(&b)).unwrap();
+            $spadd_fn(T::one(), &mut result, $factor, Op::NoOp(&b)).unwrap();
            result
        });

--- a/nalgebra-sparse/src/ops/serial/pattern.rs
+++ b/nalgebra-sparse/src/ops/serial/pattern.rs
@ -125,18 +125,22 @@ fn iterate_union<'a>(
 ) -> impl Iterator<Item = usize> + 'a {
    iter::from_fn(move || {
        if let (Some(a_item), Some(b_item)) = (sorted_a.first(), sorted_b.first()) {
-            let item = if a_item < b_item {
+            let item = match a_item.cmp(b_item) {
+                std::cmp::Ordering::Less => {
                    sorted_a = &sorted_a[1..];
                    a_item
-            } else if b_item < a_item {
+                }
+                std::cmp::Ordering::Greater => {
                    sorted_b = &sorted_b[1..];
                    b_item
-            } else {
+                }
+                std::cmp::Ordering::Equal => {
                    // Both lists contain the same element, advance both slices to avoid
                    // duplicate entries in the result
                    sorted_a = &sorted_a[1..];
                    sorted_b = &sorted_b[1..];
                    a_item
+                }
            };
            Some(*item)
        } else if let Some(a_item) = sorted_a.first() {
--- a/nalgebra-sparse/src/pattern.rs
+++ b/nalgebra-sparse/src/pattern.rs
@ -80,7 +80,7 @@ impl SparsityPattern {
    #[inline]
    #[must_use]
    pub fn major_dim(&self) -> usize {
-        assert!(self.major_offsets.len() > 0);
+        assert!(!self.major_offsets.is_empty());
        self.major_offsets.len() - 1
    }

@ -162,7 +162,7 @@ impl SparsityPattern {
                // We test for in-bounds, uniqueness and monotonicity at the same time
                // to ensure that we only visit each minor index once
                let mut iter = minor_indices.iter();
-                let mut prev = None;
+                let mut prev: Option<usize> = None;

                while let Some(next) = iter.next().copied() {
                    if next >= minor_dim {
@ -170,10 +170,10 @@ impl SparsityPattern {
                    }

                    if let Some(prev) = prev {
-                        if prev > next {
-                            return Err(NonmonotonicMinorIndices);
-                        } else if prev == next {
-                            return Err(DuplicateEntry);
+                        match prev.cmp(&next) {
+                            std::cmp::Ordering::Greater => return Err(NonmonotonicMinorIndices),
+                            std::cmp::Ordering::Equal => return Err(DuplicateEntry),
+                            std::cmp::Ordering::Less => {}
                        }
                    }
                    prev = Some(next);
@ -195,6 +195,14 @@ impl SparsityPattern {
    ///
    /// 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.
+    ///
+    /// # Safety
+    ///
+    /// Assumes that the major offsets and indices adhere to the requirements of being a valid
+    /// sparsity pattern.
+    /// Specifically, that major offsets is monotonically increasing, and
+    /// `major_offsets[i]..major_offsets[i+1]` refers to a major lane in the sparsity pattern,
+    /// and `minor_indices[major_offsets[i]..major_offsets[i+1]]` is monotonically increasing.
    pub unsafe fn from_offset_and_indices_unchecked(
        major_dim: usize,
        minor_dim: usize,
--- a/src/base/array_storage.rs
+++ b/src/base/array_storage.rs
@ -42,7 +42,6 @@ use std::mem;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[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> {
--- a/src/base/construction_view.rs
+++ b/src/base/construction_view.rs
@ -97,6 +97,8 @@ macro_rules! impl_constructors(
            }

            /// Creates, without bound checking, a new matrix view from the given data array.
+            /// # Safety
+            /// `data[start..start+rstride * cstride]` must be within bounds.
            #[inline]
            pub unsafe fn from_slice_unchecked(data: &'a [T], start: usize, $($args: usize),*) -> Self {
                Self::from_slice_generic_unchecked(data, start, $($gargs),*)
@ -113,6 +115,11 @@ macro_rules! impl_constructors(
            }

            /// Creates, without bound checking, a new matrix view with the specified strides from the given data array.
+            ///
+            /// # Safety
+            ///
+            /// `start`, `rstride`, and `cstride`, with the given matrix size will not index
+            /// outside of `data`.
            #[inline]
            pub unsafe fn from_slice_with_strides_unchecked(data: &'a [T], start: usize, $($args: usize,)* rstride: usize, cstride: usize) -> Self {
                Self::from_slice_with_strides_generic_unchecked(data, start, $($gargs,)* Dyn(rstride), Dyn(cstride))
@ -257,6 +264,10 @@ macro_rules! impl_constructors_mut(
            }

            /// Creates, without bound checking, a new mutable matrix view from the given data array.
+            ///
+            /// # Safety
+            ///
+            /// `data[start..start+(R * C)]` must be within bounds.
            #[inline]
            pub unsafe fn from_slice_unchecked(data: &'a mut [T], start: usize, $($args: usize),*) -> Self {
                Self::from_slice_generic_unchecked(data, start, $($gargs),*)
@ -274,6 +285,8 @@ macro_rules! impl_constructors_mut(
            }

            /// Creates, without bound checking, a new mutable matrix view with the specified strides from the given data array.
+            /// # Safety
+            /// `data[start..start+rstride * cstride]` must be within bounds.
            #[inline]
            pub unsafe fn from_slice_with_strides_unchecked(data: &'a mut [T], start: usize, $($args: usize,)* rstride: usize, cstride: usize) -> Self {
                Self::from_slice_with_strides_generic_unchecked(
--- a/src/base/dimension.rs
+++ b/src/base/dimension.rs
@ -23,7 +23,6 @@ use serde::{Deserialize, Deserializer, Serialize, Serializer};
    feature = "rkyv-serialize",
    archive_attr(derive(bytecheck::CheckBytes))
 )]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Dyn(pub usize);

 #[deprecated(note = "use Dyn instead.")]
@ -68,6 +67,10 @@ impl IsNotStaticOne for Dyn {}

 /// Trait implemented by any type that can be used as a dimension. This includes type-level
 /// integers and `Dyn` (for dimensions not known at compile-time).
+///
+/// # Safety
+///
+/// Hoists integers to the type level, including binary operations.
 pub unsafe trait Dim: Any + Debug + Copy + PartialEq + Send + Sync {
    #[inline(always)]
    fn is<D: Dim>() -> bool {
@ -216,7 +219,6 @@ dim_ops!(
    archive(as = "Self")
 )]
 #[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.
--- a/src/base/indexing.rs
+++ b/src/base/indexing.rs
@ -519,6 +519,10 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {

    /// Produces a view of the data at the given index, without doing
    /// any bounds checking.
+    ///
+    /// # Safety
+    ///
+    /// `index` must within bounds of the array.
    #[inline]
    #[must_use]
    pub unsafe fn get_unchecked<'a, I>(&'a self, index: I) -> I::Output
@ -530,6 +534,9 @@ impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {

    /// Returns a mutable view of the data at the given index, without doing
    /// any bounds checking.
+    /// # Safety
+    ///
+    /// `index` must within bounds of the array.
    #[inline]
    #[must_use]
    pub unsafe fn get_unchecked_mut<'a, I>(&'a mut self, index: I) -> I::OutputMut
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@ -171,7 +171,6 @@ pub type MatrixCross<T, R1, C1, R2, C2> =
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[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?
    ///
@ -313,6 +312,10 @@ where
 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.
+    ///
+    /// # Safety
+    ///
+    /// The storage dimension must match the given dimensions.
    #[inline(always)]
    pub const unsafe fn from_data_statically_unchecked(data: S) -> Matrix<T, R, C, S> {
        Matrix {
@ -1194,6 +1197,10 @@ impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> Matrix<T, R, C, S> {
    }

    /// Swaps two entries without bound-checking.
+    ///
+    /// # Safety
+    ///
+    /// Both `(r, c)` must have `r < nrows(), c < ncols()`.
    #[inline]
    pub unsafe fn swap_unchecked(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize)) {
        debug_assert!(row_cols1.0 < self.nrows() && row_cols1.1 < self.ncols());
@ -1300,6 +1307,8 @@ impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> Matrix<T, R, C, S> {

 impl<T, D: Dim, S: RawStorage<T, D>> Vector<T, D, S> {
    /// Gets a reference to the i-th element of this column vector without bound checking.
+    /// # Safety
+    /// `i` must be less than `D`.
    #[inline]
    #[must_use]
    pub unsafe fn vget_unchecked(&self, i: usize) -> &T {
@ -1311,6 +1320,8 @@ impl<T, D: Dim, S: RawStorage<T, D>> Vector<T, D, S> {

 impl<T, D: Dim, S: RawStorageMut<T, D>> Vector<T, D, S> {
    /// Gets a mutable reference to the i-th element of this column vector without bound checking.
+    /// # Safety
+    /// `i` must be less than `D`.
    #[inline]
    #[must_use]
    pub unsafe fn vget_unchecked_mut(&mut self, i: usize) -> &mut T {
--- a/src/base/matrix_view.rs
+++ b/src/base/matrix_view.rs
@ -43,6 +43,10 @@ macro_rules! view_storage_impl (

        impl<'a, T, R: Dim, C: Dim, RStride: Dim, CStride: Dim> $T<'a, T, R, C, RStride, CStride> {
            /// Create a new matrix view without bounds checking and from a raw pointer.
+            ///
+            /// # Safety
+            ///
+            /// `*ptr` must point to memory that is valid `[T; R * C]`.
            #[inline]
            pub unsafe fn from_raw_parts(ptr:     $Ptr,
                                         shape:   (R, C),
@ -63,6 +67,11 @@ macro_rules! view_storage_impl (
        // Dyn is arbitrary. It's just to be able to call the constructors with `Slice::`
        impl<'a, T, R: Dim, C: Dim> $T<'a, T, R, C, Dyn, Dyn> {
            /// Create a new matrix view without bounds checking.
+            ///
+            /// # Safety
+            ///
+            /// `storage` contains sufficient elements beyond `start + R * C` such that all
+            /// accesses are within bounds.
            #[inline]
            pub unsafe fn new_unchecked<RStor, CStor, S>(storage: $SRef, start: (usize, usize), shape: (R, C))
                -> $T<'a, T, R, C, S::RStride, S::CStride>
@ -75,6 +84,10 @@ macro_rules! view_storage_impl (
            }

            /// Create a new matrix view without bounds checking.
+            ///
+            /// # Safety
+            ///
+            /// `strides` must be a valid stride indexing.
            #[inline]
            pub unsafe fn new_with_strides_unchecked<S, RStor, CStor, RStride, CStride>(storage: $SRef,
                                                                                        start:   (usize, usize),
@ -128,12 +141,7 @@ impl<'a, T: Scalar, R: Dim, C: Dim, RStride: Dim, CStride: Dim> Clone
 {
    #[inline]
    fn clone(&self) -> Self {
-        Self {
-            ptr: self.ptr,
-            shape: self.shape,
-            strides: self.strides,
-            _phantoms: PhantomData,
-        }
+        *self
    }
 }

@ -538,8 +546,8 @@ macro_rules! matrix_view_impl (
            $me.$generic_view_with_steps(start, shape, steps)
        }

-        /// Slices this matrix starting at its component `(irow, icol)` and with `(R::dim(),
-        /// CView::dim())` consecutive components.
+        /// Slices this matrix starting at its component `(irow, icol)` and with `(RVIEW, CVIEW)`
+        /// consecutive components.
        #[inline]
        #[deprecated = slice_deprecation_note!($fixed_view)]
        pub fn $fixed_slice<const RVIEW: usize, const CVIEW: usize>($me: $Me, irow: usize, icol: usize)
@ -547,8 +555,8 @@ macro_rules! matrix_view_impl (
            $me.$fixed_view(irow, icol)
        }

-        /// Return a view of this matrix starting at its component `(irow, icol)` and with `(R::dim(),
-        /// CView::dim())` consecutive components.
+        /// Return a view of this matrix starting at its component `(irow, icol)` and with
+        /// `(RVIEW, CVIEW)` consecutive components.
        #[inline]
        pub fn $fixed_view<const RVIEW: usize, const CVIEW: usize>($me: $Me, irow: usize, icol: usize)
            -> $MatrixView<'_, T, Const<RVIEW>, Const<CVIEW>, S::RStride, S::CStride> {
--- a/src/base/norm.rs
+++ b/src/base/norm.rs
@ -336,7 +336,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
    /// Sets the magnitude of this vector unless it is smaller than `min_magnitude`.
    ///
    /// If `self.magnitude()` is smaller than `min_magnitude`, it will be left unchanged.
-    /// Otherwise this is equivalent to: `*self = self.normalize() * magnitude.
+    /// Otherwise this is equivalent to: `*self = self.normalize() * magnitude`.
    #[inline]
    pub fn try_set_magnitude(&mut self, magnitude: T::RealField, min_magnitude: T::RealField)
    where
@ -525,7 +525,7 @@ where
                let (elt, basis) = vs[..i + 1].split_last_mut().unwrap();

                for basis_element in &basis[..nbasis_elements] {
-                    *elt -= &*basis_element * elt.dot(basis_element)
+                    *elt -= basis_element * elt.dot(basis_element)
                }
            }

--- a/src/base/statistics.rs
+++ b/src/base/statistics.rs
@ -339,7 +339,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> {
            let mean = self.mean();

            self.iter().cloned().fold(T::zero(), |acc, x| {
-                acc + (x.clone() - mean.clone()) * (x.clone() - mean.clone())
+                acc + (x.clone() - mean.clone()) * (x - mean.clone())
            }) / n_elements
        }
    }
--- a/src/base/storage.rs
+++ b/src/base/storage.rs
@ -32,6 +32,8 @@ pub type CStride<T, R, C = U1> =
 /// The trait shared by all matrix data storage.
 ///
 /// TODO: doc
+/// # Safety
+///
 /// In generic code, it is recommended use the `Storage` trait bound instead. The `RawStorage`
 /// trait bound is generally used by code that needs to work with storages that contains
 /// `MaybeUninit<T>` elements.
@ -57,7 +59,7 @@ pub unsafe trait RawStorage<T, R: Dim, C: Dim = U1>: Sized {

    /// The spacing between consecutive row elements and consecutive column elements.
    ///
-    /// For example this returns `(1, 5)` for a row-major matrix with 5 columns.
+    /// For example this returns `(1, 5)` for a column-major matrix with 5 columns.
    fn strides(&self) -> (Self::RStride, Self::CStride);

    /// Compute the index corresponding to the irow-th row and icol-th column of this matrix. The
@ -129,6 +131,14 @@ pub unsafe trait RawStorage<T, R: Dim, C: Dim = U1>: Sized {
 }

 /// Trait shared by all matrix data storage that don’t contain any uninitialized elements.
+///
+/// # Safety
+///
+/// Note that `Self` must always have a number of elements compatible with the matrix length (given
+/// by `R` and `C` if they are known at compile-time). For example, implementors of this trait
+/// should **not** allow the user to modify the size of the underlying buffer with safe methods
+/// (for example the `VecStorage::data_mut` method is unsafe because the user could change the
+/// vector's size so that it no longer contains enough elements: this will lead to UB.
 pub unsafe trait Storage<T, R: Dim, C: Dim = U1>: RawStorage<T, R, C> {
    /// Builds a matrix data storage that does not contain any reference.
    fn into_owned(self) -> Owned<T, R, C>
@ -143,6 +153,8 @@ pub unsafe trait Storage<T, R: Dim, C: Dim = U1>: RawStorage<T, R, C> {

 /// Trait implemented by matrix data storage that can provide a mutable access to its elements.
 ///
+/// # Safety
+///
 /// In generic code, it is recommended use the `StorageMut` trait bound instead. The
 /// `RawStorageMut` trait bound is generally used by code that needs to work with storages that
 /// contains `MaybeUninit<T>` elements.
@ -194,10 +206,28 @@ pub unsafe trait RawStorageMut<T, R: Dim, C: Dim = U1>: RawStorage<T, R, C> {
    ///
    /// # Safety
    /// If the indices are out of bounds, the method will cause undefined behavior.
+    ///
+    /// # Validity
+    /// The default implementation of this trait function is only guaranteed to be
+    /// sound if invocations of `self.ptr_mut()` and `self.get_address_unchecked_linear_mut()`
+    /// result in stable references. If any of the data pointed to by these trait methods
+    /// moves as a consequence of invoking either of these methods then this default
+    /// trait implementation may be invalid or unsound and should be overridden.
    #[inline]
    unsafe fn swap_unchecked_linear(&mut self, i1: usize, i2: usize) {
-        let a = self.get_address_unchecked_linear_mut(i1);
-        let b = self.get_address_unchecked_linear_mut(i2);
+        // we can't just use the pointers returned from `get_address_unchecked_linear_mut` because calling a
+        // method taking self mutably invalidates any existing (mutable) pointers. since `get_address_unchecked_linear_mut` can
+        // also be overriden by a custom implementation, we can't just use `wrapping_add` assuming that's what the method does.
+        // instead, we use `offset_from` to compute the re-calculate the pointers from the base pointer.
+        // this is sound as long as this trait matches the Validity preconditions
+        // (and it's the caller's responsibility to ensure the indices are in-bounds).
+        let base = self.ptr_mut();
+        let offset1 = self.get_address_unchecked_linear_mut(i1).offset_from(base);
+        let offset2 = self.get_address_unchecked_linear_mut(i2).offset_from(base);
+
+        let base = self.ptr_mut();
+        let a = base.offset(offset1);
+        let b = base.offset(offset2);

        ptr::swap(a, b);
    }
@ -226,6 +256,10 @@ pub unsafe trait RawStorageMut<T, R: Dim, C: Dim = U1>: RawStorage<T, R, C> {
 }

 /// Trait shared by all mutable matrix data storage that don’t contain any uninitialized elements.
+///
+/// # Safety
+///
+/// See safety note for `Storage`, `RawStorageMut`.
 pub unsafe trait StorageMut<T, R: Dim, C: Dim = U1>:
    Storage<T, R, C> + RawStorageMut<T, R, C>
 {
@ -241,6 +275,8 @@ where

 /// Marker trait indicating that a storage is stored contiguously in memory.
 ///
+/// # Safety
+///
 /// The storage requirement means that for any value of `i` in `[0, nrows * ncols - 1]`, the value
 /// `.get_unchecked_linear` returns one of the matrix component. This trait is unsafe because
 /// failing to comply to this may cause Undefined Behaviors.
--- a/src/base/unit.rs
+++ b/src/base/unit.rs
@ -35,7 +35,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-// #[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct Unit<T> {
    pub(crate) value: T,
 }
@ -72,16 +71,6 @@ impl<'de, T: Deserialize<'de>> Deserialize<'de> for Unit<T> {
    }
 }

-#[cfg(feature = "cuda")]
-unsafe impl<T: cust_core::DeviceCopy, R, C, S> cust_core::DeviceCopy for Unit<Matrix<T, R, C, S>>
-where
-    T: Scalar,
-    R: Dim,
-    C: Dim,
-    S: RawStorage<T, R, C> + Copy,
-{
-}
-
 impl<T, R, C, S> PartialEq for Unit<Matrix<T, R, C, S>>
 where
    T: Scalar + PartialEq,
--- a/src/geometry/dual_quaternion.rs
+++ b/src/geometry/dual_quaternion.rs
@ -55,7 +55,6 @@ use simba::scalar::{ClosedNeg, RealField};
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub struct DualQuaternion<T> {
    /// The real component of the quaternion
    pub real: Quaternion<T>,
@ -320,6 +319,7 @@ where
 }

 impl<T: RealField> DualQuaternion<T> {
+    #[allow(clippy::wrong_self_convention)]
    fn to_vector(&self) -> OVector<T, U8> {
        self.as_ref().clone().into()
    }
--- a/src/geometry/isometry.rs
+++ b/src/geometry/isometry.rs
@ -54,7 +54,6 @@ use rkyv::bytecheck;
 ///
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
    feature = "serde-serialize-no-std",
--- a/src/geometry/orthographic.rs
+++ b/src/geometry/orthographic.rs
@ -34,7 +34,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Orthographic3<T> {
    matrix: Matrix4<T>,
--- a/src/geometry/perspective.rs
+++ b/src/geometry/perspective.rs
@ -35,7 +35,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Perspective3<T> {
    matrix: Matrix4<T>,
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@ -86,14 +86,6 @@ where
 {
 }

-#[cfg(feature = "cuda")]
-unsafe impl<T: Scalar + cust_core::DeviceCopy, D: DimName> cust_core::DeviceCopy for OPoint<T, D>
-where
-    DefaultAllocator: Allocator<T, D>,
-    OVector<T, D>: cust_core::DeviceCopy,
-{
-}
-
 #[cfg(feature = "bytemuck")]
 unsafe impl<T: Scalar, D: DimName> bytemuck::Zeroable for OPoint<T, D>
 where
@ -317,6 +309,10 @@ where
    }

    /// Gets a reference to i-th element of this point without bound-checking.
+    ///
+    /// # Safety
+    ///
+    /// `i` must be less than `self.len()`.
    #[inline]
    #[must_use]
    pub unsafe fn get_unchecked(&self, i: usize) -> &T {
@ -344,6 +340,10 @@ where
    }

    /// Gets a mutable reference to i-th element of this point without bound-checking.
+    ///
+    /// # Safety
+    ///
+    /// `i` must be less than `self.len()`.
    #[inline]
    #[must_use]
    pub unsafe fn get_unchecked_mut(&mut self, i: usize) -> &mut T {
@ -351,6 +351,10 @@ where
    }

    /// Swaps two entries without bound-checking.
+    ///
+    /// # Safety
+    ///
+    /// `i1` and `i2` must be less than `self.len()`.
    #[inline]
    pub unsafe fn swap_unchecked(&mut self, i1: usize, i2: usize) {
        self.coords.swap_unchecked((i1, 0), (i2, 0))
@ -499,10 +503,11 @@ where

        let mut it = self.coords.iter();

-        write!(f, "{}", *it.next().unwrap())?;
+        <T as fmt::Display>::fmt(it.next().unwrap(), f)?;

        for comp in it {
-            write!(f, ", {}", *comp)?;
+            write!(f, ", ")?;
+            <T as fmt::Display>::fmt(comp, f)?;
        }

        write!(f, "}}")
--- a/src/geometry/point_construction.rs
+++ b/src/geometry/point_construction.rs
@ -202,29 +202,11 @@ impl<T: Scalar> Point1<T> {
    /// assert_eq!(p.x, 1.0);
    /// ```
    #[inline]
-    #[cfg(not(feature = "cuda"))]
    pub const fn new(x: T) -> Self {
        Point {
            coords: Vector1::new(x),
        }
    }
-
-    /// Initializes this point from its components.
-    ///
-    /// # Example
-    ///
-    /// ```
-    /// # use nalgebra::Point1;
-    /// let p = Point1::new(1.0);
-    /// assert_eq!(p.x, 1.0);
-    /// ```
-    #[inline]
-    #[cfg(feature = "cuda")]
-    pub fn new(x: T) -> Self {
-        Point {
-            coords: Vector1::new(x),
-        }
-    }
 }
 macro_rules! componentwise_constructors_impl(
    ($($doc: expr; $Point: ident, $Vector: ident, $($args: ident:$irow: expr),*);* $(;)*) => {$(
@ -234,22 +216,9 @@ macro_rules! componentwise_constructors_impl(
            #[doc = $doc]
            #[doc = "```"]
            #[inline]
-            #[cfg(not(feature = "cuda"))]
            pub const fn new($($args: T),*) -> Self {
                Point { coords: $Vector::new($($args),*) }
            }
-
-            // TODO: always let new be const once CUDA updates its supported
-            //       nightly version to something more recent.
-            #[doc = "Initializes this point from its components."]
-            #[doc = "# Example\n```"]
-            #[doc = $doc]
-            #[doc = "```"]
-            #[inline]
-            #[cfg(feature = "cuda")]
-            pub fn new($($args: T),*) -> Self {
-                Point { coords: $Vector::new($($args),*) }
-            }
        }
    )*}
 );
--- a/src/geometry/quaternion.rs
+++ b/src/geometry/quaternion.rs
@ -38,7 +38,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[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>,
@ -1016,9 +1015,6 @@ 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(feature = "cuda")]
-unsafe impl<T: cust_core::DeviceCopy> cust_core::DeviceCopy for UnitQuaternion<T> {}
-
 impl<T: Scalar + ClosedNeg + PartialEq> PartialEq for UnitQuaternion<T> {
    #[inline]
    fn eq(&self, rhs: &Self) -> bool {
--- a/src/geometry/rotation.rs
+++ b/src/geometry/rotation.rs
@ -64,7 +64,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[derive(Copy, Clone)]
 pub struct Rotation<T, const D: usize> {
    matrix: SMatrix<T, D, D>,
@ -185,6 +184,10 @@ impl<T: Scalar, const D: usize> Rotation<T, D> {
    }

    /// A mutable reference to the underlying matrix representation of this rotation.
+    ///
+    /// # Safety
+    ///
+    /// Invariants of the rotation matrix should not be violated.
    #[inline]
    #[deprecated(note = "Use `.matrix_mut_unchecked()` instead.")]
    pub unsafe fn matrix_mut(&mut self) -> &mut SMatrix<T, D, D> {
--- a/src/geometry/rotation_specialization.rs
+++ b/src/geometry/rotation_specialization.rs
@ -1058,7 +1058,7 @@ impl<T: SimdRealField> Rotation3<T> {
    {
        let mut angles = [T::zero(); 3];
        let eps = T::from_subset(&1e-7);
-        let _2 = T::from_subset(&2.0);
+        let two = T::from_subset(&2.0);

        if extrinsic {
            seq.reverse();
@ -1090,7 +1090,7 @@ impl<T: SimdRealField> Rotation3<T> {
            -s1,
            c1,
        );
-        let o_t = &c * self.matrix() * (c.transpose() * r1l);
+        let o_t = c * self.matrix() * (c.transpose() * r1l);
        angles[1] = o_t.m33.acos();

        let safe1 = angles[1].abs() >= eps;
@ -1131,7 +1131,7 @@ impl<T: SimdRealField> Rotation3<T> {
        // dont adjust gimbal locked rotation
        if adjust && observable {
            angles[0] += T::pi();
-            angles[1] = _2 * lambda - angles[1];
+            angles[1] = two * lambda - angles[1];
            angles[2] -= T::pi();
        }

--- a/src/geometry/scale.rs
+++ b/src/geometry/scale.rs
@ -32,7 +32,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[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
@ -149,6 +148,10 @@ impl<T: Scalar, const D: usize> Scale<T, D> {
    ///     assert_eq!(t.inverse_unchecked() * t, Scale2::identity());
    /// }
    /// ```
+    ///
+    /// # Safety
+    ///
+    /// Should only be used if all scaling is known to be non-zero.
    #[inline]
    #[must_use]
    pub unsafe fn inverse_unchecked(&self) -> Scale<T, D>
--- a/src/geometry/scale_ops.rs
+++ b/src/geometry/scale_ops.rs
@ -83,28 +83,28 @@ add_sub_impl!(Mul, mul, ClosedMul;
    (Const<D>, U1), (Const<D>, U1) -> (Const<D>, U1)
    const D; for; where;
    self: &'a Scale<T, D>, right: &'b SVector<T, D>, Output = SVector<T, D>;
-    SVector::from(self.vector.component_mul(right));
+    self.vector.component_mul(right);
    'a, 'b);

 add_sub_impl!(Mul, mul, ClosedMul;
    (Const<D>, U1), (Const<D>, U1) -> (Const<D>, U1)
    const D; for; where;
    self: &'a Scale<T, D>, right: SVector<T, D>, Output = SVector<T, D>;
-    SVector::from(self.vector.component_mul(&right));
+    self.vector.component_mul(&right);
    'a);

 add_sub_impl!(Mul, mul, ClosedMul;
    (Const<D>, U1), (Const<D>, U1) -> (Const<D>, U1)
    const D; for; where;
    self: Scale<T, D>, right: &'b SVector<T, D>, Output = SVector<T, D>;
-    SVector::from(self.vector.component_mul(right));
+    self.vector.component_mul(right);
    'b);

 add_sub_impl!(Mul, mul, ClosedMul;
    (Const<D>, U1), (Const<D>, U1) -> (Const<D>, U1)
    const D; for; where;
    self: Scale<T, D>, right: SVector<T, D>, Output = SVector<T, D>;
-    SVector::from(self.vector.component_mul(&right)); );
+    self.vector.component_mul(&right); );

 // Scale *= Scale
 add_sub_assign_impl!(MulAssign, mul_assign, ClosedMul;
--- a/src/geometry/similarity.rs
+++ b/src/geometry/similarity.rs
@ -21,7 +21,6 @@ use rkyv::bytecheck;
 /// A similarity, i.e., an uniform scaling, followed by a rotation, followed by a translation.
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 #[cfg_attr(feature = "serde-serialize-no-std", derive(Serialize, Deserialize))]
 #[cfg_attr(
    feature = "serde-serialize-no-std",
--- a/src/geometry/transform.rs
+++ b/src/geometry/transform.rs
@ -60,17 +60,14 @@ where

 /// Tag representing the most general (not necessarily inversible) `Transform` type.
 #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub enum TGeneral {}

 /// Tag representing the most general inversible `Transform` type.
 #[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
-#[cfg_attr(feature = "cuda", derive(cust_core::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(feature = "cuda", derive(cust_core::DeviceCopy))]
 pub enum TAffine {}

 impl TCategory for TGeneral {
@ -198,16 +195,6 @@ where
 {
 }

-#[cfg(feature = "cuda")]
-unsafe impl<T: RealField + cust_core::DeviceCopy, C: TCategory, const D: usize>
-    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_core::DeviceCopy,
-{
-}
-
 impl<T: RealField, C: TCategory, const D: usize> Clone for Transform<T, C, D>
 where
    Const<D>: DimNameAdd<U1>,
@ -411,7 +398,7 @@ where
    ///                      3.0, 4.0, 0.0,
    ///                      0.0, 0.0, 1.0);
    /// let t = Transform2::from_matrix_unchecked(m);
-    /// assert_eq!(t.into_inner(), m);
+    /// assert_eq!(t.to_homogeneous(), m);
    /// ```
    #[inline]
    #[must_use]
--- a/src/geometry/translation.rs
+++ b/src/geometry/translation.rs
@ -32,7 +32,6 @@ use rkyv::bytecheck;
    )
 )]
 #[cfg_attr(feature = "rkyv-serialize", derive(bytecheck::CheckBytes))]
-#[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
--- a/src/geometry/unit_complex.rs
+++ b/src/geometry/unit_complex.rs
@ -31,9 +31,6 @@ 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(feature = "cuda")]
-unsafe impl<T: cust_core::DeviceCopy> cust_core::DeviceCopy for UnitComplex<T> {}
-
 impl<T: Scalar + PartialEq> PartialEq for UnitComplex<T> {
    #[inline]
    fn eq(&self, rhs: &Self) -> bool {
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@ -34,6 +34,17 @@ pub fn reflection_axis_mut<T: ComplexField, D: Dim, S: StorageMut<T, D>>(

    if !factor.is_zero() {
        column.unscale_mut(factor.sqrt());
+
+        // Normalize again, making sure the vector is unit-sized.
+        // If `factor` had a very small value, the first normalization
+        // (dividing by `factor.sqrt()`) might end up with a slightly
+        // non-unit vector (especially when using 32-bits float).
+        // Decompositions strongly rely on that unit-vector property,
+        // so we run a second normalization (that is much more numerically
+        // stable since the norm is close to 1) to ensure it has a unit
+        // size.
+        let _ = column.normalize_mut();
+
        (-signed_norm, true)
    } else {
        // TODO: not sure why we don't have a - sign here.
--- a/src/linalg/svd2.rs
+++ b/src/linalg/svd2.rs
@ -21,7 +21,7 @@ pub fn svd_ordered2<T: RealField>(
    // because q >= 0 and r >= 0.
    let sx = q.clone() + r.clone();
    let sy = q - r;
-    let sy_sign = if sy < T::zero() { -one.clone() } else { one };
+    let sy_sign = if sy < T::zero() { -one } else { one };
    let singular_values = Vector2::new(sx, sy * sy_sign.clone());

    if compute_u || compute_v {
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@ -360,7 +360,7 @@ mod test {
    #[test]
    fn wilkinson_shift_random() {
        for _ in 0..1000 {
-            let m = Matrix2::new_random();
+            let m = Matrix2::<f64>::new_random();
            let m = m * m.transpose();

            let expected = expected_shift(m);
--- a/src/third_party/glam/mod.rs
+++ b/src/third_party/glam/mod.rs
@ -20,3 +20,5 @@ mod v022;
 mod v023;
 #[cfg(feature = "glam024")]
 mod v024;
+#[cfg(feature = "glam025")]
+mod v025;
--- a/src/third_party/glam/v025/mod.rs
+++ b/src/third_party/glam/v025/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 glam025 as glam;
--- a/tests/geometry/point.rs
+++ b/tests/geometry/point.rs
@ -92,3 +92,11 @@ fn to_homogeneous() {

    assert_eq!(a.to_homogeneous(), expected);
 }
+
+#[test]
+fn display_fmt_respects_modifiers() {
+    let p = Point3::new(1.23, 3.45, 5.67);
+    assert_eq!(&format!("{p}"), "{1.23, 3.45, 5.67}");
+    assert_eq!(&format!("{p:.1}"), "{1.2, 3.5, 5.7}");
+    assert_eq!(&format!("{p:.0}"), "{1, 3, 6}");
+}
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@ -1,4 +1,4 @@
-use na::DMatrix;
+use na::{DMatrix, Matrix3};

 #[cfg(feature = "proptest-support")]
 mod proptest_tests {
@ -116,6 +116,31 @@ fn symmetric_eigen_singular_24x24() {
    );
 }

+// Regression test for #1368
+#[test]
+fn very_small_deviation_from_identity_issue_1368() {
+    let m = Matrix3::<f32>::new(
+        1.0,
+        3.1575704e-23,
+        8.1146196e-23,
+        3.1575704e-23,
+        1.0,
+        1.7471054e-22,
+        8.1146196e-23,
+        1.7471054e-22,
+        1.0,
+    );
+
+    for v in m
+        .try_symmetric_eigen(f32::EPSILON, 0)
+        .unwrap()
+        .eigenvalues
+        .into_iter()
+    {
+        assert_relative_eq!(*v, 1.);
+    }
+}
+
 //  #[cfg(feature = "arbitrary")]
 //  quickcheck! {
 // TODO: full eigendecomposition is not implemented yet because of its complexity when some