diff -r 981069ef919b -r 9b0a7475a786 src/linops.rs
--- a/src/linops.rs	Tue Dec 31 10:51:32 2024 -0500
+++ b/src/linops.rs	Tue Dec 31 10:57:13 2024 -0500
@@ -49,8 +49,13 @@
 
 /// Efficient in-place application for [`Linear`] operators.
 #[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> {
+pub trait GEMV<F : Num, X : Space, Y : AXPY<F> = <Self as Mapping<X>>::Codomain> : Linear<X> {
     
+    /// Computes  `αAx + βy`, where `A` is `Self`.
+    fn apply_add_mul<I : Instance<X>>(&self, y : Y, α : F, x : I, β : F) -> Y::Owned {
+        y.add_mul(α, self.apply(x), β)
+    }
+
     /// Computes  `y = αAx + βy`, where `A` is `Self`.
     fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
 
@@ -153,7 +158,7 @@
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
 where
-    Y : AXPY<F, X>,
+    Y : AXPY<F>,
     X : Clone + Space
 {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
@@ -287,6 +292,7 @@
 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
 where
     F : Num,
+    Y : AXPY<F>,
     X : Space,
     T : Linear<X>,
     S : GEMV<F, T::Codomain, Y>,
@@ -361,6 +367,7 @@
 where
     U : Space,
     V : Space,
+    Y : AXPY<F>,
     S : GEMV<F, U, Y>,
     T : GEMV<F, V, Y>,
     F : Num,
@@ -405,13 +412,16 @@
 impl<A, S, T> Linear<A> for ColOp<S, T>
 where
     A : Space,
-    S : Mapping<A>,
-    T : Mapping<A>,
+    S : Linear<A>,
+    T : Linear<A>,
 { }
 
 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
 where
     X : Space,
+    A : AXPY<F>,
+    B : AXPY<F>,
+    Pair<A, B> : AXPY<F>,
     S : GEMV<F, X, A>,
     T : GEMV<F, X, B>,
     F : Num,
@@ -548,8 +558,9 @@
 
 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
+    A : AXPY<F>,
+    B : AXPY<F>,
+    Pair<A, B> : AXPY<F>,
     U : Space,
     V : Space,
     S : GEMV<F, U, A>,