aeon_tk/element/
approx.rs

1use crate::element::{Basis, BasisFunction as _, LeastSquares, LinearOperator, Support};
2use faer::{
3    ColRef, Conj, Mat, MatRef,
4    col::generic::Col,
5    linalg::{matmul::dot, svd::SvdError},
6};
7use reborrow::{Reborrow, ReborrowMut};
8
9#[derive(Clone)]
10pub struct ApproxOperator {
11    lls: LeastSquares,
12    /// Cache for a vandermonde matrix.
13    vandermonde: Mat<f64>,
14    /// Cache for rhs of computation.
15    rhs: Mat<f64>,
16    /// Shape functions for each operation.
17    shape: Mat<f64>,
18}
19
20impl ApproxOperator {
21    pub fn build<const N: usize>(
22        &mut self,
23        support: &impl Support<N>,
24        basis: &impl Basis<N>,
25        ops: &impl LinearOperator<N>,
26    ) -> Result<(), SvdError> {
27        self.vandermonde
28            .resize_with(support.num_points(), basis.order(), |i, j| {
29                let point = support.point(i);
30                basis.func(j).value(point)
31            });
32
33        self.shape
34            .resize_with(support.num_points(), ops.num_operations(), |_, _| 0.0);
35
36        self.rhs
37            .resize_with(basis.order(), ops.num_operations(), |deg, i| {
38                ops.apply(i, &basis.func(deg))
39            });
40
41        self.lls.least_squares(
42            self.vandermonde.transpose(),
43            self.shape.rb_mut(),
44            self.rhs.rb(),
45        )?;
46
47        Ok(())
48    }
49
50    pub fn apply(&self, i: usize, values: &[f64]) -> f64 {
51        dot::inner_prod(
52            self.shape.row(i),
53            Conj::Yes,
54            Col::from_slice(values),
55            Conj::No,
56        )
57    }
58
59    pub fn weights(&self, i: usize) -> ColRef<'_, f64> {
60        self.shape.col(i)
61    }
62
63    pub fn shape(&self) -> MatRef<'_, f64> {
64        self.shape.rb()
65    }
66}
67
68impl Default for ApproxOperator {
69    fn default() -> Self {
70        Self {
71            lls: Default::default(),
72            vandermonde: Mat::zeros(0, 0),
73            rhs: Mat::zeros(0, 0),
74            shape: Mat::zeros(0, 0),
75        }
76    }
77}
aeon_tk/element/approx.rs

aeon_tk/element/
approx.rs
