devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement a generic vectorized version of Smith's algorithms for complex division.

Browse Source
This commit is contained in:
Rasmus Munk Larsen
2021-06-30 15:53:06 -07:00
parent 6035da5283
commit 9312a5bf5c
9 changed files with 46 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
Eigen/src/Core/arch/Default/GenericPacketMathFunctions.h
View File

@@ -757,6 +757,26 @@ Packet pcos_float(const Packet& x)
return psincos_float<false>(x);
}
template<typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS
EIGEN_UNUSED Packet pdiv_complex(const Packet& x, const Packet& y) {
typedef typename unpacket_traits<Packet>::as_real RealPacket;
// In the following we annotate the code for the case where the inputs
// are a pair length-2 SIMD vectors representing a single pair of complex
// numbers x = a + i*b, y = c + i*d.
const RealPacket y_abs = pabs(y.v); // |c|, |d|
const RealPacket y_abs_flip = pcplxflip(Packet(y_abs)).v; // |d|, |c|
const RealPacket y_max = pmax(y_abs, y_abs_flip); // max(|c|, |d|), max(|c|, |d|)
const RealPacket y_scaled = pdiv(y.v, y_max); // c / max(|c|, |d|), d / max(|c|, |d|)
// Compute scaled denominator.
const RealPacket y_scaled_sq = pmul(y_scaled, y_scaled); // c'**2, d'**2
const RealPacket denom = y_scaled_sq + pcplxflip(Packet(y_scaled_sq)).v;
Packet result_scaled = pmul(x, pconj(Packet(y_scaled))); // a * c' + b * d', -a * d + b * c
// Divide elementwise by denom.
result_scaled = Packet(pdiv(result_scaled.v, denom));
// Rescale result
return Packet(pdiv(result_scaled.v, y_max));
}
template<typename Packet>
EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS