--- a/src/linops.rs Sun Apr 27 20:29:43 2025 -0500
+++ b/src/linops.rs Fri May 15 14:46:30 2026 -0500
@@ -2,81 +2,143 @@
Abstract linear operators.
*/
-use numeric_literals::replace_float_literals;
-use std::marker::PhantomData;
-use serde::Serialize;
+use crate::direct_product::Pair;
+use crate::error::DynResult;
+use crate::euclidean::StaticEuclidean;
+use crate::instance::Instance;
+pub use crate::mapping::{ClosedSpace, Composition, DifferentiableImpl, Mapping, Space};
+use crate::norms::{HasDual, Linfinity, NormExponent, PairNorm, L1, L2};
use crate::types::*;
-pub use crate::mapping::{Mapping, Space, Composition};
-use crate::direct_product::Pair;
-use crate::instance::Instance;
-use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm};
+use numeric_literals::replace_float_literals;
+use serde::Serialize;
+use std::marker::PhantomData;
+use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
/// Trait for linear operators on `X`.
-pub trait Linear<X : Space> : Mapping<X>
-{ }
+pub trait Linear<X: Space>: Mapping<X> {}
+
+// impl<X: Space, A: Linear<X>> DifferentiableImpl<X> for A {
+// type Derivative = <Self as Mapping<X>>::Codomain;
+
+// /// Compute the differential of `self` at `x`, consuming the input.
+// fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
+// self.apply(x)
+// }
+// }
+
+/// Vector spaces
+#[replace_float_literals(Self::Field::cast_from(literal))]
+pub trait VectorSpace:
+ Space<Principal = Self::PrincipalV>
+ + Mul<Self::Field, Output = Self::PrincipalV>
+ + Div<Self::Field, Output = Self::PrincipalV>
+ + Add<Self, Output = Self::PrincipalV>
+ + Add<Self::PrincipalV, Output = Self::PrincipalV>
+ + Sub<Self, Output = Self::PrincipalV>
+ + Sub<Self::PrincipalV, Output = Self::PrincipalV>
+ + Neg
+ + for<'b> Add<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
+ + for<'b> Sub<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
+{
+ /// Underlying scalar field
+ type Field: Num;
+
+ /// Principal form of the space; always equal to [`Space::Principal`], but with
+ /// more traits guaranteed.
+ ///
+ /// `PrincipalV` is only assumed to be `AXPY` for itself, as [`AXPY`]
+ /// uses [`Instance`] to apply all other variants and avoid problems
+ /// of choosing multiple implementations of the trait.
+ type PrincipalV: ClosedSpace
+ + AXPY<
+ Self::PrincipalV,
+ Field = Self::Field,
+ PrincipalV = Self::PrincipalV,
+ OwnedVariant = Self::PrincipalV,
+ Principal = Self::PrincipalV,
+ >;
+
+ /// Return a similar zero as `self`.
+ fn similar_origin(&self) -> Self::PrincipalV;
+ // {
+ // self.make_origin_generator().make_origin()
+ // }
+
+ /// Return a similar zero as `x`.
+ fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::PrincipalV {
+ x.eval(|xr| xr.similar_origin())
+ }
+}
/// Efficient in-place summation.
-#[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F, X = Self> : Space + std::ops::MulAssign<F>
+#[replace_float_literals(Self::Field::cast_from(literal))]
+pub trait AXPY<X = Self>:
+ VectorSpace
+ + MulAssign<Self::Field>
+ + DivAssign<Self::Field>
+ + AddAssign<Self>
+ + AddAssign<Self::PrincipalV>
+ + SubAssign<Self>
+ + SubAssign<Self::PrincipalV>
+ + for<'b> AddAssign<&'b Self>
+ + for<'b> SubAssign<&'b Self>
where
- F : Num,
- X : Space,
+ X: Space,
{
- type Owned : AXPY<F, X>;
-
/// Computes `y = βy + αx`, where `y` is `Self`.
- fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
+ fn axpy<I: Instance<X>>(&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<X>>(&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<X>>(&mut self, α: Self::Field, x: I) {
self.axpy(α, x, 0.0)
}
- /// Return a similar zero as `self`.
- fn similar_origin(&self) -> Self::Owned;
-
/// Set self to zero.
fn set_zero(&mut self);
}
+pub trait ClosedVectorSpace: Instance<Self> + VectorSpace<PrincipalV = Self> {}
+impl<X: Instance<X> + VectorSpace<PrincipalV = Self>> ClosedVectorSpace for X {}
+
/// 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> {
+pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: 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 Y, α: F, x: I, β: F);
#[inline]
/// 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 Y, x: I) {
self.gemv(y, 1.0, x, 0.0)
}
#[inline]
/// 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 Y, 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, F = f64>: Linear<X>
where
- F : Num,
- X : Space + Norm<F, XExp>,
- XExp : NormExponent,
- CodExp : NormExponent
+ F: Num,
+ X: Space,
+ XExp: NormExponent,
+ CodExp: NormExponent,
{
/// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
/// 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;
+ ///
+ /// This may fail with an error if the bound is for some reason incalculable.
+ fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> DynResult<F>;
}
// Linear operator application into mutable target. The [`AsRef`] bound
@@ -90,18 +152,45 @@
}*/
/// Trait for forming the adjoint operator of `Self`.
-pub trait Adjointable<X, Yʹ> : Linear<X>
+pub trait Adjointable<X, Yʹ>: Linear<X>
where
- X : Space,
- Yʹ : Space,
+ X: Space,
+ Yʹ: Space,
{
- type AdjointCodomain : Space;
- type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
+ /// Codomain of the adjoint operator.
+ type AdjointCodomain: ClosedSpace;
+ /// Type of the adjoint operator.
+ type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
+ where
+ Self: 'a;
/// Form the adjoint operator of `self`.
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, Z, 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\_\*$
@@ -112,404 +201,667 @@
/// We do not make additional restrictions on `Self::Preadjoint` (in particular, it
/// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
/// have the full space as the codomain, etc.
-pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> {
- type PreadjointCodomain : Space;
- type Preadjoint<'a> : Linear<
- Ypre, Codomain=Self::PreadjointCodomain
- > where Self : 'a;
+pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
+ Linear<X>
+{
+ type PreadjointCodomain: ClosedSpace;
+ type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain>
+ where
+ Self: 'a;
/// Form the adjoint operator of `self`.
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>);
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
+pub struct IdOp<X>(PhantomData<X>);
impl<X> IdOp<X> {
- pub fn new() -> IdOp<X> { IdOp(PhantomData) }
+ pub fn new() -> IdOp<X> {
+ IdOp(PhantomData)
+ }
}
-impl<X : Clone + Space> Mapping<X> for IdOp<X> {
- type Codomain = X;
+impl<X: Space> Mapping<X> for IdOp<X> {
+ type Codomain = X::Principal;
- fn apply<I : Instance<X>>(&self, x : I) -> X {
+ fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
x.own()
}
}
-impl<X : Clone + Space> Linear<X> for IdOp<X>
-{ }
+impl<X: Space> Linear<X> for IdOp<X> {}
#[replace_float_literals(F::cast_from(literal))]
-impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
+impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
where
- Y : AXPY<F, X>,
- X : Clone + Space
+ Y: AXPY<X, Field = F>,
+ X: Space,
{
// 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 Y, α: F, x: I, β: F) {
y.axpy(α, x, β)
}
- fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
+ fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
y.copy_from(x);
}
}
impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
where
- X : Space + Clone + Norm<F, E>,
- F : Num,
- E : NormExponent
+ X: Space + Clone,
+ F: Num,
+ E: NormExponent,
{
- fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
-}
-
-impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> {
- type AdjointCodomain=X;
- type Adjoint<'a> = IdOp<X> where X : 'a;
-
- fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
+ fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> {
+ Ok(F::ONE)
+ }
}
-impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> {
- type PreadjointCodomain=X;
- type Preadjoint<'a> = IdOp<X> where X : 'a;
+impl<X: Clone + Space> Adjointable<X, X::Principal> for IdOp<X> {
+ type AdjointCodomain = X::Principal;
+ type Adjoint<'a>
+ = IdOp<X::Principal>
+ where
+ X: 'a;
- fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
+ fn adjoint(&self) -> Self::Adjoint<'_> {
+ IdOp::new()
+ }
}
+impl<X: Clone + Space> SimplyAdjointable<X, X::Principal> for IdOp<X> {
+ type AdjointCodomain = X::Principal;
+ type SimpleAdjoint = IdOp<X::Principal>;
-/// The zero operator
-#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
-pub struct ZeroOp<'a, X, XD, Y, F> {
- zero : &'a Y, // TODO: don't pass this in `new`; maybe not even store.
- dual_or_predual_zero : XD,
- _phantoms : PhantomData<(X, Y, F)>,
-}
-
-// TODO: Need to make Zero in Instance.
-
-impl<'a, F : Num, X : Space, XD, Y : Space + Clone> ZeroOp<'a, X, XD, Y, F> {
- pub fn new(zero : &'a Y, dual_or_predual_zero : XD) -> Self {
- ZeroOp{ zero, dual_or_predual_zero, _phantoms : PhantomData }
+ fn adjoint(&self) -> Self::SimpleAdjoint {
+ IdOp::new()
}
}
-impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> {
- type Codomain = Y;
+impl<X: Clone + Space> Preadjointable<X, X::Principal> for IdOp<X> {
+ type PreadjointCodomain = X::Principal;
+ type Preadjoint<'a>
+ = IdOp<X::Principal>
+ where
+ X: 'a;
- fn apply<I : Instance<X>>(&self, _x : I) -> Y {
- self.zero.clone()
+ fn preadjoint(&self) -> Self::Preadjoint<'_> {
+ IdOp::new()
}
}
-impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F>
-{ }
+/// The zero operator from a space to itself
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
+pub struct SimpleZeroOp;
+
+impl<X: VectorSpace> Mapping<X> for SimpleZeroOp {
+ type Codomain = X::PrincipalV;
+
+ fn apply<I: Instance<X>>(&self, x: I) -> X::PrincipalV {
+ X::similar_origin_inst(x)
+ }
+}
+
+impl<X: VectorSpace> Linear<X> for SimpleZeroOp {}
#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F>
+impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp
where
- F : Num,
- Y : AXPY<F, Y> + Clone,
- X : Space
+ F: Num,
+ Y: AXPY<Field = F>,
+ X: VectorSpace<Field = F> + Instance<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 Y, _α: F, _x: I, β: F) {
*y *= β;
}
- fn apply_mut<I : Instance<X>>(&self, y : &mut Y, _x : I){
+ fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
y.set_zero();
}
}
-impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F>
+impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
+where
+ F: Num,
+ X: VectorSpace<Field = F>,
+ E1: NormExponent,
+ E2: NormExponent,
+{
+ fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
+ Ok(F::ZERO)
+ }
+}
+
+impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
where
- X : Space + Norm<F, E1>,
- Y : AXPY<F> + Clone + Norm<F, E2>,
- F : Num,
- E1 : NormExponent,
- E2 : NormExponent,
+ F: Num,
+ X: VectorSpace<Field = F> + HasDual<F>,
+ X::DualSpace: ClosedVectorSpace,
{
- fn opnorm_bound(&self, _xexp : E1, _codexp : E2) -> F { F::ZERO }
+ type AdjointCodomain = X::DualSpace;
+ type Adjoint<'b>
+ = SimpleZeroOp
+ where
+ Self: 'b;
+ // () means not (pre)adjointable.
+
+ fn adjoint(&self) -> Self::Adjoint<'_> {
+ SimpleZeroOp
+ }
+}
+
+pub trait OriginGenerator<Y: VectorSpace> {
+ type Ref<'b>: OriginGenerator<Y>
+ where
+ Self: 'b;
+
+ fn origin(&self) -> Y::PrincipalV;
+ fn as_ref(&self) -> Self::Ref<'_>;
}
-impl<'a, F : Num, X, XD, Y, Yprime : Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F>
-where
- X : Space,
- Y : AXPY<F> + Clone + 'static,
- XD : AXPY<F> + Clone + 'static,
+#[derive(Copy, Clone, Debug)]
+pub struct StaticEuclideanOriginGenerator;
+
+impl<Y: StaticEuclidean<F, Field = F>, F: Float> OriginGenerator<Y>
+ for StaticEuclideanOriginGenerator
{
- type AdjointCodomain = XD;
- type Adjoint<'b> = ZeroOp<'b, Yprime, (), XD, F> where Self : 'b;
- // () means not (pre)adjointable.
+ type Ref<'b>
+ = Self
+ where
+ Self: 'b;
+
+ #[inline]
+ fn origin(&self) -> Y::PrincipalV {
+ return Y::origin();
+ }
+
+ #[inline]
+ fn as_ref(&self) -> Self::Ref<'_> {
+ *self
+ }
+}
+
+impl<Y: VectorSpace> OriginGenerator<Y> for Y {
+ type Ref<'b>
+ = &'b Y
+ where
+ Self: 'b;
+
+ #[inline]
+ fn origin(&self) -> Y::PrincipalV {
+ return self.similar_origin();
+ }
+
+ #[inline]
+ fn as_ref(&self) -> Self::Ref<'_> {
+ self
+ }
+}
- fn adjoint(&self) -> Self::Adjoint<'_> {
- ZeroOp::new(&self.dual_or_predual_zero, ())
+impl<'b, Y: VectorSpace> OriginGenerator<Y> for &'b Y {
+ type Ref<'c>
+ = Self
+ where
+ Self: 'c;
+
+ #[inline]
+ fn origin(&self) -> Y::PrincipalV {
+ return self.similar_origin();
+ }
+
+ #[inline]
+ fn as_ref(&self) -> Self::Ref<'_> {
+ self
+ }
+}
+
+/// A zero operator that can be eitherh dualised or predualised (once).
+/// This is achieved by storing an oppropriate zero.
+pub struct ZeroOp<X, Y: VectorSpace<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> {
+ codomain_origin_generator: OY,
+ other_origin_generator: O,
+ _phantoms: PhantomData<(X, Y, F)>,
+}
+
+impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F>
+where
+ OY: OriginGenerator<Y>,
+ X: VectorSpace<Field = F>,
+ Y: VectorSpace<Field = F>,
+ F: Float,
+{
+ pub fn new(y_og: OY) -> Self {
+ ZeroOp {
+ codomain_origin_generator: y_og,
+ other_origin_generator: (),
+ _phantoms: PhantomData,
+ }
}
}
-impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F>
+impl<X, Y, OY, OXprime, Xprime, Yprime, F> ZeroOp<X, Y, OY, OXprime, F>
+where
+ OY: OriginGenerator<Y>,
+ OXprime: OriginGenerator<Xprime>,
+ X: HasDual<F, DualSpace = Xprime>,
+ Y: HasDual<F, DualSpace = Yprime>,
+ F: Float,
+ Xprime: VectorSpace<Field = F, PrincipalV = Xprime>,
+ Xprime::PrincipalV: AXPY<Field = F>,
+ Yprime: Space + Instance<Yprime>,
+{
+ pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self {
+ ZeroOp {
+ codomain_origin_generator: y_og,
+ other_origin_generator: xprime_og,
+ _phantoms: PhantomData,
+ }
+ }
+}
+
+impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F>
+where
+ X: Space,
+ Y: VectorSpace<Field = F>,
+ F: Float,
+ OY: OriginGenerator<Y>,
+{
+ type Codomain = Y::PrincipalV;
+
+ fn apply<I: Instance<X>>(&self, _x: I) -> Y::PrincipalV {
+ self.codomain_origin_generator.origin()
+ }
+}
+
+impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
+where
+ X: Space,
+ Y: VectorSpace<Field = F>,
+ F: Float,
+ OY: OriginGenerator<Y>,
+{
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F>
+where
+ X: Space,
+ Y: AXPY<Field = F, PrincipalV = Y>,
+ F: Float,
+ OY: OriginGenerator<Y>,
+{
+ // Computes `y = αAx + βy`, where `A` is `Self`.
+ fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
+ *y *= β;
+ }
+
+ fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
+ y.set_zero();
+ }
+}
+
+impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
where
- F : Num,
- X : Space,
- Y : AXPY<F> + Clone,
- Ypre : Space,
- XD : AXPY<F> + Clone + 'static,
+ X: Space + Instance<X>,
+ Y: VectorSpace<Field = F>,
+ Y::PrincipalV: Clone,
+ F: Float,
+ E1: NormExponent,
+ E2: NormExponent,
+ OY: OriginGenerator<Y>,
+{
+ fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
+ Ok(F::ZERO)
+ }
+}
+
+impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> Adjointable<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>,
{
- type PreadjointCodomain = XD;
- type Preadjoint<'b> = ZeroOp<'b, Ypre, (), XD, F> where Self : 'b;
- // () means not (pre)adjointable.
+ type AdjointCodomain = Xprime;
+ type Adjoint<'c>
+ = ZeroOp<Yprime, Xprime, OXprime::Ref<'c>, (), F>
+ where
+ Self: 'c;
+ // () means not (pre)adjointable.
+
+ fn adjoint(&self) -> Self::Adjoint<'_> {
+ ZeroOp {
+ codomain_origin_generator: self.other_origin_generator.as_ref(),
+ other_origin_generator: (),
+ _phantoms: PhantomData,
+ }
+ }
+}
- fn preadjoint(&self) -> Self::Preadjoint<'_> {
- ZeroOp::new(&self.dual_or_predual_zero, ())
+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,
- T : Linear<X>,
- S : Linear<T::Codomain>
-{ }
+ X: Space,
+ T: Linear<X>,
+ S: Linear<T::Codomain>,
+{
+}
impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
where
- F : Num,
- X : Space,
- T : Linear<X>,
- S : GEMV<F, T::Codomain, Y>,
+ F: Num,
+ X: Space,
+ T: Linear<X>,
+ S: GEMV<F, T::Codomain, Y>,
{
- fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
+ fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
self.outer.gemv(y, α, self.inner.apply(x), β)
}
/// 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 Y, x: I) {
self.outer.apply_mut(y, self.inner.apply(x))
}
/// 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 Y, x: I) {
self.outer.apply_add(y, self.inner.apply(x))
}
}
impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
where
- F : Num,
- X : Space + Norm<F, Xexp>,
- Z : Space + Norm<F, Zexp>,
- Xexp : NormExponent,
- Yexp : NormExponent,
- Zexp : NormExponent,
- T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>,
- S : BoundedLinear<Z, Zexp, Yexp, F>,
+ F: Num,
+ X: Space,
+ Z: Space,
+ Xexp: NormExponent,
+ Yexp: NormExponent,
+ Zexp: NormExponent,
+ T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>,
+ S: BoundedLinear<Z, Zexp, Yexp, F>,
{
- fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F {
+ fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> DynResult<F> {
let zexp = self.intermediate_norm_exponent;
- self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp)
+ Ok(self.outer.opnorm_bound(zexp, yexp)? * self.inner.opnorm_bound(xexp, zexp)?)
}
}
/// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct RowOp<S, T>(pub S, pub T);
-use std::ops::Add;
-
impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
where
- A : Space,
- B : Space,
- S : Mapping<A>,
- T : Mapping<B>,
- S::Codomain : Add<T::Codomain>,
- <S::Codomain as Add<T::Codomain>>::Output : Space,
-
+ A: Space,
+ B: Space,
+ S: Mapping<A>,
+ T: Mapping<B>,
+ S::Codomain: Add<T::Codomain>,
+ <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
{
type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
- fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
- let Pair(a, b) = x.decompose();
- self.0.apply(a) + self.1.apply(b)
+ fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
+ x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b))
}
}
impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
where
- A : Space,
- B : Space,
- S : Linear<A>,
- T : Linear<B>,
- S::Codomain : Add<T::Codomain>,
- <S::Codomain as Add<T::Codomain>>::Output : Space,
-{ }
-
+ A: Space,
+ B: Space,
+ S: Linear<A>,
+ T: Linear<B>,
+ S::Codomain: Add<T::Codomain>,
+ <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
+{
+}
impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> 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: Space,
+ V: Space,
+ S: GEMV<F, U, Y>,
+ T: GEMV<F, V, Y>,
+ F: Num,
+ Self: Linear<Pair<U, V>, Codomain = Y>,
{
- fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
- let Pair(u, v) = x.decompose();
- self.0.gemv(y, α, u, β);
- self.1.gemv(y, α, v, F::ONE);
+ fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Y, α: F, x: I, β: F) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.gemv(y, α, u, β);
+ self.1.gemv(y, α, v, F::ONE);
+ })
}
- fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
- let Pair(u, v) = x.decompose();
- self.0.apply_mut(y, u);
- self.1.apply_add(y, v);
+ fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.apply_mut(y, u);
+ self.1.apply_add(y, v);
+ })
}
/// Computes `y += Ax`, where `A` is `Self`
- fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
- let Pair(u, v) = x.decompose();
- self.0.apply_add(y, u);
- self.1.apply_add(y, v);
+ fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.apply_add(y, u);
+ self.1.apply_add(y, v);
+ })
}
}
/// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct ColOp<S, T>(pub S, pub T);
impl<A, S, T> Mapping<A> for ColOp<S, T>
where
- A : Space,
- S : Mapping<A>,
- T : Mapping<A>,
+ A: Space,
+ S: Mapping<A>,
+ T: Mapping<A>,
{
type Codomain = Pair<S::Codomain, T::Codomain>;
- fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
- Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
+ fn apply<I: Instance<A>>(&self, a: I) -> Self::Codomain {
+ Pair(a.eval_ref(|r| self.0.apply(r)), self.1.apply(a))
}
}
impl<A, S, T> Linear<A> for ColOp<S, T>
where
- A : Space,
- S : Mapping<A>,
- T : Mapping<A>,
-{ }
+ A: Space,
+ S: Mapping<A>,
+ T: Mapping<A>,
+{
+}
impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> 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>>
+ X: Space,
+ S: GEMV<F, X, A>,
+ T: GEMV<F, X, B>,
+ F: Num,
+ Self: Linear<X, Codomain = Pair<A, B>>,
{
- fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
- self.0.gemv(&mut y.0, α, x.ref_instance(), β);
+ fn gemv<I: Instance<X>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
+ x.eval_ref(|r| self.0.gemv(&mut y.0, α, r, β));
self.1.gemv(&mut y.1, α, x, β);
}
- fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
- self.0.apply_mut(&mut y.0, x.ref_instance());
+ fn apply_mut<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
+ x.eval_ref(|r| self.0.apply_mut(&mut y.0, r));
self.1.apply_mut(&mut y.1, x);
}
/// Computes `y += Ax`, where `A` is `Self`
- fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
- self.0.apply_add(&mut y.0, x.ref_instance());
+ fn apply_add<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
+ x.eval_ref(|r| self.0.apply_add(&mut y.0, r));
self.1.apply_add(&mut y.1, x);
}
}
-
-impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
+impl<A, B, Yʹ, S, T> Adjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
where
- A : Space,
- B : Space,
- Yʹ : Space,
- S : Adjointable<A, Yʹ>,
- T : Adjointable<B, Yʹ>,
- Self : Linear<Pair<A, B>>,
+ A: Space,
+ B: Space,
+ Yʹ: Space,
+ 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=Pair<S::AdjointCodomain, T::AdjointCodomain>
// >,
{
type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
- type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+ 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ʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
+impl<A, B, Yʹ, S, T> SimplyAdjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
where
- A : Space,
- B : Space,
- Yʹ : Space,
- S : Preadjointable<A, Yʹ>,
- T : Preadjointable<B, Yʹ>,
- Self : Linear<Pair<A, B>>,
- for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear<
- Yʹ, Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>,
- >,
+ 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,
+ B: Space,
+ Yʹ: Space,
+ S: Preadjointable<A, Yʹ>,
+ T: Preadjointable<B, Yʹ>,
+ Self: Linear<Pair<A, B>>,
+ for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>:
+ Linear<Yʹ, Codomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>>,
{
type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
- type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+ 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>
+impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
where
- A : Space,
- Xʹ : Space,
- Yʹ : Space,
- R : Space + ClosedAdd,
- S : Adjointable<A, Xʹ, AdjointCodomain = R>,
- T : Adjointable<A, Yʹ, AdjointCodomain = R>,
- Self : Linear<A>,
+ A: Space,
+ Xʹ: Space,
+ Yʹ: Space,
+ R: ClosedSpace + ClosedAdd,
+ S: Adjointable<A, Xʹ, AdjointCodomain = R>,
+ T: Adjointable<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 Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
+ 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>
+impl<A, Xʹ, Yʹ, R, S, T> SimplyAdjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
where
- A : Space,
- Xʹ : Space,
- Yʹ : Space,
- R : Space + ClosedAdd,
- S : Preadjointable<A, Xʹ, PreadjointCodomain = R>,
- T : Preadjointable<A, Yʹ, PreadjointCodomain = R>,
- Self : Linear<A>,
- for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear<
- Pair<Xʹ,Yʹ>, Codomain = R,
- >,
+ 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,
+ Xʹ: Space,
+ Yʹ: Space,
+ R: ClosedSpace + ClosedAdd,
+ S: Preadjointable<A, Xʹ, PreadjointCodomain = R>,
+ T: Preadjointable<A, Yʹ, PreadjointCodomain = R>,
+ Self: Linear<A>,
+ for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
type PreadjointCodomain = R;
- type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+ 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())
@@ -517,100 +869,126 @@
}
/// Diagonal operator
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct DiagOp<S, T>(pub S, pub T);
impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
where
- A : Space,
- B : Space,
- S : Mapping<A>,
- T : Mapping<B>,
+ A: Space,
+ B: Space,
+ S: Mapping<A>,
+ T: Mapping<B>,
{
type Codomain = Pair<S::Codomain, T::Codomain>;
- fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
- let Pair(a, b) = x.decompose();
- Pair(self.0.apply(a), self.1.apply(b))
+ fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
+ x.eval_decompose(|Pair(a, b)| Pair(self.0.apply(a), self.1.apply(b)))
}
}
impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
where
- A : Space,
- B : Space,
- S : Linear<A>,
- T : Linear<B>,
-{ }
+ 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>
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>>,
+ 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>>,
{
- fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
- let Pair(u, v) = x.decompose();
- self.0.gemv(&mut y.0, α, u, β);
- self.1.gemv(&mut y.1, α, v, β);
+ fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.gemv(&mut y.0, α, u, β);
+ self.1.gemv(&mut y.1, α, v, β);
+ })
}
- fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
- let Pair(u, v) = x.decompose();
- self.0.apply_mut(&mut y.0, u);
- self.1.apply_mut(&mut y.1, v);
+ fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.apply_mut(&mut y.0, u);
+ self.1.apply_mut(&mut y.1, v);
+ })
}
/// Computes `y += Ax`, where `A` is `Self`
- fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
- let Pair(u, v) = x.decompose();
- self.0.apply_add(&mut y.0, u);
- self.1.apply_add(&mut y.1, v);
+ fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
+ x.eval_decompose(|Pair(u, v)| {
+ self.0.apply_add(&mut y.0, u);
+ self.1.apply_add(&mut y.1, v);
+ })
}
}
-impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
+impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
where
- A : Space,
- B : Space,
+ A: Space,
+ B: Space,
Xʹ: Space,
- Yʹ : Space,
- R : Space,
- 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,
- >,
+ Yʹ: Space,
+ R: ClosedSpace,
+ 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;
+ 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>
+impl<A, B, Xʹ, Yʹ, R, S, T> SimplyAdjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
where
- A : Space,
- B : Space,
+ A: Space,
+ B: Space,
Xʹ: Space,
- Yʹ : Space,
- R : Space,
- S : Preadjointable<A, Xʹ>,
- T : Preadjointable<B, Yʹ>,
- Self : Linear<Pair<A, B>>,
- for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear<
- Pair<Xʹ,Yʹ>, Codomain=R,
- >,
+ 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,
+ B: Space,
+ Xʹ: Space,
+ Yʹ: Space,
+ R: ClosedSpace,
+ S: Preadjointable<A, Xʹ>,
+ T: Preadjointable<B, Yʹ>,
+ Self: Linear<Pair<A, B>>,
+ for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
type PreadjointCodomain = R;
- type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
+ 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())
@@ -620,65 +998,90 @@
/// Block operator
pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
-
macro_rules! pairnorm {
($expj:ty) => {
impl<F, A, B, S, T, ExpA, ExpB, ExpR>
- BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F>
- for RowOp<S, T>
+ BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T>
where
- F : Float,
- A : Space + Norm<F, ExpA>,
- B : Space + 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,
- ExpA : NormExponent,
- ExpB : NormExponent,
- ExpR : NormExponent,
+ F: Float,
+ A: Space,
+ B: Space,
+ S: BoundedLinear<A, ExpA, ExpR, F>,
+ T: BoundedLinear<B, ExpB, ExpR, F>,
+ S::Codomain: Add<T::Codomain>,
+ <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
+ ExpA: NormExponent,
+ ExpB: NormExponent,
+ ExpR: NormExponent,
{
fn opnorm_bound(
&self,
- PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
- expr : ExpR
- ) -> F {
+ PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>,
+ expr: ExpR,
+ ) -> DynResult<F> {
// 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);
- let nb = self.1.opnorm_bound(expb, expr);
- na.max(nb)
+ let na = self.0.opnorm_bound(expa, expr)?;
+ let nb = self.1.opnorm_bound(expb, expr)?;
+ Ok(na.max(nb))
}
}
-
- impl<F, A, S, T, ExpA, ExpS, ExpT>
- BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
- for ColOp<S, T>
+
+ impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
+ for ColOp<S, T>
where
- F : Float,
- A : Space + Norm<F, ExpA>,
- S : BoundedLinear<A, ExpA, ExpS, F>,
- T : BoundedLinear<A, ExpA, ExpT, F>,
- ExpA : NormExponent,
- ExpS : NormExponent,
- ExpT : NormExponent,
+ F: Float,
+ A: Space,
+ S: BoundedLinear<A, ExpA, ExpS, F>,
+ T: BoundedLinear<A, ExpA, ExpT, F>,
+ ExpA: NormExponent,
+ ExpS: NormExponent,
+ ExpT: NormExponent,
{
fn opnorm_bound(
&self,
- expa : ExpA,
- PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj>
- ) -> F {
+ expa: ExpA,
+ PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>,
+ ) -> DynResult<F> {
// 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);
- let nt = self.1.opnorm_bound(expa, expt);
- ns.max(nt)
+ let ns = self.0.opnorm_bound(expa, exps)?;
+ let nt = self.1.opnorm_bound(expa, expt)?;
+ Ok(ns.max(nt))
}
}
- }
+ };
}
pairnorm!(L1);
pairnorm!(L2);
pairnorm!(Linfinity);
+/// The simplest linear mapping, scaling by a scalar.
+///
+/// TODO: redefined/replace `Weighted` by composition with [`Scaled`].
+pub struct Scaled<F: Float>(pub F);
+
+impl<Domain, F> Mapping<Domain> for Scaled<F>
+where
+ F: Float,
+ Domain: Space,
+ Domain::Principal: Mul<F>,
+ <Domain::Principal as Mul<F>>::Output: ClosedSpace,
+{
+ type Codomain = <Domain::Principal as Mul<F>>::Output;
+
+ /// Compute the value of `self` at `x`.
+ fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
+ x.own() * self.0
+ }
+}
+
+impl<Domain, F> Linear<Domain> for Scaled<F>
+where
+ F: Float,
+ Domain: Space,
+ Domain::Principal: Mul<F>,
+ <Domain::Principal as Mul<F>>::Output: ClosedSpace,
+{
+}