Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

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

Mercurial Repository: alg_tools

mercurial