diff -r 21e51de02ab6 -r b071a1b484f8 src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Fri Sep 05 00:16:08 2025 -0500
+++ b/src/nalgebra_support.rs	Fri Sep 05 00:48:59 2025 -0500
@@ -359,10 +359,11 @@
     }
 }
 
-impl<SM, SV1, M, N, E> AXPY<Matrix<E, M, N, SV1>> for Matrix<E, M, N, SM>
+// 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
-    SM: StorageMut<E, M, N>,
-    SV1: Storage<E, M, N>,
+    S: StorageMut<E, M, N>,
     M: Dim,
     N: Dim,
     E: Scalar + Zero + One + Float,
@@ -370,7 +371,7 @@
     ShapeConstraint: StridesOk<E, M, N>,
 {
     #[inline]
-    fn axpy<I: Instance<Matrix<E, M, N, SV1>>>(&mut self, α: E, x: I, β: E) {
+    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
@@ -382,7 +383,7 @@
     }
 
     #[inline]
-    fn copy_from<I: Instance<Matrix<E, M, N, SV1>>>(&mut self, y: I) {
+    fn copy_from<I: Instance<OMatrix<E, M, N>>>(&mut self, y: I) {
         y.eval_ref(|ỹ| Matrix::copy_from(self, &ỹ))
     }