Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 57
1b3b1687b9ed
parent 13
465fa2121ccb
child 59
9226980e45a7
equal deleted inserted replaced
56:5e3c1874797d 57:1b3b1687b9ed
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 crate::types::*; 7 use crate::types::*;
8 use serde::Serialize; 8 use serde::Serialize;
9 pub use crate::mapping::Apply; 9 pub use crate::mapping::Apply;
10 use crate::direct_product::Pair;
10 11
11 /// Trait for linear operators on `X`. 12 /// Trait for linear operators on `X`.
12 pub trait Linear<X> : Apply<X, Output=Self::Codomain> 13 pub trait Linear<X> : Apply<X, Output=Self::Codomain>
13 + for<'a> Apply<&'a X, Output=Self::Codomain> { 14 + for<'a> Apply<&'a X, Output=Self::Codomain> {
14 type Codomain; 15 type Codomain;
94 95
95 /// Form the preadjoint operator of `self`. 96 /// Form the preadjoint operator of `self`.
96 fn preadjoint(&self) -> Self::Preadjoint<'_>; 97 fn preadjoint(&self) -> Self::Preadjoint<'_>;
97 } 98 }
98 99
99 /// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$. 100 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
100 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {} 101 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
101 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {} 102 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
102 103
103 /// The identity operator 104 /// The identity operator
104 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 105 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
150 type AdjointCodomain=X; 151 type AdjointCodomain=X;
151 type Adjoint<'a> = IdOp<X> where X : 'a; 152 type Adjoint<'a> = IdOp<X> where X : 'a;
152 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 153 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
153 } 154 }
154 155
156 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
157 pub struct RowOp<S, T>(pub S, pub T);
158
159 use std::ops::Add;
160
161 impl<A, B, S, T> Apply<Pair<A, B>> for RowOp<S, T>
162 where
163 S : Apply<A>,
164 T : Apply<B>,
165 S::Output : Add<T::Output>
166 {
167 type Output = <S::Output as Add<T::Output>>::Output;
168
169 fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
170 self.0.apply(a) + self.1.apply(b)
171 }
172 }
173
174 impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for RowOp<S, T>
175 where
176 S : Apply<&'a A>,
177 T : Apply<&'a B>,
178 S::Output : Add<T::Output>
179 {
180 type Output = <S::Output as Add<T::Output>>::Output;
181
182 fn apply(&self, Pair(ref a, ref b) : &'a Pair<A, B>) -> Self::Output {
183 self.0.apply(a) + self.1.apply(b)
184 }
185 }
186
187 impl<A, B, S, T, D> Linear<Pair<A, B>> for RowOp<S, T>
188 where
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>,
197 T : GEMV<F, V, Y>,
198 F : Num,
199 Self : Linear<Pair<U, V>, Codomain=Y>
200 {
201 fn gemv(&self, y : &mut Y, α : F, x : &Pair<U, V>, β : F) {
202 self.0.gemv(y, α, &x.0, β);
203 self.1.gemv(y, α, &x.1, β);
204 }
205
206 fn apply_mut(&self, y : &mut Y, x : &Pair<U, V>){
207 self.0.apply_mut(y, &x.0);
208 self.1.apply_mut(y, &x.1);
209 }
210
211 /// Computes `y += Ax`, where `A` is `Self`
212 fn apply_add(&self, y : &mut Y, x : &Pair<U, V>){
213 self.0.apply_add(y, &x.0);
214 self.1.apply_add(y, &x.1);
215 }
216 }
217
218
219 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
220 pub struct ColOp<S, T>(pub S, pub T);
221
222 impl<A, S, T, O> Apply<A> for ColOp<S, T>
223 where
224 S : for<'a> Apply<&'a A, Output=O>,
225 T : Apply<A>,
226 A : std::borrow::Borrow<A>,
227 {
228 type Output = Pair<O, T::Output>;
229
230 fn apply(&self, a : A) -> Self::Output {
231 Pair(self.0.apply(a.borrow()), self.1.apply(a))
232 }
233 }
234
235 impl<A, S, T, D> Linear<A> for ColOp<S, T>
236 where
237 ColOp<S, T> : Apply<A, Output=D> + for<'a> Apply<&'a A, Output=D>,
238 {
239 type Codomain = D;
240 }
241
242 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
243 where
244 S : GEMV<F, X, A>,
245 T : GEMV<F, X, B>,
246 F : Num,
247 Self : Linear<X, Codomain=Pair<A, B>>
248 {
249 fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &X, β : F) {
250 self.0.gemv(&mut y.0, α, x, β);
251 self.1.gemv(&mut y.1, α, x, β);
252 }
253
254 fn apply_mut(&self, y : &mut Pair<A, B>, x : &X){
255 self.0.apply_mut(&mut y.0, x);
256 self.1.apply_mut(&mut y.1, x);
257 }
258
259 /// Computes `y += Ax`, where `A` is `Self`
260 fn apply_add(&self, y : &mut Pair<A, B>, x : &X){
261 self.0.apply_add(&mut y.0, x);
262 self.1.apply_add(&mut y.1, x);
263 }
264 }
265
266
267 impl<A, B, Yʹ, R, S, T> Adjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
268 where
269 S : Adjointable<A, Yʹ>,
270 T : Adjointable<B, Yʹ>,
271 Self : Linear<Pair<A, B>>,
272 for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Yʹ, Codomain=R>,
273 {
274 type AdjointCodomain = R;
275 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
276
277 fn adjoint(&self) -> Self::Adjoint<'_> {
278 ColOp(self.0.adjoint(), self.1.adjoint())
279 }
280 }
281
282 impl<A, B, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Yʹ> for RowOp<S, T>
283 where
284 S : Preadjointable<A, Yʹ>,
285 T : Preadjointable<B, Yʹ>,
286 Self : Linear<Pair<A, B>>,
287 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, Pair<A,B>, Codomain=R>,
288 {
289 type PreadjointCodomain = R;
290 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
291
292 fn preadjoint(&self) -> Self::Preadjoint<'_> {
293 ColOp(self.0.preadjoint(), self.1.preadjoint())
294 }
295 }
296
297
298 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
299 where
300 S : Adjointable<A, Xʹ>,
301 T : Adjointable<A, Yʹ>,
302 Self : Linear<A>,
303 for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
304 {
305 type AdjointCodomain = R;
306 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
307
308 fn adjoint(&self) -> Self::Adjoint<'_> {
309 RowOp(self.0.adjoint(), self.1.adjoint())
310 }
311 }
312
313 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
314 where
315 S : Preadjointable<A, Xʹ>,
316 T : Preadjointable<A, Yʹ>,
317 Self : Linear<A>,
318 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, A, Codomain=R>,
319 {
320 type PreadjointCodomain = R;
321 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
322
323 fn preadjoint(&self) -> Self::Preadjoint<'_> {
324 RowOp(self.0.preadjoint(), self.1.preadjoint())
325 }
326 }
327
328 /// Diagonal operator
329 pub struct DiagOp<S, T>(pub S, pub T);
330
331 impl<A, B, S, T> Apply<Pair<A, B>> for DiagOp<S, T>
332 where
333 S : Apply<A>,
334 T : Apply<B>,
335 {
336 type Output = Pair<S::Output, T::Output>;
337
338 fn apply(&self, Pair(a, b) : Pair<A, B>) -> Self::Output {
339 Pair(self.0.apply(a), self.1.apply(b))
340 }
341 }
342
343 impl<'a, A, B, S, T> Apply<&'a Pair<A, B>> for DiagOp<S, T>
344 where
345 S : Apply<&'a A>,
346 T : Apply<&'a B>,
347 {
348 type Output = Pair<S::Output, T::Output>;
349
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
362 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
363 where
364 S : GEMV<F, U, A>,
365 T : GEMV<F, V, B>,
366 F : Num,
367 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>
368 {
369 fn gemv(&self, y : &mut Pair<A, B>, α : F, x : &Pair<U, V>, β : F) {
370 self.0.gemv(&mut y.0, α, &x.0, β);
371 self.1.gemv(&mut y.1, α, &x.1, β);
372 }
373
374 fn apply_mut(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
375 self.0.apply_mut(&mut y.0, &x.0);
376 self.1.apply_mut(&mut y.1, &x.1);
377 }
378
379 /// Computes `y += Ax`, where `A` is `Self`
380 fn apply_add(&self, y : &mut Pair<A, B>, x : &Pair<U, V>){
381 self.0.apply_add(&mut y.0, &x.0);
382 self.1.apply_add(&mut y.1, &x.1);
383 }
384 }
385
386 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
387 where
388 S : Adjointable<A, Xʹ>,
389 T : Adjointable<B, Yʹ>,
390 Self : Linear<Pair<A, B>>,
391 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<Pair<Xʹ,Yʹ>, Codomain=R>,
392 {
393 type AdjointCodomain = R;
394 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
395
396 fn adjoint(&self) -> Self::Adjoint<'_> {
397 DiagOp(self.0.adjoint(), self.1.adjoint())
398 }
399 }
400
401 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>,Pair<Xʹ,Yʹ>> for DiagOp<S, T>
402 where
403 S : Preadjointable<A, Xʹ>,
404 T : Preadjointable<B, Yʹ>,
405 Self : Linear<Pair<A, B>>,
406 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, Pair<A, B>, Codomain=R>,
407 {
408 type PreadjointCodomain = R;
409 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
410
411 fn preadjoint(&self) -> Self::Preadjoint<'_> {
412 DiagOp(self.0.preadjoint(), self.1.preadjoint())
413 }
414 }
415
416 /// Block operator
417 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
418

Mercurial Repository: alg_tools

mercurial