diff -r e6829fbe2737 -r afe04e6b4a5b src/discrete_gradient.rs
--- a/src/discrete_gradient.rs	Sun Sep 07 09:51:26 2025 -0500
+++ b/src/discrete_gradient.rs	Sun Sep 07 10:00:14 2025 -0500
@@ -3,7 +3,7 @@
  */
 use crate::error::DynResult;
 use crate::instance::Instance;
-use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, GEMV};
+use crate::linops::{Adjointable, BoundedLinear, Linear, Mapping, SimplyAdjointable, GEMV};
 use crate::norms::{Norm, L2};
 use crate::types::Float;
 use nalgebra::{DVector, Dyn, Matrix, Storage, StorageMut, U1};
@@ -207,14 +207,7 @@
     fn new(dims: &'a [usize; N], d: usize) -> Self {
         let d_stride = dims[0..d].iter().product::<usize>();
         let len = dims.iter().product::<usize>();
-        Iter {
-            dims,
-            d,
-            d_stride,
-            i: [0; N],
-            k: 0,
-            len,
-        }
+        Iter { dims, d, d_stride, i: [0; N], k: 0, len }
     }
 }
 
@@ -336,11 +329,7 @@
     /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
     pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
         if u.len() == dims.iter().product::<usize>() {
-            Some(Grad {
-                dims,
-                h,
-                discretisation,
-            })
+            Some(Grad { dims, h, discretisation })
         } else {
             None
         }
@@ -359,11 +348,7 @@
     /// Creates a new discrete gradient operator for the vector `u`, verifying dimensions.
     pub fn new_for(u: &DVector<F>, h: F, dims: [usize; N], discretisation: B) -> Option<Self> {
         if u.len() == dims.iter().product::<usize>() * N {
-            Some(Div {
-                dims,
-                h,
-                discretisation,
-            })
+            Some(Div { dims, h, discretisation })
         } else {
             None
         }
@@ -425,13 +410,21 @@
     where
         Self: 'a;
 
-    /// Form the adjoint operator of `self`.
     fn adjoint(&self) -> Self::Adjoint<'_> {
-        Div {
-            dims: self.dims,
-            h: -self.h,
-            discretisation: self.discretisation.opposite(),
-        }
+        Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
+    }
+}
+
+impl<F, B, const N: usize> SimplyAdjointable<DVector<F>, DVector<F>> for Grad<F, B, N>
+where
+    B: Discretisation<F>,
+    F: Float + nalgebra::RealField,
+{
+    type AdjointCodomain = DVector<F>;
+    type SimpleAdjoint = Div<F, B::Opposite, N>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        Div { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
     }
 }
 
@@ -446,13 +439,21 @@
     where
         Self: 'a;
 
-    /// Form the adjoint operator of `self`.
     fn adjoint(&self) -> Self::Adjoint<'_> {
-        Grad {
-            dims: self.dims,
-            h: -self.h,
-            discretisation: self.discretisation.opposite(),
-        }
+        Grad { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
+    }
+}
+
+impl<F, B, const N: usize> SimplyAdjointable<DVector<F>, DVector<F>> for Div<F, B, N>
+where
+    B: Discretisation<F>,
+    F: Float + nalgebra::RealField,
+{
+    type AdjointCodomain = DVector<F>;
+    type SimpleAdjoint = Grad<F, B::Opposite, N>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        Grad { dims: self.dims, h: -self.h, discretisation: self.discretisation.opposite() }
     }
 }