Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 79
d63e40672dd6
parent 78
cebedc4a8331
equal deleted inserted replaced
78:cebedc4a8331 79:d63e40672dd6
10 use crate::direct_product::Pair; 10 use crate::direct_product::Pair;
11 use crate::instance::Instance; 11 use crate::instance::Instance;
12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm, HasDual}; 12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm, HasDual};
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 : AXPY> : Mapping<X, Codomain = Self::LinCodomain> {
16 { } 16 type LinCodomain : AXPY; // = Self::Codomain, but with AXPY enforced without trait bound bloat.
17 }
17 18
18 /// Efficient in-place summation. 19 /// Efficient in-place summation.
19 #[replace_float_literals(F::cast_from(literal))] 20 #[replace_float_literals(Self::Field::cast_from(literal))]
20 pub trait AXPY<F, X = Self> : Space 21 pub trait AXPY : Space {
22 type Field : Num;
23
24 /// Computes `y = βy + αx`, where `y` is `Self`.
25 fn axpy<I : Instance<Self>>(&mut self, α : Self::Field, x : I, β : Self::Field);
26
27 /// Copies `x` to `self`.
28 fn copy_from<I : Instance<Self>>(&mut self, x : I) {
29 self.axpy(1.0, x, 0.0)
30 }
31
32 /// Computes `y = αx`, where `y` is `Self`.
33 fn scale_from<I : Instance<Self>>(&mut self, α : Self::Field, x : I) {
34 self.axpy(α, x, 0.0)
35 }
36 }
37
38 macro_rules! impl_axpy {
39 ($($type:ty)*) => { $(
40 impl AXPY for $type {
41 type Field = $type;
42 fn axpy<I : Instance<Self>>(&mut self, α : $type, x : I, β : $type) {
43 *self = *self * α + x.own() * β;
44 }
45 }
46 )* };
47 }
48
49 impl_axpy!(i8 i16 i32 i64 i128 isize f32 f64);
50
51
52 /// Efficient in-place application for [`Linear`] operators.
53 #[replace_float_literals(X::Field::cast_from(literal))]
54 pub trait GEMV<X : AXPY> : Linear<X> {
55 /// Computes `y = αAx + βy`, where `A` is `Self`.
56 fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : <Self::LinCodomain as AXPY>::Field, x : I, β : <Self::LinCodomain as AXPY>::Field);
57
58 /// Computes `y = Ax`, where `A` is `Self`
59 fn apply_mut<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){
60 self.gemv(y, 1.0, x, 0.0)
61 }
62
63 /// Computes `y += Ax`, where `A` is `Self`
64 fn apply_add<I : Instance<X>>(&self, y : &mut Self::LinCodomain, x : I){
65 self.gemv(y, 1.0, x, 1.0)
66 }
67 }
68
69
70 /// Bounded linear operators
71 pub trait BoundedLinear<X, XExp, CodExp> : Linear<X>
21 where 72 where
22 F : Num, 73 F : Num,
23 X : Space, 74 X : AXP + Norm<F, XExp>,
24 {
25 /// Computes `y = βy + αx`, where `y` is `Self`.
26 fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F);
27
28 /// Copies `x` to `self`.
29 fn copy_from<I : Instance<X>>(&mut self, x : I) {
30 self.axpy(1.0, x, 0.0)
31 }
32
33 /// Computes `y = αx`, where `y` is `Self`.
34 fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) {
35 self.axpy(α, x, 0.0)
36 }
37 }
38
39 /// Efficient in-place application for [`Linear`] operators.
40 #[replace_float_literals(F::cast_from(literal))]
41 pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> {
42 /// Computes `y = αAx + βy`, where `A` is `Self`.
43 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F);
44
45 /// Computes `y = Ax`, where `A` is `Self`
46 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){
47 self.gemv(y, 1.0, x, 0.0)
48 }
49
50 /// Computes `y += Ax`, where `A` is `Self`
51 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){
52 self.gemv(y, 1.0, x, 1.0)
53 }
54 }
55
56
57 /// Bounded linear operators
58 pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X>
59 where
60 F : Num,
61 X : Space + Norm<F, XExp>,
62 XExp : NormExponent, 75 XExp : NormExponent,
63 CodExp : NormExponent 76 CodExp : NormExponent
64 { 77 {
65 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. 78 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
66 /// This is not expected to be the norm, just any bound on it that can be 79 /// This is not expected to be the norm, just any bound on it that can be
67 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm 80 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm
68 /// in `X`, and `codexp` in the codomain. 81 /// in `X`, and `codexp` in the codomain.
69 fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F; 82 fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> X::Field; // TODO: Should pick real subfield
70 } 83 }
71 84
72 // Linear operator application into mutable target. The [`AsRef`] bound 85 // Linear operator application into mutable target. The [`AsRef`] bound
73 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; 86 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
74 // the former is assumed to be e.g. a view into the latter. 87 // the former is assumed to be e.g. a view into the latter.
103 /// domain of the adjointed operator. `Self::Preadjoint` should be 116 /// domain of the adjointed operator. `Self::Preadjoint` should be
104 /// [`Adjointable`]`<'a,Ypre,X>`. 117 /// [`Adjointable`]`<'a,Ypre,X>`.
105 /// 118 ///
106 /// We do not set further restrictions on the spacds, to allow preadjointing when `X` 119 /// We do not set further restrictions on the spacds, to allow preadjointing when `X`
107 /// is on the dual space of `Ypre`, but a subset of it. 120 /// is on the dual space of `Ypre`, but a subset of it.
108 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> 121 pub trait Preadjointable<X : AXPY, Ypre : AXPY> : Linear<X>
109 //where 122 //where
110 // Ypre : HasDual<F, DualSpace=Self::Codomain>, 123 // Ypre : HasDual<F, DualSpace=Self::Codomain>,
111 { 124 {
112 type PreadjointCodomain : Space; 125 type PreadjointCodomain : AXPY;
113 type Preadjoint<'a> : Linear<Ypre, Codomain=Self::PreadjointCodomain> where Self : 'a; 126 type Preadjoint<'a> : Linear<Ypre, Codomain=Self::PreadjointCodomain> where Self : 'a;
114 127
115 /// Form the adjoint operator of `self`. 128 /// Form the adjoint operator of `self`.
116 fn preadjoint(&self) -> Self::Preadjoint<'_>; 129 fn preadjoint(&self) -> Self::Preadjoint<'_>;
117 } 130 }
118 131
119
120 /// The identity operator 132 /// The identity operator
121 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 133 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
122 pub struct IdOp<X> (PhantomData<X>); 134 pub struct IdOp<X : AXPY> (PhantomData<X>);
123 135
124 impl<X> IdOp<X> { 136 impl<X : AXPY> IdOp<X> {
125 pub fn new() -> IdOp<X> { IdOp(PhantomData) } 137 pub fn new() -> IdOp<X> { IdOp(PhantomData) }
126 } 138 }
127 139
128 impl<X : Clone + Space> Mapping<X> for IdOp<X> { 140 impl<X : Clone + AXPY> Mapping<X> for IdOp<X> {
129 type Codomain = X; 141 type Codomain = X;
130 142
131 fn apply<I : Instance<X>>(&self, x : I) -> X { 143 fn apply<I : Instance<X>>(&self, x : I) -> X {
132 x.own() 144 x.own()
133 } 145 }
134 } 146 }
135 147
136 impl<X : Clone + Space> Linear<X> for IdOp<X> 148 impl<X : Clone + AXPY> Linear<X> for IdOp<X> {
137 { } 149 type LinCodomain = X;
150 }
138 151
139 #[replace_float_literals(F::cast_from(literal))] 152 #[replace_float_literals(F::cast_from(literal))]
140 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> 153 impl<X : Clone + AXPY> GEMV<X> for IdOp<X> {
141 where
142 Y : AXPY<F, X>,
143 X : Clone + Space
144 {
145 // Computes `y = αAx + βy`, where `A` is `Self`. 154 // Computes `y = αAx + βy`, where `A` is `Self`.
146 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { 155 fn gemv<I : Instance<X>>(&self, y : &mut Self::LinCodomain, α : X::Field, x : I, β : X::Field) {
147 y.axpy(α, x, β) 156 y.axpy(α, x, β)
148 } 157 }
149 158
150 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 159 fn apply_mut<I : Instance<X>>(&self, y : &mut X, x : I){
151 y.copy_from(x); 160 y.copy_from(x);
152 } 161 }
153 } 162 }
154 163
155 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> 164 impl<X, E> BoundedLinear<X, E, E> for IdOp<X>
156 where 165 where
157 X : Space + Clone + Norm<F, E>, 166 X : AXPY + Clone + Norm<F, E>,
158 F : Num, 167 F : Num + AXPY,
159 E : NormExponent 168 E : NormExponent
160 { 169 {
161 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } 170 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> X::Field { F::ONE }
162 } 171 }
163 172
164 impl<X : Clone + HasDual<F, DualSpace=X>, F : Float> Adjointable<X, F> for IdOp<X> { 173 impl<X : Clone + HasDual<F, DualSpace=X>, F : Float> Adjointable<X, F> for IdOp<X> {
165 type AdjointCodomain = X; 174 type AdjointCodomain = X;
166 type Adjoint<'a> = IdOp<X::DualSpace> where X : 'a; 175 type Adjoint<'a> = IdOp<X::DualSpace> where X : 'a;
167 176
168 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 177 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
169 } 178 }
170 179
171 impl<X : Clone + Space> Preadjointable<X, X> for IdOp<X> { 180 impl<X : Clone + AXPY> Preadjointable<X, X> for IdOp<X> {
172 type PreadjointCodomain = X; 181 type PreadjointCodomain = X;
173 type Preadjoint<'a> = IdOp<X> where X : 'a; 182 type Preadjoint<'a> = IdOp<X> where X : 'a;
174 183
175 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() } 184 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
176 } 185 }
177 186
178 187
179 impl<S, T, E, X> Linear<X> for Composition<S, T, E> 188 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
180 where 189 where
181 X : Space, 190 X : AXPY,
182 T : Linear<X>, 191 T : Linear<X>,
183 S : Linear<T::Codomain> 192 S : Linear<T::LinCodomain>
184 { } 193 {
194 type LinCodomain = S::LinCodomain;
195 }
185 196
186 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E> 197 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
187 where 198 where
188 F : Num, 199 F : Num,
189 X : Space, 200 X : Space,
225 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. 236 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
226 pub struct RowOp<S, T>(pub S, pub T); 237 pub struct RowOp<S, T>(pub S, pub T);
227 238
228 use std::ops::Add; 239 use std::ops::Add;
229 240
230 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T> 241 impl<A, B, S, T, Y> Mapping<Pair<A, B>> for RowOp<S, T>
231 where 242 where
232 A : Space, 243 A : Space,
233 B : Space, 244 B : Space,
245 Y : Space,
234 S : Mapping<A>, 246 S : Mapping<A>,
235 T : Mapping<B>, 247 T : Mapping<B>,
236 S::Codomain : Add<T::Codomain>, 248 S::Codomain : Add<T::Codomain, Output = Y>,
237 <S::Codomain as Add<T::Codomain>>::Output : Space, 249 {
238 250 type Codomain = Y;
239 {
240 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
241 251
242 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { 252 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain {
243 let Pair(a, b) = x.decompose(); 253 let Pair(a, b) = x.decompose();
244 self.0.apply(a) + self.1.apply(b) 254 self.0.apply(a) + self.1.apply(b)
245 } 255 }
246 } 256 }
247 257
248 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> 258 impl<F, A, B, S, T, Y> Linear<Pair<A, B>> for RowOp<S, T>
249 where 259 where
250 A : Space, 260 F : Num + AXPY,
251 B : Space, 261 A : AXPY<Field=F>,
262 B : AXPY<Field=F>,
263 Y : AXPY,
252 S : Linear<A>, 264 S : Linear<A>,
253 T : Linear<B>, 265 T : Linear<B>,
254 S::Codomain : Add<T::Codomain>, 266 S::LinCodomain : Add<T::LinCodomain, Output=Y>,
255 <S::Codomain as Add<T::Codomain>>::Output : Space, 267 {
256 { } 268 type LinCodomain = Y;
257 269 }
258 270
259 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> 271
260 where 272 impl<'b, F, S, T, Y, G, U, V> GEMV<Pair<U, V>> for RowOp<S, T>
261 U : Space, 273 where
262 V : Space, 274 U : AXPY<Field=G>,
263 S : GEMV<F, U, Y>, 275 V : AXPY<Field=G>,
264 T : GEMV<F, V, Y>, 276 Y : AXPY<Field=F>,
265 F : Num, 277 S : GEMV<U>,
266 Self : Linear<Pair<U, V>, Codomain=Y> 278 T : GEMV<V>,
279 F : Num + AXPY,
280 G : Num + AXPY,
281 S::LinCodomain : Add<T::LinCodomain, Output=Y>,
267 { 282 {
268 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) { 283 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) {
269 let Pair(u, v) = x.decompose(); 284 let Pair(u, v) = x.decompose();
270 self.0.gemv(y, α, u, β); 285 self.0.gemv(y, α, u, β);
271 self.1.gemv(y, α, v, F::ONE); 286 self.1.gemv(y, α, v, F::ONE);
286 } 301 }
287 302
288 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. 303 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
289 pub struct ColOp<S, T>(pub S, pub T); 304 pub struct ColOp<S, T>(pub S, pub T);
290 305
291 impl<A, S, T> Mapping<A> for ColOp<S, T> 306 impl<X, S, T> Mapping<X> for ColOp<S, T>
292 where 307 where
293 A : Space, 308 X : AXPY,
294 S : Mapping<A>, 309 S : Mapping<X>,
295 T : Mapping<A>, 310 T : Mapping<X>,
296 { 311 {
297 type Codomain = Pair<S::Codomain, T::Codomain>; 312 type Codomain = Pair<S::Codomain, T::Codomain>;
298 313
299 fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain { 314 fn apply<I : Instance<X>>(&self, a : I) -> Self::Codomain {
300 Pair(self.0.apply(a.ref_instance()), self.1.apply(a)) 315 Pair(self.0.apply(a.ref_instance()), self.1.apply(a))
301 } 316 }
302 } 317 }
303 318
304 impl<A, S, T> Linear<A> for ColOp<S, T> 319 impl<X, S, T, A, B, F> Linear<X> for ColOp<S, T>
305 where 320 where
306 A : Space, 321 F : Num + AXPY,
307 S : Mapping<A>, 322 X : AXPY,
308 T : Mapping<A>, 323 S : Linear<X, LinCodomain=A>,
309 { } 324 T : Linear<X, LinCodomain=B>,
310 325 A : AXPY<Field=F>,
311 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> 326 B : AXPY<Field=F>,
312 where 327 {
313 X : Space, 328 type LinCodomain = Pair<S::LinCodomain, T::LinCodomain>;
314 S : GEMV<F, X, A>, 329 }
315 T : GEMV<F, X, B>, 330
316 F : Num, 331 impl<F, S, T, A, B, X> GEMV<X> for ColOp<S, T>
317 Self : Linear<X, Codomain=Pair<A, B>> 332 where
333 F : Num + AXPY,
334 X : AXPY,
335 S : GEMV<X, LinCodomain=A>,
336 T : GEMV<X, LinCodomain=B>,
337 A : AXPY<Field=F>,
338 B : AXPY<Field=F>,
318 { 339 {
319 fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 340 fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) {
320 self.0.gemv(&mut y.0, α, x.ref_instance(), β); 341 self.0.gemv(&mut y.0, α, x.ref_instance(), β);
321 self.1.gemv(&mut y.1, α, x, β); 342 self.1.gemv(&mut y.1, α, x, β);
322 } 343 }
339 F : Float, 360 F : Float,
340 A : HasDual<F>, 361 A : HasDual<F>,
341 B : HasDual<F>, 362 B : HasDual<F>,
342 S : Adjointable<A, F>, 363 S : Adjointable<A, F>,
343 T : Adjointable<B, F>, 364 T : Adjointable<B, F>,
344 Yʹ : Space, 365 Yʹ : AXPY,
345 S :: Codomain : HasDual<F, DualSpace=Yʹ>, 366 S :: Codomain : HasDual<F, DualSpace=Yʹ>,
346 T :: Codomain : HasDual<F, DualSpace=Yʹ>, 367 T :: Codomain : HasDual<F, DualSpace=Yʹ>,
347 S::Codomain : Add<T::Codomain>, 368 S::Codomain : Add<T::Codomain>,
348 <S::Codomain as Add<T::Codomain>>::Output : HasDual<F, DualSpace=Yʹ>, 369 <S::Codomain as Add<T::Codomain>>::Output : HasDual<F, DualSpace=Yʹ>,
349 Self::Codomain : HasDual<F, DualSpace=Yʹ>, 370 Self::Codomain : HasDual<F, DualSpace=Yʹ>,
359 fn adjoint(&self) -> Self::Adjoint<'_> { 380 fn adjoint(&self) -> Self::Adjoint<'_> {
360 ColOp(self.0.adjoint(), self.1.adjoint()) 381 ColOp(self.0.adjoint(), self.1.adjoint())
361 } 382 }
362 } 383 }
363 384
364 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 385 impl<F, A, B, Y, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T>
365 where 386 wherea
366 A : Space, 387 A : AXPY,
367 B : Space, 388 B : AXPY,
368 Yʹ : Space, 389 Yʹ : AXPY,
369 S : Preadjointable<A, Yʹ>, 390 S : Preadjointable<A, Yʹ>,
370 T : Preadjointable<B, Yʹ>, 391 T : Preadjointable<B, Yʹ>,
371 S::Codomain : Add<T::Codomain>, 392 S::Codomain : Add<T::Codomain>,
372 <S::Codomain as Add<T::Codomain>>::Output : Space, 393 <S::Codomain as Add<T::Codomain>>::Output : AXPY,
373 //Self : Linear<Pair<A, B>, Codomain=Y>, 394 //Self : Linear<Pair<A, B>, Codomain=Y>,
374 //for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, F>, 395 //for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, F>,
375 { 396 {
376 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; 397 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
377 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 398 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
405 } 426 }
406 } 427 }
407 428
408 impl<A, Aʹ, Xʹ, Yʹ, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 429 impl<A, Aʹ, Xʹ, Yʹ, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
409 where 430 where
410 A : Space, 431 A : AXPY,
411 Aʹ : Space, 432 Aʹ : AXPY,
412 Xʹ : Space, 433 Xʹ : AXPY,
413 Yʹ : Space, 434 Yʹ : AXPY,
414 S : Preadjointable<A, Xʹ, PreadjointCodomain=Aʹ>, 435 S : Preadjointable<A, Xʹ, PreadjointCodomain=Aʹ>,
415 T : Preadjointable<A, Yʹ, PreadjointCodomain=Aʹ>, 436 T : Preadjointable<A, Yʹ, PreadjointCodomain=Aʹ>,
416 Aʹ : ClosedAdd, 437 Aʹ : ClosedAdd,
417 //for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, F>, 438 //for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, F>,
418 { 439 {
427 /// Diagonal operator 448 /// Diagonal operator
428 pub struct DiagOp<S, T>(pub S, pub T); 449 pub struct DiagOp<S, T>(pub S, pub T);
429 450
430 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T> 451 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
431 where 452 where
432 A : Space, 453 A : AXPY,
433 B : Space, 454 B : AXPY,
434 S : Mapping<A>, 455 S : Mapping<A>,
435 T : Mapping<B>, 456 T : Mapping<B>,
436 { 457 {
437 type Codomain = Pair<S::Codomain, T::Codomain>; 458 type Codomain = Pair<S::Codomain, T::Codomain>;
438 459
440 let Pair(a, b) = x.decompose(); 461 let Pair(a, b) = x.decompose();
441 Pair(self.0.apply(a), self.1.apply(b)) 462 Pair(self.0.apply(a), self.1.apply(b))
442 } 463 }
443 } 464 }
444 465
445 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> 466 impl<F, A, B, S, T, G, U, V> Linear<Pair<A, B>> for DiagOp<S, T>
446 where 467 where
447 A : Space, 468 F : Num + AXPY,
448 B : Space, 469 G : Num + AXPY,
449 S : Linear<A>, 470 A : AXPY<Field=F>,
450 T : Linear<B>, 471 B : AXPY<Field=F>,
451 { } 472 U : AXPY<Field=G>,
452 473 V : AXPY<Field=G>,
453 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> 474 S : Linear<A, LinCodomain=U>,
454 where 475 T : Linear<B, LinCodomain=V>,
455 A : Space, 476 {
456 B : Space, 477 type LinCodomain = Pair<U, V>;
457 U : Space, 478 }
458 V : Space, 479
459 S : GEMV<F, U, A>, 480 impl<F, A, B, S, T, G, U, V> GEMV<Pair<U, V>> for DiagOp<S, T>
460 T : GEMV<F, V, B>, 481 where
461 F : Num, 482 F : Num + AXPY,
462 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>, 483 G : Num + AXPY,
463 { 484 A : AXPY<Field=F>,
464 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 485 B : AXPY<Field=F>,
486 U : AXPY<Field=G>,
487 V : AXPY<Field=G>,
488 S : GEMV<A, LinCodomain=U>,
489 T : GEMV<B, LinCodomain=V>,
490 {
491 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : G, x : I, β : G) {
465 let Pair(u, v) = x.decompose(); 492 let Pair(u, v) = x.decompose();
466 self.0.gemv(&mut y.0, α, u, β); 493 self.0.gemv(&mut y.0, α, u, β);
467 self.1.gemv(&mut y.1, α, v, β); 494 self.1.gemv(&mut y.1, α, v, β);
468 } 495 }
469 496
499 } 526 }
500 } 527 }
501 528
502 impl<A, B, Xʹ, Yʹ, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 529 impl<A, B, Xʹ, Yʹ, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
503 where 530 where
504 A : Space, 531 A : AXPY,
505 B : Space, 532 B : AXPY,
506 Xʹ : Space, 533 Xʹ : AXPY,
507 Yʹ : Space, 534 Yʹ : AXPY,
508 S : Preadjointable<A, Xʹ>, 535 S : Preadjointable<A, Xʹ>,
509 T : Preadjointable<B, Yʹ>, 536 T : Preadjointable<B, Yʹ>,
510 { 537 {
511 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; 538 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
512 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 539 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
521 548
522 549
523 macro_rules! pairnorm { 550 macro_rules! pairnorm {
524 ($expj:ty) => { 551 ($expj:ty) => {
525 impl<F, A, B, S, T, ExpA, ExpB, ExpR> 552 impl<F, A, B, S, T, ExpA, ExpB, ExpR>
526 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> 553 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR>
527 for RowOp<S, T> 554 for RowOp<S, T>
528 where 555 where
529 F : Float, 556 F : Float,
530 A : Space + Norm<F, ExpA>, 557 A : AXPY + Norm<F, ExpA>,
531 B : Space + Norm<F, ExpB>, 558 B : AXPY + Norm<F, ExpB>,
532 S : BoundedLinear<A, ExpA, ExpR, F>, 559 S : BoundedLinear<A, ExpA, ExpR, F>,
533 T : BoundedLinear<B, ExpB, ExpR, F>, 560 T : BoundedLinear<B, ExpB, ExpR, F>,
534 S::Codomain : Add<T::Codomain>, 561 S::Codomain : Add<T::Codomain>,
535 <S::Codomain as Add<T::Codomain>>::Output : Space, 562 <S::Codomain as Add<T::Codomain>>::Output : AXPY,
536 ExpA : NormExponent, 563 ExpA : NormExponent,
537 ExpB : NormExponent, 564 ExpB : NormExponent,
538 ExpR : NormExponent, 565 ExpR : NormExponent,
539 { 566 {
540 fn opnorm_bound( 567 fn opnorm_bound(
541 &self, 568 &self,
542 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>, 569 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>,
543 expr : ExpR 570 expr : ExpR
544 ) -> F { 571 ) -> <Pair<A, B> as AXPY>::Field {
545 // An application of the triangle inequality bounds the norm by the maximum 572 // An application of the triangle inequality bounds the norm by the maximum
546 // of the individual norms. A simple observation shows this to be exact. 573 // of the individual norms. A simple observation shows this to be exact.
547 let na = self.0.opnorm_bound(expa, expr); 574 let na = self.0.opnorm_bound(expa, expr);
548 let nb = self.1.opnorm_bound(expb, expr); 575 let nb = self.1.opnorm_bound(expb, expr);
549 na.max(nb) 576 na.max(nb)
550 } 577 }
551 } 578 }
552 579
553 impl<F, A, S, T, ExpA, ExpS, ExpT> 580 impl<A, S, T, ExpA, ExpS, ExpT>
554 BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> 581 BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>>
555 for ColOp<S, T> 582 for ColOp<S, T>
556 where 583 where
557 F : Float, 584 F : Float,
558 A : Space + Norm<F, ExpA>, 585 A : AXPY + Norm<F, ExpA>,
559 S : BoundedLinear<A, ExpA, ExpS, F>, 586 S : BoundedLinear<A, ExpA, ExpS, F>,
560 T : BoundedLinear<A, ExpA, ExpT, F>, 587 T : BoundedLinear<A, ExpA, ExpT, F>,
561 ExpA : NormExponent, 588 ExpA : NormExponent,
562 ExpS : NormExponent, 589 ExpS : NormExponent,
563 ExpT : NormExponent, 590 ExpT : NormExponent,
564 { 591 {
565 fn opnorm_bound( 592 fn opnorm_bound(
566 &self, 593 &self,
567 expa : ExpA, 594 expa : ExpA,
568 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj> 595 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj>
569 ) -> F { 596 ) -> A::Field {
570 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, 597 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
571 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} 598 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
572 let ns = self.0.opnorm_bound(expa, exps); 599 let ns = self.0.opnorm_bound(expa, exps);
573 let nt = self.1.opnorm_bound(expa, expt); 600 let nt = self.1.opnorm_bound(expa, expt);
574 ns.max(nt) 601 ns.max(nt)

Mercurial Repository: alg_tools

mercurial