faer_sparse/
lib.rs

1#![cfg_attr(not(feature = "std"), no_std)]
2#![allow(clippy::missing_safety_doc)]
3#![allow(clippy::type_complexity)]
4#![allow(clippy::too_many_arguments)]
5#![forbid(elided_lifetimes_in_paths)]
6#![allow(non_snake_case)]
7
8use bytemuck::Pod;
9use dyn_stack::{PodStack, SizeOverflow, StackReq};
10use faer_core::{
11    group_helpers::*,
12    permutation::PermutationRef,
13    sparse::{windows2, *},
14    Entity, Side,
15};
16
17const CHOLESKY_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
18const QR_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
19const LU_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
20
21#[derive(Copy, Clone, Debug)]
22pub struct SymbolicSupernodalParams<'a> {
23    pub relax: Option<&'a [(usize, f64)]>,
24}
25
26impl Default for SymbolicSupernodalParams<'_> {
27    #[inline]
28    fn default() -> Self {
29        Self {
30            relax: Some(&[(4, 1.0), (16, 0.8), (48, 0.1), (usize::MAX, 0.05)]),
31        }
32    }
33}
34
35#[derive(Copy, Clone, Debug, PartialEq)]
36pub struct SupernodalThreshold(pub f64);
37
38impl Default for SupernodalThreshold {
39    #[inline]
40    fn default() -> Self {
41        Self(1.0)
42    }
43}
44
45impl SupernodalThreshold {
46    pub const FORCE_SIMPLICIAL: Self = Self(f64::INFINITY);
47    pub const FORCE_SUPERNODAL: Self = Self(0.0);
48    pub const AUTO: Self = Self(1.0);
49}
50
51use faer_core::sparse::util::{
52    ghost_adjoint_symbolic, ghost_permute_hermitian_unsorted,
53    ghost_permute_hermitian_unsorted_symbolic, transpose,
54};
55
56pub use faer_core::{
57    permutation::{Index, SignedIndex},
58    FaerError,
59};
60
61#[allow(unused_macros)]
62macro_rules! shadow {
63    ($id: ident) => {
64        let $id = {
65            struct Shadowed;
66            impl ::core::ops::Drop for Shadowed {
67                fn drop(&mut self) {}
68            }
69            Shadowed
70        };
71        ::core::mem::drop($id);
72    };
73}
74
75macro_rules! impl_copy {
76    (< $($lt_param: lifetime),* >< $($ty_param: ident $(: $tt: tt)?),* > <$ty: ty>) => {
77        impl<$($lt_param,)* $($ty_param $(: $tt)?,)*> Copy for $ty {}
78        impl<$($lt_param,)* $($ty_param $(: $tt)?,)*> Clone for $ty {
79            #[inline(always)]
80            fn clone(&self) -> Self {
81                *self
82            }
83        }
84    };
85}
86
87#[inline]
88#[track_caller]
89fn try_zeroed<I: Pod>(n: usize) -> Result<alloc::vec::Vec<I>, FaerError> {
90    let mut v = alloc::vec::Vec::new();
91    v.try_reserve_exact(n).map_err(nomem)?;
92    unsafe {
93        core::ptr::write_bytes::<I>(v.as_mut_ptr(), 0u8, n);
94        v.set_len(n);
95    }
96    Ok(v)
97}
98
99#[inline]
100#[track_caller]
101fn try_collect<I: IntoIterator>(iter: I) -> Result<alloc::vec::Vec<I::Item>, FaerError> {
102    let iter = iter.into_iter();
103    let mut v = alloc::vec::Vec::new();
104    v.try_reserve_exact(iter.size_hint().0).map_err(nomem)?;
105    v.extend(iter);
106    Ok(v)
107}
108
109#[inline]
110fn nomem<T>(_: T) -> FaerError {
111    FaerError::OutOfMemory
112}
113
114fn make_raw_req<E: Entity>(size: usize) -> Result<StackReq, SizeOverflow> {
115    let req = Ok(StackReq::empty());
116    let additional = StackReq::try_new::<E::Unit>(size)?;
117    let (req, _) = E::faer_map_with_context(req, E::UNIT, &mut {
118        #[inline(always)]
119        |req, ()| {
120            let req = match req {
121                Ok(req) => req.try_and(additional),
122                _ => Err(SizeOverflow),
123            };
124            (req, ())
125        }
126    });
127    req
128}
129
130fn make_raw<E: Entity>(size: usize, stack: PodStack<'_>) -> (SliceGroupMut<'_, E>, PodStack<'_>) {
131    let (stack, array) = E::faer_map_with_context(stack, E::UNIT, &mut {
132        #[inline(always)]
133        |stack, ()| {
134            let (alloc, stack) = stack.make_raw::<E::Unit>(size);
135            (stack, alloc)
136        }
137    });
138    (SliceGroupMut::new(array), stack)
139}
140
141#[cfg(test)]
142macro_rules! monomorphize_test {
143    ($name: ident) => {
144        monomorphize_test!($name, u32);
145        monomorphize_test!($name, u64);
146    };
147
148    ($name: ident, $ty: ident) => {
149        paste::paste! {
150            #[test]
151            fn [<$name _ $ty>]() {
152                $name::<$ty>();
153            }
154        }
155    };
156}
157
158extern crate alloc;
159
160pub mod triangular_solve;
161
162pub mod amd;
163pub mod colamd;
164
165pub mod cholesky;
166
167#[doc(hidden)]
168pub mod lu;
169#[doc(hidden)]
170pub mod qr;
171
172#[doc(hidden)]
173pub mod superlu;
174
175mod ghost;
176
177mod mem;
178
179#[cfg(test)]
180pub(crate) mod qd {
181    // https://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf
182    // https://gitlab.com/hodge_star/mantis
183
184    use bytemuck::{Pod, Zeroable};
185    use faer_entity::*;
186    use pulp::{Scalar, Simd};
187
188    /// Value representing the implicit sum of two floating point terms, such that the absolute
189    /// value of the second term is less half a ULP of the first term.
190    #[derive(Copy, Clone, Debug, PartialEq, Eq, PartialOrd, Ord)]
191    #[repr(C)]
192    pub struct Double<T>(pub T, pub T);
193
194    unsafe impl<T: Zeroable> Zeroable for Double<T> {}
195    unsafe impl<T: Pod> Pod for Double<T> {}
196
197    impl<I: Iterator> Iterator for Double<I> {
198        type Item = Double<I::Item>;
199
200        #[inline(always)]
201        fn next(&mut self) -> Option<Self::Item> {
202            let x0 = self.0.next()?;
203            let x1 = self.1.next()?;
204            Some(Double(x0, x1))
205        }
206    }
207
208    #[inline(always)]
209    fn quick_two_sum<S: Simd>(simd: S, a: S::f64s, b: S::f64s) -> (S::f64s, S::f64s) {
210        let s = simd.f64s_add(a, b);
211        let err = simd.f64s_sub(b, simd.f64s_sub(s, a));
212        (s, err)
213    }
214
215    #[inline(always)]
216    fn two_sum<S: Simd>(simd: S, a: S::f64s, b: S::f64s) -> (S::f64s, S::f64s) {
217        let s = simd.f64s_add(a, b);
218        let bb = simd.f64s_sub(s, a);
219
220        // (a - (s - bb)) + (b - bb)
221        let err = simd.f64s_add(simd.f64s_sub(a, simd.f64s_sub(s, bb)), simd.f64s_sub(b, bb));
222        (s, err)
223    }
224
225    #[inline(always)]
226    #[allow(dead_code)]
227    fn quick_two_diff<S: Simd>(simd: S, a: S::f64s, b: S::f64s) -> (S::f64s, S::f64s) {
228        let s = simd.f64s_sub(a, b);
229        let err = simd.f64s_sub(simd.f64s_sub(a, s), b);
230        (s, err)
231    }
232
233    #[inline(always)]
234    fn two_diff<S: Simd>(simd: S, a: S::f64s, b: S::f64s) -> (S::f64s, S::f64s) {
235        let s = simd.f64s_sub(a, b);
236        let bb = simd.f64s_sub(s, a);
237
238        // (a - (s - bb)) - (b + bb)
239        let err = simd.f64s_sub(simd.f64s_sub(a, simd.f64s_sub(s, bb)), simd.f64s_add(b, bb));
240        (s, err)
241    }
242
243    #[inline(always)]
244    fn two_prod<S: Simd>(simd: S, a: S::f64s, b: S::f64s) -> (S::f64s, S::f64s) {
245        let p = simd.f64s_mul(a, b);
246        let err = simd.f64s_mul_add(a, b, simd.f64s_neg(p));
247
248        (p, err)
249    }
250
251    pub mod double {
252        use super::*;
253
254        #[inline(always)]
255        pub fn simd_add<S: Simd>(
256            simd: S,
257            a: Double<S::f64s>,
258            b: Double<S::f64s>,
259        ) -> Double<S::f64s> {
260            let (s, e) = two_sum(simd, a.0, b.0);
261            let e = simd.f64s_add(e, simd.f64s_add(a.1, b.1));
262            let (s, e) = quick_two_sum(simd, s, e);
263            Double(s, e)
264        }
265
266        #[inline(always)]
267        pub fn simd_sub<S: Simd>(
268            simd: S,
269            a: Double<S::f64s>,
270            b: Double<S::f64s>,
271        ) -> Double<S::f64s> {
272            let (s, e) = two_diff(simd, a.0, b.0);
273            let e = simd.f64s_add(e, a.1);
274            let e = simd.f64s_sub(e, b.1);
275            let (s, e) = quick_two_sum(simd, s, e);
276            Double(s, e)
277        }
278
279        #[inline(always)]
280        pub fn simd_neg<S: Simd>(simd: S, a: Double<S::f64s>) -> Double<S::f64s> {
281            Double(simd.f64s_neg(a.0), simd.f64s_neg(a.1))
282        }
283
284        #[inline(always)]
285        pub fn simd_mul<S: Simd>(
286            simd: S,
287            a: Double<S::f64s>,
288            b: Double<S::f64s>,
289        ) -> Double<S::f64s> {
290            let (p1, p2) = two_prod(simd, a.0, b.0);
291            let p2 = simd.f64s_add(
292                p2,
293                simd.f64s_add(simd.f64s_mul(a.0, b.1), simd.f64s_mul(a.1, b.0)),
294            );
295            let (p1, p2) = quick_two_sum(simd, p1, p2);
296            Double(p1, p2)
297        }
298
299        #[inline(always)]
300        fn simd_mul_f64<S: Simd>(simd: S, a: Double<S::f64s>, b: S::f64s) -> Double<S::f64s> {
301            let (p1, p2) = two_prod(simd, a.0, b);
302            let p2 = simd.f64s_add(p2, simd.f64s_mul(a.1, b));
303            let (p1, p2) = quick_two_sum(simd, p1, p2);
304            Double(p1, p2)
305        }
306
307        pub fn simd_select<S: Simd>(
308            simd: S,
309            mask: S::m64s,
310            if_true: Double<S::f64s>,
311            if_false: Double<S::f64s>,
312        ) -> Double<S::f64s> {
313            Double(
314                simd.m64s_select_f64s(mask, if_true.0, if_false.0),
315                simd.m64s_select_f64s(mask, if_true.1, if_false.1),
316            )
317        }
318
319        #[inline]
320        pub fn simd_div<S: Simd>(
321            simd: S,
322            a: Double<S::f64s>,
323            b: Double<S::f64s>,
324        ) -> Double<S::f64s> {
325            simd.vectorize(
326                #[inline(always)]
327                || {
328                    let pos_zero = simd.f64s_splat(0.0);
329                    let pos_infty = simd.f64s_splat(f64::INFINITY);
330                    let sign_bit = simd.f64s_splat(-0.0);
331
332                    let a_sign = simd.f64s_and(a.0, sign_bit);
333                    let b_sign = simd.f64s_and(b.0, sign_bit);
334
335                    let combined_sign = simd.f64s_xor(a_sign, b_sign);
336
337                    let a_is_zero = simd_eq(simd, a, Double(pos_zero, pos_zero));
338                    let b_is_zero = simd_eq(simd, b, Double(pos_zero, pos_zero));
339                    let a_is_infty = simd_eq(
340                        simd,
341                        Double(simd.f64s_abs(a.0), simd.f64s_abs(a.1)),
342                        Double(pos_infty, pos_infty),
343                    );
344                    let b_is_infty = simd_eq(
345                        simd,
346                        Double(simd.f64s_abs(b.0), simd.f64s_abs(b.1)),
347                        Double(pos_infty, pos_infty),
348                    );
349
350                    let q1 = simd.f64s_div(a.0, b.0);
351                    let r = simd_mul_f64(simd, b, q1);
352
353                    let (s1, s2) = two_diff(simd, a.0, r.0);
354                    let s2 = simd.f64s_sub(s2, r.1);
355                    let s2 = simd.f64s_add(s2, a.1);
356
357                    let q2 = simd.f64s_div(simd.f64s_add(s1, s2), b.0);
358                    let (q0, q1) = quick_two_sum(simd, q1, q2);
359
360                    simd_select(
361                        simd,
362                        simd.m64s_and(b_is_zero, simd.m64s_not(a_is_zero)),
363                        Double(
364                            simd.f64s_or(combined_sign, pos_infty),
365                            simd.f64s_or(combined_sign, pos_infty),
366                        ),
367                        simd_select(
368                            simd,
369                            simd.m64s_and(b_is_infty, simd.m64s_not(a_is_infty)),
370                            Double(
371                                simd.f64s_or(combined_sign, pos_zero),
372                                simd.f64s_or(combined_sign, pos_zero),
373                            ),
374                            Double(q0, q1),
375                        ),
376                    )
377                },
378            )
379        }
380
381        #[inline(always)]
382        pub fn simd_abs<S: Simd>(simd: S, a: Double<S::f64s>) -> Double<S::f64s> {
383            let is_negative = simd.f64s_less_than(a.0, simd.f64s_splat(0.0));
384            Double(
385                simd.m64s_select_f64s(is_negative, simd.f64s_neg(a.0), a.0),
386                simd.m64s_select_f64s(is_negative, simd.f64s_neg(a.1), a.1),
387            )
388        }
389
390        #[inline(always)]
391        pub fn simd_less_than<S: Simd>(simd: S, a: Double<S::f64s>, b: Double<S::f64s>) -> S::m64s {
392            let lt0 = simd.f64s_less_than(a.0, b.0);
393            let eq0 = simd.f64s_equal(a.0, b.0);
394            let lt1 = simd.f64s_less_than(a.1, b.1);
395            simd.m64s_or(lt0, simd.m64s_and(eq0, lt1))
396        }
397
398        #[inline(always)]
399        pub fn simd_less_than_or_equal<S: Simd>(
400            simd: S,
401            a: Double<S::f64s>,
402            b: Double<S::f64s>,
403        ) -> S::m64s {
404            let lt0 = simd.f64s_less_than(a.0, b.0);
405            let eq0 = simd.f64s_equal(a.0, b.0);
406            let lt1 = simd.f64s_less_than_or_equal(a.1, b.1);
407            simd.m64s_or(lt0, simd.m64s_and(eq0, lt1))
408        }
409
410        #[inline(always)]
411        pub fn simd_greater_than<S: Simd>(
412            simd: S,
413            a: Double<S::f64s>,
414            b: Double<S::f64s>,
415        ) -> S::m64s {
416            let lt0 = simd.f64s_greater_than(a.0, b.0);
417            let eq0 = simd.f64s_equal(a.0, b.0);
418            let lt1 = simd.f64s_greater_than(a.1, b.1);
419            simd.m64s_or(lt0, simd.m64s_and(eq0, lt1))
420        }
421
422        #[inline(always)]
423        pub fn simd_greater_than_or_equal<S: Simd>(
424            simd: S,
425            a: Double<S::f64s>,
426            b: Double<S::f64s>,
427        ) -> S::m64s {
428            let lt0 = simd.f64s_greater_than(a.0, b.0);
429            let eq0 = simd.f64s_equal(a.0, b.0);
430            let lt1 = simd.f64s_greater_than_or_equal(a.1, b.1);
431            simd.m64s_or(lt0, simd.m64s_and(eq0, lt1))
432        }
433
434        #[inline(always)]
435        pub fn simd_eq<S: Simd>(simd: S, a: Double<S::f64s>, b: Double<S::f64s>) -> S::m64s {
436            let eq0 = simd.f64s_equal(a.0, b.0);
437            let eq1 = simd.f64s_equal(a.1, b.1);
438            simd.m64s_and(eq0, eq1)
439        }
440    }
441
442    impl core::ops::Add for Double<f64> {
443        type Output = Self;
444
445        #[inline(always)]
446        fn add(self, rhs: Self) -> Self::Output {
447            double::simd_add(Scalar::new(), self, rhs)
448        }
449    }
450
451    impl core::ops::Sub for Double<f64> {
452        type Output = Self;
453
454        #[inline(always)]
455        fn sub(self, rhs: Self) -> Self::Output {
456            double::simd_sub(Scalar::new(), self, rhs)
457        }
458    }
459
460    impl core::ops::Mul for Double<f64> {
461        type Output = Self;
462
463        #[inline(always)]
464        fn mul(self, rhs: Self) -> Self::Output {
465            double::simd_mul(Scalar::new(), self, rhs)
466        }
467    }
468
469    impl core::ops::Div for Double<f64> {
470        type Output = Self;
471
472        #[inline(always)]
473        fn div(self, rhs: Self) -> Self::Output {
474            double::simd_div(Scalar::new(), self, rhs)
475        }
476    }
477
478    impl core::ops::AddAssign for Double<f64> {
479        #[inline(always)]
480        fn add_assign(&mut self, rhs: Self) {
481            *self = *self + rhs;
482        }
483    }
484
485    impl core::ops::SubAssign for Double<f64> {
486        #[inline(always)]
487        fn sub_assign(&mut self, rhs: Self) {
488            *self = *self - rhs;
489        }
490    }
491
492    impl core::ops::MulAssign for Double<f64> {
493        #[inline(always)]
494        fn mul_assign(&mut self, rhs: Self) {
495            *self = *self * rhs;
496        }
497    }
498
499    impl core::ops::DivAssign for Double<f64> {
500        #[inline(always)]
501        fn div_assign(&mut self, rhs: Self) {
502            *self = *self / rhs;
503        }
504    }
505
506    impl core::ops::Neg for Double<f64> {
507        type Output = Self;
508
509        #[inline(always)]
510        fn neg(self) -> Self::Output {
511            Self(-self.0, -self.1)
512        }
513    }
514
515    impl Double<f64> {
516        /// 2.0^{-100}
517        pub const EPSILON: Self = Self(7.888609052210118e-31, 0.0);
518        /// 2.0^{-970}: precision below this value begins to degrade.
519        pub const MIN_POSITIVE: Self = Self(1.0020841800044864e-292, 0.0);
520
521        pub const ZERO: Self = Self(0.0, 0.0);
522        pub const NAN: Self = Self(f64::NAN, f64::NAN);
523        pub const INFINITY: Self = Self(f64::INFINITY, f64::INFINITY);
524
525        #[inline(always)]
526        pub fn abs(self) -> Self {
527            double::simd_abs(Scalar::new(), self)
528        }
529
530        #[inline(always)]
531        pub fn recip(self) -> Self {
532            double::simd_div(Scalar::new(), Self(1.0, 0.0), self)
533        }
534
535        #[inline]
536        pub fn sqrt(self) -> Self {
537            if self == Self::ZERO {
538                Self::ZERO
539            } else if self < Self::ZERO {
540                Self::NAN
541            } else if self == Self::INFINITY {
542                Self::INFINITY
543            } else {
544                let a = self;
545                let x = a.0.sqrt().recip();
546                let ax = Self(a.0 * x, 0.0);
547
548                ax + (a - ax * ax) * Double(x * 0.5, 0.0)
549            }
550        }
551    }
552
553    pub struct DoubleGroup {
554        __private: (),
555    }
556
557    impl ForType for DoubleGroup {
558        type FaerOf<T> = Double<T>;
559    }
560    impl ForCopyType for DoubleGroup {
561        type FaerOfCopy<T: Copy> = Double<T>;
562    }
563    impl ForDebugType for DoubleGroup {
564        type FaerOfDebug<T: core::fmt::Debug> = Double<T>;
565    }
566
567    mod faer_impl {
568        use super::*;
569        use faer_core::{ComplexField, Conjugate, Entity, RealField};
570
571        unsafe impl Entity for Double<f64> {
572            type Unit = f64;
573            type Index = u64;
574
575            type SimdUnit<S: Simd> = S::f64s;
576            type SimdMask<S: Simd> = S::m64s;
577            type SimdIndex<S: Simd> = S::u64s;
578
579            type Group = DoubleGroup;
580            type Iter<I: Iterator> = Double<I>;
581
582            type PrefixUnit<'a, S: Simd> = pulp::Prefix<'a, f64, S, S::m64s>;
583            type SuffixUnit<'a, S: Simd> = pulp::Suffix<'a, f64, S, S::m64s>;
584            type PrefixMutUnit<'a, S: Simd> = pulp::PrefixMut<'a, f64, S, S::m64s>;
585            type SuffixMutUnit<'a, S: Simd> = pulp::SuffixMut<'a, f64, S, S::m64s>;
586
587            const N_COMPONENTS: usize = 2;
588            const UNIT: GroupCopyFor<Self, ()> = Double((), ());
589
590            #[inline(always)]
591            fn faer_first<T>(group: GroupFor<Self, T>) -> T {
592                group.0
593            }
594
595            #[inline(always)]
596            fn faer_from_units(group: GroupFor<Self, Self::Unit>) -> Self {
597                group
598            }
599
600            #[inline(always)]
601            fn faer_into_units(self) -> GroupFor<Self, Self::Unit> {
602                self
603            }
604
605            #[inline(always)]
606            fn faer_as_ref<T>(group: &GroupFor<Self, T>) -> GroupFor<Self, &T> {
607                Double(&group.0, &group.1)
608            }
609
610            #[inline(always)]
611            fn faer_as_mut<T>(group: &mut GroupFor<Self, T>) -> GroupFor<Self, &mut T> {
612                Double(&mut group.0, &mut group.1)
613            }
614
615            #[inline(always)]
616            fn faer_as_ptr<T>(group: *mut GroupFor<Self, T>) -> GroupFor<Self, *mut T> {
617                unsafe {
618                    Double(
619                        core::ptr::addr_of_mut!((*group).0),
620                        core::ptr::addr_of_mut!((*group).1),
621                    )
622                }
623            }
624
625            #[inline(always)]
626            fn faer_map_impl<T, U>(
627                group: GroupFor<Self, T>,
628                f: &mut impl FnMut(T) -> U,
629            ) -> GroupFor<Self, U> {
630                Double((*f)(group.0), (*f)(group.1))
631            }
632
633            #[inline(always)]
634            fn faer_zip<T, U>(
635                first: GroupFor<Self, T>,
636                second: GroupFor<Self, U>,
637            ) -> GroupFor<Self, (T, U)> {
638                Double((first.0, second.0), (first.1, second.1))
639            }
640
641            #[inline(always)]
642            fn faer_unzip<T, U>(
643                zipped: GroupFor<Self, (T, U)>,
644            ) -> (GroupFor<Self, T>, GroupFor<Self, U>) {
645                (
646                    Double(zipped.0 .0, zipped.1 .0),
647                    Double(zipped.0 .1, zipped.1 .1),
648                )
649            }
650
651            #[inline(always)]
652            fn faer_map_with_context<Ctx, T, U>(
653                ctx: Ctx,
654                group: GroupFor<Self, T>,
655                f: &mut impl FnMut(Ctx, T) -> (Ctx, U),
656            ) -> (Ctx, GroupFor<Self, U>) {
657                let (ctx, x0) = (*f)(ctx, group.0);
658                let (ctx, x1) = (*f)(ctx, group.1);
659                (ctx, Double(x0, x1))
660            }
661
662            #[inline(always)]
663            fn faer_into_iter<I: IntoIterator>(iter: GroupFor<Self, I>) -> Self::Iter<I::IntoIter> {
664                Double(iter.0.into_iter(), iter.1.into_iter())
665            }
666        }
667
668        unsafe impl Conjugate for Double<f64> {
669            type Conj = Double<f64>;
670            type Canonical = Double<f64>;
671            #[inline(always)]
672            fn canonicalize(self) -> Self::Canonical {
673                self
674            }
675        }
676
677        impl RealField for Double<f64> {
678            #[inline(always)]
679            fn faer_epsilon() -> Option<Self> {
680                Some(Self::EPSILON)
681            }
682            #[inline(always)]
683            fn faer_zero_threshold() -> Option<Self> {
684                Some(Self::MIN_POSITIVE)
685            }
686
687            #[inline(always)]
688            fn faer_div(self, rhs: Self) -> Self {
689                self / rhs
690            }
691
692            #[inline(always)]
693            fn faer_usize_to_index(a: usize) -> Self::Index {
694                a as _
695            }
696
697            #[inline(always)]
698            fn faer_index_to_usize(a: Self::Index) -> usize {
699                a as _
700            }
701
702            #[inline(always)]
703            fn faer_max_index() -> Self::Index {
704                Self::Index::MAX
705            }
706
707            #[inline(always)]
708            fn faer_simd_less_than<S: Simd>(
709                simd: S,
710                a: SimdGroupFor<Self, S>,
711                b: SimdGroupFor<Self, S>,
712            ) -> Self::SimdMask<S> {
713                double::simd_less_than(simd, a, b)
714            }
715
716            #[inline(always)]
717            fn faer_simd_less_than_or_equal<S: Simd>(
718                simd: S,
719                a: SimdGroupFor<Self, S>,
720                b: SimdGroupFor<Self, S>,
721            ) -> Self::SimdMask<S> {
722                double::simd_less_than_or_equal(simd, a, b)
723            }
724
725            #[inline(always)]
726            fn faer_simd_greater_than<S: Simd>(
727                simd: S,
728                a: SimdGroupFor<Self, S>,
729                b: SimdGroupFor<Self, S>,
730            ) -> Self::SimdMask<S> {
731                double::simd_greater_than(simd, a, b)
732            }
733
734            #[inline(always)]
735            fn faer_simd_greater_than_or_equal<S: Simd>(
736                simd: S,
737                a: SimdGroupFor<Self, S>,
738                b: SimdGroupFor<Self, S>,
739            ) -> Self::SimdMask<S> {
740                double::simd_greater_than_or_equal(simd, a, b)
741            }
742
743            #[inline(always)]
744            fn faer_simd_select<S: Simd>(
745                simd: S,
746                mask: Self::SimdMask<S>,
747                if_true: SimdGroupFor<Self, S>,
748                if_false: SimdGroupFor<Self, S>,
749            ) -> SimdGroupFor<Self, S> {
750                double::simd_select(simd, mask, if_true, if_false)
751            }
752
753            #[inline(always)]
754            fn faer_simd_index_select<S: Simd>(
755                simd: S,
756                mask: Self::SimdMask<S>,
757                if_true: Self::SimdIndex<S>,
758                if_false: Self::SimdIndex<S>,
759            ) -> Self::SimdIndex<S> {
760                simd.m64s_select_u64s(mask, if_true, if_false)
761            }
762
763            #[inline(always)]
764            fn faer_simd_index_seq<S: Simd>(simd: S) -> Self::SimdIndex<S> {
765                let _ = simd;
766                pulp::cast_lossy([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15_u64])
767            }
768
769            #[inline(always)]
770            fn faer_simd_index_splat<S: Simd>(simd: S, value: Self::Index) -> Self::SimdIndex<S> {
771                simd.u64s_splat(value)
772            }
773
774            #[inline(always)]
775            fn faer_simd_index_add<S: Simd>(
776                simd: S,
777                a: Self::SimdIndex<S>,
778                b: Self::SimdIndex<S>,
779            ) -> Self::SimdIndex<S> {
780                simd.u64s_add(a, b)
781            }
782
783            #[inline(always)]
784            fn faer_min_positive() -> Self {
785                Self::MIN_POSITIVE
786            }
787
788            #[inline(always)]
789            fn faer_min_positive_inv() -> Self {
790                Self::MIN_POSITIVE.recip()
791            }
792
793            #[inline(always)]
794            fn faer_min_positive_sqrt() -> Self {
795                Self::MIN_POSITIVE.sqrt()
796            }
797
798            #[inline(always)]
799            fn faer_min_positive_sqrt_inv() -> Self {
800                Self::MIN_POSITIVE.sqrt().recip()
801            }
802
803            #[inline(always)]
804            fn faer_simd_index_rotate_left<S: Simd>(
805                simd: S,
806                values: SimdIndexFor<Self, S>,
807                amount: usize,
808            ) -> SimdIndexFor<Self, S> {
809                simd.u64s_rotate_left(values, amount)
810            }
811
812            #[inline(always)]
813            fn faer_simd_abs<S: Simd>(
814                simd: S,
815                values: SimdGroupFor<Self, S>,
816            ) -> SimdGroupFor<Self, S> {
817                double::simd_abs(simd, values)
818            }
819        }
820
821        impl ComplexField for Double<f64> {
822            type Real = Double<f64>;
823            type Simd = pulp::Arch;
824            type ScalarSimd = pulp::Arch;
825            type PortableSimd = pulp::Arch;
826
827            #[inline(always)]
828            fn faer_sqrt(self) -> Self {
829                self.sqrt()
830            }
831
832            #[inline(always)]
833            fn faer_from_f64(value: f64) -> Self {
834                Self(value, 0.0)
835            }
836
837            #[inline(always)]
838            fn faer_add(self, rhs: Self) -> Self {
839                self + rhs
840            }
841
842            #[inline(always)]
843            fn faer_sub(self, rhs: Self) -> Self {
844                self - rhs
845            }
846
847            #[inline(always)]
848            fn faer_mul(self, rhs: Self) -> Self {
849                self * rhs
850            }
851
852            #[inline(always)]
853            fn faer_neg(self) -> Self {
854                -self
855            }
856
857            #[inline(always)]
858            fn faer_inv(self) -> Self {
859                self.recip()
860            }
861
862            #[inline(always)]
863            fn faer_conj(self) -> Self {
864                self
865            }
866
867            #[inline(always)]
868            fn faer_scale_real(self, rhs: Self::Real) -> Self {
869                self * rhs
870            }
871
872            #[inline(always)]
873            fn faer_scale_power_of_two(self, rhs: Self::Real) -> Self {
874                Self(self.0 * rhs.0, self.1 * rhs.0)
875            }
876
877            #[inline(always)]
878            fn faer_score(self) -> Self::Real {
879                self.abs()
880            }
881
882            #[inline(always)]
883            fn faer_abs(self) -> Self::Real {
884                self.abs()
885            }
886
887            #[inline(always)]
888            fn faer_abs2(self) -> Self::Real {
889                self * self
890            }
891
892            #[inline(always)]
893            fn faer_nan() -> Self {
894                Self::NAN
895            }
896
897            #[inline(always)]
898            fn faer_from_real(real: Self::Real) -> Self {
899                real
900            }
901
902            #[inline(always)]
903            fn faer_real(self) -> Self::Real {
904                self
905            }
906
907            #[inline(always)]
908            fn faer_imag(self) -> Self::Real {
909                Self::ZERO
910            }
911
912            #[inline(always)]
913            fn faer_zero() -> Self {
914                Self::ZERO
915            }
916
917            #[inline(always)]
918            fn faer_one() -> Self {
919                Self(1.0, 0.0)
920            }
921
922            #[inline(always)]
923            fn faer_slice_as_simd<S: Simd>(
924                slice: &[Self::Unit],
925            ) -> (&[Self::SimdUnit<S>], &[Self::Unit]) {
926                S::f64s_as_simd(slice)
927            }
928
929            #[inline(always)]
930            fn faer_slice_as_simd_mut<S: Simd>(
931                slice: &mut [Self::Unit],
932            ) -> (&mut [Self::SimdUnit<S>], &mut [Self::Unit]) {
933                S::f64s_as_mut_simd(slice)
934            }
935
936            #[inline(always)]
937            fn faer_partial_load_unit<S: Simd>(simd: S, slice: &[Self::Unit]) -> Self::SimdUnit<S> {
938                simd.f64s_partial_load(slice)
939            }
940
941            #[inline(always)]
942            fn faer_partial_store_unit<S: Simd>(
943                simd: S,
944                slice: &mut [Self::Unit],
945                values: Self::SimdUnit<S>,
946            ) {
947                simd.f64s_partial_store(slice, values)
948            }
949
950            #[inline(always)]
951            fn faer_partial_load_last_unit<S: Simd>(
952                simd: S,
953                slice: &[Self::Unit],
954            ) -> Self::SimdUnit<S> {
955                simd.f64s_partial_load_last(slice)
956            }
957
958            #[inline(always)]
959            fn faer_partial_store_last_unit<S: Simd>(
960                simd: S,
961                slice: &mut [Self::Unit],
962                values: Self::SimdUnit<S>,
963            ) {
964                simd.f64s_partial_store_last(slice, values)
965            }
966
967            #[inline(always)]
968            fn faer_simd_splat_unit<S: Simd>(simd: S, unit: Self::Unit) -> Self::SimdUnit<S> {
969                simd.f64s_splat(unit)
970            }
971
972            #[inline(always)]
973            fn faer_simd_neg<S: Simd>(
974                simd: S,
975                values: SimdGroupFor<Self, S>,
976            ) -> SimdGroupFor<Self, S> {
977                double::simd_neg(simd, values)
978            }
979
980            #[inline(always)]
981            fn faer_simd_conj<S: Simd>(
982                simd: S,
983                values: SimdGroupFor<Self, S>,
984            ) -> SimdGroupFor<Self, S> {
985                let _ = simd;
986                values
987            }
988
989            #[inline(always)]
990            fn faer_simd_add<S: Simd>(
991                simd: S,
992                lhs: SimdGroupFor<Self, S>,
993                rhs: SimdGroupFor<Self, S>,
994            ) -> SimdGroupFor<Self, S> {
995                double::simd_add(simd, lhs, rhs)
996            }
997
998            #[inline(always)]
999            fn faer_simd_sub<S: Simd>(
1000                simd: S,
1001                lhs: SimdGroupFor<Self, S>,
1002                rhs: SimdGroupFor<Self, S>,
1003            ) -> SimdGroupFor<Self, S> {
1004                double::simd_sub(simd, lhs, rhs)
1005            }
1006
1007            #[inline(always)]
1008            fn faer_simd_mul<S: Simd>(
1009                simd: S,
1010                lhs: SimdGroupFor<Self, S>,
1011                rhs: SimdGroupFor<Self, S>,
1012            ) -> SimdGroupFor<Self, S> {
1013                double::simd_mul(simd, lhs, rhs)
1014            }
1015
1016            #[inline(always)]
1017            fn faer_simd_scale_real<S: Simd>(
1018                simd: S,
1019                lhs: SimdGroupFor<Self, S>,
1020                rhs: SimdGroupFor<Self, S>,
1021            ) -> SimdGroupFor<Self, S> {
1022                double::simd_mul(simd, lhs, rhs)
1023            }
1024
1025            #[inline(always)]
1026            fn faer_simd_conj_mul<S: Simd>(
1027                simd: S,
1028                lhs: SimdGroupFor<Self, S>,
1029                rhs: SimdGroupFor<Self, S>,
1030            ) -> SimdGroupFor<Self, S> {
1031                double::simd_mul(simd, lhs, rhs)
1032            }
1033
1034            #[inline(always)]
1035            fn faer_simd_mul_adde<S: Simd>(
1036                simd: S,
1037                lhs: SimdGroupFor<Self, S>,
1038                rhs: SimdGroupFor<Self, S>,
1039                acc: SimdGroupFor<Self, S>,
1040            ) -> SimdGroupFor<Self, S> {
1041                double::simd_add(simd, acc, double::simd_mul(simd, lhs, rhs))
1042            }
1043
1044            #[inline(always)]
1045            fn faer_simd_conj_mul_adde<S: Simd>(
1046                simd: S,
1047                lhs: SimdGroupFor<Self, S>,
1048                rhs: SimdGroupFor<Self, S>,
1049                acc: SimdGroupFor<Self, S>,
1050            ) -> SimdGroupFor<Self, S> {
1051                double::simd_add(simd, acc, double::simd_mul(simd, lhs, rhs))
1052            }
1053
1054            #[inline(always)]
1055            fn faer_simd_score<S: Simd>(
1056                simd: S,
1057                values: SimdGroupFor<Self, S>,
1058            ) -> SimdGroupFor<Self::Real, S> {
1059                double::simd_abs(simd, values)
1060            }
1061
1062            #[inline(always)]
1063            fn faer_simd_abs2_adde<S: Simd>(
1064                simd: S,
1065                values: SimdGroupFor<Self, S>,
1066                acc: SimdGroupFor<Self::Real, S>,
1067            ) -> SimdGroupFor<Self::Real, S> {
1068                Self::faer_simd_add(simd, acc, Self::faer_simd_mul(simd, values, values))
1069            }
1070
1071            #[inline(always)]
1072            fn faer_simd_abs2<S: Simd>(
1073                simd: S,
1074                values: SimdGroupFor<Self, S>,
1075            ) -> SimdGroupFor<Self::Real, S> {
1076                Self::faer_simd_mul(simd, values, values)
1077            }
1078
1079            #[inline(always)]
1080            fn faer_simd_scalar_mul<S: Simd>(simd: S, lhs: Self, rhs: Self) -> Self {
1081                let _ = simd;
1082                lhs * rhs
1083            }
1084
1085            #[inline(always)]
1086            fn faer_simd_scalar_conj_mul<S: Simd>(simd: S, lhs: Self, rhs: Self) -> Self {
1087                let _ = simd;
1088                lhs * rhs
1089            }
1090
1091            #[inline(always)]
1092            fn faer_simd_scalar_mul_adde<S: Simd>(
1093                simd: S,
1094                lhs: Self,
1095                rhs: Self,
1096                acc: Self,
1097            ) -> Self {
1098                let _ = simd;
1099                lhs * rhs + acc
1100            }
1101
1102            #[inline(always)]
1103            fn faer_simd_scalar_conj_mul_adde<S: Simd>(
1104                simd: S,
1105                lhs: Self,
1106                rhs: Self,
1107                acc: Self,
1108            ) -> Self {
1109                let _ = simd;
1110                lhs * rhs + acc
1111            }
1112
1113            #[inline(always)]
1114            fn faer_slice_as_aligned_simd<S: Simd>(
1115                simd: S,
1116                slice: &[UnitFor<Self>],
1117                offset: pulp::Offset<SimdMaskFor<Self, S>>,
1118            ) -> (
1119                pulp::Prefix<'_, UnitFor<Self>, S, SimdMaskFor<Self, S>>,
1120                &[SimdUnitFor<Self, S>],
1121                pulp::Suffix<'_, UnitFor<Self>, S, SimdMaskFor<Self, S>>,
1122            ) {
1123                simd.f64s_as_aligned_simd(slice, offset)
1124            }
1125            #[inline(always)]
1126            fn faer_slice_as_aligned_simd_mut<S: Simd>(
1127                simd: S,
1128                slice: &mut [UnitFor<Self>],
1129                offset: pulp::Offset<SimdMaskFor<Self, S>>,
1130            ) -> (
1131                pulp::PrefixMut<'_, UnitFor<Self>, S, SimdMaskFor<Self, S>>,
1132                &mut [SimdUnitFor<Self, S>],
1133                pulp::SuffixMut<'_, UnitFor<Self>, S, SimdMaskFor<Self, S>>,
1134            ) {
1135                simd.f64s_as_aligned_mut_simd(slice, offset)
1136            }
1137
1138            #[inline(always)]
1139            fn faer_simd_rotate_left<S: Simd>(
1140                simd: S,
1141                values: SimdGroupFor<Self, S>,
1142                amount: usize,
1143            ) -> SimdGroupFor<Self, S> {
1144                Double(
1145                    simd.f64s_rotate_left(values.0, amount),
1146                    simd.f64s_rotate_left(values.1, amount),
1147                )
1148            }
1149
1150            #[inline(always)]
1151            fn faer_align_offset<S: Simd>(
1152                simd: S,
1153                ptr: *const UnitFor<Self>,
1154                len: usize,
1155            ) -> pulp::Offset<SimdMaskFor<Self, S>> {
1156                simd.f64s_align_offset(ptr, len)
1157            }
1158        }
1159    }
1160}
faer_sparse/lib.rs

faer_sparse/
lib.rs
