1use faer::{
2 Accum, Conj, Par,
3 linalg::matmul::{dot::inner_prod, matmul},
4};
5use faer_traits::ComplexField;
6use mdarray::{Array, Dim, Layout, Shape, Slice};
7use mdarray_linalg::matvec::{MatVec, MatVecBuilder, Outer, OuterBuilder, VecOps};
8use num_complex::{Complex32, Complex64, ComplexFloat};
9use num_traits::{One, Zero};
10
11use crate::{Faer, into_faer, into_faer_col, into_faer_col_mut, into_faer_mut, into_faer_row};
12
13trait FaerVectorScalar: ComplexFloat + ComplexField {
14 fn from_faer_real(real: <Self as ComplexField>::Real) -> <Self as ComplexFloat>::Real;
15}
16
17impl FaerVectorScalar for f32 {
18 fn from_faer_real(real: <Self as ComplexField>::Real) -> <Self as ComplexFloat>::Real {
19 real
20 }
21}
22
23impl FaerVectorScalar for f64 {
24 fn from_faer_real(real: <Self as ComplexField>::Real) -> <Self as ComplexFloat>::Real {
25 real
26 }
27}
28
29impl FaerVectorScalar for Complex32 {
30 fn from_faer_real(real: <Self as ComplexField>::Real) -> <Self as ComplexFloat>::Real {
31 real
32 }
33}
34
35impl FaerVectorScalar for Complex64 {
36 fn from_faer_real(real: <Self as ComplexField>::Real) -> <Self as ComplexFloat>::Real {
37 real
38 }
39}
40
41struct FaerMatVecBuilder<'a, T, La, Lx, D0, D1>
42where
43 La: Layout,
44 Lx: Layout,
45 D0: Dim,
46 D1: Dim,
47{
48 alpha: T,
49 a: &'a Slice<T, (D0, D1), La>,
50 x: &'a Slice<T, (D1,), Lx>,
51 par: Par,
52}
53
54struct FaerOuterBuilder<'a, T, Lx, Ly, Dx, Dy>
55where
56 Lx: Layout,
57 Ly: Layout,
58 Dx: Dim,
59 Dy: Dim,
60{
61 alpha: T,
62 x: &'a Slice<T, (Dx,), Lx>,
63 y: &'a Slice<T, (Dy,), Ly>,
64 par: Par,
65}
66
67impl<'a, T, La, Lx, D0, D1> MatVecBuilder<'a, T, La, Lx, D0, D1>
68 for FaerMatVecBuilder<'a, T, La, Lx, D0, D1>
69where
70 La: Layout,
71 Lx: Layout,
72 D0: Dim,
73 D1: Dim,
74 T: FaerVectorScalar + One,
75{
76 fn scale(mut self, alpha: T) -> Self {
77 self.alpha *= alpha;
78 self
79 }
80
81 fn eval(self) -> Array<T, (D0,)> {
82 let m = self.a.shape().dim(0);
83 let n = self.a.shape().dim(1);
84 let x_len = self.x.shape().dim(0);
85 assert_eq!(n, x_len, "Matrix columns must match vector length");
86
87 let shape = <(D0,) as Shape>::from_dims(&[m]);
88 let mut y = Array::<T, (D0,)>::from_elem(shape, T::zero());
89 self.write(&mut y);
90 y
91 }
92
93 fn write<Ly: Layout>(self, y: &mut Slice<T, (D0,), Ly>) {
94 let m = self.a.shape().dim(0);
95 let n = self.a.shape().dim(1);
96 let x_len = self.x.shape().dim(0);
97 let y_len = y.shape().dim(0);
98 assert_eq!(n, x_len, "Matrix columns must match vector length");
99 assert_eq!(m, y_len, "Matrix rows must match output vector length");
100
101 matmul(
102 into_faer_col_mut(y),
103 Accum::Replace,
104 into_faer(self.a),
105 into_faer_col(self.x),
106 self.alpha,
107 self.par,
108 );
109 }
110
111 fn add_to_vec<Ly: Layout>(self, y: &mut Slice<T, (D0,), Ly>) {
112 let m = self.a.shape().dim(0);
113 let n = self.a.shape().dim(1);
114 let x_len = self.x.shape().dim(0);
115 let y_len = y.shape().dim(0);
116 assert_eq!(n, x_len, "Matrix columns must match vector length");
117 assert_eq!(m, y_len, "Matrix rows must match output vector length");
118
119 matmul(
120 into_faer_col_mut(y),
121 Accum::Add,
122 into_faer(self.a),
123 into_faer_col(self.x),
124 self.alpha,
125 self.par,
126 );
127 }
128
129 fn add_to_scaled_vec<Ly: Layout>(self, y: &mut Slice<T, (D0,), Ly>, beta: T) {
130 let m = self.a.shape().dim(0);
131 let n = self.a.shape().dim(1);
132 let x_len = self.x.shape().dim(0);
133 let y_len = y.shape().dim(0);
134 assert_eq!(n, x_len, "Matrix columns must match vector length");
135 assert_eq!(m, y_len, "Matrix rows must match output vector length");
136
137 for yi in y.iter_mut() {
138 *yi = beta * *yi;
139 }
140
141 matmul(
142 into_faer_col_mut(y),
143 Accum::Add,
144 into_faer(self.a),
145 into_faer_col(self.x),
146 self.alpha,
147 self.par,
148 );
149 }
150}
151
152impl<T, D0: Dim, D1: Dim> MatVec<T, D0, D1> for Faer
153where
154 T: FaerVectorScalar + One,
155{
156 fn matvec<'a, La, Lx>(
157 &self,
158 a: &'a Slice<T, (D0, D1), La>,
159 x: &'a Slice<T, (D1,), Lx>,
160 ) -> impl MatVecBuilder<'a, T, La, Lx, D0, D1>
161 where
162 La: Layout,
163 Lx: Layout,
164 {
165 FaerMatVecBuilder {
166 alpha: T::one(),
167 a,
168 x,
169 par: faer::get_global_parallelism(),
170 }
171 }
172}
173
174impl<T, D: Dim> VecOps<T, D> for Faer
175where
176 T: FaerVectorScalar + Zero + Copy,
177{
178 type Real = <T as ComplexFloat>::Real;
179
180 fn add_to_scaled<Lx: Layout, Ly: Layout>(
181 &self,
182 _alpha: T,
183 _x: &Slice<T, (D,), Lx>,
184 _y: &mut Slice<T, (D,), Ly>,
185 ) {
186 unimplemented!();
187 }
188
189 fn dot<Lx: Layout, Ly: Layout>(&self, x: &Slice<T, (D,), Lx>, y: &Slice<T, (D,), Ly>) -> T {
190 assert_eq!(
191 x.shape().dim(0),
192 y.shape().dim(0),
193 "Vectors must have same length"
194 );
195 inner_prod(into_faer_row(x), Conj::No, into_faer_col(y), Conj::No)
196 }
197
198 fn dotc<Lx: Layout, Ly: Layout>(&self, x: &Slice<T, (D,), Lx>, y: &Slice<T, (D,), Ly>) -> T {
199 assert_eq!(
200 x.shape().dim(0),
201 y.shape().dim(0),
202 "Vectors must have same length"
203 );
204 inner_prod(into_faer_row(x), Conj::Yes, into_faer_col(y), Conj::No)
205 }
206
207 fn norm2<Lx: Layout>(&self, x: &Slice<T, (D,), Lx>) -> Self::Real {
208 T::from_faer_real(into_faer_col(x).norm_l2())
209 }
210
211 fn norm1<Lx: Layout>(&self, x: &Slice<T, (D,), Lx>) -> Self::Real {
212 T::from_faer_real(into_faer_col(x).norm_l1())
213 }
214
215 fn rot<Lx: Layout, Ly: Layout>(
216 &self,
217 _x: &mut Slice<T, (D,), Lx>,
218 _y: &mut Slice<T, (D,), Ly>,
219 _c: Self::Real,
220 _s: T,
221 ) {
222 unimplemented!();
223 }
224}
225
226impl<T, Dx, Dy> Outer<T, Dx, Dy> for Faer
227where
228 T: FaerVectorScalar + One,
229 Dx: Dim,
230 Dy: Dim,
231{
232 fn outer<'a, Lx, Ly>(
233 &self,
234 x: &'a Slice<T, (Dx,), Lx>,
235 y: &'a Slice<T, (Dy,), Ly>,
236 ) -> impl OuterBuilder<'a, T, Lx, Ly, Dx, Dy>
237 where
238 Lx: Layout,
239 Ly: Layout,
240 {
241 FaerOuterBuilder {
242 alpha: T::one(),
243 x,
244 y,
245 par: faer::get_global_parallelism(),
246 }
247 }
248}
249
250impl<'a, T, Lx, Ly, Dx, Dy> OuterBuilder<'a, T, Lx, Ly, Dx, Dy>
251 for FaerOuterBuilder<'a, T, Lx, Ly, Dx, Dy>
252where
253 Lx: Layout,
254 Ly: Layout,
255 Dx: Dim,
256 Dy: Dim,
257 T: FaerVectorScalar,
258{
259 fn scale(mut self, alpha: T) -> Self {
260 self.alpha *= alpha;
261 self
262 }
263
264 fn eval(self) -> Array<T, (Dx, Dy)> {
265 let m = self.x.shape().dim(0);
266 let n = self.y.shape().dim(0);
267 let shape = <(Dx, Dy) as Shape>::from_dims(&[m, n]);
268 let mut a = Array::<T, (Dx, Dy)>::from_elem(shape, T::zero());
269 self.write(&mut a);
270 a
271 }
272
273 fn write<La: Layout>(self, a: &mut Slice<T, (Dx, Dy), La>) {
274 let m = a.shape().dim(0);
275 let n = a.shape().dim(1);
276 assert_eq!(m, self.x.shape().dim(0), "Output rows must match x length");
277 assert_eq!(n, self.y.shape().dim(0), "Output cols must match y length");
278
279 matmul(
280 into_faer_mut(a),
281 Accum::Replace,
282 into_faer_col(self.x),
283 into_faer_row(self.y),
284 self.alpha,
285 self.par,
286 );
287 }
288
289 fn add_to<La: Layout>(self, a: &mut Slice<T, (Dx, Dy), La>) {
290 let m = a.shape().dim(0);
291 let n = a.shape().dim(1);
292 assert_eq!(m, self.x.shape().dim(0), "Output rows must match x length");
293 assert_eq!(n, self.y.shape().dim(0), "Output cols must match y length");
294
295 matmul(
296 into_faer_mut(a),
297 Accum::Add,
298 into_faer_col(self.x),
299 into_faer_row(self.y),
300 self.alpha,
301 self.par,
302 );
303 }
304}