diff -r 102421d462d1 -r 53ab61a41d70 src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Wed Sep 03 21:03:47 2025 -0500
+++ b/src/nalgebra_support.rs	Wed Sep 03 21:31:13 2025 -0500
@@ -389,7 +389,7 @@
 
 impl<SM, M, E> Projection<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: StorageMut<E, M> + Clone,
+    SM: StorageMut<E, M>,
     M: Dim,
     E: Scalar + Zero + One + Float + RealField,
     DefaultAllocator: Allocator<M>,
@@ -405,7 +405,7 @@
 
 impl<SM, M, E> ProjectionMut<E, Linfinity> for Vector<E, M, SM>
 where
-    SM: StorageMut<E, M> + Clone,
+    SM: StorageMut<E, M>,
     M: Dim,
     E: Scalar + Zero + One + Copy + Float + RealField,
     DefaultAllocator: Allocator<M>,
@@ -464,26 +464,25 @@
     })
 }
 
-// TODO: should allow different input storages in `Euclidean`.
-
-impl<E, M, S> Euclidean<E> for Vector<E, M, S>
+impl<E, M, N, S> Euclidean<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
-    type PrincipalE = OVector<E, M>;
+    type PrincipalE = OMatrix<E, M, N>;
 
     #[inline]
     fn dot<I: Instance<Self>>(&self, other: I) -> E {
-        other.eval_ref(|ref r| Vector::<E, M, S>::dot(self, r))
+        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]
@@ -507,13 +506,14 @@
 }
 
 /// The default norm for `Vector` is [`L2`].
-impl<E, M, S> Normed<E> for Vector<E, M, S>
+impl<E, M, N, S> Normed<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     type NormExp = L2;
 
@@ -524,22 +524,24 @@
 
     #[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>
+impl<E, M, N, S> HasDual<E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
-    type DualSpace = OVector<E, M>;
+    type DualSpace = OMatrix<E, M, N>;
 
-    fn dual_origin(&self) -> OVector<E, M> {
-        OVector::zeros_generic(M::from_usize(self.len()), Const)
+    fn dual_origin(&self) -> OMatrix<E, M, N> {
+        let (m, n) = self.shape_generic();
+        OMatrix::zeros_generic(m, n)
     }
 }
 
@@ -560,7 +562,7 @@
 impl<E, M, S> Dist<L1, E> for Vector<E, M, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    S: Storage<E, M>,
     E: Float + Scalar + Zero + One + RealField,
     DefaultAllocator: Allocator<M>,
     ShapeConstraint: StridesOk<E, M, U1, S>,
@@ -571,13 +573,14 @@
     }
 }
 
-impl<E, M, S> Norm<L2, E> for Vector<E, M, S>
+impl<E, M, N, S> Norm<L2, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn norm(&self, _: L2) -> E {
@@ -585,13 +588,14 @@
     }
 }
 
-impl<E, M, S> Dist<L2, E> for Vector<E, M, S>
+impl<E, M, N, S> Dist<L2, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: L2) -> E {
@@ -599,13 +603,14 @@
     }
 }
 
-impl<E, M, S> Norm<Linfinity, E> for Vector<E, M, S>
+impl<E, M, N, S> Norm<Linfinity, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M>,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn norm(&self, _: Linfinity) -> E {
@@ -613,13 +618,14 @@
     }
 }
 
-impl<E, M, S> Dist<Linfinity, E> for Vector<E, M, S>
+impl<E, M, N, S> Dist<Linfinity, E> for Matrix<E, M, N, S>
 where
     M: Dim,
-    S: Storage<E, M> + Clone,
+    N: Dim,
+    S: Storage<E, M, N>,
     E: Float + Scalar + Zero + One + RealField,
-    DefaultAllocator: Allocator<M>,
-    ShapeConstraint: StridesOk<E, M, U1, S>,
+    DefaultAllocator: Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, M, N, S>,
 {
     #[inline]
     fn dist<I: Instance<Self>>(&self, other: I, _: Linfinity) -> E {