Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 59
9226980e45a7
parent 57
1b3b1687b9ed
child 61
05089fbc0310
equal deleted inserted replaced
58:1a38447a89fa 59:9226980e45a7
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};
9 pub use crate::mapping::Apply;
10 use crate::direct_product::Pair; 10 use crate::direct_product::Pair;
11 use crate::instance::Instance;
12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity};
11 13
12 /// Trait for linear operators on `X`. 14 /// Trait for linear operators on `X`.
13 pub trait Linear<X> : Apply<X, Output=Self::Codomain> 15 pub trait Linear<X : Space> : Mapping<X>
14 + for<'a> Apply<&'a X, Output=Self::Codomain> { 16 { }
15 type Codomain;
16 }
17 17
18 /// Efficient in-place summation. 18 /// Efficient in-place summation.
19 #[replace_float_literals(F::cast_from(literal))] 19 #[replace_float_literals(F::cast_from(literal))]
20 pub trait AXPY<F : Num, X = Self> { 20 pub trait AXPY<F, X = Self> : Space
21 where
22 F : Num,
23 X : Space,
24 {
21 /// Computes `y = βy + αx`, where `y` is `Self`. 25 /// Computes `y = βy + αx`, where `y` is `Self`.
22 fn axpy(&mut self, α : F, x : &X, β : F); 26 fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
23 27
24 /// Copies `x` to `self`. 28 /// Copies `x` to `self`.
25 fn copy_from(&mut self, x : &X) { 29 fn copy_from<I : Instance<X>>(&mut self, x : I) {
26 self.axpy(1.0, x, 0.0) 30 self.axpy(1.0, x, 0.0)
27 } 31 }
28 32
29 /// Computes `y = αx`, where `y` is `Self`. 33 /// Computes `y = αx`, where `y` is `Self`.
30 fn scale_from(&mut self, α : F, x : &X) { 34 fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) {
31 self.axpy(α, x, 0.0) 35 self.axpy(α, x, 0.0)
32 } 36 }
33 } 37 }
34 38
35 /// Efficient in-place application for [`Linear`] operators. 39 /// Efficient in-place application for [`Linear`] operators.
36 #[replace_float_literals(F::cast_from(literal))] 40 #[replace_float_literals(F::cast_from(literal))]
37 pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> { 41 pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> {
38 /// Computes `y = αAx + βy`, where `A` is `Self`. 42 /// Computes `y = αAx + βy`, where `A` is `Self`.
39 fn gemv(&self, y : &mut Y, α : F, x : &X, β : F); 43 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
40 44
41 /// Computes `y = Ax`, where `A` is `Self` 45 /// Computes `y = Ax`, where `A` is `Self`
42 fn apply_mut(&self, y : &mut Y, x : &X){ 46 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
43 self.gemv(y, 1.0, x, 0.0) 47 self.gemv(y, 1.0, x, 0.0)
44 } 48 }
45 49
46 /// Computes `y += Ax`, where `A` is `Self` 50 /// Computes `y += Ax`, where `A` is `Self`
47 fn apply_add(&self, y : &mut Y, x : &X){ 51 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
48 self.gemv(y, 1.0, x, 1.0) 52 self.gemv(y, 1.0, x, 1.0)
49 } 53 }
50 } 54 }
51 55
52 56
53 /// Bounded linear operators 57 /// Bounded linear operators
54 pub trait BoundedLinear<X> : Linear<X> { 58 pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
55 type FloatType : Float; 59 where
60 F : Num,
61 X : Space,
62 XExp : NormExponent,
63 CodExp : NormExponent
64 {
56 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. 65 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
57 /// This is not expected to be the norm, just any bound on it that can be 66 /// This is not expected to be the norm, just any bound on it that can be
58 /// reasonably implemented. 67 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm
59 fn opnorm_bound(&self) -> Self::FloatType; 68 /// in `X`, and `codexp` in the codomain.
69 fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F;
60 } 70 }
61 71
62 // Linear operator application into mutable target. The [`AsRef`] bound 72 // Linear operator application into mutable target. The [`AsRef`] bound
63 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; 73 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
64 // the former is assumed to be e.g. a view into the latter. 74 // the former is assumed to be e.g. a view into the latter.
68 self.apply(x) 78 self.apply(x)
69 } 79 }
70 }*/ 80 }*/
71 81
72 /// Trait for forming the adjoint operator of `Self`. 82 /// Trait for forming the adjoint operator of `Self`.
73 pub trait Adjointable<X,Yʹ> : Linear<X> { 83 pub trait Adjointable<X, Yʹ> : Linear<X>
74 type AdjointCodomain; 84 where
85 X : Space,
86 Yʹ : Space,
87 {
88 type AdjointCodomain : Space;
75 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; 89 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
76 90
77 /// Form the adjoint operator of `self`. 91 /// Form the adjoint operator of `self`.
78 fn adjoint(&self) -> Self::Adjoint<'_>; 92 fn adjoint(&self) -> Self::Adjoint<'_>;
79
80 /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain {
81 self.adjoint().apply(y)
82 }*/
83 } 93 }
84 94
85 /// Trait for forming a preadjoint of an operator. 95 /// Trait for forming a preadjoint of an operator.
86 /// 96 ///
87 /// For an operator $A$ this is an operator $A\_\*$ 97 /// For an operator $A$ this is an operator $A\_\*$
88 /// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self` 98 /// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self`
89 /// operator. The space `Ypre` is the predual of its codomain, and should be the 99 /// operator. The space `Ypre` is the predual of its codomain, and should be the
90 /// domain of the adjointed operator. `Self::Preadjoint` should be 100 /// domain of the adjointed operator. `Self::Preadjoint` should be
91 /// [`Adjointable`]`<'a,Ypre,X>`. 101 /// [`Adjointable`]`<'a,Ypre,X>`.
92 pub trait Preadjointable<X,Ypre> : Linear<X> { 102
93 type PreadjointCodomain; 103 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> {
94 type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a; 104 type PreadjointCodomain : Space;
95 105 type Preadjoint<'a> : Adjointable<
96 /// Form the preadjoint operator of `self`. 106 Ypre, X, AdjointCodomain=Self::Codomain, Codomain=Self::PreadjointCodomain
107 > where Self : 'a;
108
109 /// Form the adjoint operator of `self`.
97 fn preadjoint(&self) -> Self::Preadjoint<'_>; 110 fn preadjoint(&self) -> Self::Preadjoint<'_>;
98 } 111 }
99 112
100 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$. 113 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
101 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {} 114 pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
102 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {} 115 impl<'a,X : Space, T> SimplyAdjointable<X> for T
116 where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
103 117
104 /// The identity operator 118 /// The identity operator
105 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 119 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
106 pub struct IdOp<X> (PhantomData<X>); 120 pub struct IdOp<X> (PhantomData<X>);
107 121
108 impl<X> IdOp<X> { 122 impl<X> IdOp<X> {
109 fn new() -> IdOp<X> { IdOp(PhantomData) } 123 pub fn new() -> IdOp<X> { IdOp(PhantomData) }
110 } 124 }
111 125
112 impl<X> Apply<X> for IdOp<X> { 126 impl<X : Clone + Space> Mapping<X> for IdOp<X> {
113 type Output = X;
114
115 fn apply(&self, x : X) -> X {
116 x
117 }
118 }
119
120 impl<'a, X> Apply<&'a X> for IdOp<X> where X : Clone {
121 type Output = X;
122
123 fn apply(&self, x : &'a X) -> X {
124 x.clone()
125 }
126 }
127
128 impl<X> Linear<X> for IdOp<X> where X : Clone {
129 type Codomain = X; 127 type Codomain = X;
130 } 128
131 129 fn apply<I : Instance<X>>(&self, x : I) -> X {
130 x.own()
131 }
132 }
133
134 impl<X : Clone + Space> Linear<X> for IdOp<X>
135 { }
132 136
133 #[replace_float_literals(F::cast_from(literal))] 137 #[replace_float_literals(F::cast_from(literal))]
134 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone { 138 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X>
139 where
140 Y : AXPY<F, X>,
141 X : Clone + Space
142 {
135 // Computes `y = αAx + βy`, where `A` is `Self`. 143 // Computes `y = αAx + βy`, where `A` is `Self`.
136 fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) { 144 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) {
137 y.axpy(α, x, β) 145 y.axpy(α, x, β)
138 } 146 }
139 147
140 fn apply_mut(&self, y : &mut Y, x : &X){ 148 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
141 y.copy_from(x); 149 y.copy_from(x);
142 } 150 }
143 } 151 }
144 152
145 impl<X> BoundedLinear<X> for IdOp<X> where X : Clone { 153 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
146 type FloatType = float; 154 where
147 fn opnorm_bound(&self) -> float { 1.0 } 155 X : Space + Clone,
148 } 156 F : Num,
149 157 E : NormExponent
150 impl<X> Adjointable<X,X> for IdOp<X> where X : Clone { 158 {
159 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
160 }
161
162 impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> {
151 type AdjointCodomain=X; 163 type AdjointCodomain=X;
152 type Adjoint<'a> = IdOp<X> where X : 'a; 164 type Adjoint<'a> = IdOp<X> where X : 'a;
165
153 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 166 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
167 }
168
169 impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> {
170 type PreadjointCodomain=X;
171 type Preadjoint<'a> = IdOp<X> where X : 'a;
172
173 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
154 } 174 }
155 175
156 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. 176 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
157 pub struct RowOp<S, T>(pub S, pub T); 177 pub struct RowOp<S, T>(pub S, pub T);
158 178
159 use std::ops::Add; 179 use std::ops::Add;
160 180
161 impl<A, B, S, T> Apply<Pair<A, B>> for RowOp<S, T> 181 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
162 where 182 where
163 S : Apply<A>, 183 A : Space,
164 T : Apply<B>, 184 B : Space,
165 S::Output : Add<T::Output> 185 S : Mapping<A>,
166 { 186 T : Mapping<B>,
167 type Output = <S::Output as Add<T::Output>>::Output; 187 S::Codomain : Add<T::Codomain>,
168 188 <S::Codomain as Add<T::Codomain>>::Output : Space,
169 fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output { 189
190 {
191 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
192
193 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
194 let Pair(a, b) = x.decompose();
170 self.0.apply(a) + self.1.apply(b) 195 self.0.apply(a) + self.1.apply(b)
171 } 196 }
172 } 197 }
173 198
174 impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for RowOp<S, T> 199 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
175 where 200 where
176 S : Apply<&'a A>, 201 A : Space,
177 T : Apply<&'a B>, 202 B : Space,
178 S::Output : Add<T::Output> 203 S : Linear<A>,
179 { 204 T : Linear<B>,
180 type Output = <S::Output as Add<T::Output>>::Output; 205 S::Codomain : Add<T::Codomain>,
181 206 <S::Codomain as Add<T::Codomain>>::Output : Space,
182 fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output { 207 { }
183 self.0.apply(a) + self.1.apply(b) 208
184 } 209
185 } 210 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
186 211 where
187 impl<A, B, S, T, D> Linear<Pair<A, B>> for RowOp<S, T> 212 U : Space,
188 where 213 V : Space,
189 RowOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a> Apply<&'a Pair<A, B>, Output=D>,
190 {
191 type Codomain = D;
192 }
193
194 impl<F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
195 where
196 S : GEMV<F, U, Y>, 214 S : GEMV<F, U, Y>,
197 T : GEMV<F, V, Y>, 215 T : GEMV<F, V, Y>,
198 F : Num, 216 F : Num,
199 Self : Linear<Pair<U, V>, Codomain=Y> 217 Self : Linear<Pair<U, V>, Codomain=Y>
200 { 218 {
201 fn gemv(&self, y : &mut Y, α : F, x : &Pair<U, V>, β : F) { 219 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
202 self.0.gemv(y, α, &x.0, β); 220 let Pair(u, v) = x.decompose();
203 self.1.gemv(y, α, &x.1, β); 221 self.0.gemv(y, α, u, β);
204 } 222 self.1.gemv(y, α, v, F::ONE);
205 223 }
206 fn apply_mut(&self, y : &mut Y, x : &Pair<U, V>){ 224
207 self.0.apply_mut(y, &x.0); 225 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
208 self.1.apply_mut(y, &x.1); 226 let Pair(u, v) = x.decompose();
227 self.0.apply_mut(y, u);
228 self.1.apply_mut(y, v);
209 } 229 }
210 230
211 /// Computes `y += Ax`, where `A` is `Self` 231 /// Computes `y += Ax`, where `A` is `Self`
212 fn apply_add(&self, y : &mut Y, x : &Pair<U, V>){ 232 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) {
213 self.0.apply_add(y, &x.0); 233 let Pair(u, v) = x.decompose();
214 self.1.apply_add(y, &x.1); 234 self.0.apply_add(y, u);
215 } 235 self.1.apply_add(y, v);
216 } 236 }
217 237 }
218 238
219 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. 239 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
220 pub struct ColOp<S, T>(pub S, pub T); 240 pub struct ColOp<S, T>(pub S, pub T);
221 241
222 impl<A, S, T, O> Apply<A> for ColOp<S, T> 242 impl<A, S, T> Mapping<A> for ColOp<S, T>
223 where 243 where
224 S : for<'a> Apply<&'a A, Output=O>, 244 A : Space,
225 T : Apply<A>, 245 S : Mapping<A>,
226 A : std::borrow::Borrow<A>, 246 T : Mapping<A>,
227 { 247 {
228 type Output = Pair<O, T::Output>; 248 type Codomain = Pair<S::Codomain, T::Codomain>;
229 249
230 fn apply(&self, a : A) -> Self::Output { 250 fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain {
231 Pair(self.0.apply(a.borrow()), self.1.apply(a)) 251 Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
232 } 252 }
233 } 253 }
234 254
235 impl<A, S, T, D> Linear<A> for ColOp<S, T> 255 impl<A, S, T> Linear<A> for ColOp<S, T>
236 where 256 where
237 ColOp<S, T> : Apply<A, Output=D> + for<'a> Apply<&'a A, Output=D>, 257 A : Space,
238 { 258 S : Mapping<A>,
239 type Codomain = D; 259 T : Mapping<A>,
240 } 260 { }
241 261
242 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> 262 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
243 where 263 where
264 X : Space,
244 S : GEMV<F, X, A>, 265 S : GEMV<F, X, A>,
245 T : GEMV<F, X, B>, 266 T : GEMV<F, X, B>,
246 F : Num, 267 F : Num,
247 Self : Linear<X, Codomain=Pair<A, B>> 268 Self : Linear<X, Codomain=Pair<A, B>>
248 { 269 {
249 fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &X, β : F) { 270 fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
250 self.0.gemv(&mut y.0, α, x, β); 271 self.0.gemv(&mut y.0, α, x.ref_instance(), β);
251 self.1.gemv(&mut y.1, α, x, β); 272 self.1.gemv(&mut y.1, α, x, β);
252 } 273 }
253 274
254 fn apply_mut(&self, y : &mut Pair<A, B>, x : &X){ 275 fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
255 self.0.apply_mut(&mut y.0, x); 276 self.0.apply_mut(&mut y.0, x.ref_instance());
256 self.1.apply_mut(&mut y.1, x); 277 self.1.apply_mut(&mut y.1, x);
257 } 278 }
258 279
259 /// Computes `y += Ax`, where `A` is `Self` 280 /// Computes `y += Ax`, where `A` is `Self`
260 fn apply_add(&self, y : &mut Pair<A, B>, x : &X){ 281 fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){
261 self.0.apply_add(&mut y.0, x); 282 self.0.apply_add(&mut y.0, x.ref_instance());
262 self.1.apply_add(&mut y.1, x); 283 self.1.apply_add(&mut y.1, x);
263 } 284 }
264 } 285 }
265 286
266 287
267 impl<A, B, Yʹ, R, S, T> Adjointable<Pair<A,B>,Yʹ> for RowOp<S, T> 288 impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
268 where 289 where
290 A : Space,
291 B : Space,
292 Yʹ : Space,
269 S : Adjointable<A, Yʹ>, 293 S : Adjointable<A, Yʹ>,
270 T : Adjointable<B, Yʹ>, 294 T : Adjointable<B, Yʹ>,
271 Self : Linear<Pair<A, B>>, 295 Self : Linear<Pair<A, B>>,
272 for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Yʹ, Codomain=R>, 296 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
273 { 297 // Yʹ,
274 type AdjointCodomain = R; 298 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
299 // >,
300 {
301 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
275 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 302 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
276 303
277 fn adjoint(&self) -> Self::Adjoint<'_> { 304 fn adjoint(&self) -> Self::Adjoint<'_> {
278 ColOp(self.0.adjoint(), self.1.adjoint()) 305 ColOp(self.0.adjoint(), self.1.adjoint())
279 } 306 }
280 } 307 }
281 308
282 impl<A, B, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Yʹ> for RowOp<S, T> 309 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
283 where 310 where
311 A : Space,
312 B : Space,
313 Yʹ : Space,
284 S : Preadjointable<A, Yʹ>, 314 S : Preadjointable<A, Yʹ>,
285 T : Preadjointable<B, Yʹ>, 315 T : Preadjointable<B, Yʹ>,
286 Self : Linear<Pair<A, B>>, 316 Self : Linear<Pair<A, B>>,
287 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, Pair<A,B>, Codomain=R>, 317 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
288 { 318 Yʹ, Pair<A,B>,
289 type PreadjointCodomain = R; 319 Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>,
320 AdjointCodomain = Self::Codomain,
321 >,
322 {
323 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
290 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 324 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
291 325
292 fn preadjoint(&self) -> Self::Preadjoint<'_> { 326 fn preadjoint(&self) -> Self::Preadjoint<'_> {
293 ColOp(self.0.preadjoint(), self.1.preadjoint()) 327 ColOp(self.0.preadjoint(), self.1.preadjoint())
294 } 328 }
295 } 329 }
296 330
297 331
298 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 332 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
299 where 333 where
300 S : Adjointable<A, Xʹ>, 334 A : Space,
301 T : Adjointable<A, Yʹ>, 335 Xʹ : Space,
336 Yʹ : Space,
337 R : Space + ClosedAdd,
338 S : Adjointable<A, Xʹ, AdjointCodomain = R>,
339 T : Adjointable<A, Yʹ, AdjointCodomain = R>,
302 Self : Linear<A>, 340 Self : Linear<A>,
303 for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>, 341 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
342 // Pair<Xʹ,Yʹ>,
343 // Codomain=R,
344 // >,
304 { 345 {
305 type AdjointCodomain = R; 346 type AdjointCodomain = R;
306 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 347 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
307 348
308 fn adjoint(&self) -> Self::Adjoint<'_> { 349 fn adjoint(&self) -> Self::Adjoint<'_> {
310 } 351 }
311 } 352 }
312 353
313 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 354 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
314 where 355 where
315 S : Preadjointable<A, Xʹ>, 356 A : Space,
316 T : Preadjointable<A, Yʹ>, 357 Xʹ : Space,
358 Yʹ : Space,
359 R : Space + ClosedAdd,
360 S : Preadjointable<A, Xʹ, PreadjointCodomain = R>,
361 T : Preadjointable<A, Yʹ, PreadjointCodomain = R>,
317 Self : Linear<A>, 362 Self : Linear<A>,
318 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, A, Codomain=R>, 363 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
364 Pair<Xʹ,Yʹ>, A,
365 Codomain = R,
366 AdjointCodomain = Self::Codomain,
367 >,
319 { 368 {
320 type PreadjointCodomain = R; 369 type PreadjointCodomain = R;
321 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 370 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
322 371
323 fn preadjoint(&self) -> Self::Preadjoint<'_> { 372 fn preadjoint(&self) -> Self::Preadjoint<'_> {
326 } 375 }
327 376
328 /// Diagonal operator 377 /// Diagonal operator
329 pub struct DiagOp<S, T>(pub S, pub T); 378 pub struct DiagOp<S, T>(pub S, pub T);
330 379
331 impl<A, B, S, T> Apply<Pair<A, B>> for DiagOp<S, T> 380 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
332 where 381 where
333 S : Apply<A>, 382 A : Space,
334 T : Apply<B>, 383 B : Space,
335 { 384 S : Mapping<A>,
336 type Output = Pair<S::Output, T::Output>; 385 T : Mapping<B>,
337 386 {
338 fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output { 387 type Codomain = Pair<S::Codomain, T::Codomain>;
388
389 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
390 let Pair(a, b) = x.decompose();
339 Pair(self.0.apply(a), self.1.apply(b)) 391 Pair(self.0.apply(a), self.1.apply(b))
340 } 392 }
341 } 393 }
342 394
343 impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for DiagOp<S, T> 395 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
344 where 396 where
345 S : Apply<&'a A>, 397 A : Space,
346 T : Apply<&'a B>, 398 B : Space,
347 { 399 S : Linear<A>,
348 type Output = Pair<S::Output, T::Output>; 400 T : Linear<B>,
349 401 { }
350 fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
351 Pair(self.0.apply(a), self.1.apply(b))
352 }
353 }
354
355 impl<A, B, S, T, D> Linear<Pair<A, B>> for DiagOp<S, T>
356 where
357 DiagOp<S, T> : Apply<Pair<A, B>, Output=D> + for<'a> Apply<&'a Pair<A, B>, Output=D>,
358 {
359 type Codomain = D;
360 }
361 402
362 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> 403 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
363 where 404 where
405 A : Space,
406 B : Space,
407 U : Space,
408 V : Space,
364 S : GEMV<F, U, A>, 409 S : GEMV<F, U, A>,
365 T : GEMV<F, V, B>, 410 T : GEMV<F, V, B>,
366 F : Num, 411 F : Num,
367 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>> 412 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>,
368 { 413 {
369 fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &Pair<U, V>, β : F) { 414 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
370 self.0.gemv(&mut y.0, α, &x.0, β); 415 let Pair(u, v) = x.decompose();
371 self.1.gemv(&mut y.1, α, &x.1, β); 416 self.0.gemv(&mut y.0, α, u, β);
372 } 417 self.1.gemv(&mut y.1, α, v, β);
373 418 }
374 fn apply_mut(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){ 419
375 self.0.apply_mut(&mut y.0, &x.0); 420 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
376 self.1.apply_mut(&mut y.1, &x.1); 421 let Pair(u, v) = x.decompose();
422 self.0.apply_mut(&mut y.0, u);
423 self.1.apply_mut(&mut y.1, v);
377 } 424 }
378 425
379 /// Computes `y += Ax`, where `A` is `Self` 426 /// Computes `y += Ax`, where `A` is `Self`
380 fn apply_add(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){ 427 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){
381 self.0.apply_add(&mut y.0, &x.0); 428 let Pair(u, v) = x.decompose();
382 self.1.apply_add(&mut y.1, &x.1); 429 self.0.apply_add(&mut y.0, u);
383 } 430 self.1.apply_add(&mut y.1, v);
384 } 431 }
385 432 }
386 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T> 433
387 where 434 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
435 where
436 A : Space,
437 B : Space,
438 Xʹ: Space,
439 Yʹ : Space,
440 R : Space,
388 S : Adjointable<A, Xʹ>, 441 S : Adjointable<A, Xʹ>,
389 T : Adjointable<B, Yʹ>, 442 T : Adjointable<B, Yʹ>,
390 Self : Linear<Pair<A, B>>, 443 Self : Linear<Pair<A, B>>,
391 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>, 444 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
445 Pair<Xʹ,Yʹ>, Codomain=R,
446 >,
392 { 447 {
393 type AdjointCodomain = R; 448 type AdjointCodomain = R;
394 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 449 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
395 450
396 fn adjoint(&self) -> Self::Adjoint<'_> { 451 fn adjoint(&self) -> Self::Adjoint<'_> {
397 DiagOp(self.0.adjoint(), self.1.adjoint()) 452 DiagOp(self.0.adjoint(), self.1.adjoint())
398 } 453 }
399 } 454 }
400 455
401 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T> 456 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
402 where 457 where
458 A : Space,
459 B : Space,
460 Xʹ: Space,
461 Yʹ : Space,
462 R : Space,
403 S : Preadjointable<A, Xʹ>, 463 S : Preadjointable<A, Xʹ>,
404 T : Preadjointable<B, Yʹ>, 464 T : Preadjointable<B, Yʹ>,
405 Self : Linear<Pair<A, B>>, 465 Self : Linear<Pair<A, B>>,
406 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, Pair<A, B>, Codomain=R>, 466 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
467 Pair<Xʹ,Yʹ>, Pair<A, B>,
468 Codomain=R,
469 AdjointCodomain = Self::Codomain,
470 >,
407 { 471 {
408 type PreadjointCodomain = R; 472 type PreadjointCodomain = R;
409 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 473 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
410 474
411 fn preadjoint(&self) -> Self::Preadjoint<'_> { 475 fn preadjoint(&self) -> Self::Preadjoint<'_> {
414 } 478 }
415 479
416 /// Block operator 480 /// Block operator
417 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>; 481 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
418 482
483
484 macro_rules! pairnorm {
485 ($expj:ty) => {
486 impl<F, A, B, S, T, ExpA, ExpB, ExpR>
487 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F>
488 for RowOp<S, T>
489 where
490 F : Float,
491 A : Space,
492 B : Space,
493 S : BoundedLinear<A, ExpA, ExpR, F>,
494 T : BoundedLinear<B, ExpB, ExpR, F>,
495 S::Codomain : Add<T::Codomain>,
496 <S::Codomain as Add<T::Codomain>>::Output : Space,
497 ExpA : NormExponent,
498 ExpB : NormExponent,
499 ExpR : NormExponent,
500 {
501 fn opnorm_bound(
502 &self,
503 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
504 expr : ExpR
505 ) -> F {
506 // An application of the triangle inequality bounds the norm by the maximum
507 // of the individual norms. A simple observation shows this to be exact.
508 let na = self.0.opnorm_bound(expa, expr);
509 let nb = self.1.opnorm_bound(expb, expr);
510 na.max(nb)
511 }
512 }
513
514 impl<F, A, S, T, ExpA, ExpS, ExpT>
515 BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
516 for ColOp<S, T>
517 where
518 F : Float,
519 A : Space,
520 S : BoundedLinear<A, ExpA, ExpS, F>,
521 T : BoundedLinear<A, ExpA, ExpT, F>,
522 ExpA : NormExponent,
523 ExpS : NormExponent,
524 ExpT : NormExponent,
525 {
526 fn opnorm_bound(
527 &self,
528 expa : ExpA,
529 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj>
530 ) -> F {
531 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
532 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
533 let ns = self.0.opnorm_bound(expa, exps);
534 let nt = self.1.opnorm_bound(expa, expt);
535 ns.max(nt)
536 }
537 }
538 }
539 }
540
541 pairnorm!(L1);
542 pairnorm!(L2);
543 pairnorm!(Linfinity);
544

Mercurial Repository: alg_tools

mercurial