--- a/src/linops.rs Sat Dec 14 09:31:27 2024 -0500
+++ b/src/linops.rs Tue Dec 31 08:30:02 2024 -0500
@@ -4,59 +4,69 @@
use numeric_literals::replace_float_literals;
use std::marker::PhantomData;
+use serde::Serialize;
use crate::types::*;
-use serde::Serialize;
-pub use crate::mapping::Apply;
+pub use crate::mapping::{Mapping, Space};
use crate::direct_product::Pair;
+use crate::instance::Instance;
+use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity};
/// Trait for linear operators on `X`.
-pub trait Linear<X> : Apply<X, Output=Self::Codomain>
- + for<'a> Apply<&'a X, Output=Self::Codomain> {
- type Codomain;
-}
+pub trait Linear<X : Space> : Mapping<X>
+{ }
/// Efficient in-place summation.
#[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F : Num, X = Self> {
+pub trait AXPY<F, X = Self> : Space
+where
+ F : Num,
+ X : Space,
+{
/// Computes `y = βy + αx`, where `y` is `Self`.
- fn axpy(&mut self, α : F, x : &X, β : F);
+ fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
/// Copies `x` to `self`.
- fn copy_from(&mut self, x : &X) {
+ 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(&mut self, α : F, x : &X) {
+ fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) {
self.axpy(α, x, 0.0)
}
}
/// Efficient in-place application for [`Linear`] operators.
#[replace_float_literals(F::cast_from(literal))]
-pub trait GEMV<F : Num, X, Y = <Self as Linear<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(&self, y : &mut Y, α : F, x : &X, β : F);
+ fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
/// Computes `y = Ax`, where `A` is `Self`
- fn apply_mut(&self, y : &mut Y, x : &X){
+ fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
self.gemv(y, 1.0, x, 0.0)
}
/// Computes `y += Ax`, where `A` is `Self`
- fn apply_add(&self, y : &mut Y, x : &X){
+ 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> : Linear<X> {
- type FloatType : Float;
+pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
+where
+ 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.
- fn opnorm_bound(&self) -> Self::FloatType;
+ /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm
+ /// in `X`, and `codexp` in the codomain.
+ fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F;
}
// Linear operator application into mutable target. The [`AsRef`] bound
@@ -70,16 +80,16 @@
}*/
/// Trait for forming the adjoint operator of `Self`.
-pub trait Adjointable<X,Yʹ> : Linear<X> {
- type AdjointCodomain;
+pub trait Adjointable<X, Yʹ> : Linear<X>
+where
+ X : Space,
+ Yʹ : Space,
+{
+ type AdjointCodomain : Space;
type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
/// Form the adjoint operator of `self`.
fn adjoint(&self) -> Self::Adjoint<'_>;
-
- /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain {
- self.adjoint().apply(y)
- }*/
}
/// Trait for forming a preadjoint of an operator.
@@ -89,189 +99,206 @@
/// operator. The space `Ypre` is the predual of its codomain, and should be the
/// domain of the adjointed operator. `Self::Preadjoint` should be
/// [`Adjointable`]`<'a,Ypre,X>`.
-pub trait Preadjointable<X,Ypre> : Linear<X> {
- type PreadjointCodomain;
- type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;
- /// Form the preadjoint operator of `self`.
+pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> {
+ type PreadjointCodomain : Space;
+ type Preadjoint<'a> : Adjointable<
+ Ypre, X, AdjointCodomain=Self::Codomain, 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> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
-impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
+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>);
impl<X> IdOp<X> {
- fn new() -> IdOp<X> { IdOp(PhantomData) }
+ pub fn new() -> IdOp<X> { IdOp(PhantomData) }
}
-impl<X> Apply<X> for IdOp<X> {
- type Output = X;
+impl<X : Clone + Space> Mapping<X> for IdOp<X> {
+ type Codomain = X;
- fn apply(&self, x : X) -> X {
- x
+ fn apply<I : Instance<X>>(&self, x : I) -> X {
+ x.own()
}
}
-impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
- type Output = X;
-
- fn apply(&self, x : &'a X) -> X {
- x.clone()
- }
-}
-
-impl<X> Linear<X> for IdOp<X> where X : Clone {
- type Codomain = X;
-}
-
+impl<X : Clone + 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> where Y : AXPY<F, X>, X : Clone {
+impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
+where
+ Y : AXPY<F, X>,
+ X : Clone + Space
+{
// Computes `y = αAx + βy`, where `A` is `Self`.
- fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) {
+ fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
y.axpy(α, x, β)
}
- fn apply_mut(&self, y : &mut Y, x : &X){
+ fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
y.copy_from(x);
}
}
-impl<X> BoundedLinear<X> for IdOp<X> where X : Clone {
- type FloatType = float;
- fn opnorm_bound(&self) -> float { 1.0 }
+impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
+where
+ X : Space + Clone,
+ F : Num,
+ E : NormExponent
+{
+ fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
}
-impl<X> Adjointable<X,X> for IdOp<X> where X : Clone {
+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() }
}
+impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> {
+ type PreadjointCodomain=X;
+ type Preadjoint<'a> = IdOp<X> where X : 'a;
+
+ fn preadjoint(&self) -> Self::Preadjoint<'_> { 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>
+impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
where
- S : Apply<A>,
- T : Apply<B>,
- S::Output : Add<T::Output>
+ A : Space,
+ B : Space,
+ S : Mapping<A>,
+ T : Mapping<B>,
+ S::Codomain : Add<T::Codomain>,
+ <S::Codomain as Add<T::Codomain>>::Output : Space,
+
{
- type Output = <S::Output as Add<T::Output>>::Output;
+ type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
- fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::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)
}
}
-impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for RowOp<S, T>
+impl<A, B, S, T> Linear<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;
+ A : Space,
+ B : Space,
+ S : Linear<A>,
+ T : Linear<B>,
+ S::Codomain : Add<T::Codomain>,
+ <S::Codomain as Add<T::Codomain>>::Output : Space,
+{ }
- 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>
+impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> 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
+ 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(&self, y : &mut Y, α : F, x : &Pair<U, V>, β : F) {
- self.0.gemv(y, α, &x.0, β);
- self.1.gemv(y, α, &x.1, β);
+ 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 apply_mut(&self, y : &mut Y, x : &Pair<U, V>){
- self.0.apply_mut(y, &x.0);
- self.1.apply_mut(y, &x.1);
+ 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_mut(y, v);
}
/// 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);
+ 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);
}
}
-
/// “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>
+impl<A, S, T> Mapping<A> for ColOp<S, T>
where
- S : for<'a> Apply<&'a A, Output=O>,
- T : Apply<A>,
- A : std::borrow::Borrow<A>,
+ A : Space,
+ S : Mapping<A>,
+ T : Mapping<A>,
{
- type Output = Pair<O, T::Output>;
+ type Codomain = Pair<S::Codomain, T::Codomain>;
- fn apply(&self, a : A) -> Self::Output {
- Pair(self.0.apply(a.borrow()), self.1.apply(a))
+ fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
+ Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
}
}
-impl<A, S, T, D> Linear<A> for ColOp<S, T>
+impl<A, S, T> Linear<A> for ColOp<S, T>
where
- ColOp<S, T> : Apply<A, Output=D> + for<'a> Apply<&'a A, Output=D>,
-{
- type Codomain = D;
-}
+ 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>>
{
- fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &X, β : F) {
- self.0.gemv(&mut y.0, α, x, β);
+ fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
+ self.0.gemv(&mut y.0, α, x.ref_instance(), β);
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);
+ fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+ self.0.apply_mut(&mut y.0, x.ref_instance());
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);
+ fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
+ self.0.apply_add(&mut y.0, x.ref_instance());
self.1.apply_add(&mut y.1, x);
}
}
-impl<A, B, Yʹ, R, 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>>,
- for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Yʹ, Codomain=R>,
+ // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
+ // Yʹ,
+ // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
+ // >,
{
- type AdjointCodomain = R;
+ type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
fn adjoint(&self) -> Self::Adjoint<'_> {
@@ -279,14 +306,21 @@
}
}
-impl<A, B, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
+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>> : Adjointable<Yʹ, Pair<A,B>, Codomain=R>,
+ for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+ Yʹ, Pair<A,B>,
+ Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>,
+ AdjointCodomain = Self::Codomain,
+ >,
{
- type PreadjointCodomain = R;
+ type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
fn preadjoint(&self) -> Self::Preadjoint<'_> {
@@ -297,10 +331,17 @@
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ʹ>,
+ A : Space,
+ Xʹ : Space,
+ Yʹ : Space,
+ R : Space + 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>,
+ // 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;
@@ -312,10 +353,18 @@
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ʹ>,
+ 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>> : Adjointable<Pair<Xʹ,Yʹ>, A, Codomain=R>,
+ for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+ Pair<Xʹ,Yʹ>, A,
+ Codomain = R,
+ AdjointCodomain = Self::Codomain,
+ >,
{
type PreadjointCodomain = R;
type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
@@ -328,67 +377,73 @@
/// Diagonal operator
pub struct DiagOp<S, T>(pub S, pub T);
-impl<A, B, S, T> Apply<Pair<A, B>> for DiagOp<S, T>
+impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
where
- S : Apply<A>,
- T : Apply<B>,
+ A : Space,
+ B : Space,
+ S : Mapping<A>,
+ T : Mapping<B>,
{
- type Output = Pair<S::Output, T::Output>;
+ type Codomain = Pair<S::Codomain, T::Codomain>;
- 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 {
+ 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))
}
}
-impl<A, B, S, T, D> Linear<Pair<A, B>> for DiagOp<S, T>
+impl<A, B, S, T> 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;
-}
+ 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>>
+ 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 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 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);
+ 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);
}
/// 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);
+ 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);
}
}
-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,
+ 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>,
+ 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;
@@ -398,12 +453,21 @@
}
}
-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> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
where
+ 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>> : Adjointable<Pair<Xʹ,Yʹ>, Pair<A, B>, Codomain=R>,
+ for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
+ Pair<Xʹ,Yʹ>, Pair<A, B>,
+ Codomain=R,
+ AdjointCodomain = Self::Codomain,
+ >,
{
type PreadjointCodomain = R;
type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
@@ -416,3 +480,65 @@
/// 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>
+ where
+ 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 : Space,
+ ExpA : NormExponent,
+ ExpB : NormExponent,
+ ExpR : NormExponent,
+ {
+ fn opnorm_bound(
+ &self,
+ PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
+ expr : ExpR
+ ) -> 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)
+ }
+ }
+
+ 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,
+ 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 {
+ // 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)
+ }
+ }
+ }
+}
+
+pairnorm!(L1);
+pairnorm!(L2);
+pairnorm!(Linfinity);
+