diff -r 1f19c6bbf07b -r 3868555d135c src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Sun Apr 27 20:29:43 2025 -0500
+++ b/src/nalgebra_support.rs	Fri May 15 14:46:30 2026 -0500
@@ -8,107 +8,390 @@
 [`num_traits`] does.
 */
 
-use nalgebra::{
-    Matrix, Storage, StorageMut, OMatrix, Dim, DefaultAllocator, Scalar,
-    ClosedAddAssign, ClosedMulAssign, SimdComplexField, Vector, OVector, RealField,
-    LpNorm, UniformNorm
-};
-use nalgebra::base::constraint::{
-    ShapeConstraint, SameNumberOfRows, SameNumberOfColumns
-};
-use nalgebra::base::dimension::*;
-use nalgebra::base::allocator::Allocator;
-use std::ops::Mul;
-use num_traits::identities::{Zero, One};
+use crate::euclidean::*;
+use crate::instance::{Decomposition, Instance, MyCow, Ownable, Space};
 use crate::linops::*;
-use crate::euclidean::*;
-use crate::mapping::{Space, BasicDecomposition};
+use crate::norms::*;
 use crate::types::Float;
-use crate::norms::*;
-use crate::instance::Instance;
+use nalgebra::base::allocator::Allocator;
+use nalgebra::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint};
+use nalgebra::base::dimension::*;
+use nalgebra::{
+    ArrayStorage, ClosedAddAssign, ClosedMulAssign, DefaultAllocator, Dim, LpNorm, Matrix,
+    MatrixView, OMatrix, OVector, RawStorage, RealField, SimdComplexField, Storage, StorageMut,
+    UniformNorm, VecStorage, Vector, ViewStorage, ViewStorageMut, U1,
+};
+use num_traits::identities::Zero;
+use std::ops::Mul;
+
+impl<S, M, N, E> Ownable for Matrix<E, M, N, S>
+where
+    S: Storage<E, M, N>,
+    M: Dim,
+    N: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, N>,
+{
+    type OwnedVariant = OMatrix<E, M, N>;
+
+    #[inline]
+    fn into_owned(self) -> Self::OwnedVariant {
+        Matrix::into_owned(self)
+    }
+
+    /// Returns an owned instance of a reference.
+    fn clone_owned(&self) -> Self::OwnedVariant {
+        Matrix::clone_owned(self)
+    }
 
-impl<SM,N,M,E> Space for Matrix<E,N,M,SM>
+    fn cow_owned<'b>(self) -> MyCow<'b, Self::OwnedVariant>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.into_owned())
+    }
+
+    fn ref_cow_owned<'b>(&'b self) -> MyCow<'b, Self::OwnedVariant>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.clone_owned())
+    }
+}
+
+trait StridesOk<E, N, M = U1, S = <DefaultAllocator as Allocator<N, M>>::Buffer<E>>:
+    DimEq<Dyn, S::RStride> + DimEq<Dyn, S::CStride>
 where
-    SM: Storage<E,N,M> + Clone,
-    N : Dim, M : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
-    DefaultAllocator : Allocator<N,M>,
+    S: RawStorage<E, N, M>,
+    E: Float,
+    N: Dim,
+    M: Dim,
+    DefaultAllocator: Allocator<N, M>,
 {
-    type Decomp = BasicDecomposition;
+}
+
+impl<E, M> StridesOk<E, Dyn, M, VecStorage<E, Dyn, M>> for ShapeConstraint
+where
+    M: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<Dyn, M>,
+{
 }
 
-impl<SM,SV,N,M,K,E> Mapping<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
-        N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
-        DefaultAllocator : Allocator<N,K>,
-        DefaultAllocator : Allocator<M,K>,
-        DefaultAllocator : Allocator<N,M>,
-        DefaultAllocator : Allocator<M,N> {
-    type Codomain = OMatrix<E,N,K>;
+impl<E, const N: usize, const M: usize> StridesOk<E, Const<N>, Const<M>, ArrayStorage<E, N, M>>
+    for ShapeConstraint
+where
+    E: Float,
+{
+}
+
+macro_rules! strides_ok {
+    ($R:ty, $C:ty where $($qual:tt)*) => {
+        impl<'a, E, N, M, $($qual)*> StridesOk<E, N, M, ViewStorage<'a, E, N, M, $R, $C>> for ShapeConstraint
+        where
+            N: Dim,
+            M: Dim,
+            E: Float,
+            DefaultAllocator: Allocator<N, M>,
+        {
+        }
+        impl<'a, E, N, M, $($qual)*> StridesOk<E, N, M, ViewStorageMut<'a, E, N, M, $R, $C>> for ShapeConstraint
+        where
+            N: Dim,
+            M: Dim,
+            E: Float,
+            DefaultAllocator: Allocator<N, M>,
+        {
+        }
+    };
+}
+
+strides_ok!(Dyn, Dyn where );
+strides_ok!(Dyn, Const<C> where const C : usize);
+strides_ok!(Const<R>, Dyn where const R : usize);
+strides_ok!(Const<R>, Const<C> where const R : usize, const C : usize);
+
+impl<SM, N, M, E> Space for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<N, M>,
+    ShapeConstraint: StridesOk<E, N, M>,
+{
+    type Principal = OMatrix<E, N, M>;
+    type Decomp = MatrixDecomposition;
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct MatrixDecomposition;
+
+impl<E, M, K, S> Decomposition<Matrix<E, M, K, S>> for MatrixDecomposition
+where
+    S: Storage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K>,
+{
+    type Decomposition<'b>
+        = MyCow<'b, OMatrix<E, M, K>>
+    where
+        Matrix<E, M, K, S>: 'b;
+    type Reference<'b>
+        = MatrixView<'b, E, M, K, Dyn, Dyn>
+    where
+        Matrix<E, M, K, S>: 'b;
 
     #[inline]
-    fn apply<I : Instance<Matrix<E,M,K,SV>>>(
-        &self, x : I
-    ) -> Self::Codomain {
-        x.either(|owned| self.mul(owned), |refr| self.mul(refr))
+    fn lift<'b>(r: Self::Reference<'b>) -> Self::Decomposition<'b>
+    where
+        S: 'b,
+    {
+        MyCow::Owned(r.into_owned())
     }
 }
 
+impl<S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition> for Matrix<E, M, K, S2>
+where
+    S1: Storage<E, M, K>,
+    S2: Storage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, S2> + StridesOk<E, M, K>,
+{
+    #[inline]
+    fn eval_ref<'b, R>(
+        &'b self,
+        f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
+    ) -> R
+    where
+        Self: 'b,
+        Matrix<E, M, K, S1>: 'b,
+    {
+        f(self.as_view::<M, K, Dyn, Dyn>())
+    }
 
-impl<'a, SM,SV,N,M,K,E> Linear<Matrix<E,M,K,SV>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV: Storage<E,M,K> + Clone,
-        N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + ClosedAddAssign + ClosedMulAssign,
-        DefaultAllocator : Allocator<N,K>,
-        DefaultAllocator : Allocator<M,K>,
-        DefaultAllocator : Allocator<N,M>,
-        DefaultAllocator : Allocator<M,N> {
+    #[inline]
+    fn own(self) -> OMatrix<E, M, K> {
+        self.into_owned()
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
+
+    #[inline]
+    fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
+}
+
+impl<'a, S1, S2, M, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
+    for &'a Matrix<E, M, K, S2>
+where
+    S1: Storage<E, M, K>,
+    S2: Storage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, S2> + StridesOk<E, M, K>,
+{
+    fn eval_ref<'b, R>(
+        &'b self,
+        f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
+    ) -> R
+    where
+        Self: 'b,
+        Matrix<E, M, K, S1>: 'b,
+    {
+        f((*self).as_view::<M, K, Dyn, Dyn>())
+    }
+
+    #[inline]
+    fn own(self) -> OMatrix<E, M, K> {
+        self.into_owned()
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
+
+    #[inline]
+    fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
 }
 
-impl<SM,SV1,SV2,N,M,K,E> GEMV<E, Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>> for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: StorageMut<E,N,K>,
-      N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + Float,
-      DefaultAllocator : Allocator<N,K>,
-      DefaultAllocator : Allocator<M,K>,
-      DefaultAllocator : Allocator<N,M>,
-      DefaultAllocator : Allocator<M,N> {
+impl<'a, S1, M, SM, K, E> Instance<Matrix<E, M, K, S1>, MatrixDecomposition>
+    for MyCow<'a, Matrix<E, M, K, SM>>
+where
+    S1: Storage<E, M, K>,
+    SM: Storage<E, M, K>,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K, SM> + StridesOk<E, M, K>,
+{
+    #[inline]
+    fn eval_ref<'b, R>(
+        &'b self,
+        f: impl FnOnce(<MatrixDecomposition as Decomposition<Matrix<E, M, K, S1>>>::Reference<'b>) -> R,
+    ) -> R
+    where
+        Self: 'b,
+        Matrix<E, M, K, S1>: 'b,
+    {
+        f(self.as_view::<M, K, Dyn, Dyn>())
+    }
+
+    #[inline]
+    fn own(self) -> OMatrix<E, M, K> {
+        self.into_owned()
+    }
+
+    #[inline]
+    fn cow<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
 
     #[inline]
-    fn gemv<I : Instance<Matrix<E,M,K,SV1>>>(
-        &self, y : &mut Matrix<E,N,K,SV2>, α : E, x : I, β : E
+    fn decompose<'b>(self) -> MyCow<'b, OMatrix<E, M, K>>
+    where
+        Self: 'b,
+    {
+        self.cow_owned()
+    }
+}
+
+impl<SM, N, M, K, E> Mapping<OMatrix<E, M, K>> for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, K> + Allocator<N, K>,
+    ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
+{
+    type Codomain = OMatrix<E, N, K>;
+
+    #[inline]
+    fn apply<I: Instance<OMatrix<E, M, K>>>(&self, x: I) -> Self::Codomain {
+        x.eval_ref(|refr| self.mul(refr))
+    }
+}
+
+impl<'a, SM, N, M, K, E> Linear<OMatrix<E, M, K>> for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float + ClosedMulAssign + ClosedAddAssign,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
+{
+}
+
+impl<SM, SV2, N, M, K, E> GEMV<E, OMatrix<E, M, K>, Matrix<E, N, K, SV2>> for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    SV2: StorageMut<E, N, K>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K>,
+    ShapeConstraint: StridesOk<E, M, K> + StridesOk<E, N, K>,
+{
+    #[inline]
+    fn gemv<I: Instance<OMatrix<E, M, K>>>(
+        &self, y: &mut Matrix<E, N, K, SV2>, α: E, x: I, β: E
     ) {
         x.eval(|x̃| Matrix::gemm(y, α, self, x̃, β))
     }
 
     #[inline]
-    fn apply_mut<'a, I : Instance<Matrix<E,M,K,SV1>>>(&self, y : &mut Matrix<E,N,K,SV2>, x : I) {
+    fn apply_mut<'a, I: Instance<OMatrix<E, M, K>>>(&self, y: &mut Matrix<E, N, K, SV2>, x: I) {
         x.eval(|x̃| self.mul_to(x̃, y))
     }
 }
 
-impl<SM,SV1,M,E> AXPY<E, Vector<E,M,SV1>> for Vector<E,M,SM>
-where SM: StorageMut<E,M> + Clone, SV1: Storage<E,M> + Clone,
-      M : Dim, E : Scalar + Zero + One + Float,
-      DefaultAllocator : Allocator<M> {
-    type Owned = OVector<E, M>;
+impl<S, M, N, E> VectorSpace for Matrix<E, M, N, S>
+where
+    S: Storage<E, M, N>,
+    M: Dim,
+    N: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
+    type Field = E;
+    type PrincipalV = OMatrix<E, M, N>;
 
     #[inline]
-    fn axpy<I : Instance<Vector<E,M,SV1>>>(&mut self, α : E, x : I, β : E) {
-        x.eval(|x̃| Matrix::axpy(self, α, x̃, β))
+    fn similar_origin(&self) -> Self::PrincipalV {
+        let (n, m) = self.shape_generic();
+        OMatrix::zeros_generic(n, m)
+    }
+}
+
+// This can only be implemented for the “principal” OMatrix as parameter, as otherwise
+// we run into problems of multiple implementations when calling the methods.
+impl<M, N, E, S> AXPY<OMatrix<E, M, N>> for Matrix<E, M, N, S>
+where
+    S: StorageMut<E, M, N>,
+    M: Dim,
+    N: Dim,
+    E: Float,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
+    #[inline]
+    fn axpy<I: Instance<OMatrix<E, M, N>>>(&mut self, α: E, x: I, β: E) {
+        x.eval(|x̃| {
+            assert_eq!(self.ncols(), x̃.ncols());
+            // nalgebra does not implement axpy for matrices, and flattenining
+            // also seems difficult, so loop over columns.
+            for (mut y, ỹ) in self.column_iter_mut().zip(x̃.column_iter()) {
+                Vector::axpy(&mut y, α, &ỹ, β)
+            }
+        })
     }
 
     #[inline]
-    fn copy_from<I : Instance<Vector<E,M,SV1>>>(&mut self, y : I) {
-        y.eval(|ỹ| Matrix::copy_from(self, ỹ))
+    fn copy_from<I: Instance<OMatrix<E, M, N>>>(&mut self, y: I) {
+        y.eval_ref(|ỹ| Matrix::copy_from(self, &ỹ))
     }
 
     #[inline]
     fn set_zero(&mut self) {
         self.iter_mut().for_each(|e| *e = E::ZERO);
     }
-
-    #[inline]
-    fn similar_origin(&self) -> Self::Owned {
-        OVector::zeros_generic(M::from_usize(self.len()), Const)
-    }
 }
 
 /* Implemented automatically as Euclidean.
@@ -125,26 +408,52 @@
     }
 }*/
 
-impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
-where SM: StorageMut<E,M> + Clone,
-      M : Dim, E : Scalar + Zero + One + Float + RealField,
-      DefaultAllocator : Allocator<M> {
+impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM>
+where
+    SM: StorageMut<E, M>,
+    M: Dim,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M>,
+{
     #[inline]
-    fn proj_ball_mut(&mut self, ρ : E, _ : Linfinity) {
-        self.iter_mut().for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
+    fn proj_ball(self, ρ: E, exp: Linfinity) -> <Self as Space>::Principal {
+        let mut owned = self.into_owned();
+        owned.proj_ball_mut(ρ, exp);
+        owned
     }
 }
 
-impl<'own,SV1,SV2,SM,N,M,K,E> Adjointable<Matrix<E,M,K,SV1>, Matrix<E,N,K,SV2>>
-for Matrix<E,N,M,SM>
-where SM: Storage<E,N,M>, SV1: Storage<E,M,K> + Clone, SV2: Storage<E,N,K> + Clone,
-      N : Dim, M : Dim, K : Dim, E : Scalar + Zero + One + SimdComplexField,
-      DefaultAllocator : Allocator<N,K>,
-      DefaultAllocator : Allocator<M,K>,
-      DefaultAllocator : Allocator<N,M>,
-      DefaultAllocator : Allocator<M,N> {
-    type AdjointCodomain = OMatrix<E,M,K>;
-    type Adjoint<'a> = OMatrix<E,M,N> where SM : 'a;
+impl<SM, M, E> ProjectionMut<E, Linfinity> for Vector<E, M, SM>
+where
+    SM: StorageMut<E, M>,
+    M: Dim,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M>,
+{
+    #[inline]
+    fn proj_ball_mut(&mut self, ρ: E, _: Linfinity) {
+        self.iter_mut()
+            .for_each(|v| *v = num_traits::clamp(*v, -ρ, ρ))
+    }
+}
+
+impl<'own, SM, N, M, K, E> Adjointable<OMatrix<E, M, K>, OMatrix<E, N, K>> for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, K> + StridesOk<E, M, K>,
+{
+    type AdjointCodomain = OMatrix<E, M, K>;
+    type Adjoint<'a>
+        = OMatrix<E, M, N>
+    where
+        SM: 'a;
 
     #[inline]
     fn adjoint(&self) -> Self::Adjoint<'_> {
@@ -152,6 +461,25 @@
     }
 }
 
+impl<'own, SM, N, M, K, E> SimplyAdjointable<OMatrix<E, M, K>, OMatrix<E, N, K>>
+    for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, K> + StridesOk<E, M, K>,
+{
+    type AdjointCodomain = OMatrix<E, M, K>;
+    type SimpleAdjoint = OMatrix<E, M, N>;
+
+    #[inline]
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        Matrix::adjoint(self)
+    }
+}
 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
 #[inline]
 fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>(
@@ -160,7 +488,7 @@
     m2: &Matrix<T, R2, C2, S2>,
 ) -> T::SimdRealField
 where
-    T:  SimdComplexField,
+    T: SimdComplexField,
     R1: Dim,
     C1: Dim,
     S1: Storage<T, R1, C1>,
@@ -175,40 +503,41 @@
     })
 }
 
-// TODO: should allow different input storages in `Euclidean`.
-
-impl<E,M,S> Euclidean<E>
-for Vector<E,M,S>
-where M : Dim,
-      S : StorageMut<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
-    type Output = OVector<E, M>;
+impl<E, M, N, S> Euclidean<E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
+    type PrincipalE = OMatrix<E, M, N>;
 
     #[inline]
-    fn dot<I : Instance<Self>>(&self, other : I) -> E {
-        Vector::<E,M,S>::dot(self, other.ref_instance())
+    fn dot<I: Instance<Self>>(&self, other: I) -> E {
+        other.eval_ref(|ref r| Matrix::<E, M, N, S>::dot(self, r))
     }
 
     #[inline]
     fn norm2_squared(&self) -> E {
-        Vector::<E,M,S>::norm_squared(self)
+        Matrix::<E, M, N, S>::norm_squared(self)
     }
 
     #[inline]
-    fn dist2_squared<I : Instance<Self>>(&self, other : I) -> E {
-        metric_distance_squared(self, other.ref_instance())
+    fn dist2_squared<I: Instance<Self>>(&self, other: I) -> E {
+        other.eval_ref(|ref r| metric_distance_squared(self, r))
     }
 }
 
-impl<E,M,S> StaticEuclidean<E>
-for Vector<E,M,S>
-where M : DimName,
-      S : StorageMut<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
+impl<E, M, S> StaticEuclidean<E> for Vector<E, M, S>
+where
+    M: DimName,
+    S: Storage<E, M>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M>,
+{
     #[inline]
     fn origin() -> OVector<E, M> {
         OVector::zeros()
@@ -216,13 +545,15 @@
 }
 
 /// The default norm for `Vector` is [`L2`].
-impl<E,M,S> Normed<E>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
+impl<E, M, N, S> Normed<E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
     type NormExp = L2;
 
     #[inline]
@@ -232,92 +563,106 @@
 
     #[inline]
     fn is_zero(&self) -> bool {
-        Vector::<E,M,S>::norm_squared(self) == E::ZERO
+        Matrix::<E, M, N, S>::norm_squared(self) == E::ZERO
     }
 }
 
-impl<E,M,S> HasDual<E>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-    // TODO: Doesn't work with different storage formats.
-    type DualSpace = Self;
+impl<E, M, N, S> HasDual<E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
+    type DualSpace = OMatrix<E, M, N>;
+
+    fn dual_origin(&self) -> OMatrix<E, M, N> {
+        let (m, n) = self.shape_generic();
+        OMatrix::zeros_generic(m, n)
+    }
 }
 
-impl<E,M,S> Norm<E, L1>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M>,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
+impl<E, M, S> Norm<L1, E> for Vector<E, M, S>
+where
+    M: Dim,
+    S: Storage<E, M>,
+    E: Float + RealField,
+{
     #[inline]
-    fn norm(&self, _ : L1) -> E {
+    fn norm(&self, _: L1) -> E {
         nalgebra::Norm::norm(&LpNorm(1), self)
     }
 }
 
-impl<E,M,S> Dist<E, L1>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
+impl<E, M, S> Dist<L1, E> for Vector<E, M, S>
+where
+    M: Dim,
+    S: Storage<E, M>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M>,
+    ShapeConstraint: StridesOk<E, M>,
+{
     #[inline]
-    fn dist<I : Instance<Self>>(&self, other : I, _ : L1) -> E {
-        nalgebra::Norm::metric_distance(&LpNorm(1), self, other.ref_instance())
+    fn dist<I: Instance<Self>>(&self, other: I, _: L1) -> E {
+        other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&LpNorm(1), self, r))
     }
 }
 
-impl<E,M,S> Norm<E, L2>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M>,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
+impl<E, M, N, S> Norm<L2, E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+{
     #[inline]
-    fn norm(&self, _ : L2) -> E {
+    fn norm(&self, _: L2) -> E {
         nalgebra::Norm::norm(&LpNorm(2), self)
     }
 }
 
-impl<E,M,S> Dist<E, L2>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
+impl<E, M, N, S> Dist<L2, E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
     #[inline]
-    fn dist<I : Instance<Self>>(&self, other : I, _ : L2) -> E {
-        nalgebra::Norm::metric_distance(&LpNorm(2), self, other.ref_instance())
+    fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E {
+        other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&LpNorm(2), self, r))
     }
 }
 
-impl<E,M,S> Norm<E, Linfinity>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M>,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
-
+impl<E, M, N, S> Norm<Linfinity, E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+{
     #[inline]
-    fn norm(&self, _ : Linfinity) -> E {
+    fn norm(&self, _: Linfinity) -> E {
         nalgebra::Norm::norm(&UniformNorm, self)
     }
 }
 
-impl<E,M,S> Dist<E, Linfinity>
-for Vector<E,M,S>
-where M : Dim,
-      S : Storage<E,M> + Clone,
-      E : Float + Scalar + Zero + One + RealField,
-      DefaultAllocator : Allocator<M> {
+impl<E, M, N, S> Dist<Linfinity, E> for Matrix<E, M, N, S>
+where
+    M: Dim,
+    N: Dim,
+    S: Storage<E, M, N>,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N>,
+{
     #[inline]
-    fn dist<I : Instance<Self>>(&self, other : I, _ : Linfinity) -> E {
-        nalgebra::Norm::metric_distance(&UniformNorm, self, other.ref_instance())
+    fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E {
+        other.eval_ref(|ref r| nalgebra::Norm::metric_distance(&UniformNorm, self, r))
     }
 }
 
@@ -329,15 +674,15 @@
 /// from [`nalgebra`] conflicting with them. Only when absolutely necessary to work with
 /// nalgebra, one can convert to the nalgebra view of the same type using the methods of
 /// this trait.
-pub trait ToNalgebraRealField : Float {
+pub trait ToNalgebraRealField: Float {
     /// The nalgebra type corresponding to this type. Usually same as `Self`.
     ///
     /// This type only carries `nalgebra` traits.
-    type NalgebraType : RealField;
+    type NalgebraType: RealField;
     /// The “mixed” type corresponding to this type. Usually same as `Self`.
     ///
     /// This type carries both `num_traits` and `nalgebra` traits.
-    type MixedType : RealField + Float;
+    type MixedType: RealField + Float;
 
     /// Convert to the nalgebra view of `self`.
     fn to_nalgebra(self) -> Self::NalgebraType;
@@ -346,10 +691,10 @@
     fn to_nalgebra_mixed(self) -> Self::MixedType;
 
     /// Convert from the nalgebra view of `self`.
-    fn from_nalgebra(t : Self::NalgebraType) -> Self;
+    fn from_nalgebra(t: Self::NalgebraType) -> Self;
 
     /// Convert from the mixed (nalgebra and num_traits) view to `self`.
-    fn from_nalgebra_mixed(t : Self::MixedType) -> Self;
+    fn from_nalgebra_mixed(t: Self::MixedType) -> Self;
 }
 
 impl ToNalgebraRealField for f32 {
@@ -357,17 +702,24 @@
     type MixedType = f32;
 
     #[inline]
-    fn to_nalgebra(self) -> Self::NalgebraType { self }
-
-    #[inline]
-    fn to_nalgebra_mixed(self) -> Self::MixedType { self }
+    fn to_nalgebra(self) -> Self::NalgebraType {
+        self
+    }
 
     #[inline]
-    fn from_nalgebra(t : Self::NalgebraType) -> Self { t }
+    fn to_nalgebra_mixed(self) -> Self::MixedType {
+        self
+    }
 
     #[inline]
-    fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t }
+    fn from_nalgebra(t: Self::NalgebraType) -> Self {
+        t
+    }
 
+    #[inline]
+    fn from_nalgebra_mixed(t: Self::MixedType) -> Self {
+        t
+    }
 }
 
 impl ToNalgebraRealField for f64 {
@@ -375,15 +727,22 @@
     type MixedType = f64;
 
     #[inline]
-    fn to_nalgebra(self) -> Self::NalgebraType { self }
+    fn to_nalgebra(self) -> Self::NalgebraType {
+        self
+    }
 
     #[inline]
-    fn to_nalgebra_mixed(self) -> Self::MixedType { self }
+    fn to_nalgebra_mixed(self) -> Self::MixedType {
+        self
+    }
 
     #[inline]
-    fn from_nalgebra(t : Self::NalgebraType) -> Self { t }
+    fn from_nalgebra(t: Self::NalgebraType) -> Self {
+        t
+    }
 
     #[inline]
-    fn from_nalgebra_mixed(t : Self::MixedType) -> Self { t }
+    fn from_nalgebra_mixed(t: Self::MixedType) -> Self {
+        t
+    }
 }
-