Commit 584d94e5 authored by Michael Tesch's avatar Michael Tesch
Browse files

Merge branch 'feature/neon-fma-optimization' into 'master' 

Optimize NEON packet math with AArch64 FMA, fix benchmark naming

See merge request !52
parents 6e56e992 5414a1ed

CMakeLists.txt

+1 −1
Original line number Diff line number Diff line
@@ -3,7 +3,7 @@ 
#

cmake_minimum_required (VERSION 3.14)

project (cppduals

  VERSION 0.8.2

  VERSION 0.8.3

  LANGUAGES C CXX

  )

include (GNUInstallDirs)

README.md

+75 −0
Original line number Diff line number Diff line
@@ -227,6 +227,73 @@ is slower than K separate `2N` tiles for K=3 because `((K+1)N)^3` 
grows cubically while K separate `(2N)^3` calls grow linearly in K.

It also uses 4x more memory than tiling for K=3.



### Eigen NEON vectorization



Eigen matrix operations benefit from SIMD vectorization of `dual<T>`

packet arithmetic.  The `duals/arch/NEON/` headers provide ARM NEON

specializations: `Packet2df` packs two `dual<float>` into a 128-bit

register, and `Packet1dd` packs one `dual<double>`.  The key

optimization is a fused multiply-add (`pmadd`) using AArch64 FMA

instructions, which is the inner loop of Eigen's GEMM kernel.



Benchmarks comparing NEON-vectorized Eigen vs scalar (no Eigen

vectorization, `-DEIGEN_DONT_VECTORIZE`) on Apple M4 Pro.  Times in

nanoseconds; **bold** speedups indicate NEON is faster.



**Single precision — Eigen matrix operations:**



| Type | Operation | Scalar | NEON | Speedup |

|------|-----------|-------:|-----:|--------:|

| `float` | MatVec (N=4096) | 25,987,000 | 4,739,000 | **5.5x** |

| `float` | MatMat (N=1024) | 119,841,000 | 24,146,000 | **5.0x** |

| `complex<float>` | VecVec Mul (N=2048) | 1,010 | 460 | **2.2x** |

| `complex<float>` | MatVec (N=2048) | 6,056,000 | 4,017,000 | **1.5x** |

| `complex<float>` | MatMat (N=512) | 39,198,000 | 17,593,000 | **2.2x** |

| `dual<float>` | VecVec Mul (N=2048) | 364 | 394 | 0.92x |

| `dual<float>` | MatVec (N=2048) | 6,462,000 | 3,421,000 | **1.9x** |

| `dual<float>` | MatMat (N=512) | 24,153,000 | 20,508,000 | **1.18x** |

| `dual<float>` | MatDiv (N=512) | 36,592,000 | 32,778,000 | **1.12x** |

| `complex<dual<float>>` | VecVec Mul (N=1024) | 649 | 789 | 0.82x |

| `complex<dual<float>>` | MatVec (N=1024) | 2,677,000 | 2,527,000 | **1.06x** |

| `complex<dual<float>>` | MatMat (N=256) | 14,162,000 | 18,073,000 | 0.78x |



**Double precision — Eigen matrix operations:**



| Type | Operation | Scalar | NEON | Speedup |

|------|-----------|-------:|-----:|--------:|

| `double` | MatVec (N=4096) | 25,430,000 | 9,531,000 | **2.7x** |

| `double` | MatMat (N=1024) | 162,100,000 | 49,418,000 | **3.3x** |

| `complex<double>` | VecVec Mul (N=2048) | 1,107 | 697 | **1.6x** |

| `complex<double>` | MatVec (N=2048) | 7,618,000 | 7,567,000 | 1.0x |

| `complex<double>` | MatMat (N=512) | 48,192,000 | 35,369,000 | **1.4x** |

| `dual<double>` | VecVec Mul (N=2048) | 799 | 591 | **1.35x** |

| `dual<double>` | MatVec (N=2048) | 7,549,000 | 6,660,000 | **1.13x** |

| `dual<double>` | MatMat (N=512) | 32,114,000 | 32,354,000 | 1.0x |

| `dual<double>` | MatDiv (N=512) | 51,311,000 | 48,674,000 | **1.05x** |

| `complex<dual<double>>` | VecVec Mul (N=1024) | 1,335 | 1,334 | 1.0x |

| `complex<dual<double>>` | MatVec (N=1024) | 3,390,000 | 3,185,000 | 1.06x |

| `complex<dual<double>>` | MatMat (N=256) | 21,846,000 | 21,853,000 | 1.0x |

| `complex<dual<double>>` | MatDiv (N=256) | 32,096,000 | 32,163,000 | 1.0x |



For matrix operations (MatMat, MatVec, MatDiv), NEON vectorization

provides a clear speedup for all scalar and dual types.  The dual

speedup is smaller than the scalar baseline because `dual<float>`

packs only 2 elements per 128-bit register (vs 4 for float) and

dual multiplication requires more arithmetic (FMA + bit-select vs

simple multiply); `dual<double>` fits only 1 per register, so gains

come entirely from fused multiply-add instructions.

`complex<dual<float>>` fits a single value per 128-bit register

(packet size 1), limiting GEMM speedup — MatVec still benefits but

MatMat sees overhead from the complex dual packet arithmetic.

`complex<dual<double>>` requires 4 doubles (256 bits) which exceeds

NEON's 128-bit registers, so there is no NEON packet type for it —

Eigen falls back to scalar, giving 1.0x across the board.  AVX

(256-bit) can pack one `complex<dual<double>>` per register and does

provide vectorized packet math for this type.

Element-wise operations on small vectors sometimes show overhead from

Eigen's packet dispatch; the compiler's auto-vectorization of the

scalar path is competitive for trivial loops.



Contributions

=============
@@ -258,6 +325,14 @@ also licensed under MPL-2.0. 
ChangeLog

=========



v0.8.3

------



- optimize ARM NEON packet math for AArch64: use FMA intrinsics for dual<float> and dual<double>.

- add pmadd specializations for dual and complex<dual<float>> NEON packets, improving Eigen GEMM performance.

- add NEON benchmark comparison tables to README.

- fix inverted _novec benchmark binary naming in CMakeLists.



v0.8.2

------

duals/arch/NEON/ComplexDual.h

+74 −8
Original line number Diff line number Diff line
@@ -70,12 +70,47 @@ template<> EIGEN_STRONG_INLINE Packet1cdf pconj(const Packet1cdf& a) 


template<> EIGEN_STRONG_INLINE Packet1cdf pmul<Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)

{

  // can probably do this more efficiently...

  //= a0.b0         - a2.b2

  //= a0.b1 + a1.b0 - (a2.b3 + a3.b2)

  //= a0.b2 +         a2.b0

  //= a0.b3 + a1.b2 + a2.b1 + a3.b0

  // complex<dual<float>> multiply:  a = [ar.r, ar.d, ai.r, ai.d]

  //                                 b = [br.r, br.d, bi.r, bi.d]

  // result_re = ar*br - ai*bi  (dual multiply & subtract)

  // result_im = ar*bi + ai*br  (dual multiply & add)

  //

  // r[0] = a0*b0 - a2*b2

  // r[1] = a0*b1 + a1*b0 - a2*b3 - a3*b2

  // r[2] = a0*b2 + a2*b0

  // r[3] = a0*b3 + a1*b2 + a2*b1 + a3*b0

#if EIGEN_ARCH_ARM64

  float32x4_t av = a.v.v, bv = b.v.v;



  // Broadcast each scalar element of a

  float32x4_t a0 = vdupq_laneq_f32(av, 0);   // ar.r

  float32x4_t a2 = vdupq_laneq_f32(av, 2);   // ai.r

  float32x4_t a1 = vdupq_laneq_f32(av, 1);   // ar.d

  float32x4_t a3 = vdupq_laneq_f32(av, 3);   // ai.d



  // Swap complex halves: [bi.r, bi.d, br.r, br.d]

  float32x4_t b_swap = vextq_f32(bv, bv, 2);



  // Negate real pair for complex structure: [-bi.r, -bi.d, br.r, br.d]

  const uint32x4_t neg_re = {0x80000000, 0x80000000, 0, 0};

  float32x4_t b_swap_neg = vreinterpretq_f32_u32(

    veorq_u32(vreinterpretq_u32_f32(b_swap), neg_re));



  // Real-part contributions (even lanes are final):

  // r = a0*b + a2*b_swap_neg

  float32x4_t r = vfmaq_f32(vmulq_f32(a0, bv), a2, b_swap_neg);



  // Dual-part contributions (odd lanes only):

  // d = a1*b_reals + a3*b_swap_reals_neg  (same complex structure on real parts)

  float32x4_t b_r = vtrn1q_f32(bv, bv);              // [b0, b0, b2, b2]

  float32x4_t b_swap_r_neg = vreinterpretq_f32_u32(

    veorq_u32(vreinterpretq_u32_f32(vtrn1q_f32(b_swap, b_swap)), neg_re));

  float32x4_t d = vfmaq_f32(vmulq_f32(a1, b_r), a3, b_swap_r_neg);



  // Blend: even lanes from r, odd lanes from r+d

  const uint32x4_t odd_mask = {0, 0xffffffff, 0, 0xffffffff};

  return Packet1cdf(vbslq_f32(odd_mask, vaddq_f32(r, d), r));

#else

  uint32x4_t mask = {0x00000000, 0xffffffff, 0x00000000, 0xffffffff};

  uint32x4_t nega = {0x80000000, 0x80000000, 0x00000000, 0x00000000};

  return Packet1cdf(vaddq_f32(vaddq_f32(vmulq_f32(vdupq_lane_f32(vget_low_f32(a.v.v), 0),
@@ -98,6 +133,37 @@ template<> EIGEN_STRONG_INLINE Packet1cdf pmul<Packet1cdf>(const Packet1cdf& a, 
                                                             (float32x4_t)vcombine_u64(vsli_n_u64(vget_low_u64((uint64x2_t)b.v.v),

                                                                                     vget_high_u64((uint64x2_t)b.v.v), 32),

                                                                          (uint64x1_t)vdup_lane_u32(vget_low_u32((uint32x4_t)b.v.v), 0))))))));

#endif

}



template<> EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& a, const Packet1cdf& b, const Packet1cdf& c)

{

#if EIGEN_ARCH_ARM64

  float32x4_t av = a.v.v, bv = b.v.v;



  float32x4_t a0 = vdupq_laneq_f32(av, 0);

  float32x4_t a2 = vdupq_laneq_f32(av, 2);

  float32x4_t a1 = vdupq_laneq_f32(av, 1);

  float32x4_t a3 = vdupq_laneq_f32(av, 3);



  float32x4_t b_swap = vextq_f32(bv, bv, 2);

  const uint32x4_t neg_re = {0x80000000, 0x80000000, 0, 0};

  float32x4_t b_swap_neg = vreinterpretq_f32_u32(

    veorq_u32(vreinterpretq_u32_f32(b_swap), neg_re));



  // r = c + a0*b + a2*b_swap_neg

  float32x4_t r = vfmaq_f32(vfmaq_f32(c.v.v, a0, bv), a2, b_swap_neg);



  float32x4_t b_r = vtrn1q_f32(bv, bv);

  float32x4_t b_swap_r_neg = vreinterpretq_f32_u32(

    veorq_u32(vreinterpretq_u32_f32(vtrn1q_f32(b_swap, b_swap)), neg_re));

  float32x4_t d = vfmaq_f32(vmulq_f32(a1, b_r), a3, b_swap_r_neg);



  const uint32x4_t odd_mask = {0, 0xffffffff, 0, 0xffffffff};

  return Packet1cdf(vbslq_f32(odd_mask, vaddq_f32(r, d), r));

#else

  return padd(pmul(a, b), c);

#endif

}



template<> EIGEN_STRONG_INLINE Packet1cdf pand   <Packet1cdf>(const Packet1cdf& a, const Packet1cdf& b)
@@ -169,7 +235,7 @@ struct palign_impl<Offset,Packet1cdf> 
template<> struct conj_helper<Packet1cdf, Packet1cdf, false,true>

{

  EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(x, pconj(y), c); }



  EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const

  {
@@ -180,7 +246,7 @@ template<> struct conj_helper<Packet1cdf, Packet1cdf, false,true> 
template<> struct conj_helper<Packet1cdf, Packet1cdf, true,false>

{

  EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(pconj(x), y, c); }



  EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const

  {
@@ -191,7 +257,7 @@ template<> struct conj_helper<Packet1cdf, Packet1cdf, true,false> 
template<> struct conj_helper<Packet1cdf, Packet1cdf, true,true>

{

  EIGEN_STRONG_INLINE Packet1cdf pmadd(const Packet1cdf& x, const Packet1cdf& y, const Packet1cdf& c) const

  { return padd(pmul(x,y),c); }

  { return pconj(internal::pmadd(x, y, pconj(c))); }



  EIGEN_STRONG_INLINE Packet1cdf pmul(const Packet1cdf& a, const Packet1cdf& b) const

  {

duals/arch/NEON/Dual.h

+68 −27
Original line number Diff line number Diff line
@@ -87,6 +87,15 @@ template<> EIGEN_STRONG_INLINE Packet2df pconj(const Packet2df& a) { return a; } 


template<> EIGEN_STRONG_INLINE Packet2df pmul<Packet2df>(const Packet2df& a, const Packet2df& b)

{

#if EIGEN_ARCH_ARM64

  // Dual multiply: (a.r*b.r, a.r*b.d + a.d*b.r) using AArch64 FMA

  float32x4_t ar = vtrn1q_f32(a.v, a.v);         // [a0.r, a0.r, a1.r, a1.r]

  float32x4_t br = vtrn1q_f32(b.v, b.v);         // [b0.r, b0.r, b1.r, b1.r]

  float32x4_t t1 = vmulq_f32(ar, b.v);           // [a.r*b.r, a.r*b.d, ...]

  float32x4_t t2 = vfmaq_f32(t1, a.v, br);       // [2*a.r*b.r, a.r*b.d+a.d*b.r, ...]

  const uint32x4_t odd_mask = {0, 0xffffffff, 0, 0xffffffff};

  return Packet2df(vbslq_f32(odd_mask, t2, t1));  // even from t1, odd from t2

#else

  Packet4f v1, v2;

  uint32x4_t mask = {0x00000000, 0xffffffff, 0x00000000, 0xffffffff};

  v1 = vcombine_f32(vdup_lane_f32(vget_low_f32(b.v), 0), vdup_lane_f32(vget_high_f32(b.v), 0));
@@ -94,6 +103,22 @@ template<> EIGEN_STRONG_INLINE Packet2df pmul<Packet2df>(const Packet2df& a, con 
  return Packet2df(vaddq_f32(vreinterpretq_f32_u32(vandq_u32(mask,

                                                             vreinterpretq_u32_f32(vmulq_f32(a.v, v1)))),

                             vmulq_f32(b.v, v2)));

#endif

}



template<> EIGEN_STRONG_INLINE Packet2df pmadd(const Packet2df& a, const Packet2df& b, const Packet2df& c)

{

#if EIGEN_ARCH_ARM64

  // Fused dual multiply-add: a*b + c using two FMAs and a bit-select

  float32x4_t ar = vtrn1q_f32(a.v, a.v);

  float32x4_t br = vtrn1q_f32(b.v, b.v);

  float32x4_t t1 = vfmaq_f32(c.v, ar, b.v);      // [c.r+a.r*b.r, c.d+a.r*b.d, ...]

  float32x4_t t2 = vfmaq_f32(t1, a.v, br);        // [c.r+2*a.r*b.r, c.d+a.r*b.d+a.d*b.r, ...]

  const uint32x4_t odd_mask = {0, 0xffffffff, 0, 0xffffffff};

  return Packet2df(vbslq_f32(odd_mask, t2, t1));

#else

  return padd(pmul(a, b), c);

#endif

}



template<> EIGEN_STRONG_INLINE Packet2df pand   <Packet2df>(const Packet2df& a, const Packet2df& b)
@@ -214,7 +239,7 @@ struct palign_impl<Offset,Packet2df> 
template<> struct conj_helper<Packet2df, Packet2df, false,true>

{

  EIGEN_STRONG_INLINE Packet2df pmadd(const Packet2df& x, const Packet2df& y, const Packet2df& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(x, pconj(y), c); }



  EIGEN_STRONG_INLINE Packet2df pmul(const Packet2df& a, const Packet2df& b) const

  {
@@ -225,7 +250,7 @@ template<> struct conj_helper<Packet2df, Packet2df, false,true> 
template<> struct conj_helper<Packet2df, Packet2df, true,false>

{

  EIGEN_STRONG_INLINE Packet2df pmadd(const Packet2df& x, const Packet2df& y, const Packet2df& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(pconj(x), y, c); }



  EIGEN_STRONG_INLINE Packet2df pmul(const Packet2df& a, const Packet2df& b) const

  {
@@ -236,7 +261,7 @@ template<> struct conj_helper<Packet2df, Packet2df, true,false> 
template<> struct conj_helper<Packet2df, Packet2df, true,true>

{

  EIGEN_STRONG_INLINE Packet2df pmadd(const Packet2df& x, const Packet2df& y, const Packet2df& c) const

  { return padd(pmul(x,y),c); }

  { return pconj(internal::pmadd(x, y, pconj(c))); }



  EIGEN_STRONG_INLINE Packet2df pmul(const Packet2df& a, const Packet2df& b) const

  {
@@ -248,6 +273,17 @@ EIGEN_MAKE_CONJ_HELPER_CPLX_REAL(Packet2df,Packet4f) 


template<> EIGEN_STRONG_INLINE Packet2df pdiv<Packet2df>(const Packet2df& a, const Packet2df& b)

{

#if EIGEN_ARCH_ARM64

  // Dual division: (a.r/b.r, (a.d*b.r - a.r*b.d) / b.r^2)

  float32x4_t ar = vtrn1q_f32(a.v, a.v);         // [a.r, a.r, ...]

  float32x4_t br = vtrn1q_f32(b.v, b.v);         // [b.r, b.r, ...]

  // cross = a*br - ar*b = [0, a.d*b.r - a.r*b.d, ...]

  float32x4_t cross = vfmsq_f32(vmulq_f32(a.v, br), ar, b.v);

  const uint32x4_t odd_mask = {0, 0xffffffff, 0, 0xffffffff};

  float32x4_t num = vbslq_f32(odd_mask, cross, a.v);            // [a.r, cross, ...]

  float32x4_t den = vbslq_f32(odd_mask, vmulq_f32(br, br), br); // [b.r, b.r^2, ...]

  return Packet2df(vdivq_f32(num, den));

#else

  uint32x4_t mask = {0xffffffff, 0x00000000, 0xffffffff, 0x00000000};

  float32x4_t c = vcombine_f32(vdup_lane_f32(vget_low_f32(b.v), 0), vdup_lane_f32(vget_high_f32(b.v), 0));

  float32x4_t d = vcombine_f32(vdup_lane_f32(vget_low_f32(a.v), 0), vdup_lane_f32(vget_high_f32(a.v), 0));
@@ -263,6 +299,7 @@ template<> EIGEN_STRONG_INLINE Packet2df pdiv<Packet2df>(const Packet2df& a, con 
                                                              vreinterpretq_u32_f32(b.v))),

                              vreinterpretq_f32_u32(vandq_u32(vmvnq_u32(mask),

                                                              vreinterpretq_u32_f32(vmulq_f32(c, c)))))));

#endif

}



EIGEN_DEVICE_FUNC inline void
@@ -330,15 +367,24 @@ template<> EIGEN_STRONG_INLINE Packet1dd pconj(const Packet1dd& a) { return a; } 


template<> EIGEN_STRONG_INLINE Packet1dd pmul<Packet1dd>(const Packet1dd& a, const Packet1dd& b)

{

  const uint64x2_t mask = {0, 0xffffffffffffffff};

  return Packet1dd

    (vaddq_f64(vreinterpretq_f64_u64

               (vandq_u64(mask,

                          vreinterpretq_u64_f64(vmulq_f64(a.v,

                                                          vdupq_lane_f64(vget_low_f64(b.v), 0))))),

               vmulq_f64(vdupq_lane_f64(vget_low_f64(a.v), 0),

                         b.v)));

  // Dual multiply: (a.r*b.r, a.r*b.d + a.d*b.r)

  float64x2_t ar = vdupq_laneq_f64(a.v, 0);     // [a.r, a.r]

  float64x2_t br = vdupq_laneq_f64(b.v, 0);     // [b.r, b.r]

  float64x2_t t1 = vmulq_f64(ar, b.v);           // [a.r*b.r, a.r*b.d]

  float64x2_t t2 = vfmaq_f64(t1, a.v, br);       // [2*a.r*b.r, a.r*b.d+a.d*b.r]

  const uint64x2_t odd_mask = {0, 0xffffffffffffffff};

  return Packet1dd(vbslq_f64(odd_mask, t2, t1));

}



template<> EIGEN_STRONG_INLINE Packet1dd pmadd(const Packet1dd& a, const Packet1dd& b, const Packet1dd& c)

{

  // Fused dual multiply-add: a*b + c

  float64x2_t ar = vdupq_laneq_f64(a.v, 0);

  float64x2_t br = vdupq_laneq_f64(b.v, 0);

  float64x2_t t1 = vfmaq_f64(c.v, ar, b.v);      // [c.r+a.r*b.r, c.d+a.r*b.d]

  float64x2_t t2 = vfmaq_f64(t1, a.v, br);        // [c.r+2*a.r*b.r, c.d+a.r*b.d+a.d*b.r]

  const uint64x2_t odd_mask = {0, 0xffffffffffffffff};

  return Packet1dd(vbslq_f64(odd_mask, t2, t1));

}



template<> EIGEN_STRONG_INLINE Packet1dd pand   <Packet1dd>(const Packet1dd& a, const Packet1dd& b)
@@ -408,7 +454,7 @@ struct palign_impl<Offset,Packet1dd> 
template<> struct conj_helper<Packet1dd, Packet1dd, false,true>

{

  EIGEN_STRONG_INLINE Packet1dd pmadd(const Packet1dd& x, const Packet1dd& y, const Packet1dd& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(x, pconj(y), c); }



  EIGEN_STRONG_INLINE Packet1dd pmul(const Packet1dd& a, const Packet1dd& b) const

  {
@@ -419,7 +465,7 @@ template<> struct conj_helper<Packet1dd, Packet1dd, false,true> 
template<> struct conj_helper<Packet1dd, Packet1dd, true,false>

{

  EIGEN_STRONG_INLINE Packet1dd pmadd(const Packet1dd& x, const Packet1dd& y, const Packet1dd& c) const

  { return padd(pmul(x,y),c); }

  { return internal::pmadd(pconj(x), y, c); }



  EIGEN_STRONG_INLINE Packet1dd pmul(const Packet1dd& a, const Packet1dd& b) const

  {
@@ -430,7 +476,7 @@ template<> struct conj_helper<Packet1dd, Packet1dd, true,false> 
template<> struct conj_helper<Packet1dd, Packet1dd, true,true>

{

  EIGEN_STRONG_INLINE Packet1dd pmadd(const Packet1dd& x, const Packet1dd& y, const Packet1dd& c) const

  { return padd(pmul(x,y),c); }

  { return pconj(internal::pmadd(x, y, pconj(c))); }



  EIGEN_STRONG_INLINE Packet1dd pmul(const Packet1dd& a, const Packet1dd& b) const

  {
@@ -442,19 +488,14 @@ EIGEN_MAKE_CONJ_HELPER_CPLX_REAL(Packet1dd,Packet2d) 


template<> EIGEN_STRONG_INLINE Packet1dd pdiv<Packet1dd>(const Packet1dd& a, const Packet1dd& b)

{

  const uint64x2_t mask = {0xffffffffffffffff, 0};

  const uint64x2_t imask = {0, 0xffffffffffffffff};

  return Packet1dd(vdivq_f64

                   (vaddq_f64(vreinterpretq_f64_u64(vandq_u64(mask, vreinterpretq_u64_f64(a.v))),

                              vreinterpretq_f64_u64(vandq_u64(imask,

                                        vreinterpretq_u64_f64(vsubq_f64(vmulq_f64(a.v,

                                                            vdupq_lane_f64(vget_low_f64(b.v), 0)),

                                                  vmulq_f64(vdupq_lane_f64(vget_low_f64(a.v), 0),

                                                            b.v)))))),

                    vaddq_f64(vreinterpretq_f64_u64(vandq_u64(mask, vreinterpretq_u64_f64(b.v))),

                              vreinterpretq_f64_u64(vandq_u64(imask,

                                        vreinterpretq_u64_f64(vmulq_f64(vdupq_lane_f64(vget_low_f64(b.v), 0),

                                                                        vdupq_lane_f64(vget_low_f64(b.v), 0))))))));

  // Dual division: (a.r/b.r, (a.d*b.r - a.r*b.d) / b.r^2)

  float64x2_t ar = vdupq_laneq_f64(a.v, 0);     // [a.r, a.r]

  float64x2_t br = vdupq_laneq_f64(b.v, 0);     // [b.r, b.r]

  float64x2_t cross = vfmsq_f64(vmulq_f64(a.v, br), ar, b.v); // [0, a.d*b.r - a.r*b.d]

  const uint64x2_t odd_mask = {0, 0xffffffffffffffff};

  float64x2_t num = vbslq_f64(odd_mask, cross, a.v);            // [a.r, cross]

  float64x2_t den = vbslq_f64(odd_mask, vmulq_f64(br, br), br); // [b.r, b.r^2]

  return Packet1dd(vdivq_f64(num, den));

}



EIGEN_STRONG_INLINE Packet1dd pcplxflip/*<Packet1dd>*/(const Packet1dd& x)

tests/CMakeLists.txt

+3 −3
Original line number Diff line number Diff line
@@ -219,7 +219,7 @@ if (CPPDUALS_BENCHMARK) 


  foreach (VECTORIZE YES NO)

    foreach (BENCH  bench_dual bench_eigen bench_exp bench_gemm bench_example bench_fmt)

      if (NOT VECTORIZE)

      if (VECTORIZE)

        set (BENCHE ${BENCH})

      else ()

        set (BENCHE ${BENCH}_novec)
@@ -244,7 +244,7 @@ if (CPPDUALS_BENCHMARK) 
  # Compiled 4 ways: old/new × vectorized/novec

  foreach (VECTORIZE YES NO)

    foreach (IMPL old new)

      if (NOT VECTORIZE)

      if (VECTORIZE)

        set (BMT bench_expm_deriv_${IMPL})

      else ()

        set (BMT bench_expm_deriv_${IMPL}_novec)
@@ -267,7 +267,7 @@ if (CPPDUALS_BENCHMARK) 
  # bench_multidual: compiled 4 ways (old/new × vectorized/novec)

  foreach (VECTORIZE YES NO)

    foreach (IMPL old new)

      if (NOT VECTORIZE)

      if (VECTORIZE)

        set (BMT bench_multidual_${IMPL})

      else ()

        set (BMT bench_multidual_${IMPL}_novec)