Block-Structured AMR Software Framework
AMReX_GpuComplex.H
Go to the documentation of this file.
1 #ifndef AMREX_GPUCOMPLEX_H_
2 #define AMREX_GPUCOMPLEX_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_Algorithm.H>
6 #include <AMReX_Math.H>
7 #include <cmath>
8 #include <iostream>
9 
10 namespace amrex {
11 
12 template <typename T> struct GpuComplex;
13 
14 template <typename T>
16 T norm (const GpuComplex<T>& a_z) noexcept;
17 
27 template <typename T>
28 struct alignas(2*sizeof(T)) GpuComplex
29 {
30  using value_type = T;
31 
33 
38  constexpr GpuComplex (const T& a_r = T(), const T& a_i = T()) noexcept
39  : m_real(a_r), m_imag(a_i) {}
40 
45  constexpr T real () const noexcept { return m_real; }
46 
51  constexpr T imag () const noexcept { return m_imag; }
52 
57  GpuComplex<T>& operator+= (const T& a_t) noexcept
58  {
59  m_real += a_t;
60  return *this;
61  }
62 
67  GpuComplex<T>& operator-= (const T& a_t) noexcept
68  {
69  m_real -= a_t;
70  return *this;
71  }
72 
77  GpuComplex<T>& operator*= (const T& a_t) noexcept
78  {
79  m_real *= a_t;
80  m_imag *= a_t;
81  return *this;
82  }
83 
88  GpuComplex<T>& operator/= (const T& a_t) noexcept
89  {
90  m_real /= a_t;
91  m_imag /= a_t;
92  return *this;
93  }
94 
98  template <typename U>
101  {
102  m_real += a_z.real();
103  m_imag += a_z.imag();
104  return *this;
105  }
106 
110  template <typename U>
113  {
114  m_real -= a_z.real();
115  m_imag -= a_z.imag();
116  return *this;
117  }
118 
122  template <typename U>
125  {
126  const T a_r = m_real * a_z.real() - m_imag * a_z.imag();
127  m_imag = m_real * a_z.imag() + m_imag * a_z.real();
128  m_real = a_r;
129  return *this;
130  }
131 
135  template <typename U>
138  {
139  const T r = m_real * a_z.real() + m_imag * a_z.imag();
140  const T n = amrex::norm(a_z);
141  m_imag = (m_imag * a_z.real() - m_real * a_z.imag()) / n;
142  m_real = r / n;
143  return *this;
144  }
145 
149  template <typename U>
150  friend std::ostream & operator<< (std::ostream &out, const GpuComplex<U>& c);
151 };
152 
153 template <typename U>
154 std::ostream& operator<< (std::ostream& out, const GpuComplex<U>& c)
155 {
156  out << c.m_real;
157  out << " + " << c.m_imag << "i";
158  return out;
159 }
160 
164 template<typename T>
166 GpuComplex<T> operator+ (const GpuComplex<T>& a_x) { return a_x; }
167 
171 template<typename T>
173 GpuComplex<T> operator- (const GpuComplex<T>& a_x) { return GpuComplex<T>(-a_x.real(), -a_x.imag()); }
174 
178 template <typename T>
180 GpuComplex<T> operator- (const GpuComplex<T>& a_x, const GpuComplex<T>& a_y) noexcept
181 {
182  GpuComplex<T> r = a_x;
183  r -= a_y;
184  return r;
185 }
186 
190 template <typename T>
192 GpuComplex<T> operator- (const GpuComplex<T>& a_x, const T& a_y) noexcept
193 {
194  GpuComplex<T> r = a_x;
195  r -= a_y;
196  return r;
197 }
198 
202 template <typename T>
204 GpuComplex<T> operator- (const T& a_x, const GpuComplex<T>& a_y) noexcept
205 {
206  GpuComplex<T> r = a_y;
207  r -= a_x;
208  return -r;
209 }
210 
214 template <typename T>
216 GpuComplex<T> operator+ (const GpuComplex<T>& a_x, const GpuComplex<T>& a_y) noexcept
217 {
218  GpuComplex<T> r = a_x;
219  r += a_y;
220  return r;
221 }
222 
226 template <typename T>
228 GpuComplex<T> operator+ (const GpuComplex<T>& a_x, const T& a_y) noexcept
229 {
230  GpuComplex<T> r = a_x;
231  r += a_y;
232  return r;
233 }
234 
238 template <typename T>
240 GpuComplex<T> operator+ (const T& a_x, const GpuComplex<T>& a_y) noexcept
241 {
242  GpuComplex<T> r = a_y;
243  r += a_x;
244  return r;
245 }
246 
250 template <typename T>
252 GpuComplex<T> operator* (const GpuComplex<T>& a_x, const GpuComplex<T>& a_y) noexcept
253 {
254  GpuComplex<T> r = a_x;
255  r *= a_y;
256  return r;
257 }
258 
262 template <typename T>
264 GpuComplex<T> operator* (const GpuComplex<T>& a_x, const T& a_y) noexcept
265 {
266  GpuComplex<T> r = a_x;
267  r *= a_y;
268  return r;
269 }
270 
274 template <typename T>
276 GpuComplex<T> operator* (const T& a_x, const GpuComplex<T>& a_y) noexcept
277 {
278  GpuComplex<T> r = a_y;
279  r *= a_x;
280  return r;
281 }
282 
286 template <typename T>
288 GpuComplex<T> operator/ (const GpuComplex<T>& a_x, const GpuComplex<T>& a_y) noexcept
289 {
290  GpuComplex<T> r = a_x;
291  r /= a_y;
292  return r;
293 }
294 
298 template <typename T>
300 GpuComplex<T> operator/ (const GpuComplex<T>& a_x, const T& a_y) noexcept
301 {
302  GpuComplex<T> r = a_x;
303  r /= a_y;
304  return r;
305 }
306 
310 template <typename T>
312 GpuComplex<T> operator/ (const T& a_x, const GpuComplex<T>& a_y) noexcept
313 {
314  GpuComplex<T> r = a_x;
315  r /= a_y;
316  return r;
317 }
318 
322 template <typename T>
324 GpuComplex<T> polar (const T& a_r, const T& a_theta) noexcept
325 {
326  return GpuComplex<T>(a_r * std::cos(a_theta), a_r * std::sin(a_theta));
327 }
328 
332 template <typename T>
334 GpuComplex<T> exp (const GpuComplex<T>& a_z) noexcept
335 {
336  return polar<T>(std::exp(a_z.real()), a_z.imag());
337 }
338 
342 template <typename T>
344 T norm (const GpuComplex<T>& a_z) noexcept
345 {
346  const T x = a_z.real();
347  const T y = a_z.imag();
348  return x * x + y * y;
349 }
350 
354 template <typename T>
356 T abs (const GpuComplex<T>& a_z) noexcept
357 {
358  T x = a_z.real();
359  T y = a_z.imag();
360 
361  const T s = amrex::max(std::abs(x), std::abs(y));
362  if (s == T()) { return s; }
363  x /= s;
364  y /= s;
365  return s * std::sqrt(x * x + y * y);
366 }
367 
371 template <typename T>
373 GpuComplex<T> sqrt (const GpuComplex<T>& a_z) noexcept
374 {
375  T x = a_z.real();
376  T y = a_z.imag();
377 
378  if (x == T())
379  {
380  T t = std::sqrt(std::abs(y) / 2);
381  return GpuComplex<T>(t, y < T() ? -t : t);
382  }
383  else
384  {
385  T t = std::sqrt(2 * (amrex::abs(a_z) + std::abs(x)));
386  T u = t / 2;
387  return x > T()
388  ? GpuComplex<T>(u, y / t)
389  : GpuComplex<T>(std::abs(y) / t, y < T() ? -u : u);
390  }
391 }
392 
396 template <typename T>
398 T arg (const GpuComplex<T>& a_z) noexcept
399 {
400  return std::atan2(a_z.imag(), a_z.real());
401 }
402 
406 template <typename T>
408 GpuComplex<T> log (const GpuComplex<T>& a_z) noexcept
409 {
410  return GpuComplex<T>(std::log(amrex::abs(a_z)), amrex::arg(a_z));
411 }
412 
416 template <typename T>
418 GpuComplex<T> pow (const GpuComplex<T>& a_z, const T& a_y) noexcept
419 {
420  if (a_z.imag() == T() && a_z.real() > T()) {
421  return std::pow(a_z.real(), a_y);
422  }
423 
424  GpuComplex<T> t = amrex::log(a_z);
425  return amrex::polar<T>(std::exp(a_y * t.real()), a_y * t.imag());
426 }
427 
428 namespace detail
429 {
430  template <typename T>
433  {
434  GpuComplex<T> y = a_n % 2 ? a_z : GpuComplex<T>(1);
435 
436  while (a_n >>= 1)
437  {
438  a_z *= a_z;
439  if (a_n % 2) {
440  y *= a_z;
441  }
442  }
443 
444  return y;
445  }
446 }
447 
451 template <typename T>
453 GpuComplex<T> pow (const GpuComplex<T>& a_z, int a_n) noexcept
454 {
455  return a_n < 0
456  ? GpuComplex<T>(1) / detail::complex_pow_unsigned(a_z, (unsigned)(-a_n))
457  : detail::complex_pow_unsigned(a_z, a_n);
458 }
459 
460 }
461 
462 #endif
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > complex_pow_unsigned(GpuComplex< T > a_z, unsigned a_n)
Definition: AMReX_GpuComplex.H:432
Definition: AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > operator/(const GpuComplex< T > &a_x, const GpuComplex< T > &a_y) noexcept
Divide a complex number by another one.
Definition: AMReX_GpuComplex.H:288
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > log(const GpuComplex< T > &a_z) noexcept
Complex natural logarithm function.
Definition: AMReX_GpuComplex.H:408
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > polar(const T &a_r, const T &a_theta) noexcept
Return a complex number given its polar representation.
Definition: AMReX_GpuComplex.H:324
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T abs(const GpuComplex< T > &a_z) noexcept
Return the absolute value of a complex number.
Definition: AMReX_GpuComplex.H:356
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, int a_n) noexcept
Raise a complex number to an integer power.
Definition: AMReX_GpuComplex.H:453
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T norm(const GpuComplex< T > &a_z) noexcept
Return the norm (magnitude squared) of a complex number.
Definition: AMReX_GpuComplex.H:344
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T arg(const GpuComplex< T > &a_z) noexcept
Return the angle of a complex number's polar representation.
Definition: AMReX_GpuComplex.H:398
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > pow(const GpuComplex< T > &a_z, const T &a_y) noexcept
Raise a complex number to a (real) power.
Definition: AMReX_GpuComplex.H:418
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > operator-(const GpuComplex< T > &a_x)
Negate a complex number.
Definition: AMReX_GpuComplex.H:173
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > exp(const GpuComplex< T > &a_z) noexcept
Complex expotential function.
Definition: AMReX_GpuComplex.H:334
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > operator+(const GpuComplex< T > &a_x)
Identity operation on a complex number.
Definition: AMReX_GpuComplex.H:166
std::ostream & operator<<(std::ostream &os, AmrMesh const &amr_mesh)
Definition: AMReX_AmrMesh.cpp:1236
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > operator*(const GpuComplex< T > &a_x, const GpuComplex< T > &a_y) noexcept
Multiply two complex numbers.
Definition: AMReX_GpuComplex.H:252
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > sqrt(const GpuComplex< T > &a_z) noexcept
Return the square root of a complex number.
Definition: AMReX_GpuComplex.H:373
Definition: AMReX_FabArrayCommI.H:841
A host / device complex number type, because std::complex doesn't work in device code with Cuda yet.
Definition: AMReX_GpuComplex.H:29
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE GpuComplex(const T &a_r=T(), const T &a_i=T()) noexcept
Construct a complex number given the real and imaginary part.
Definition: AMReX_GpuComplex.H:38
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T real() const noexcept
Return the real part.
Definition: AMReX_GpuComplex.H:45
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > & operator-=(const T &a_t) noexcept
Subtract a real number from this complex number.
Definition: AMReX_GpuComplex.H:67
friend std::ostream & operator<<(std::ostream &out, const GpuComplex< U > &c)
Print this complex number to an ostream.
Definition: AMReX_GpuComplex.H:154
T m_imag
Definition: AMReX_GpuComplex.H:32
T value_type
Definition: AMReX_GpuComplex.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > & operator/=(const T &a_t) noexcept
Divide this complex number by a real.
Definition: AMReX_GpuComplex.H:88
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > & operator*=(const T &a_t) noexcept
Multiply this complex number by a real.
Definition: AMReX_GpuComplex.H:77
T m_real
Definition: AMReX_GpuComplex.H:32
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE T imag() const noexcept
Return the imaginary part.
Definition: AMReX_GpuComplex.H:51
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE GpuComplex< T > & operator+=(const T &a_t) noexcept
Add a real number to this complex number.
Definition: AMReX_GpuComplex.H:57