diff -r 848ecc05becf -r 05089fbc0310 src/linops.rs
--- a/src/linops.rs	Tue Dec 31 09:02:55 2024 -0500
+++ b/src/linops.rs	Tue Dec 31 08:30:43 2024 -0500
@@ -6,10 +6,10 @@
 use std::marker::PhantomData;
 use serde::Serialize;
 use crate::types::*;
-pub use crate::mapping::{Mapping, Space};
+pub use crate::mapping::{Mapping, Space, Composition};
 use crate::direct_product::Pair;
 use crate::instance::Instance;
-use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity};
+use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm};
 
 /// Trait for linear operators on `X`.
 pub trait Linear<X : Space> : Mapping<X>
@@ -58,7 +58,7 @@
 pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
 where
     F : Num,
-    X : Space,
+    X : Space + Norm<F, XExp>,
     XExp : NormExponent,
     CodExp : NormExponent
 {
@@ -152,7 +152,7 @@
 
 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
 where
-    X : Space + Clone,
+    X : Space + Clone + Norm<F, E>,
     F : Num,
     E : NormExponent
 {
@@ -173,6 +173,53 @@
     fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
 }
 
+
+impl<S, T, E, X> Linear<X> for Composition<S, T, E>
+where
+    X : Space,
+    T : Linear<X>,
+    S : Linear<T::Codomain>
+{ }
+
+impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
+where
+    F : Num,
+    X : Space,
+    T : Linear<X>,
+    S : GEMV<F, T::Codomain, Y>,
+{
+    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
+        self.outer.gemv(y, α, self.inner.apply(x), β)
+    }
+
+    /// Computes `y = Ax`, where `A` is `Self`
+    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+        self.outer.apply_mut(y, self.inner.apply(x))
+    }
+
+    /// Computes `y += Ax`, where `A` is `Self`
+    fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
+        self.outer.apply_add(y, self.inner.apply(x))
+    }
+}
+
+impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
+where
+    F : Num,
+    X : Space + Norm<F, Xexp>,
+    Z : Space + Norm<F, Zexp>,
+    Xexp : NormExponent,
+    Yexp : NormExponent,
+    Zexp : NormExponent,
+    T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>,
+    S : BoundedLinear<Z, Zexp, Yexp, F>,
+{
+    fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F {
+        let zexp = self.intermediate_norm_exponent;
+        self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp)
+    }
+}
+
 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
 pub struct RowOp<S, T>(pub S, pub T);
 
@@ -488,8 +535,8 @@
         for RowOp<S, T>
         where
             F : Float,
-            A : Space,
-            B : Space,
+            A : Space + Norm<F, ExpA>,
+            B : Space + Norm<F, ExpB>,
             S : BoundedLinear<A, ExpA, ExpR, F>,
             T : BoundedLinear<B, ExpB, ExpR, F>,
             S::Codomain : Add<T::Codomain>,
@@ -516,7 +563,7 @@
         for ColOp<S, T>
         where
             F : Float,
-            A : Space,
+            A : Space + Norm<F, ExpA>,
             S : BoundedLinear<A, ExpA, ExpS, F>,
             T : BoundedLinear<A, ExpA, ExpT, F>,
             ExpA : NormExponent,