src/linops.rs

Wed, 02 Nov 2022 00:11:49 +0200

author
Tuomo Valkonen <tuomov@iki.fi>
date
Wed, 02 Nov 2022 00:11:49 +0200
changeset 11
fc50a2d39053
parent 8
4e09b7829b51
child 13
465fa2121ccb
permissions
-rw-r--r--

work-units-sketching

0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
1 /*!
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
2 Abstract linear operators.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
3 */
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
4
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
5 use numeric_literals::replace_float_literals;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
6 use std::marker::PhantomData;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
7 use crate::types::*;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
8 use serde::Serialize;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
9
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
10 /// Trait for linear operators on `X`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
11 pub trait Linear<X> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
12 /// The range space of the operator.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
13 type Codomain;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
14 /// Apply the linear operator to `x`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
15 fn apply(&self, x : &X) -> Self::Codomain;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
16 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
17
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
18 /// Efficient in-place summation.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
19 #[replace_float_literals(F::cast_from(literal))]
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
20 pub trait AXPY<F : Num, X = Self> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
21 /// Computes `y = βy + αx`, where `y` is `Self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
22 fn axpy(&mut self, α : F, x : &X, β : F);
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
23
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
24 /// Copies `x` to `self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
25 fn copy_from(&mut self, x : &X) {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
26 self.axpy(1.0, x, 0.0)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
27 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
28
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
29 /// Computes `y = αx`, where `y` is `Self`.
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
30 fn scale_from(&mut self, α : F, x : &X) {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
31 self.axpy(α, x, 0.0)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
32 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
33 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
34
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
35 /// Efficient in-place application for [`Linear`] operators.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
36 #[replace_float_literals(F::cast_from(literal))]
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
37 pub trait GEMV<F : Num, X, Y = <Self as Linear<X>>::Codomain> : Linear<X> {
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
38 /// Computes `y = αAx + βy`, where `A` is `Self`.
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
39 fn gemv(&self, y : &mut Y, α : F, x : &X, β : F);
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
40
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
41 /// Computes `y = Ax`, where `A` is `Self`
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
42 fn apply_mut(&self, y : &mut Y, x : &X){
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
43 self.gemv(y, 1.0, x, 0.0)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
44 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
45
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
46 /// Computes `y += Ax`, where `A` is `Self`
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
47 fn apply_add(&self, y : &mut Y, x : &X){
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
48 self.gemv(y, 1.0, x, 1.0)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
49 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
50 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
51
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
52
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
53 /// Bounded linear operators
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
54 pub trait BoundedLinear<X> : Linear<X> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
55 type FloatType : Float;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
56 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
57 /// This is not expected to be the norm, just any bound on it that can be
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
58 /// reasonably implemented.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
59 fn opnorm_bound(&self) -> Self::FloatType;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
60 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
61
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
62 // Linear operator application into mutable target. The [`AsRef`] bound
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
63 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
64 // the former is assumed to be e.g. a view into the latter.
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
65
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
66 /*impl<X,Y,T> Fn(&X) -> Y for T where T : Linear<X,Codomain=Y> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
67 fn call(&self, x : &X) -> Y {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
68 self.apply(x)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
69 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
70 }*/
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
71
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
72 /// Trait for forming the adjoint operator of `Self`.
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
73 pub trait Adjointable<X,Yʹ> : Linear<X> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
74 type AdjointCodomain;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
75 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
76
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
77 /// Form the adjoint operator of `self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
78 fn adjoint(&self) -> Self::Adjoint<'_>;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
79
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
80 /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
81 self.adjoint().apply(y)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
82 }*/
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
83 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
84
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
85 /// Trait for forming a preadjoint of an operator.
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
86 ///
8
4e09b7829b51 Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents: 5
diff changeset
87 /// For an operator $A$ this is an operator $A\_\*$
4e09b7829b51 Multithreaded bisection tree operations
Tuomo Valkonen <tuomov@iki.fi>
parents: 5
diff changeset
88 /// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self`
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
89 /// operator. The space `Ypre` is the predual of its codomain, and should be the
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
90 /// domain of the adjointed operator. `Self::Preadjoint` should be
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
91 /// [`Adjointable`]`<'a,Ypre,X>`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
92 pub trait Preadjointable<X,Ypre> : Linear<X> {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
93 type PreadjointCodomain;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
94 type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
95
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
96 /// Form the preadjoint operator of `self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
97 fn preadjoint(&self) -> Self::Preadjoint<'_>;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
98 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
99
5
59dc4c5883f4 Improve documentation
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
100 /// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$.
0
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
101 pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
102 impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
103
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
104 /// The identity operator
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
105 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
106 pub struct IdOp<X : Clone> (PhantomData<X>);
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
107
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
108 impl<X> IdOp<X> where X : Clone {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
109 fn new() -> IdOp<X> { IdOp(PhantomData) }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
110 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
111
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
112 impl<X> Linear<X> for IdOp<X> where X : Clone {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
113 type Codomain = X;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
114 fn apply(&self, x : &X) -> X {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
115 x.clone()
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
116 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
117 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
118
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
119 #[replace_float_literals(F::cast_from(literal))]
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
120 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
121 // Computes `y = αAx + βy`, where `A` is `Self`.
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
122 fn gemv(&self, y : &mut Y, α : F, x : &X, β : F) {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
123 y.axpy(α, x, β)
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
124 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
125
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
126 fn apply_mut(&self, y : &mut Y, x : &X){
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
127 y.copy_from(x);
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
128 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
129 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
130
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
131 impl<X> BoundedLinear<X> for IdOp<X> where X : Clone {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
132 type FloatType = float;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
133 fn opnorm_bound(&self) -> float { 1.0 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
134 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
135
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
136 impl<X> Adjointable<X,X> for IdOp<X> where X : Clone {
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
137 type AdjointCodomain=X;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
138 type Adjoint<'a> = IdOp<X> where X : 'a;
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
139 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
140 }
9f27689eb130 Initialise new clean repository
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
141

mercurial