Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 99
9e5b9fc81c52
parent 74
2c76df38d02b
child 101
997961aa6eee
equal deleted inserted replaced
98:2e4517b55442 99:9e5b9fc81c52
1 /*! 1 /*!
2 Abstract linear operators. 2 Abstract linear operators.
3 */ 3 */
4 4
5 use numeric_literals::replace_float_literals;
6 use std::marker::PhantomData;
7 use serde::Serialize;
8 use crate::types::*;
9 pub use crate::mapping::{Mapping, Space, Composition};
10 use crate::direct_product::Pair; 5 use crate::direct_product::Pair;
11 use crate::instance::Instance; 6 use crate::instance::Instance;
12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm}; 7 pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space};
8 use crate::norms::{Linfinity, Norm, NormExponent, PairNorm, L1, L2};
9 use crate::types::*;
10 use numeric_literals::replace_float_literals;
11 use serde::Serialize;
12 use std::marker::PhantomData;
13 13
14 /// Trait for linear operators on `X`. 14 /// Trait for linear operators on `X`.
15 pub trait Linear<X : Space> : Mapping<X> 15 pub trait Linear<X: Space>: Mapping<X> {}
16 { } 16
17 // impl<X: Space, A: Linear<X>> DifferentiableImpl<X> for A {
18 // type Derivative = <Self as Mapping<X>>::Codomain;
19
20 // /// Compute the differential of `self` at `x`, consuming the input.
21 // fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
22 // self.apply(x)
23 // }
24 // }
17 25
18 /// Efficient in-place summation. 26 /// Efficient in-place summation.
19 #[replace_float_literals(F::cast_from(literal))] 27 #[replace_float_literals(F::cast_from(literal))]
20 pub trait AXPY<F, X = Self> : Space + std::ops::MulAssign<F> 28 pub trait AXPY<F, X = Self>: Space + std::ops::MulAssign<F>
21 where 29 where
22 F : Num, 30 F: Num,
23 X : Space, 31 X: Space,
24 { 32 {
25 type Owned : AXPY<F, X>; 33 type Owned: AXPY<F, X>;
26 34
27 /// Computes `y = βy + αx`, where `y` is `Self`. 35 /// Computes `y = βy + αx`, where `y` is `Self`.
28 fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F); 36 fn axpy<I: Instance<X>>(&mut self, α: F, x: I, β: F);
29 37
30 /// Copies `x` to `self`. 38 /// Copies `x` to `self`.
31 fn copy_from<I : Instance<X>>(&mut self, x : I) { 39 fn copy_from<I: Instance<X>>(&mut self, x: I) {
32 self.axpy(1.0, x, 0.0) 40 self.axpy(1.0, x, 0.0)
33 } 41 }
34 42
35 /// Computes `y = αx`, where `y` is `Self`. 43 /// Computes `y = αx`, where `y` is `Self`.
36 fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) { 44 fn scale_from<I: Instance<X>>(&mut self, α: F, x: I) {
37 self.axpy(α, x, 0.0) 45 self.axpy(α, x, 0.0)
38 } 46 }
39 47
40 /// Return a similar zero as `self`. 48 /// Return a similar zero as `self`.
41 fn similar_origin(&self) -> Self::Owned; 49 fn similar_origin(&self) -> Self::Owned;
44 fn set_zero(&mut self); 52 fn set_zero(&mut self);
45 } 53 }
46 54
47 /// Efficient in-place application for [`Linear`] operators. 55 /// Efficient in-place application for [`Linear`] operators.
48 #[replace_float_literals(F::cast_from(literal))] 56 #[replace_float_literals(F::cast_from(literal))]
49 pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> { 57 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> {
50 /// Computes `y = αAx + βy`, where `A` is `Self`. 58 /// Computes `y = αAx + βy`, where `A` is `Self`.
51 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F); 59 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F);
52 60
53 #[inline] 61 #[inline]
54 /// Computes `y = Ax`, where `A` is `Self` 62 /// Computes `y = Ax`, where `A` is `Self`
55 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 63 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
56 self.gemv(y, 1.0, x, 0.0) 64 self.gemv(y, 1.0, x, 0.0)
57 } 65 }
58 66
59 #[inline] 67 #[inline]
60 /// Computes `y += Ax`, where `A` is `Self` 68 /// Computes `y += Ax`, where `A` is `Self`
61 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ 69 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
62 self.gemv(y, 1.0, x, 1.0) 70 self.gemv(y, 1.0, x, 1.0)
63 } 71 }
64 } 72 }
65 73
66
67 /// Bounded linear operators 74 /// Bounded linear operators
68 pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X> 75 pub trait BoundedLinear<X, XExp, CodExp, F = f64>: Linear<X>
69 where 76 where
70 F : Num, 77 F: Num,
71 X : Space + Norm<F, XExp>, 78 X: Space + Norm<F, XExp>,
72 XExp : NormExponent, 79 XExp: NormExponent,
73 CodExp : NormExponent 80 CodExp: NormExponent,
74 { 81 {
75 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. 82 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
76 /// This is not expected to be the norm, just any bound on it that can be 83 /// This is not expected to be the norm, just any bound on it that can be
77 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm 84 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm
78 /// in `X`, and `codexp` in the codomain. 85 /// in `X`, and `codexp` in the codomain.
79 fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F; 86 fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> F;
80 } 87 }
81 88
82 // Linear operator application into mutable target. The [`AsRef`] bound 89 // Linear operator application into mutable target. The [`AsRef`] bound
83 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; 90 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
84 // the former is assumed to be e.g. a view into the latter. 91 // the former is assumed to be e.g. a view into the latter.
88 self.apply(x) 95 self.apply(x)
89 } 96 }
90 }*/ 97 }*/
91 98
92 /// Trait for forming the adjoint operator of `Self`. 99 /// Trait for forming the adjoint operator of `Self`.
93 pub trait Adjointable<X, Yʹ> : Linear<X> 100 pub trait Adjointable<X, Yʹ>: Linear<X>
94 where 101 where
95 X : Space, 102 X: Space,
96 Yʹ : Space, 103 Yʹ: Space,
97 { 104 {
98 type AdjointCodomain : Space; 105 type AdjointCodomain: Space;
99 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; 106 type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
107 where
108 Self: 'a;
100 109
101 /// Form the adjoint operator of `self`. 110 /// Form the adjoint operator of `self`.
102 fn adjoint(&self) -> Self::Adjoint<'_>; 111 fn adjoint(&self) -> Self::Adjoint<'_>;
103 } 112 }
104 113
110 /// domain of the adjointed operator. `Self::Preadjoint` should be 119 /// domain of the adjointed operator. `Self::Preadjoint` should be
111 /// [`Adjointable`]`<'a,Ypre,X>`. 120 /// [`Adjointable`]`<'a,Ypre,X>`.
112 /// We do not make additional restrictions on `Self::Preadjoint` (in particular, it 121 /// We do not make additional restrictions on `Self::Preadjoint` (in particular, it
113 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint 122 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
114 /// have the full space as the codomain, etc. 123 /// have the full space as the codomain, etc.
115 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> { 124 pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
116 type PreadjointCodomain : Space; 125 Linear<X>
117 type Preadjoint<'a> : Linear< 126 {
118 Ypre, Codomain=Self::PreadjointCodomain 127 type PreadjointCodomain: Space;
119 > where Self : 'a; 128 type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain>
129 where
130 Self: 'a;
120 131
121 /// Form the adjoint operator of `self`. 132 /// Form the adjoint operator of `self`.
122 fn preadjoint(&self) -> Self::Preadjoint<'_>; 133 fn preadjoint(&self) -> Self::Preadjoint<'_>;
123 } 134 }
124 135
125 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$. 136 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
126 pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {} 137 pub trait SimplyAdjointable<X: Space>: Adjointable<X, <Self as Mapping<X>>::Codomain> {}
127 impl<'a,X : Space, T> SimplyAdjointable<X> for T 138 impl<'a, X: Space, T> SimplyAdjointable<X> for T where
128 where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {} 139 T: Adjointable<X, <Self as Mapping<X>>::Codomain>
140 {
141 }
129 142
130 /// The identity operator 143 /// The identity operator
131 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 144 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
132 pub struct IdOp<X> (PhantomData<X>); 145 pub struct IdOp<X>(PhantomData<X>);
133 146
134 impl<X> IdOp<X> { 147 impl<X> IdOp<X> {
135 pub fn new() -> IdOp<X> { IdOp(PhantomData) } 148 pub fn new() -> IdOp<X> {
136 } 149 IdOp(PhantomData)
137 150 }
138 impl<X : Clone + Space> Mapping<X> for IdOp<X> { 151 }
152
153 impl<X: Clone + Space> Mapping<X> for IdOp<X> {
139 type Codomain = X; 154 type Codomain = X;
140 155
141 fn apply<I : Instance<X>>(&self, x : I) -> X { 156 fn apply<I: Instance<X>>(&self, x: I) -> X {
142 x.own() 157 x.own()
143 } 158 }
144 } 159 }
145 160
146 impl<X : Clone + Space> Linear<X> for IdOp<X> 161 impl<X: Clone + Space> Linear<X> for IdOp<X> {}
147 { }
148 162
149 #[replace_float_literals(F::cast_from(literal))] 163 #[replace_float_literals(F::cast_from(literal))]
150 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> 164 impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
151 where 165 where
152 Y : AXPY<F, X>, 166 Y: AXPY<F, X>,
153 X : Clone + Space 167 X: Clone + Space,
154 { 168 {
155 // Computes `y = αAx + βy`, where `A` is `Self`. 169 // Computes `y = αAx + βy`, where `A` is `Self`.
156 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { 170 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
157 y.axpy(α, x, β) 171 y.axpy(α, x, β)
158 } 172 }
159 173
160 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 174 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
161 y.copy_from(x); 175 y.copy_from(x);
162 } 176 }
163 } 177 }
164 178
165 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> 179 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
166 where 180 where
167 X : Space + Clone + Norm<F, E>, 181 X: Space + Clone + Norm<F, E>,
168 F : Num, 182 F: Num,
169 E : NormExponent 183 E: NormExponent,
170 { 184 {
171 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } 185 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> F {
172 } 186 F::ONE
173 187 }
174 impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> { 188 }
175 type AdjointCodomain=X; 189
176 type Adjoint<'a> = IdOp<X> where X : 'a; 190 impl<X: Clone + Space> Adjointable<X, X> for IdOp<X> {
177 191 type AdjointCodomain = X;
178 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 192 type Adjoint<'a>
179 } 193 = IdOp<X>
180 194 where
181 impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> { 195 X: 'a;
182 type PreadjointCodomain=X; 196
183 type Preadjoint<'a> = IdOp<X> where X : 'a; 197 fn adjoint(&self) -> Self::Adjoint<'_> {
184 198 IdOp::new()
185 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() } 199 }
186 } 200 }
187 201
202 impl<X: Clone + Space> Preadjointable<X, X> for IdOp<X> {
203 type PreadjointCodomain = X;
204 type Preadjoint<'a>
205 = IdOp<X>
206 where
207 X: 'a;
208
209 fn preadjoint(&self) -> Self::Preadjoint<'_> {
210 IdOp::new()
211 }
212 }
188 213
189 /// The zero operator 214 /// The zero operator
190 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 215 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
191 pub struct ZeroOp<'a, X, XD, Y, F> { 216 pub struct ZeroOp<'a, X, XD, Y, F> {
192 zero : &'a Y, // TODO: don't pass this in `new`; maybe not even store. 217 zero: &'a Y, // TODO: don't pass this in `new`; maybe not even store.
193 dual_or_predual_zero : XD, 218 dual_or_predual_zero: XD,
194 _phantoms : PhantomData<(X, Y, F)>, 219 _phantoms: PhantomData<(X, Y, F)>,
195 } 220 }
196 221
197 // TODO: Need to make Zero in Instance. 222 // TODO: Need to make Zero in Instance.
198 223
199 impl<'a, F : Num, X : Space, XD, Y : Space + Clone> ZeroOp<'a, X, XD, Y, F> { 224 impl<'a, F: Num, X: Space, XD, Y: Space + Clone> ZeroOp<'a, X, XD, Y, F> {
200 pub fn new(zero : &'a Y, dual_or_predual_zero : XD) -> Self { 225 pub fn new(zero: &'a Y, dual_or_predual_zero: XD) -> Self {
201 ZeroOp{ zero, dual_or_predual_zero, _phantoms : PhantomData } 226 ZeroOp {
202 } 227 zero,
203 } 228 dual_or_predual_zero,
204 229 _phantoms: PhantomData,
205 impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> { 230 }
231 }
232 }
233
234 impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> {
206 type Codomain = Y; 235 type Codomain = Y;
207 236
208 fn apply<I : Instance<X>>(&self, _x : I) -> Y { 237 fn apply<I: Instance<X>>(&self, _x: I) -> Y {
209 self.zero.clone() 238 self.zero.clone()
210 } 239 }
211 } 240 }
212 241
213 impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> 242 impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> {}
214 { }
215 243
216 #[replace_float_literals(F::cast_from(literal))] 244 #[replace_float_literals(F::cast_from(literal))]
217 impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F> 245 impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F>
218 where 246 where
219 F : Num, 247 F: Num,
220 Y : AXPY<F, Y> + Clone, 248 Y: AXPY<F, Y> + Clone,
221 X : Space 249 X: Space,
222 { 250 {
223 // Computes `y = αAx + βy`, where `A` is `Self`. 251 // Computes `y = αAx + βy`, where `A` is `Self`.
224 fn gemv<I : Instance<X>>(&self, y : &mut Y, _α : F, _x : I, β : F) { 252 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
225 *y *= β; 253 *y *= β;
226 } 254 }
227 255
228 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, _x : I){ 256 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
229 y.set_zero(); 257 y.set_zero();
230 } 258 }
231 } 259 }
232 260
233 impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F> 261 impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F>
234 where 262 where
235 X : Space + Norm<F, E1>, 263 X: Space + Norm<F, E1>,
236 Y : AXPY<F> + Clone + Norm<F, E2>, 264 Y: AXPY<F> + Clone + Norm<F, E2>,
237 F : Num, 265 F: Num,
238 E1 : NormExponent, 266 E1: NormExponent,
239 E2 : NormExponent, 267 E2: NormExponent,
240 { 268 {
241 fn opnorm_bound(&self, _xexp : E1, _codexp : E2) -> F { F::ZERO } 269 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> F {
242 } 270 F::ZERO
243 271 }
244 impl<'a, F : Num, X, XD, Y, Yprime : Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F> 272 }
245 where 273
246 X : Space, 274 impl<'a, F: Num, X, XD, Y, Yprime: Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F>
247 Y : AXPY<F> + Clone + 'static, 275 where
248 XD : AXPY<F> + Clone + 'static, 276 X: Space,
277 Y: AXPY<F> + Clone + 'static,
278 XD: AXPY<F> + Clone + 'static,
249 { 279 {
250 type AdjointCodomain = XD; 280 type AdjointCodomain = XD;
251 type Adjoint<'b> = ZeroOp<'b, Yprime, (), XD, F> where Self : 'b; 281 type Adjoint<'b>
252 // () means not (pre)adjointable. 282 = ZeroOp<'b, Yprime, (), XD, F>
283 where
284 Self: 'b;
285 // () means not (pre)adjointable.
253 286
254 fn adjoint(&self) -> Self::Adjoint<'_> { 287 fn adjoint(&self) -> Self::Adjoint<'_> {
255 ZeroOp::new(&self.dual_or_predual_zero, ()) 288 ZeroOp::new(&self.dual_or_predual_zero, ())
256 } 289 }
257 } 290 }
258 291
259 impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F> 292 impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F>
260 where 293 where
261 F : Num, 294 F: Num,
262 X : Space, 295 X: Space,
263 Y : AXPY<F> + Clone, 296 Y: AXPY<F> + Clone,
264 Ypre : Space, 297 Ypre: Space,
265 XD : AXPY<F> + Clone + 'static, 298 XD: AXPY<F> + Clone + 'static,
266 { 299 {
267 type PreadjointCodomain = XD; 300 type PreadjointCodomain = XD;
268 type Preadjoint<'b> = ZeroOp<'b, Ypre, (), XD, F> where Self : 'b; 301 type Preadjoint<'b>
269 // () means not (pre)adjointable. 302 = ZeroOp<'b, Ypre, (), XD, F>
303 where
304 Self: 'b;
305 // () means not (pre)adjointable.
270 306
271 fn preadjoint(&self) -> Self::Preadjoint<'_> { 307 fn preadjoint(&self) -> Self::Preadjoint<'_> {
272 ZeroOp::new(&self.dual_or_predual_zero, ()) 308 ZeroOp::new(&self.dual_or_predual_zero, ())
273 } 309 }
274 } 310 }
275 311
276 impl<S, T, E, X> Linear<X> for Composition<S, T, E> 312 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
277 where 313 where
278 X : Space, 314 X: Space,
279 T : Linear<X>, 315 T: Linear<X>,
280 S : Linear<T::Codomain> 316 S: Linear<T::Codomain>,
281 { } 317 {
318 }
282 319
283 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E> 320 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
284 where 321 where
285 F : Num, 322 F: Num,
286 X : Space, 323 X: Space,
287 T : Linear<X>, 324 T: Linear<X>,
288 S : GEMV<F, T::Codomain, Y>, 325 S: GEMV<F, T::Codomain, Y>,
289 { 326 {
290 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { 327 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
291 self.outer.gemv(y, α, self.inner.apply(x), β) 328 self.outer.gemv(y, α, self.inner.apply(x), β)
292 } 329 }
293 330
294 /// Computes `y = Ax`, where `A` is `Self` 331 /// Computes `y = Ax`, where `A` is `Self`
295 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 332 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
296 self.outer.apply_mut(y, self.inner.apply(x)) 333 self.outer.apply_mut(y, self.inner.apply(x))
297 } 334 }
298 335
299 /// Computes `y += Ax`, where `A` is `Self` 336 /// Computes `y += Ax`, where `A` is `Self`
300 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ 337 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
301 self.outer.apply_add(y, self.inner.apply(x)) 338 self.outer.apply_add(y, self.inner.apply(x))
302 } 339 }
303 } 340 }
304 341
305 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp> 342 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
306 where 343 where
307 F : Num, 344 F: Num,
308 X : Space + Norm<F, Xexp>, 345 X: Space + Norm<F, Xexp>,
309 Z : Space + Norm<F, Zexp>, 346 Z: Space + Norm<F, Zexp>,
310 Xexp : NormExponent, 347 Xexp: NormExponent,
311 Yexp : NormExponent, 348 Yexp: NormExponent,
312 Zexp : NormExponent, 349 Zexp: NormExponent,
313 T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>, 350 T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>,
314 S : BoundedLinear<Z, Zexp, Yexp, F>, 351 S: BoundedLinear<Z, Zexp, Yexp, F>,
315 { 352 {
316 fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F { 353 fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> F {
317 let zexp = self.intermediate_norm_exponent; 354 let zexp = self.intermediate_norm_exponent;
318 self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp) 355 self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp)
319 } 356 }
320 } 357 }
321 358
324 361
325 use std::ops::Add; 362 use std::ops::Add;
326 363
327 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T> 364 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
328 where 365 where
329 A : Space, 366 A: Space,
330 B : Space, 367 B: Space,
331 S : Mapping<A>, 368 S: Mapping<A>,
332 T : Mapping<B>, 369 T: Mapping<B>,
333 S::Codomain : Add<T::Codomain>, 370 S::Codomain: Add<T::Codomain>,
334 <S::Codomain as Add<T::Codomain>>::Output : Space, 371 <S::Codomain as Add<T::Codomain>>::Output: Space,
335
336 { 372 {
337 type Codomain = <S::Codomain as Add<T::Codomain>>::Output; 373 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
338 374
339 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { 375 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
340 let Pair(a, b) = x.decompose(); 376 let Pair(a, b) = x.decompose();
341 self.0.apply(a) + self.1.apply(b) 377 self.0.apply(a) + self.1.apply(b)
342 } 378 }
343 } 379 }
344 380
345 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> 381 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
346 where 382 where
347 A : Space, 383 A: Space,
348 B : Space, 384 B: Space,
349 S : Linear<A>, 385 S: Linear<A>,
350 T : Linear<B>, 386 T: Linear<B>,
351 S::Codomain : Add<T::Codomain>, 387 S::Codomain: Add<T::Codomain>,
352 <S::Codomain as Add<T::Codomain>>::Output : Space, 388 <S::Codomain as Add<T::Codomain>>::Output: Space,
353 { } 389 {
354 390 }
355 391
356 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> 392 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
357 where 393 where
358 U : Space, 394 U: Space,
359 V : Space, 395 V: Space,
360 S : GEMV<F, U, Y>, 396 S: GEMV<F, U, Y>,
361 T : GEMV<F, V, Y>, 397 T: GEMV<F, V, Y>,
362 F : Num, 398 F: Num,
363 Self : Linear<Pair<U, V>, Codomain=Y> 399 Self: Linear<Pair<U, V>, Codomain = Y>,
364 { 400 {
365 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) { 401 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Y, α: F, x: I, β: F) {
366 let Pair(u, v) = x.decompose(); 402 let Pair(u, v) = x.decompose();
367 self.0.gemv(y, α, u, β); 403 self.0.gemv(y, α, u, β);
368 self.1.gemv(y, α, v, F::ONE); 404 self.1.gemv(y, α, v, F::ONE);
369 } 405 }
370 406
371 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { 407 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
372 let Pair(u, v) = x.decompose(); 408 let Pair(u, v) = x.decompose();
373 self.0.apply_mut(y, u); 409 self.0.apply_mut(y, u);
374 self.1.apply_add(y, v); 410 self.1.apply_add(y, v);
375 } 411 }
376 412
377 /// Computes `y += Ax`, where `A` is `Self` 413 /// Computes `y += Ax`, where `A` is `Self`
378 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { 414 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
379 let Pair(u, v) = x.decompose(); 415 let Pair(u, v) = x.decompose();
380 self.0.apply_add(y, u); 416 self.0.apply_add(y, u);
381 self.1.apply_add(y, v); 417 self.1.apply_add(y, v);
382 } 418 }
383 } 419 }
385 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. 421 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
386 pub struct ColOp<S, T>(pub S, pub T); 422 pub struct ColOp<S, T>(pub S, pub T);
387 423
388 impl<A, S, T> Mapping<A> for ColOp<S, T> 424 impl<A, S, T> Mapping<A> for ColOp<S, T>
389 where 425 where
390 A : Space, 426 A: Space,
391 S : Mapping<A>, 427 S: Mapping<A>,
392 T : Mapping<A>, 428 T: Mapping<A>,
393 { 429 {
394 type Codomain = Pair<S::Codomain, T::Codomain>; 430 type Codomain = Pair<S::Codomain, T::Codomain>;
395 431
396 fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain { 432 fn apply<I: Instance<A>>(&self, a: I) -> Self::Codomain {
397 Pair(self.0.apply(a.ref_instance()), self.1.apply(a)) 433 Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
398 } 434 }
399 } 435 }
400 436
401 impl<A, S, T> Linear<A> for ColOp<S, T> 437 impl<A, S, T> Linear<A> for ColOp<S, T>
402 where 438 where
403 A : Space, 439 A: Space,
404 S : Mapping<A>, 440 S: Mapping<A>,
405 T : Mapping<A>, 441 T: Mapping<A>,
406 { } 442 {
443 }
407 444
408 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> 445 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
409 where 446 where
410 X : Space, 447 X: Space,
411 S : GEMV<F, X, A>, 448 S: GEMV<F, X, A>,
412 T : GEMV<F, X, B>, 449 T: GEMV<F, X, B>,
413 F : Num, 450 F: Num,
414 Self : Linear<X, Codomain=Pair<A, B>> 451 Self: Linear<X, Codomain = Pair<A, B>>,
415 { 452 {
416 fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 453 fn gemv<I: Instance<X>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
417 self.0.gemv(&mut y.0, α, x.ref_instance(), β); 454 self.0.gemv(&mut y.0, α, x.ref_instance(), β);
418 self.1.gemv(&mut y.1, α, x, β); 455 self.1.gemv(&mut y.1, α, x, β);
419 } 456 }
420 457
421 fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ 458 fn apply_mut<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
422 self.0.apply_mut(&mut y.0, x.ref_instance()); 459 self.0.apply_mut(&mut y.0, x.ref_instance());
423 self.1.apply_mut(&mut y.1, x); 460 self.1.apply_mut(&mut y.1, x);
424 } 461 }
425 462
426 /// Computes `y += Ax`, where `A` is `Self` 463 /// Computes `y += Ax`, where `A` is `Self`
427 fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ 464 fn apply_add<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
428 self.0.apply_add(&mut y.0, x.ref_instance()); 465 self.0.apply_add(&mut y.0, x.ref_instance());
429 self.1.apply_add(&mut y.1, x); 466 self.1.apply_add(&mut y.1, x);
430 } 467 }
431 } 468 }
432 469
433 470 impl<A, B, Yʹ, S, T> Adjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
434 impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 471 where
435 where 472 A: Space,
436 A : Space, 473 B: Space,
437 B : Space, 474 Yʹ: Space,
438 Yʹ : Space, 475 S: Adjointable<A, Yʹ>,
439 S : Adjointable<A, Yʹ>, 476 T: Adjointable<B, Yʹ>,
440 T : Adjointable<B, Yʹ>, 477 Self: Linear<Pair<A, B>>,
441 Self : Linear<Pair<A, B>>,
442 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 478 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
443 // Yʹ, 479 // Yʹ,
444 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain> 480 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
445 // >, 481 // >,
446 { 482 {
447 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>; 483 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
448 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 484 type Adjoint<'a>
485 = ColOp<S::Adjoint<'a>, T::Adjoint<'a>>
486 where
487 Self: 'a;
449 488
450 fn adjoint(&self) -> Self::Adjoint<'_> { 489 fn adjoint(&self) -> Self::Adjoint<'_> {
451 ColOp(self.0.adjoint(), self.1.adjoint()) 490 ColOp(self.0.adjoint(), self.1.adjoint())
452 } 491 }
453 } 492 }
454 493
455 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 494 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
456 where 495 where
457 A : Space, 496 A: Space,
458 B : Space, 497 B: Space,
459 Yʹ : Space, 498 Yʹ: Space,
460 S : Preadjointable<A, Yʹ>, 499 S: Preadjointable<A, Yʹ>,
461 T : Preadjointable<B, Yʹ>, 500 T: Preadjointable<B, Yʹ>,
462 Self : Linear<Pair<A, B>>, 501 Self: Linear<Pair<A, B>>,
463 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 502 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>:
464 Yʹ, Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>, 503 Linear<Yʹ, Codomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>>,
465 >,
466 { 504 {
467 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; 505 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
468 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 506 type Preadjoint<'a>
507 = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
508 where
509 Self: 'a;
469 510
470 fn preadjoint(&self) -> Self::Preadjoint<'_> { 511 fn preadjoint(&self) -> Self::Preadjoint<'_> {
471 ColOp(self.0.preadjoint(), self.1.preadjoint()) 512 ColOp(self.0.preadjoint(), self.1.preadjoint())
472 } 513 }
473 } 514 }
474 515
475 516 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
476 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 517 where
477 where 518 A: Space,
478 A : Space, 519 Xʹ: Space,
479 Xʹ : Space, 520 Yʹ: Space,
480 Yʹ : Space, 521 R: Space + ClosedAdd,
481 R : Space + ClosedAdd, 522 S: Adjointable<A, Xʹ, AdjointCodomain = R>,
482 S : Adjointable<A, Xʹ, AdjointCodomain = R>, 523 T: Adjointable<A, Yʹ, AdjointCodomain = R>,
483 T : Adjointable<A, Yʹ, AdjointCodomain = R>, 524 Self: Linear<A>,
484 Self : Linear<A>,
485 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 525 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
486 // Pair<Xʹ,Yʹ>, 526 // Pair<Xʹ,Yʹ>,
487 // Codomain=R, 527 // Codomain=R,
488 // >, 528 // >,
489 { 529 {
490 type AdjointCodomain = R; 530 type AdjointCodomain = R;
491 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 531 type Adjoint<'a>
532 = RowOp<S::Adjoint<'a>, T::Adjoint<'a>>
533 where
534 Self: 'a;
492 535
493 fn adjoint(&self) -> Self::Adjoint<'_> { 536 fn adjoint(&self) -> Self::Adjoint<'_> {
494 RowOp(self.0.adjoint(), self.1.adjoint()) 537 RowOp(self.0.adjoint(), self.1.adjoint())
495 } 538 }
496 } 539 }
497 540
498 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 541 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
499 where 542 where
500 A : Space, 543 A: Space,
501 Xʹ : Space, 544 Xʹ: Space,
502 Yʹ : Space, 545 Yʹ: Space,
503 R : Space + ClosedAdd, 546 R: Space + ClosedAdd,
504 S : Preadjointable<A, Xʹ, PreadjointCodomain = R>, 547 S: Preadjointable<A, Xʹ, PreadjointCodomain = R>,
505 T : Preadjointable<A, Yʹ, PreadjointCodomain = R>, 548 T: Preadjointable<A, Yʹ, PreadjointCodomain = R>,
506 Self : Linear<A>, 549 Self: Linear<A>,
507 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 550 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
508 Pair<Xʹ,Yʹ>, Codomain = R,
509 >,
510 { 551 {
511 type PreadjointCodomain = R; 552 type PreadjointCodomain = R;
512 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 553 type Preadjoint<'a>
554 = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
555 where
556 Self: 'a;
513 557
514 fn preadjoint(&self) -> Self::Preadjoint<'_> { 558 fn preadjoint(&self) -> Self::Preadjoint<'_> {
515 RowOp(self.0.preadjoint(), self.1.preadjoint()) 559 RowOp(self.0.preadjoint(), self.1.preadjoint())
516 } 560 }
517 } 561 }
519 /// Diagonal operator 563 /// Diagonal operator
520 pub struct DiagOp<S, T>(pub S, pub T); 564 pub struct DiagOp<S, T>(pub S, pub T);
521 565
522 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T> 566 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
523 where 567 where
524 A : Space, 568 A: Space,
525 B : Space, 569 B: Space,
526 S : Mapping<A>, 570 S: Mapping<A>,
527 T : Mapping<B>, 571 T: Mapping<B>,
528 { 572 {
529 type Codomain = Pair<S::Codomain, T::Codomain>; 573 type Codomain = Pair<S::Codomain, T::Codomain>;
530 574
531 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { 575 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
532 let Pair(a, b) = x.decompose(); 576 let Pair(a, b) = x.decompose();
533 Pair(self.0.apply(a), self.1.apply(b)) 577 Pair(self.0.apply(a), self.1.apply(b))
534 } 578 }
535 } 579 }
536 580
537 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> 581 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
538 where 582 where
539 A : Space, 583 A: Space,
540 B : Space, 584 B: Space,
541 S : Linear<A>, 585 S: Linear<A>,
542 T : Linear<B>, 586 T: Linear<B>,
543 { } 587 {
588 }
544 589
545 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> 590 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
546 where 591 where
547 A : Space, 592 A: Space,
548 B : Space, 593 B: Space,
549 U : Space, 594 U: Space,
550 V : Space, 595 V: Space,
551 S : GEMV<F, U, A>, 596 S: GEMV<F, U, A>,
552 T : GEMV<F, V, B>, 597 T: GEMV<F, V, B>,
553 F : Num, 598 F: Num,
554 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>, 599 Self: Linear<Pair<U, V>, Codomain = Pair<A, B>>,
555 { 600 {
556 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 601 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
557 let Pair(u, v) = x.decompose(); 602 let Pair(u, v) = x.decompose();
558 self.0.gemv(&mut y.0, α, u, β); 603 self.0.gemv(&mut y.0, α, u, β);
559 self.1.gemv(&mut y.1, α, v, β); 604 self.1.gemv(&mut y.1, α, v, β);
560 } 605 }
561 606
562 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ 607 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
563 let Pair(u, v) = x.decompose(); 608 let Pair(u, v) = x.decompose();
564 self.0.apply_mut(&mut y.0, u); 609 self.0.apply_mut(&mut y.0, u);
565 self.1.apply_mut(&mut y.1, v); 610 self.1.apply_mut(&mut y.1, v);
566 } 611 }
567 612
568 /// Computes `y += Ax`, where `A` is `Self` 613 /// Computes `y += Ax`, where `A` is `Self`
569 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ 614 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
570 let Pair(u, v) = x.decompose(); 615 let Pair(u, v) = x.decompose();
571 self.0.apply_add(&mut y.0, u); 616 self.0.apply_add(&mut y.0, u);
572 self.1.apply_add(&mut y.1, v); 617 self.1.apply_add(&mut y.1, v);
573 } 618 }
574 } 619 }
575 620
576 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 621 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
577 where 622 where
578 A : Space, 623 A: Space,
579 B : Space, 624 B: Space,
580 Xʹ: Space, 625 Xʹ: Space,
581 Yʹ : Space, 626 Yʹ: Space,
582 R : Space, 627 R: Space,
583 S : Adjointable<A, Xʹ>, 628 S: Adjointable<A, Xʹ>,
584 T : Adjointable<B, Yʹ>, 629 T: Adjointable<B, Yʹ>,
585 Self : Linear<Pair<A, B>>, 630 Self: Linear<Pair<A, B>>,
586 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 631 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
587 Pair<Xʹ,Yʹ>, Codomain=R,
588 >,
589 { 632 {
590 type AdjointCodomain = R; 633 type AdjointCodomain = R;
591 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 634 type Adjoint<'a>
635 = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>
636 where
637 Self: 'a;
592 638
593 fn adjoint(&self) -> Self::Adjoint<'_> { 639 fn adjoint(&self) -> Self::Adjoint<'_> {
594 DiagOp(self.0.adjoint(), self.1.adjoint()) 640 DiagOp(self.0.adjoint(), self.1.adjoint())
595 } 641 }
596 } 642 }
597 643
598 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 644 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
599 where 645 where
600 A : Space, 646 A: Space,
601 B : Space, 647 B: Space,
602 Xʹ: Space, 648 Xʹ: Space,
603 Yʹ : Space, 649 Yʹ: Space,
604 R : Space, 650 R: Space,
605 S : Preadjointable<A, Xʹ>, 651 S: Preadjointable<A, Xʹ>,
606 T : Preadjointable<B, Yʹ>, 652 T: Preadjointable<B, Yʹ>,
607 Self : Linear<Pair<A, B>>, 653 Self: Linear<Pair<A, B>>,
608 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 654 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
609 Pair<Xʹ,Yʹ>, Codomain=R,
610 >,
611 { 655 {
612 type PreadjointCodomain = R; 656 type PreadjointCodomain = R;
613 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 657 type Preadjoint<'a>
658 = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
659 where
660 Self: 'a;
614 661
615 fn preadjoint(&self) -> Self::Preadjoint<'_> { 662 fn preadjoint(&self) -> Self::Preadjoint<'_> {
616 DiagOp(self.0.preadjoint(), self.1.preadjoint()) 663 DiagOp(self.0.preadjoint(), self.1.preadjoint())
617 } 664 }
618 } 665 }
619 666
620 /// Block operator 667 /// Block operator
621 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>; 668 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
622
623 669
624 macro_rules! pairnorm { 670 macro_rules! pairnorm {
625 ($expj:ty) => { 671 ($expj:ty) => {
626 impl<F, A, B, S, T, ExpA, ExpB, ExpR> 672 impl<F, A, B, S, T, ExpA, ExpB, ExpR>
627 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> 673 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T>
628 for RowOp<S, T>
629 where 674 where
630 F : Float, 675 F: Float,
631 A : Space + Norm<F, ExpA>, 676 A: Space + Norm<F, ExpA>,
632 B : Space + Norm<F, ExpB>, 677 B: Space + Norm<F, ExpB>,
633 S : BoundedLinear<A, ExpA, ExpR, F>, 678 S: BoundedLinear<A, ExpA, ExpR, F>,
634 T : BoundedLinear<B, ExpB, ExpR, F>, 679 T: BoundedLinear<B, ExpB, ExpR, F>,
635 S::Codomain : Add<T::Codomain>, 680 S::Codomain: Add<T::Codomain>,
636 <S::Codomain as Add<T::Codomain>>::Output : Space, 681 <S::Codomain as Add<T::Codomain>>::Output: Space,
637 ExpA : NormExponent, 682 ExpA: NormExponent,
638 ExpB : NormExponent, 683 ExpB: NormExponent,
639 ExpR : NormExponent, 684 ExpR: NormExponent,
640 { 685 {
641 fn opnorm_bound( 686 fn opnorm_bound(
642 &self, 687 &self,
643 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>, 688 PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>,
644 expr : ExpR 689 expr: ExpR,
645 ) -> F { 690 ) -> F {
646 // An application of the triangle inequality bounds the norm by the maximum 691 // An application of the triangle inequality bounds the norm by the maximum
647 // of the individual norms. A simple observation shows this to be exact. 692 // of the individual norms. A simple observation shows this to be exact.
648 let na = self.0.opnorm_bound(expa, expr); 693 let na = self.0.opnorm_bound(expa, expr);
649 let nb = self.1.opnorm_bound(expb, expr); 694 let nb = self.1.opnorm_bound(expb, expr);
650 na.max(nb) 695 na.max(nb)
651 } 696 }
652 } 697 }
653 698
654 impl<F, A, S, T, ExpA, ExpS, ExpT> 699 impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
655 BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> 700 for ColOp<S, T>
656 for ColOp<S, T>
657 where 701 where
658 F : Float, 702 F: Float,
659 A : Space + Norm<F, ExpA>, 703 A: Space + Norm<F, ExpA>,
660 S : BoundedLinear<A, ExpA, ExpS, F>, 704 S: BoundedLinear<A, ExpA, ExpS, F>,
661 T : BoundedLinear<A, ExpA, ExpT, F>, 705 T: BoundedLinear<A, ExpA, ExpT, F>,
662 ExpA : NormExponent, 706 ExpA: NormExponent,
663 ExpS : NormExponent, 707 ExpS: NormExponent,
664 ExpT : NormExponent, 708 ExpT: NormExponent,
665 { 709 {
666 fn opnorm_bound( 710 fn opnorm_bound(
667 &self, 711 &self,
668 expa : ExpA, 712 expa: ExpA,
669 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj> 713 PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>,
670 ) -> F { 714 ) -> F {
671 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, 715 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
672 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} 716 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
673 let ns = self.0.opnorm_bound(expa, exps); 717 let ns = self.0.opnorm_bound(expa, exps);
674 let nt = self.1.opnorm_bound(expa, expt); 718 let nt = self.1.opnorm_bound(expa, expt);
675 ns.max(nt) 719 ns.max(nt)
676 } 720 }
677 } 721 }
678 } 722 };
679 } 723 }
680 724
681 pairnorm!(L1); 725 pairnorm!(L1);
682 pairnorm!(L2); 726 pairnorm!(L2);
683 pairnorm!(Linfinity); 727 pairnorm!(Linfinity);
684

Mercurial Repository: alg_tools

mercurial