diff -r cebedc4a8331 -r d63e40672dd6 src/linops.rs
--- a/src/linops.rs	Sat Dec 21 14:27:14 2024 -0500
+++ b/src/linops.rs	Fri Dec 20 16:14:17 2024 -0500
@@ -12,53 +12,66 @@
 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm, HasDual};
 
 /// Trait for linear operators on `X`.
-pub trait Linear<X : Space> : Mapping<X>
-{ }
+pub trait Linear<X : AXPY> : Mapping<X, Codomain = Self::LinCodomain> {
+    type LinCodomain : AXPY; // = Self::Codomain, but with AXPY enforced without trait bound bloat.
+}
 
 /// Efficient in-place summation.
-#[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F, X = Self> : Space
-where
-    F : Num,
-    X : Space,
-{
+#[replace_float_literals(Self::Field::cast_from(literal))]
+pub trait AXPY : Space {
+    type Field : Num;
+
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
+    fn axpy<I : Instance<Self>>(&mut self, α : Self::Field, x : I, β : Self::Field);
 
     /// Copies `x` to `self`.
-    fn copy_from<I : Instance<X>>(&mut self, x : I) {
+    fn copy_from<I : Instance<Self>>(&mut self, x : I) {
         self.axpy(1.0, x, 0.0)
     }
 
     /// Computes  `y = αx`, where `y` is `Self`.
-    fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) {
+    fn scale_from<I : Instance<Self>>(&mut self, α : Self::Field, x : I) {
         self.axpy(α, x, 0.0)
     }
 }
 
+macro_rules! impl_axpy {
+    ($($type:ty)*) => { $(
+        impl AXPY for $type {
+            type Field = $type;
+            fn axpy<I : Instance<Self>>(&mut self, α : $type, x : I, β : $type) {
+                *self = *self * α + x.own() * β;
+            }
+        }
+    )* };
+}
+
+impl_axpy!(i8 i16 i32 i64 i128 isize f32 f64);
+
+
 /// 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> {
+#[replace_float_literals(X::Field::cast_from(literal))]
+pub trait GEMV<X : AXPY> : Linear<X> {
     /// Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
+    fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : <Self::LinCodomain as AXPY>::Field, x : I, β : <Self::LinCodomain as AXPY>::Field);
 
     /// Computes `y = Ax`, where `A` is `Self`
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_mut<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){
         self.gemv(y, 1.0, x, 0.0)
     }
 
     /// Computes `y += Ax`, where `A` is `Self`
-    fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_add<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){
         self.gemv(y, 1.0, x, 1.0)
     }
 }
 
 
 /// Bounded linear operators
-pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
+pub trait BoundedLinear<X, XExp, CodExp> : Linear<X>
 where
     F : Num,
-    X : Space + Norm<F, XExp>,
+    X : AXP + Norm<F, XExp>,
     XExp : NormExponent,
     CodExp : NormExponent
 {
@@ -66,7 +79,7 @@
     /// This is not expected to be the norm, just any bound on it that can be
     /// reasonably implemented. The [`NormExponent`] `xexp`  indicates the norm
     /// in `X`, and `codexp` in the codomain.
-    fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F;
+    fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> X::Field; // TODO: Should pick real subfield
 }
 
 // Linear operator application into mutable target. The [`AsRef`] bound
@@ -105,27 +118,26 @@
 ///
 /// We do not set further restrictions on the spacds, to allow preadjointing when `X`
 /// is on the dual space of `Ypre`, but a subset of it.
-pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X>
+pub trait Preadjointable<X : AXPY, Ypre : AXPY> : Linear<X>
 //where
 //    Ypre : HasDual<F, DualSpace=Self::Codomain>,
 {
-    type PreadjointCodomain : Space;
+    type PreadjointCodomain : AXPY;
     type Preadjoint<'a> : Linear<Ypre, Codomain=Self::PreadjointCodomain> where Self : 'a;
 
     /// Form the adjoint operator of `self`.
     fn preadjoint(&self) -> Self::Preadjoint<'_>;
 }
 
-
 /// The identity operator
 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct IdOp<X> (PhantomData<X>);
+pub struct IdOp<X : AXPY> (PhantomData<X>);
 
-impl<X> IdOp<X> {
+impl<X : AXPY> IdOp<X> {
     pub fn new() -> IdOp<X> { IdOp(PhantomData) }
 }
 
-impl<X : Clone + Space> Mapping<X> for IdOp<X> {
+impl<X : Clone + AXPY> Mapping<X> for IdOp<X> {
     type Codomain = X;
 
     fn apply<I : Instance<X>>(&self, x : I) -> X {
@@ -133,32 +145,29 @@
     }
 }
 
-impl<X : Clone + Space> Linear<X> for IdOp<X>
-{ }
+impl<X : Clone + AXPY> Linear<X> for IdOp<X> {
+    type LinCodomain = X;
+}
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
-where
-    Y : AXPY<F, X>,
-    X : Clone + Space
-{
+impl<X : Clone + AXPY> GEMV<X> for IdOp<X> {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
-    fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
+    fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : X::Field, x : I, β : X::Field) {
         y.axpy(α, x, β)
     }
 
-    fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+    fn apply_mut<I : Instance<X>>(&self, y : &mut X, x : I){
         y.copy_from(x);
     }
 }
 
-impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
+impl<X, E> BoundedLinear<X, E, E> for IdOp<X>
 where
-    X : Space + Clone + Norm<F, E>,
-    F : Num,
+    X : AXPY + Clone + Norm<F, E>,
+    F : Num + AXPY,
     E : NormExponent
 {
-    fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
+    fn opnorm_bound(&self, _xexp : E, _codexp : E) -> X::Field { F::ONE }
 }
 
 impl<X : Clone + HasDual<F, DualSpace=X>, F : Float> Adjointable<X, F> for IdOp<X> {
@@ -168,7 +177,7 @@
     fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
 }
 
-impl<X : Clone + Space> Preadjointable<X, X> for IdOp<X> {
+impl<X : Clone + AXPY> Preadjointable<X, X> for IdOp<X> {
     type PreadjointCodomain = X;
     type Preadjoint<'a> = IdOp<X> where X : 'a;
 
@@ -178,10 +187,12 @@
 
 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
 where
-    X : Space,
+    X : AXPY,
     T : Linear<X>,
-    S : Linear<T::Codomain>
-{ }
+    S : Linear<T::LinCodomain>
+{
+    type LinCodomain = S::LinCodomain;
+}
 
 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
 where
@@ -227,17 +238,16 @@
 
 use std::ops::Add;
 
-impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
+impl<A, B, S, T, Y> Mapping<Pair<A, B>> for RowOp<S, T>
 where
     A : Space,
     B : Space,
+    Y : Space,
     S : Mapping<A>,
     T : Mapping<B>,
-    S::Codomain : Add<T::Codomain>,
-    <S::Codomain as Add<T::Codomain>>::Output : Space,
-
+    S::Codomain : Add<T::Codomain, Output = Y>,
 {
-    type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
+    type Codomain = Y;
 
     fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
         let Pair(a, b) = x.decompose();
@@ -245,25 +255,30 @@
     }
 }
 
-impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
+impl<F, A, B, S, T, Y> Linear<Pair<A, B>> for RowOp<S, T>
 where
-    A : Space,
-    B : Space,
+    F : Num + AXPY,
+    A : AXPY<Field=F>,
+    B : AXPY<Field=F>,
+    Y : AXPY,
     S : Linear<A>,
     T : Linear<B>,
-    S::Codomain : Add<T::Codomain>,
-    <S::Codomain as Add<T::Codomain>>::Output : Space,
-{ }
+    S::LinCodomain : Add<T::LinCodomain, Output=Y>,
+{
+    type LinCodomain = Y;
+}
 
 
-impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
+impl<'b, F, S, T, Y, G, U, V> GEMV<Pair<U, V>> for RowOp<S, T>
 where
-    U : Space,
-    V : Space,
-    S : GEMV<F, U, Y>,
-    T : GEMV<F, V, Y>,
-    F : Num,
-    Self : Linear<Pair<U, V>, Codomain=Y>
+    U : AXPY<Field=G>,
+    V : AXPY<Field=G>,
+    Y : AXPY<Field=F>,
+    S : GEMV<U>,
+    T : GEMV<V>,
+    F : Num + AXPY,
+    G : Num + AXPY,
+    S::LinCodomain : Add<T::LinCodomain, Output=Y>,
 {
     fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
         let Pair(u, v) = x.decompose();
@@ -288,33 +303,39 @@
 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
 pub struct ColOp<S, T>(pub S, pub T);
 
-impl<A, S, T> Mapping<A> for ColOp<S, T>
+impl<X, S, T> Mapping<X> for ColOp<S, T>
 where
-    A : Space,
-    S : Mapping<A>,
-    T : Mapping<A>,
+    X : AXPY,
+    S : Mapping<X>,
+    T : Mapping<X>,
 {
     type Codomain = Pair<S::Codomain, T::Codomain>;
 
-    fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
+    fn apply<I : Instance<X>>(&self, a : I) -> Self::Codomain {
         Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
     }
 }
 
-impl<A, S, T> Linear<A> for ColOp<S, T>
+impl<X, S, T, A, B, F> Linear<X> for ColOp<S, T>
 where
-    A : Space,
-    S : Mapping<A>,
-    T : Mapping<A>,
-{ }
+    F : Num + AXPY,
+    X : AXPY,
+    S : Linear<X, LinCodomain=A>,
+    T : Linear<X, LinCodomain=B>,
+    A : AXPY<Field=F>,
+    B : AXPY<Field=F>,
+{
+    type LinCodomain = Pair<S::LinCodomain, T::LinCodomain>;
+}
 
-impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
+impl<F, S, T, A, B, X> GEMV<X> for ColOp<S, T>
 where
-    X : Space,
-    S : GEMV<F, X, A>,
-    T : GEMV<F, X, B>,
-    F : Num,
-    Self : Linear<X, Codomain=Pair<A, B>>
+    F : Num + AXPY,
+    X : AXPY,
+    S : GEMV<X, LinCodomain=A>,
+    T : GEMV<X, LinCodomain=B>,
+    A : AXPY<Field=F>,
+    B : AXPY<Field=F>,
 {
     fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
         self.0.gemv(&mut y.0, α, x.ref_instance(), β);
@@ -341,7 +362,7 @@
     B : HasDual<F>,
     S : Adjointable<A, F>,
     T : Adjointable<B, F>,
-    Yʹ : Space,
+    Yʹ : AXPY,
     S :: Codomain : HasDual<F, DualSpace=Yʹ>,
     T :: Codomain : HasDual<F, DualSpace=Yʹ>,
     S::Codomain : Add<T::Codomain>,
@@ -361,15 +382,15 @@
     }
 }
 
-impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
-where
-    A : Space,
-    B : Space,
-    Yʹ : Space,
+impl<F, A, B, Y, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
+wherea
+    A : AXPY,
+    B : AXPY,
+    Yʹ : AXPY,
     S : Preadjointable<A, Yʹ>,
     T : Preadjointable<B, Yʹ>,
     S::Codomain : Add<T::Codomain>,
-    <S::Codomain as Add<T::Codomain>>::Output : Space,
+    <S::Codomain as Add<T::Codomain>>::Output : AXPY,
     //Self : Linear<Pair<A, B>, Codomain=Y>,
     //for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, F>,
 {
@@ -407,10 +428,10 @@
 
 impl<A, Aʹ, Xʹ, Yʹ, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
 where
-    A : Space,
-    Aʹ : Space,
-    Xʹ : Space,
-    Yʹ : Space,
+    A : AXPY,
+    Aʹ : AXPY,
+    Xʹ : AXPY,
+    Yʹ : AXPY,
     S : Preadjointable<A, Xʹ, PreadjointCodomain=Aʹ>,
     T : Preadjointable<A, Yʹ, PreadjointCodomain=Aʹ>,
     Aʹ : ClosedAdd,
@@ -429,8 +450,8 @@
 
 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
+    A : AXPY,
+    B : AXPY,
     S : Mapping<A>,
     T : Mapping<B>,
 {
@@ -442,26 +463,32 @@
     }
 }
 
-impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
-where
-    A : Space,
-    B : Space,
-    S : Linear<A>,
-    T : Linear<B>,
-{ }
-
-impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
+impl<F, A, B, S, T, G, U, V> Linear<Pair<A, B>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
-    U : Space,
-    V : Space,
-    S : GEMV<F, U, A>,
-    T : GEMV<F, V, B>,
-    F : Num,
-    Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>,
+    F : Num + AXPY,
+    G : Num + AXPY,
+    A : AXPY<Field=F>,
+    B : AXPY<Field=F>,
+    U : AXPY<Field=G>,
+    V : AXPY<Field=G>,
+    S : Linear<A, LinCodomain=U>,
+    T : Linear<B, LinCodomain=V>,
 {
-    fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+    type LinCodomain = Pair<U, V>;
+}
+
+impl<F, A, B, S, T, G, U, V> GEMV<Pair<U, V>> for DiagOp<S, T>
+where
+    F : Num + AXPY,
+    G : Num + AXPY,
+    A : AXPY<Field=F>,
+    B : AXPY<Field=F>,
+    U : AXPY<Field=G>,
+    V : AXPY<Field=G>,
+    S : GEMV<A, LinCodomain=U>,
+    T : GEMV<B, LinCodomain=V>,
+{
+    fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : G, x : I, β : G) {
         let Pair(u, v) = x.decompose();
         self.0.gemv(&mut y.0, α, u, β);
         self.1.gemv(&mut y.1, α, v, β);
@@ -501,10 +528,10 @@
 
 impl<A, B, Xʹ, Yʹ, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
 where
-    A : Space,
-    B : Space,
-    Xʹ : Space,
-    Yʹ : Space,
+    A : AXPY,
+    B : AXPY,
+    Xʹ : AXPY,
+    Yʹ : AXPY,
     S : Preadjointable<A, Xʹ>,
     T : Preadjointable<B, Yʹ>,
 {
@@ -523,16 +550,16 @@
 macro_rules! pairnorm {
     ($expj:ty) => {
         impl<F, A, B, S, T, ExpA, ExpB, ExpR>
-        BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F>
+        BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR>
         for RowOp<S, T>
         where
             F : Float,
-            A : Space + Norm<F, ExpA>,
-            B : Space + Norm<F, ExpB>,
+            A : AXPY + Norm<F, ExpA>,
+            B : AXPY + Norm<F, ExpB>,
             S : BoundedLinear<A, ExpA, ExpR, F>,
             T : BoundedLinear<B, ExpB, ExpR, F>,
             S::Codomain : Add<T::Codomain>,
-            <S::Codomain as Add<T::Codomain>>::Output : Space,
+            <S::Codomain as Add<T::Codomain>>::Output : AXPY,
             ExpA : NormExponent,
             ExpB : NormExponent,
             ExpR : NormExponent,
@@ -541,7 +568,7 @@
                 &self,
                 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
                 expr : ExpR
-            ) -> F {
+            ) -> <Pair<A, B> as AXPY>::Field {
                 // An application of the triangle inequality bounds the norm by the maximum
                 // of the individual norms. A simple observation shows this to be exact.
                 let na = self.0.opnorm_bound(expa, expr);
@@ -550,12 +577,12 @@
             }
         }
         
-        impl<F, A, S, T, ExpA, ExpS, ExpT>
-        BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
+        impl<A, S, T, ExpA, ExpS, ExpT>
+        BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>>
         for ColOp<S, T>
         where
             F : Float,
-            A : Space + Norm<F, ExpA>,
+            A : AXPY + Norm<F, ExpA>,
             S : BoundedLinear<A, ExpA, ExpS, F>,
             T : BoundedLinear<A, ExpA, ExpT, F>,
             ExpA : NormExponent,
@@ -566,7 +593,7 @@
                 &self,
                 expa : ExpA,
                 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj>
-            ) -> F {
+            ) -> A::Field {
                 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
                 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
                 let ns = self.0.opnorm_bound(expa, exps);