diff -r 5e3c1874797d -r 1b3b1687b9ed src/linops.rs
--- a/src/linops.rs	Wed Dec 11 20:45:17 2024 -0500
+++ b/src/linops.rs	Fri Dec 13 22:37:12 2024 -0500
@@ -7,6 +7,7 @@
 use crate::types::*;
 use serde::Serialize;
 pub use crate::mapping::Apply;
+use crate::direct_product::Pair;
 
 /// Trait for linear operators on `X`.
 pub trait Linear<X> : Apply<X, Output=Self::Codomain>
@@ -96,7 +97,7 @@
     fn preadjoint(&self) -> Self::Preadjoint<'_>;
 }
 
-/// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$.
+/// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
 
@@ -152,3 +153,266 @@
     fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
 }
 
+/// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
+pub struct RowOp<S, T>(pub S, pub T);
+
+use std::ops::Add;
+
+impl<A, B, S, T> Apply<Pair<A, B>> for RowOp<S, T>
+where
+    S : Apply<A>,
+    T : Apply<B>,
+    S::Output : Add<T::Output>
+{
+    type Output = <S::Output as Add<T::Output>>::Output;
+
+    fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
+        self.0.apply(a) + self.1.apply(b)
+    }
+}
+
+impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for RowOp<S, T>
+where
+    S : Apply<&'a A>,
+    T : Apply<&'a B>,
+    S::Output : Add<T::Output>
+{
+    type Output = <S::Output as Add<T::Output>>::Output;
+
+    fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
+        self.0.apply(a) + self.1.apply(b)
+    }
+}
+
+impl<A, B, S, T, D> Linear<Pair<A, B>> for RowOp<S, T>
+where
+    RowOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a>  Apply<&'a Pair<A, B>, Output=D>,
+{
+    type Codomain = D;
+}
+
+impl<F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
+where
+    S : GEMV<F, U, Y>,
+    T : GEMV<F, V, Y>,
+    F : Num,
+    Self : Linear<Pair<U, V>, Codomain=Y>
+{
+    fn gemv(&self, y : &mut Y, α : F, x : &Pair<U, V>, β : F) {
+        self.0.gemv(y, α, &x.0, β);
+        self.1.gemv(y, α, &x.1, β);
+    }
+
+    fn apply_mut(&self, y : &mut Y, x : &Pair<U, V>){
+        self.0.apply_mut(y, &x.0);
+        self.1.apply_mut(y, &x.1);
+    }
+
+    /// Computes `y += Ax`, where `A` is `Self`
+    fn apply_add(&self, y : &mut Y, x : &Pair<U, V>){
+        self.0.apply_add(y, &x.0);
+        self.1.apply_add(y, &x.1);
+    }
+}
+
+
+/// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
+pub struct ColOp<S, T>(pub S, pub T);
+
+impl<A, S, T, O> Apply<A> for ColOp<S, T>
+where
+    S : for<'a> Apply<&'a A, Output=O>,
+    T : Apply<A>,
+    A : std::borrow::Borrow<A>,
+{
+    type Output = Pair<O, T::Output>;
+
+    fn apply(&self, a : A) -> Self::Output {
+        Pair(self.0.apply(a.borrow()), self.1.apply(a))
+    }
+}
+
+impl<A, S, T, D> Linear<A> for ColOp<S, T>
+where
+    ColOp<S, T> : Apply<A, Output=D> + for<'a>  Apply<&'a A, Output=D>,
+{
+    type Codomain = D;
+}
+
+impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
+where
+    S : GEMV<F, X, A>,
+    T : GEMV<F, X, B>,
+    F : Num,
+    Self : Linear<X, Codomain=Pair<A, B>>
+{
+    fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &X, β : F) {
+        self.0.gemv(&mut y.0, α, x, β);
+        self.1.gemv(&mut y.1, α, x, β);
+    }
+
+    fn apply_mut(&self, y : &mut Pair<A, B>, x : &X){
+        self.0.apply_mut(&mut y.0, x);
+        self.1.apply_mut(&mut y.1, x);
+    }
+
+    /// Computes `y += Ax`, where `A` is `Self`
+    fn apply_add(&self, y : &mut Pair<A, B>, x : &X){
+        self.0.apply_add(&mut y.0, x);
+        self.1.apply_add(&mut y.1, x);
+    }
+}
+
+
+impl<A, B, Yʹ, R, S, T> Adjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
+where
+    S : Adjointable<A, Yʹ>,
+    T : Adjointable<B, Yʹ>,
+    Self : Linear<Pair<A, B>>,
+    for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Yʹ, Codomain=R>,
+{
+    type AdjointCodomain = R;
+    type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+
+    fn adjoint(&self) -> Self::Adjoint<'_> {
+        ColOp(self.0.adjoint(), self.1.adjoint())
+    }
+}
+
+impl<A, B, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
+where
+    S : Preadjointable<A, Yʹ>,
+    T : Preadjointable<B, Yʹ>,
+    Self : Linear<Pair<A, B>>,
+    for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, Pair<A,B>, Codomain=R>,
+{
+    type PreadjointCodomain = R;
+    type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+
+    fn preadjoint(&self) -> Self::Preadjoint<'_> {
+        ColOp(self.0.preadjoint(), self.1.preadjoint())
+    }
+}
+
+
+impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
+where
+    S : Adjointable<A, Xʹ>,
+    T : Adjointable<A, Yʹ>,
+    Self : Linear<A>,
+    for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
+{
+    type AdjointCodomain = R;
+    type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+
+    fn adjoint(&self) -> Self::Adjoint<'_> {
+        RowOp(self.0.adjoint(), self.1.adjoint())
+    }
+}
+
+impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
+where
+    S : Preadjointable<A, Xʹ>,
+    T : Preadjointable<A, Yʹ>,
+    Self : Linear<A>,
+    for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, A, Codomain=R>,
+{
+    type PreadjointCodomain = R;
+    type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+
+    fn preadjoint(&self) -> Self::Preadjoint<'_> {
+        RowOp(self.0.preadjoint(), self.1.preadjoint())
+    }
+}
+
+/// Diagonal operator
+pub struct DiagOp<S, T>(pub S, pub T);
+
+impl<A, B, S, T> Apply<Pair<A, B>> for DiagOp<S, T>
+where
+    S : Apply<A>,
+    T : Apply<B>,
+{
+    type Output = Pair<S::Output, T::Output>;
+
+    fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
+        Pair(self.0.apply(a), self.1.apply(b))
+    }
+}
+
+impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for DiagOp<S, T>
+where
+    S : Apply<&'a A>,
+    T : Apply<&'a B>,
+{
+    type Output = Pair<S::Output, T::Output>;
+
+    fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
+        Pair(self.0.apply(a), self.1.apply(b))
+    }
+}
+
+impl<A, B, S, T, D> Linear<Pair<A, B>> for DiagOp<S, T>
+where
+    DiagOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a>  Apply<&'a Pair<A, B>, Output=D>,
+{
+    type Codomain = D;
+}
+
+impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
+where
+    S : GEMV<F, U, A>,
+    T : GEMV<F, V, B>,
+    F : Num,
+    Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>
+{
+    fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &Pair<U, V>, β : F) {
+        self.0.gemv(&mut y.0, α, &x.0, β);
+        self.1.gemv(&mut y.1, α, &x.1, β);
+    }
+
+    fn apply_mut(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
+        self.0.apply_mut(&mut y.0, &x.0);
+        self.1.apply_mut(&mut y.1, &x.1);
+    }
+
+    /// Computes `y += Ax`, where `A` is `Self`
+    fn apply_add(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
+        self.0.apply_add(&mut y.0, &x.0);
+        self.1.apply_add(&mut y.1, &x.1);
+    }
+}
+
+impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
+where
+    S : Adjointable<A, Xʹ>,
+    T : Adjointable<B, Yʹ>,
+    Self : Linear<Pair<A, B>>,
+    for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
+{
+    type AdjointCodomain = R;
+    type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+
+    fn adjoint(&self) -> Self::Adjoint<'_> {
+        DiagOp(self.0.adjoint(), self.1.adjoint())
+    }
+}
+
+impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
+where
+    S : Preadjointable<A, Xʹ>,
+    T : Preadjointable<B, Yʹ>,
+    Self : Linear<Pair<A, B>>,
+    for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, Pair<A, B>, Codomain=R>,
+{
+    type PreadjointCodomain = R;
+    type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+
+    fn preadjoint(&self) -> Self::Preadjoint<'_> {
+        DiagOp(self.0.preadjoint(), self.1.preadjoint())
+    }
+}
+
+/// Block operator
+pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
+