Eigen  5.0.1
 
Loading...
Searching...
No Matches
Complex.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2016 Konstantinos Margaritis <markos@freevec.org>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10// SPDX-License-Identifier: MPL-2.0
11
12#ifndef EIGEN_COMPLEX32_ZVECTOR_H
13#define EIGEN_COMPLEX32_ZVECTOR_H
14
15// IWYU pragma: private
16#include "../../InternalHeaderCheck.h"
17
18namespace Eigen {
19
20namespace internal {
21
22EIGEN_GCC_FAST_MATH_COMPLEX_VECTORIZE_WORKAROUND_PUSH
23
24inline Packet4ui p4ui_CONJ_XOR() {
25 return Packet4ui{0x00000000, 0x80000000, 0x00000000,
26 0x80000000}; // vec_mergeh((Packet4ui)p4i_ZERO, (Packet4ui)p4f_MZERO);
27}
28
29inline Packet2ul p2ul_CONJ_XOR1() {
30 return (Packet2ul)vec_sld((Packet4ui)p2d_ZERO_, (Packet4ui)p2l_ZERO,
31 8); //{ 0x8000000000000000, 0x0000000000000000 };
32}
33inline Packet2ul p2ul_CONJ_XOR2() {
34 return (Packet2ul)vec_sld((Packet4ui)p2l_ZERO, (Packet4ui)p2d_ZERO_,
35 8); //{ 0x8000000000000000, 0x0000000000000000 };
36}
37
38struct Packet1cd {
39 EIGEN_STRONG_INLINE Packet1cd() {}
40 EIGEN_STRONG_INLINE explicit Packet1cd(const Packet2d& a) : v(a) {}
41 Packet2d v;
42};
43
44struct Packet2cf {
45 EIGEN_STRONG_INLINE Packet2cf() {}
46 EIGEN_STRONG_INLINE explicit Packet2cf(const Packet4f& a) : v(a) {}
47 Packet4f v;
48};
49
50template <>
51struct packet_traits<std::complex<float> > : default_packet_traits {
52 typedef Packet2cf type;
53 typedef Packet2cf half;
54 enum {
55 Vectorizable = 1,
56 AlignedOnScalar = 1,
57 size = 2,
58
59 HasAdd = 1,
60 HasSub = 1,
61 HasMul = 1,
62 HasDiv = 1,
63 HasLog = 1,
64 HasExp = 1,
65 HasNegate = 1,
66 HasAbs = 0,
67 HasAbs2 = 0,
68 HasMin = 0,
69 HasMax = 0,
70 HasSetLinear = 0
71 };
72};
73
74template <>
75struct packet_traits<std::complex<double> > : default_packet_traits {
76 typedef Packet1cd type;
77 typedef Packet1cd half;
78 enum {
79 Vectorizable = 1,
80 AlignedOnScalar = 1,
81 size = 1,
82
83 HasAdd = 1,
84 HasSub = 1,
85 HasMul = 1,
86 HasDiv = 1,
87 HasLog = 1,
88 HasNegate = 1,
89 HasAbs = 0,
90 HasAbs2 = 0,
91 HasMin = 0,
92 HasMax = 0,
93 HasSetLinear = 0
94 };
95};
96
97template <>
98struct unpacket_traits<Packet2cf> {
99 typedef std::complex<float> type;
100 enum {
101 size = 2,
102 alignment = Aligned16,
103 vectorizable = true,
104 masked_load_available = false,
105 masked_store_available = false
106 };
107 typedef Packet2cf half;
108 typedef Packet4f as_real;
109};
110template <>
111struct unpacket_traits<Packet1cd> {
112 typedef std::complex<double> type;
113 enum {
114 size = 1,
115 alignment = Aligned16,
116 vectorizable = true,
117 masked_load_available = false,
118 masked_store_available = false
119 };
120 typedef Packet1cd half;
121 typedef Packet2d as_real;
122};
123
124/* Forward declaration */
125EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf, 2>& kernel);
126
127/* complex<double> first */
128template <>
129EIGEN_STRONG_INLINE Packet1cd pload<Packet1cd>(const std::complex<double>* from) {
130 EIGEN_DEBUG_ALIGNED_LOAD return Packet1cd(pload<Packet2d>((const double*)from));
131}
132template <>
133EIGEN_STRONG_INLINE Packet1cd ploadu<Packet1cd>(const std::complex<double>* from) {
134 EIGEN_DEBUG_UNALIGNED_LOAD return Packet1cd(ploadu<Packet2d>((const double*)from));
135}
136template <>
137EIGEN_STRONG_INLINE void pstore<std::complex<double> >(std::complex<double>* to, const Packet1cd& from) {
138 EIGEN_DEBUG_ALIGNED_STORE pstore((double*)to, from.v);
139}
140template <>
141EIGEN_STRONG_INLINE void pstoreu<std::complex<double> >(std::complex<double>* to, const Packet1cd& from) {
142 EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, from.v);
143}
144
145template <>
146EIGEN_STRONG_INLINE Packet1cd
147pset1<Packet1cd>(const std::complex<double>& from) { /* here we really have to use unaligned loads :( */
148 return ploadu<Packet1cd>(&from);
149}
150
151template <>
152EIGEN_DEVICE_FUNC inline Packet1cd pgather<std::complex<double>, Packet1cd>(const std::complex<double>* from,
153 Index stride EIGEN_UNUSED) {
154 return pload<Packet1cd>(from);
155}
156template <>
157EIGEN_DEVICE_FUNC inline void pscatter<std::complex<double>, Packet1cd>(std::complex<double>* to, const Packet1cd& from,
158 Index stride EIGEN_UNUSED) {
159 pstore<std::complex<double> >(to, from);
160}
161template <>
162EIGEN_STRONG_INLINE Packet1cd padd<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
163 return Packet1cd(a.v + b.v);
164}
165template <>
166EIGEN_STRONG_INLINE Packet1cd psub<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
167 return Packet1cd(a.v - b.v);
168}
169template <>
170EIGEN_STRONG_INLINE Packet1cd pnegate(const Packet1cd& a) {
171 return Packet1cd(pnegate(Packet2d(a.v)));
172}
173template <>
174EIGEN_STRONG_INLINE Packet1cd pconj(const Packet1cd& a) {
175 return Packet1cd((Packet2d)vec_xor((Packet2d)a.v, (Packet2d)p2ul_CONJ_XOR2()));
176}
177template <>
178EIGEN_STRONG_INLINE Packet1cd pmul<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
179 Packet2d a_re, a_im, v1, v2;
180
181 // Permute and multiply the real parts of a and b
182 a_re = vec_perm(a.v, a.v, p16uc_PSET64_HI);
183 // Get the imaginary parts of a
184 a_im = vec_perm(a.v, a.v, p16uc_PSET64_LO);
185 // multiply a_re * b
186 v1 = vec_madd(a_re, b.v, p2d_ZERO);
187 // multiply a_im * b and get the conjugate result
188 v2 = vec_madd(a_im, b.v, p2d_ZERO);
189 v2 = (Packet2d)vec_sld((Packet4ui)v2, (Packet4ui)v2, 8);
190 v2 = (Packet2d)vec_xor((Packet2d)v2, (Packet2d)p2ul_CONJ_XOR1());
191
192 return Packet1cd(v1 + v2);
193}
194template <>
195EIGEN_STRONG_INLINE Packet1cd pand<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
196 return Packet1cd(vec_and(a.v, b.v));
197}
198template <>
199EIGEN_STRONG_INLINE Packet1cd por<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
200 return Packet1cd(vec_or(a.v, b.v));
201}
202template <>
203EIGEN_STRONG_INLINE Packet1cd pxor<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
204 return Packet1cd(vec_xor(a.v, b.v));
205}
206template <>
207EIGEN_STRONG_INLINE Packet1cd pandnot<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
208 return Packet1cd(vec_and(a.v, vec_nor(b.v, b.v)));
209}
210template <>
211EIGEN_STRONG_INLINE Packet1cd pcmp_eq(const Packet1cd& a, const Packet1cd& b) {
212 Packet2d eq = vec_cmpeq(a.v, b.v);
213 Packet2d tmp = {eq[1], eq[0]};
214 return (Packet1cd)pand<Packet2d>(eq, tmp);
215}
216
217template <>
218EIGEN_STRONG_INLINE void prefetch<std::complex<double> >(const std::complex<double>* addr) {
219 EIGEN_ZVECTOR_PREFETCH(addr);
220}
221
222template <>
223EIGEN_STRONG_INLINE std::complex<double> pfirst<Packet1cd>(const Packet1cd& a) {
224 EIGEN_ALIGN16 std::complex<double> res;
225 pstore<std::complex<double> >(&res, a);
226
227 return res;
228}
229
230EIGEN_MAKE_CONJ_HELPER_CPLX_REAL(Packet1cd, Packet2d)
231
232template <>
233EIGEN_STRONG_INLINE Packet1cd pdiv<Packet1cd>(const Packet1cd& a, const Packet1cd& b) {
234 return pdiv_complex(a, b);
235}
236
237EIGEN_INSTANTIATE_COMPLEX_MATH_FUNCS_NO_EXP(Packet1cd)
238EIGEN_INSTANTIATE_COMPLEX_MATH_FUNCS(Packet2cf)
239
240EIGEN_STRONG_INLINE Packet1cd pcplxflip /*<Packet1cd>*/ (const Packet1cd& x) {
241 return Packet1cd(preverse(Packet2d(x.v)));
242}
243
244EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet1cd, 2>& kernel) {
245 Packet2d tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_HI);
246 kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO);
247 kernel.packet[0].v = tmp;
248}
249
250/* complex<float> follows */
251template <>
252EIGEN_STRONG_INLINE Packet2cf pload<Packet2cf>(const std::complex<float>* from) {
253 EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>((const float*)from));
254}
255template <>
256EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) {
257 EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>((const float*)from));
258}
259template <>
260EIGEN_STRONG_INLINE void pstore<std::complex<float> >(std::complex<float>* to, const Packet2cf& from) {
261 EIGEN_DEBUG_ALIGNED_STORE pstore((float*)to, from.v);
262}
263template <>
264EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float>* to, const Packet2cf& from) {
265 EIGEN_DEBUG_UNALIGNED_STORE pstoreu((float*)to, from.v);
266}
267
268template <>
269EIGEN_STRONG_INLINE std::complex<float> pfirst<Packet2cf>(const Packet2cf& a) {
270 EIGEN_ALIGN16 std::complex<float> res[2];
271 pstore<std::complex<float> >(res, a);
272
273 return res[0];
274}
275
276template <>
277EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<float>& from) {
278 Packet2cf res;
279 if ((std::ptrdiff_t(&from) % 16) == 0)
280 res.v = pload<Packet4f>((const float*)&from);
281 else
282 res.v = ploadu<Packet4f>((const float*)&from);
283 res.v = vec_perm(res.v, res.v, p16uc_PSET64_HI);
284 return res;
285}
286
287template <>
288EIGEN_DEVICE_FUNC inline Packet2cf pgather<std::complex<float>, Packet2cf>(const std::complex<float>* from,
289 Index stride) {
290 EIGEN_ALIGN16 std::complex<float> af[2];
291 af[0] = from[0 * stride];
292 af[1] = from[1 * stride];
293 return pload<Packet2cf>(af);
294}
295template <>
296EIGEN_DEVICE_FUNC inline void pscatter<std::complex<float>, Packet2cf>(std::complex<float>* to, const Packet2cf& from,
297 Index stride) {
298 EIGEN_ALIGN16 std::complex<float> af[2];
299 pstore<std::complex<float> >((std::complex<float>*)af, from);
300 to[0 * stride] = af[0];
301 to[1 * stride] = af[1];
302}
303
304template <>
305EIGEN_STRONG_INLINE Packet2cf padd<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
306 return Packet2cf(padd<Packet4f>(a.v, b.v));
307}
308template <>
309EIGEN_STRONG_INLINE Packet2cf psub<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
310 return Packet2cf(psub<Packet4f>(a.v, b.v));
311}
312template <>
313EIGEN_STRONG_INLINE Packet2cf pnegate(const Packet2cf& a) {
314 return Packet2cf(pnegate(Packet4f(a.v)));
315}
316
317template <>
318EIGEN_STRONG_INLINE Packet2cf pand<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
319 return Packet2cf(pand<Packet4f>(a.v, b.v));
320}
321template <>
322EIGEN_STRONG_INLINE Packet2cf por<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
323 return Packet2cf(por<Packet4f>(a.v, b.v));
324}
325template <>
326EIGEN_STRONG_INLINE Packet2cf pxor<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
327 return Packet2cf(pxor<Packet4f>(a.v, b.v));
328}
329template <>
330EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
331 return Packet2cf(pandnot<Packet4f>(a.v, b.v));
332}
333
334template <>
335EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>* from) {
336 return pset1<Packet2cf>(*from);
337}
338
339template <>
340EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float>* addr) {
341 EIGEN_ZVECTOR_PREFETCH(addr);
342}
343
344template <>
345EIGEN_STRONG_INLINE Packet2cf pcmp_eq(const Packet2cf& a, const Packet2cf& b) {
346 Packet4f eq = vec_cmpeq(a.v, b.v);
347 Packet4f tmp = {eq[1], eq[0], eq[3], eq[2]};
348 return (Packet2cf)pand<Packet4f>(eq, tmp);
349}
350template <>
351EIGEN_STRONG_INLINE Packet2cf pconj(const Packet2cf& a) {
352 return Packet2cf(pxor<Packet4f>(a.v, reinterpret_cast<Packet4f>(p4ui_CONJ_XOR())));
353}
354template <>
355EIGEN_STRONG_INLINE Packet2cf pmul<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
356 Packet4f a_re, a_im, prod, prod_im;
357
358 // Permute and multiply the real parts of a and b
359 a_re = vec_perm(a.v, a.v, p16uc_PSET32_WODD);
360
361 // Get the imaginary parts of a
362 a_im = vec_perm(a.v, a.v, p16uc_PSET32_WEVEN);
363
364 // multiply a_im * b and get the conjugate result
365 prod_im = a_im * b.v;
366 prod_im = pxor<Packet4f>(prod_im, reinterpret_cast<Packet4f>(p4ui_CONJ_XOR()));
367 // permute back to a proper order
368 prod_im = vec_perm(prod_im, prod_im, p16uc_COMPLEX32_REV);
369
370 // multiply a_re * b, add prod_im
371 prod = pmadd<Packet4f>(a_re, b.v, prod_im);
372
373 return Packet2cf(prod);
374}
375
376template <>
377EIGEN_STRONG_INLINE Packet2cf preverse(const Packet2cf& a) {
378 Packet4f rev_a;
379 rev_a = vec_perm(a.v, a.v, p16uc_COMPLEX32_REV2);
380 return Packet2cf(rev_a);
381}
382
383template <>
384EIGEN_STRONG_INLINE std::complex<float> predux<Packet2cf>(const Packet2cf& a) {
385 Packet4f b;
386 b = vec_sld(a.v, a.v, 8);
387 b = padd<Packet4f>(a.v, b);
388 return pfirst<Packet2cf>(Packet2cf(b));
389}
390
391template <>
392EIGEN_STRONG_INLINE std::complex<float> predux_mul<Packet2cf>(const Packet2cf& a) {
393 Packet4f b;
394 Packet2cf prod;
395 b = vec_sld(a.v, a.v, 8);
396 prod = pmul<Packet2cf>(a, Packet2cf(b));
397
398 return pfirst<Packet2cf>(prod);
399}
400
401EIGEN_MAKE_CONJ_HELPER_CPLX_REAL(Packet2cf, Packet4f)
402
403template <>
404EIGEN_STRONG_INLINE Packet2cf pdiv<Packet2cf>(const Packet2cf& a, const Packet2cf& b) {
405 return pdiv_complex(a, b);
406}
407
408template <>
409EIGEN_STRONG_INLINE Packet2cf pcplxflip<Packet2cf>(const Packet2cf& x) {
410 return Packet2cf(vec_perm(x.v, x.v, p16uc_COMPLEX32_REV));
411}
412
413EIGEN_STRONG_INLINE void ptranspose(PacketBlock<Packet2cf, 2>& kernel) {
414 Packet4f tmp = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_HI);
415 kernel.packet[1].v = vec_perm(kernel.packet[0].v, kernel.packet[1].v, p16uc_TRANSPOSE64_LO);
416 kernel.packet[0].v = tmp;
417}
418
419EIGEN_GCC_FAST_MATH_COMPLEX_VECTORIZE_WORKAROUND_POP
420
421} // end namespace internal
422
423} // end namespace Eigen
424
425#endif // EIGEN_COMPLEX32_ZVECTOR_H
Eigen::Aligned16
@ Aligned16
Definition Constants.h:238