Mercurial > repos > alg_tools / changeset

--- a/src/linops.rs	Sat Sep 06 23:29:34 2025 -0500
+++ b/src/linops.rs	Sun Sep 07 09:44:43 2025 -0500
@@ -156,7 +156,9 @@
     X: Space,
     Yʹ: Space,
 {
+    /// Codomain of the adjoint operator.
     type AdjointCodomain: ClosedSpace;
+    /// Type of the adjoint operator.
     type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
     where
         Self: 'a;
@@ -165,6 +167,29 @@
     fn adjoint(&self) -> Self::Adjoint<'_>;
 }
 
+/// Variant of [`Adjointable`] where the adjoint does not depend on a lifetime parameter.
+/// This exists due to restrictions of Rust's type system: if `A :: Adjointable`, and we make
+/// further restrictions on the adjoint operator, through, e.g.
+/// ```
+/// for<'a> A::Adjoint<'a> : GEMV<F, X, Y>,
+/// ```
+/// Then `'static` lifetime is forced on `X`. Having `A::SimpleAdjoint` not depend on `'a`
+/// avoids this, but makes it impossible for the adjoint to be just a light wrapper around the
+/// original operator.
+pub trait SimplyAdjointable<X, Yʹ>: Linear<X>
+where
+    X: Space,
+    Yʹ: Space,
+{
+    /// Codomain of the adjoint operator.
+    type AdjointCodomain: ClosedSpace;
+    /// Type of the adjoint operator.
+    type SimpleAdjoint: Linear<Yʹ, Codomain = Self::AdjointCodomain>;
+
+    /// Form the adjoint operator of `self`.
+    fn adjoint(&self) -> Self::SimpleAdjoint;
+}
+
 /// Trait for forming a preadjoint of an operator.
 ///
 /// For an operator $A$ this is an operator $A\_\*$
@@ -187,13 +212,6 @@
     fn preadjoint(&self) -> Self::Preadjoint<'_>;
 }
 
-/// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
-pub trait SimplyAdjointable<X: Space>: Adjointable<X, <Self as Mapping<X>>::Codomain> {}
-impl<'a, X: Space, T> SimplyAdjointable<X> for T where
-    T: Adjointable<X, <Self as Mapping<X>>::Codomain>
-{
-}
-
 /// The identity operator
 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
 pub struct IdOp<X>(PhantomData<X>);
@@ -253,6 +271,15 @@
     }
 }
 
+impl<X: Clone + Space> SimplyAdjointable<X, X::Principal> for IdOp<X> {
+    type AdjointCodomain = X::Principal;
+    type SimpleAdjoint = IdOp<X::Principal>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        IdOp::new()
+    }
+}
+
 impl<X: Clone + Space> Preadjointable<X, X::Principal> for IdOp<X> {
     type PreadjointCodomain = X::Principal;
     type Preadjoint<'a>
@@ -497,6 +524,31 @@
     }
 }
 
+impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> SimplyAdjointable<X, Yprime>
+    for ZeroOp<X, Y, OY, OXprime, F>
+where
+    X: HasDual<F, DualSpace = Xprime>,
+    Y: HasDual<F, DualSpace = Yprime>,
+    F: Float,
+    Xprime: ClosedVectorSpace<Field = F>,
+    //Xprime::Owned: AXPY<Field = F>,
+    Yprime: ClosedSpace,
+    OY: OriginGenerator<Y>,
+    OXprime: OriginGenerator<X::DualSpace> + Clone,
+{
+    type AdjointCodomain = Xprime;
+    type SimpleAdjoint = ZeroOp<Yprime, Xprime, OXprime, (), F>;
+    // () means not (pre)adjointable.
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        ZeroOp {
+            codomain_origin_generator: self.other_origin_generator.clone(),
+            other_origin_generator: (),
+            _phantoms: PhantomData,
+        }
+    }
+}
+
 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
 where
     X: Space,
@@ -681,6 +733,27 @@
     }
 }
 
+impl<A, B, Yʹ, S, T> SimplyAdjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
+where
+    A: Space,
+    B: Space,
+    Yʹ: Space,
+    S: SimplyAdjointable<A, Yʹ>,
+    T: SimplyAdjointable<B, Yʹ>,
+    Self: Linear<Pair<A, B>>,
+    // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+    //     Yʹ,
+    //     Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
+    // >,
+{
+    type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
+    type SimpleAdjoint = ColOp<S::SimpleAdjoint, T::SimpleAdjoint>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        ColOp(self.0.adjoint(), self.1.adjoint())
+    }
+}
+
 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
 where
     A: Space,
@@ -728,6 +801,28 @@
     }
 }
 
+impl<A, Xʹ, Yʹ, R, S, T> SimplyAdjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
+where
+    A: Space,
+    Xʹ: Space,
+    Yʹ: Space,
+    R: ClosedSpace + ClosedAdd,
+    S: SimplyAdjointable<A, Xʹ, AdjointCodomain = R>,
+    T: SimplyAdjointable<A, Yʹ, AdjointCodomain = R>,
+    Self: Linear<A>,
+    // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+    //     Pair<Xʹ,Yʹ>,
+    //     Codomain=R,
+    // >,
+{
+    type AdjointCodomain = R;
+    type SimpleAdjoint = RowOp<S::SimpleAdjoint, T::SimpleAdjoint>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        RowOp(self.0.adjoint(), self.1.adjoint())
+    }
+}
+
 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
 where
     A: Space,
@@ -834,6 +929,26 @@
     }
 }
 
+impl<A, B, Xʹ, Yʹ, R, S, T> SimplyAdjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
+where
+    A: Space,
+    B: Space,
+    Xʹ: Space,
+    Yʹ: Space,
+    R: ClosedSpace,
+    S: SimplyAdjointable<A, Xʹ>,
+    T: SimplyAdjointable<B, Yʹ>,
+    Self: Linear<Pair<A, B>>,
+    for<'a> DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
+{
+    type AdjointCodomain = R;
+    type SimpleAdjoint = DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>;
+
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        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
     A: Space,

--- a/src/nalgebra_support.rs	Sat Sep 06 23:29:34 2025 -0500
+++ b/src/nalgebra_support.rs	Sun Sep 07 09:44:43 2025 -0500
@@ -461,6 +461,25 @@
     }
 }
 
+impl<'own, SM, N, M, K, E> SimplyAdjointable<OMatrix<E, M, K>, OMatrix<E, N, K>>
+    for Matrix<E, N, M, SM>
+where
+    SM: Storage<E, N, M>,
+    N: Dim,
+    M: Dim,
+    K: Dim,
+    E: Float + RealField,
+    DefaultAllocator: Allocator<N, K> + Allocator<M, K> + Allocator<M, N>,
+    ShapeConstraint: StridesOk<E, N, K> + StridesOk<E, M, K>,
+{
+    type AdjointCodomain = OMatrix<E, M, K>;
+    type SimpleAdjoint = OMatrix<E, M, N>;
+
+    #[inline]
+    fn adjoint(&self) -> Self::SimpleAdjoint {
+        Matrix::adjoint(self)
+    }
+}
 /// This function is [`nalgebra::EuclideanNorm::metric_distance`] without the `sqrt`.
 #[inline]
 fn metric_distance_squared<T, R1, C1, S1, R2, C2, S2>(