# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1756703056 18000
# Node ID 26ef556870fd2a147835b299b27af7d1a06928fe
# Parent  d6009939e8324152a236ff9c0b5543623077b5bc
AXPY for nalgebra matrices, not just vectors

diff -r d6009939e832 -r 26ef556870fd src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Sat Aug 30 22:43:37 2025 -0500
+++ b/src/nalgebra_support.rs	Mon Sep 01 00:04:16 2025 -0500
@@ -103,24 +103,32 @@
     }
 }
 
-impl<SM, SV1, M, E> AXPY<Vector<E, M, SV1>> for Vector<E, M, SM>
+impl<SM, SV1, M, N, E> AXPY<Matrix<E, M, N, SV1>> for Matrix<E, M, N, SM>
 where
-    SM: StorageMut<E, M> + Clone,
-    SV1: Storage<E, M> + Clone,
+    SM: StorageMut<E, M, N> + Clone,
+    SV1: Storage<E, M, N> + Clone,
     M: Dim,
+    N: Dim,
     E: Scalar + Zero + One + Float,
-    DefaultAllocator: Allocator<M>,
+    DefaultAllocator: Allocator<M, N>,
 {
     type Field = E;
-    type Owned = OVector<E, M>;
+    type Owned = 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 axpy<I: Instance<Matrix<E, M, N, SV1>>>(&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) {
+    fn copy_from<I: Instance<Matrix<E, M, N, SV1>>>(&mut self, y: I) {
         y.eval(|ỹ| Matrix::copy_from(self, ỹ))
     }
 
@@ -131,7 +139,8 @@
 
     #[inline]
     fn similar_origin(&self) -> Self::Owned {
-        OVector::zeros_generic(M::from_usize(self.len()), Const)
+        let (n, m) = self.shape_generic();
+        OMatrix::zeros_generic(n, m)
     }
 }