# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1734147432 18000
# Node ID 1b3b1687b9ede104853037eaebd2a8acb49e1844
# Parent  5e3c1874797d3a73bec9e98ce8657e10c0467f38
Add direct products (Pair, RowOp, ColOp, DiagOp)

diff -r 5e3c1874797d -r 1b3b1687b9ed src/direct_product.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/direct_product.rs	Fri Dec 13 22:37:12 2024 -0500
@@ -0,0 +1,228 @@
+/*!
+Direct products of the form $A \times B$.
+
+TODO: This could be easily much more generic if `derive_more` could derive arithmetic
+operations on references.
+*/
+
+use core::ops::{Mul,MulAssign,Div,DivAssign,Add,AddAssign,Sub,SubAssign,Neg};
+use std::clone::Clone;
+use serde::{Serialize, Deserialize};
+use crate::types::{Num, Float};
+use crate::{maybe_lifetime, maybe_ref, impl_vectorspace_ops};
+use crate::euclidean::{Dot, Euclidean};
+
+#[derive(Debug,Clone,PartialEq,Eq,Serialize,Deserialize)]
+pub struct Pair<A, B> (pub A, pub B);
+
+impl<A, B> Pair<A,B> {
+    pub fn new(a : A, b : B) -> Pair<A,B> { Pair{ 0 : a, 1 : b } }
+}
+
+impl<A, B> From<(A,B)> for Pair<A,B> {
+    #[inline]
+    fn from((a, b) : (A, B)) -> Pair<A,B> { Pair{ 0 : a, 1 : b } }
+}
+
+macro_rules! impl_binop {
+    ($trait : ident, $fn : ident, $refl:ident, $refr:ident) => {
+        impl_binop!(@doit: $trait, $fn;
+                           maybe_lifetime!($refl, &'l Pair<A,B>),
+                               (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+                           maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
+                               (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+                           $refl, $refr);
+    };
+
+    (@doit: $trait:ident, $fn:ident;
+            $self:ty, ($aself:ty, $bself:ty);
+            $in:ty, ($ain:ty, $bin:ty);
+            $refl:ident, $refr:ident) => {
+        impl<'l, 'r, A, B, Ai, Bi> $trait<$in>
+        for $self
+        where $aself: $trait<$ain>,
+              $bself: $trait<$bin> {
+            type Output = Pair<<$aself as $trait<$ain>>::Output,
+                                       <$bself as $trait<$bin>>::Output>;
+
+            #[inline]
+            fn $fn(self, y : $in) -> Self::Output {
+                Pair { 0 : maybe_ref!($refl, self.0).$fn(maybe_ref!($refr, y.0)),
+                               1 : maybe_ref!($refl, self.1).$fn(maybe_ref!($refr, y.1)) }
+            }
+        }
+    };
+}
+
+macro_rules! impl_assignop {
+    ($trait : ident, $fn : ident, $refr:ident) => {
+        impl_assignop!(@doit: $trait, $fn;
+                              maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
+                                  (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+                              $refr);
+    };
+    (@doit: $trait:ident, $fn:ident;
+            $in:ty, ($ain:ty, $bin:ty);
+            $refr:ident) => {
+        impl<'r, A, B, Ai, Bi> $trait<$in>
+        for Pair<A,B>
+        where A: $trait<$ain>,
+              B: $trait<$bin> {
+            #[inline]
+            fn $fn(&mut self, y : $in) -> () {
+                self.0.$fn(maybe_ref!($refr, y.0));
+                self.1.$fn(maybe_ref!($refr, y.1));
+            }
+        }
+    }
+}
+
+macro_rules! impl_scalarop {
+    ($trait : ident, $fn : ident, $refl:ident) => {
+        impl_scalarop!(@doit: $trait, $fn;
+                              maybe_lifetime!($refl, &'l Pair<A,B>),
+                                  (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+                              $refl);
+    };
+    (@doit: $trait:ident, $fn:ident;
+            $self:ty, ($aself:ty, $bself:ty);
+            $refl:ident) => {
+        // Scalar as Rhs
+        impl<'l, F : Num, A, B> $trait<F>
+        for $self
+        where $aself: $trait<F>,
+              $bself: $trait<F> {
+            type Output = Pair<<$aself as $trait<F>>::Output,
+                                       <$bself as $trait<F>>::Output>;
+            #[inline]
+            fn $fn(self, a : F) -> Self::Output {
+                Pair{ 0 : maybe_ref!($refl, self.0).$fn(a),
+                              1 : maybe_ref!($refl, self.1).$fn(a)}
+            }
+        }
+    }
+}
+
+// Not used due to compiler overflow
+#[allow(unused_macros)]
+macro_rules! impl_scalarlhs_op {
+    ($trait:ident, $fn:ident, $refr:ident, $field:ty) => {
+        impl_scalarlhs_op!(@doit: $trait, $fn, 
+                                  maybe_lifetime!($refr, &'r Pair<Ai,Bi>),
+                                  (maybe_lifetime!($refr, &'r Ai), maybe_lifetime!($refr, &'r Bi));
+                                  $refr, $field);
+    };
+    (@doit: $trait:ident, $fn:ident, 
+            $in:ty, ($ain:ty, $bin:ty);
+            $refr:ident, $field:ty) => {
+        impl<'r, Ai, Bi> $trait<$in>
+        for $field
+        where $field : $trait<$ain>
+                     + $trait<$bin> {
+            type Output = Pair<<$field as $trait<$ain>>::Output,
+                                       <$field as $trait<$bin>>::Output>;
+            #[inline]
+            fn $fn(self, x : $in) -> Self::Output {
+                Pair{ 0 : self.$fn(maybe_ref!($refr, x.0)),
+                              1 : self.$fn(maybe_ref!($refr, x.1))}
+            }
+        }
+    };
+}
+
+macro_rules! impl_scalar_assignop {
+    ($trait : ident, $fn : ident) => {
+        impl<'r, F : Num, A, B> $trait<F>
+        for Pair<A,B>
+        where A: $trait<F>, B: $trait<F> {
+            #[inline]
+            fn $fn(&mut self, a : F) -> () {
+                self.0.$fn(a);
+                self.1.$fn(a);
+            }
+        }
+    }
+}
+
+macro_rules! impl_unaryop {
+    ($trait:ident, $fn:ident, $refl:ident) => {
+        impl_unaryop!(@doit: $trait, $fn;
+                             maybe_lifetime!($refl, &'l Pair<A,B>),
+                                 (maybe_lifetime!($refl, &'l A), maybe_lifetime!($refl, &'l B));
+                             $refl);
+    };
+    (@doit: $trait:ident, $fn:ident;
+            $self:ty, ($aself:ty, $bself:ty);
+            $refl : ident) => {
+        impl<'l, A, B> $trait
+        for $self
+        where $aself: $trait,
+              $bself: $trait {
+            type Output = Pair<<$aself as $trait>::Output,
+                                       <$bself as $trait>::Output>;
+            #[inline]
+            fn $fn(self) -> Self::Output {
+                Pair{ 0 : maybe_ref!($refl, self.0).$fn(),
+                              1 : maybe_ref!($refl, self.1).$fn()}
+            }
+        }
+    }
+}
+
+
+impl_vectorspace_ops!(impl_binop, impl_assignop, impl_scalarop, impl_scalarlhs_op,
+                      impl_scalar_assignop, impl_unaryop);
+
+
+impl<A, B, U, V, F> Dot<Pair<U, V>, F> for Pair<A, B>
+where
+    A : Dot<U, F>,
+    B : Dot<V, F>,
+    F : Num
+{
+
+    fn dot(&self, Pair(ref u, ref v) : &Pair<U, V>) -> F {
+        self.0.dot(u) + self.1.dot(v)
+    }
+}
+
+impl<A, B, F> Euclidean<F>  for Pair<A, B>
+where
+    A : Euclidean<F>,
+    B : Euclidean<F>,
+    F : Float
+{
+    type Output = Pair<<A as Euclidean<F>>::Output, <B as Euclidean<F>>::Output>;
+
+    fn similar_origin(&self) -> <Self as Euclidean<F>>::Output {
+        Pair(self.0.similar_origin(), self.1.similar_origin())
+    }
+
+    fn dist2_squared(&self, Pair(ref u, ref v) : &Self) -> F {
+        self.0.dist2_squared(u) + self.1.dist2_squared(v)
+    }
+}
+
+use crate::linops::AXPY;
+
+impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B>
+where
+    A : AXPY<F, U>,
+    B : AXPY<F, V>,
+    F : Num
+{
+    fn axpy(&mut self, α : F, x : &Pair<U, V>, β : F) {
+        self.0.axpy(α, &x.0, β);
+        self.1.axpy(α, &x.1, β);
+    }
+
+    fn copy_from(&mut self, x : &Pair<U, V>,) {
+        self.0.copy_from(&x.0);
+        self.1.copy_from(&x.1);
+    }
+
+    fn scale_from(&mut self, α : F, x : &Pair<U, V>) {
+        self.0.scale_from(α, &x.0);
+        self.1.scale_from(α, &x.1);
+    }
+}
\ No newline at end of file
diff -r 5e3c1874797d -r 1b3b1687b9ed src/lib.rs
--- a/src/lib.rs	Wed Dec 11 20:45:17 2024 -0500
+++ b/src/lib.rs	Fri Dec 13 22:37:12 2024 -0500
@@ -40,6 +40,7 @@
 pub mod fe_model;
 pub mod bisection_tree;
 pub mod nalgebra_support;
-
+pub(crate) mod metaprogramming;
+pub mod direct_product;
 
 pub use types::*;
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>>;
+
diff -r 5e3c1874797d -r 1b3b1687b9ed src/metaprogramming.rs
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/metaprogramming.rs	Fri Dec 13 22:37:12 2024 -0500
@@ -0,0 +1,103 @@
+/*!
+Metaprogramming tools
+*/
+
+/// Reference `x` if so indicated by the first parameter.
+/// Typically to be used from another macro. See the implementation of
+/// [power][crate::vectorspace::powerspace] and [product spaces][crate::vectorspace::productspace].
+///
+/// ```ignore
+/// maybe_ref!(ref, V)   // ➡ &V
+/// maybe_ref!(noref, V) // ➡ V
+/// ```
+#[macro_export]
+macro_rules! maybe_ref {
+    (ref, $x:expr) => { &$x };
+    (noref, $x:expr) => { $x };
+    (ref, $x:ty) => { &$x };
+    (noref, $x:ty) => { $x };
+}
+
+/// Choose `a` if first argument is the literal `ref`, otherwise `b`.
+#[macro_export]
+macro_rules! ifref {
+    (noref, $a:expr, $b:expr) => { $b };
+    (ref, $a:expr, $b:expr) => { $a };
+}
+
+
+/// Annotate `x` with a lifetime if the first parameter
+/// Typically to be used from another macro. See the implementation of
+/// [power][crate::vectorspace::powerspace] and [product spaces][crate::vectorspace::productspace].
+///
+/// ```ignore
+/// maybe_ref!(ref, &'a V)    // ➡ &'a V
+/// maybe_ref!(noref, &'a V)  // ➡ V
+/// ```
+#[macro_export]
+macro_rules! maybe_lifetime {
+    (ref, $x:ty) => { $x };
+    (noref, &$lt:lifetime $x:ty) => { $x };
+    (noref, &$x:ty) => { $x };
+}
+
+/// Use as
+/// ```ignore
+/// impl_vectorspace_ops!(impl_binop, impl_assignop, impl_scalarop, impl_scalar_assignop, impl_unaryop);
+/// ```
+/// with `impl_binop`, `impl_assignop`, `impl_scalarop`, and `impl_unaryop` macros provided.
+/// For example usage see the [power][crate::vectorspace::powerspace] and
+/// [product spaces][crate::vectorspace::productspace] implementations.
+#[macro_export]
+macro_rules! impl_vectorspace_ops {
+    ($impl_binop:ident, $impl_assignop:ident, $impl_scalarop:ident, $impl_scalarlhs_op:ident,
+        $impl_scalar_assignop:ident, $impl_unaryop:ident) => {
+        impl_vectorspace_ops!($impl_binop, $impl_assignop, $impl_scalarop, $impl_scalarlhs_op,
+                                $impl_scalar_assignop, $impl_unaryop,
+                                (f32, f64,
+                                num::complex::Complex<f32>,
+                                num::complex::Complex<f64>));
+    };
+    ($impl_binop:ident, $impl_assignop:ident, $impl_scalarop:ident, $impl_scalarlhs_op:ident,
+        $impl_scalar_assignop:ident, $impl_unaryop:ident, ($($field:ty),+)) => {
+        impl_vectorspace_ops!(@binary, $impl_binop, Add, add);
+        impl_vectorspace_ops!(@binary, $impl_binop, Sub, sub);
+        impl_vectorspace_ops!(@assign, $impl_assignop, AddAssign, add_assign);
+        impl_vectorspace_ops!(@assign, $impl_assignop, SubAssign, sub_assign);
+        impl_vectorspace_ops!(@scalar, $impl_scalarop, Mul, mul);
+        impl_vectorspace_ops!(@scalar, $impl_scalarop, Div, div);
+        // Compiler overflow
+        // $(
+        //     impl_vectorspace_ops!(@scalar_lhs, $impl_scalarlhs_op, Mul, mul, $field);
+        // )*
+        impl_vectorspace_ops!(@scalar_assign, $impl_scalar_assignop, MulAssign, mul_assign);
+        impl_vectorspace_ops!(@scalar_assign, $impl_scalar_assignop, DivAssign, div_assign);
+        impl_vectorspace_ops!(@unary, $impl_unaryop, Neg, neg);
+    };
+    (@binary, $impl:ident, $trait : ident, $fn : ident) => {
+        $impl!($trait, $fn, ref, ref);
+        $impl!($trait, $fn, ref, noref);
+        $impl!($trait, $fn, noref, ref);
+        $impl!($trait, $fn, noref, noref);
+    };
+    (@assign, $impl:ident, $trait : ident, $fn :ident) => {
+        $impl!($trait, $fn, ref);
+        $impl!($trait, $fn, noref);
+    };
+    (@scalar, $impl:ident, $trait : ident, $fn :ident) => {
+        $impl!($trait, $fn, ref);
+        $impl!($trait, $fn, noref);
+    };
+    (@scalar_lhs, $impl:ident, $trait : ident, $fn : ident, $field: ty) => {
+        // These operators need workarounds
+        $impl!($trait, $fn, ref, $field);
+        $impl!($trait, $fn, noref, $field);
+    };
+    (@scalar_assign, $impl:ident, $trait : ident, $fn :ident) => {
+        $impl!($trait, $fn);
+    };
+    (@unary, $impl:ident, $trait : ident, $fn :ident) => {
+        $impl!($trait, $fn, ref);
+        $impl!($trait, $fn, noref);
+    };
+}