Kokkos Core Kernels Package  Version of the Day
Kokkos_Complex.hpp
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Kokkos v. 3.0
6 // Copyright (2020) National Technology & Engineering
7 // Solutions of Sandia, LLC (NTESS).
8 //
9 // Under the terms of Contract DE-NA0003525 with NTESS,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY NTESS "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL NTESS OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Christian R. Trott (crtrott@sandia.gov)
40 //
41 // ************************************************************************
42 //@HEADER
43 */
44 #ifndef KOKKOS_COMPLEX_HPP
45 #define KOKKOS_COMPLEX_HPP
46 
47 #include <Kokkos_Atomic.hpp>
48 #include <Kokkos_MathematicalFunctions.hpp>
49 #include <Kokkos_NumericTraits.hpp>
50 #include <impl/Kokkos_Error.hpp>
51 #include <complex>
52 #include <type_traits>
53 #include <iosfwd>
54 
55 namespace Kokkos {
56 
64 template <class RealType>
65 class
66 #ifdef KOKKOS_ENABLE_COMPLEX_ALIGN
67  alignas(2 * sizeof(RealType))
68 #endif
69  complex {
70  private:
71  RealType re_{};
72  RealType im_{};
73 
74  public:
76  using value_type = RealType;
77 
79  KOKKOS_DEFAULTED_FUNCTION
80  complex() noexcept = default;
81 
83  KOKKOS_DEFAULTED_FUNCTION
84  complex(const complex&) noexcept = default;
85 
86  KOKKOS_DEFAULTED_FUNCTION
87  complex& operator=(const complex&) noexcept = default;
88 
90  template <class RType,
91  typename std::enable_if<std::is_convertible<RType, RealType>::value,
92  int>::type = 0>
93  KOKKOS_INLINE_FUNCTION complex(const complex<RType>& other) noexcept
94  // Intentionally do the conversions implicitly here so that users don't
95  // get any warnings about narrowing, etc., that they would expect to get
96  // otherwise.
97  : re_(other.real()), im_(other.imag()) {}
98 
104  KOKKOS_INLINE_FUNCTION
105  complex(const std::complex<RealType>& src) noexcept
106  // We can use this aspect of the standard to avoid calling
107  // non-device-marked functions `std::real` and `std::imag`: "For any
108  // object z of type complex<T>, reinterpret_cast<T(&)[2]>(z)[0] is the
109  // real part of z and reinterpret_cast<T(&)[2]>(z)[1] is the imaginary
110  // part of z." Now we don't have to provide a whole bunch of the overloads
111  // of things taking either Kokkos::complex or std::complex
112  : re_(reinterpret_cast<const RealType (&)[2]>(src)[0]),
113  im_(reinterpret_cast<const RealType (&)[2]>(src)[1]) {}
114 
120  // TODO: make explicit. DJS 2019-08-28
121  operator std::complex<RealType>() const noexcept {
122  return std::complex<RealType>(re_, im_);
123  }
124 
127  KOKKOS_INLINE_FUNCTION complex(const RealType& val) noexcept
128  : re_(val), im_(static_cast<RealType>(0)) {}
129 
131  KOKKOS_INLINE_FUNCTION
132  complex(const RealType& re, const RealType& im) noexcept : re_(re), im_(im) {}
133 
135  KOKKOS_INLINE_FUNCTION complex& operator=(const RealType& val) noexcept {
136  re_ = val;
137  im_ = RealType(0);
138  return *this;
139  }
140 
146  complex& operator=(const std::complex<RealType>& src) noexcept {
147  *this = complex(src);
148  return *this;
149  }
150 
152  KOKKOS_INLINE_FUNCTION
153  KOKKOS_CONSTEXPR_14 RealType& imag() noexcept { return im_; }
154 
156  KOKKOS_INLINE_FUNCTION
157  KOKKOS_CONSTEXPR_14 RealType& real() noexcept { return re_; }
158 
160  KOKKOS_INLINE_FUNCTION
161  constexpr RealType imag() const noexcept { return im_; }
162 
164  KOKKOS_INLINE_FUNCTION
165  constexpr RealType real() const noexcept { return re_; }
166 
168  KOKKOS_INLINE_FUNCTION
169  KOKKOS_CONSTEXPR_14
170  void imag(RealType v) noexcept { im_ = v; }
171 
173  KOKKOS_INLINE_FUNCTION
174  KOKKOS_CONSTEXPR_14
175  void real(RealType v) noexcept { re_ = v; }
176 
177  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator+=(
178  const complex<RealType>& src) noexcept {
179  re_ += src.re_;
180  im_ += src.im_;
181  return *this;
182  }
183 
184  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator+=(
185  const RealType& src) noexcept {
186  re_ += src;
187  return *this;
188  }
189 
190  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator-=(
191  const complex<RealType>& src) noexcept {
192  re_ -= src.re_;
193  im_ -= src.im_;
194  return *this;
195  }
196 
197  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator-=(
198  const RealType& src) noexcept {
199  re_ -= src;
200  return *this;
201  }
202 
203  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator*=(
204  const complex<RealType>& src) noexcept {
205  const RealType realPart = re_ * src.re_ - im_ * src.im_;
206  const RealType imagPart = re_ * src.im_ + im_ * src.re_;
207  re_ = realPart;
208  im_ = imagPart;
209  return *this;
210  }
211 
212  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator*=(
213  const RealType& src) noexcept {
214  re_ *= src;
215  im_ *= src;
216  return *this;
217  }
218 
219  // Conditional noexcept, just in case RType throws on divide-by-zero
220  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator/=(
221  const complex<RealType>& y) noexcept(noexcept(RealType{} / RealType{})) {
222  using Kokkos::Experimental::fabs;
223  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
224  // If the real part is +/-Inf and the imaginary part is -/+Inf,
225  // this won't change the result.
226  const RealType s = fabs(y.real()) + fabs(y.imag());
227 
228  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
229  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
230  // because y/s is NaN.
231  // TODO mark this branch unlikely
232  if (s == RealType(0)) {
233  this->re_ /= s;
234  this->im_ /= s;
235  } else {
236  const complex x_scaled(this->re_ / s, this->im_ / s);
237  const complex y_conj_scaled(y.re_ / s, -(y.im_) / s);
238  const RealType y_scaled_abs =
239  y_conj_scaled.re_ * y_conj_scaled.re_ +
240  y_conj_scaled.im_ * y_conj_scaled.im_; // abs(y) == abs(conj(y))
241  *this = x_scaled * y_conj_scaled;
242  *this /= y_scaled_abs;
243  }
244  return *this;
245  }
246 
247  KOKKOS_CONSTEXPR_14
248  KOKKOS_INLINE_FUNCTION complex& operator/=(
249  const std::complex<RealType>& y) noexcept(noexcept(RealType{} /
250  RealType{})) {
251  using Kokkos::Experimental::fabs;
252  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
253  // If the real part is +/-Inf and the imaginary part is -/+Inf,
254  // this won't change the result.
255  const RealType s = fabs(y.real()) + fabs(y.imag());
256 
257  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
258  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
259  // because y/s is NaN.
260  if (s == RealType(0)) {
261  this->re_ /= s;
262  this->im_ /= s;
263  } else {
264  const complex x_scaled(this->re_ / s, this->im_ / s);
265  const complex y_conj_scaled(y.re_ / s, -(y.im_) / s);
266  const RealType y_scaled_abs =
267  y_conj_scaled.re_ * y_conj_scaled.re_ +
268  y_conj_scaled.im_ * y_conj_scaled.im_; // abs(y) == abs(conj(y))
269  *this = x_scaled * y_conj_scaled;
270  *this /= y_scaled_abs;
271  }
272  return *this;
273  }
274 
275  KOKKOS_CONSTEXPR_14 KOKKOS_INLINE_FUNCTION complex& operator/=(
276  const RealType& src) noexcept(noexcept(RealType{} / RealType{})) {
277  re_ /= src;
278  im_ /= src;
279  return *this;
280  }
281 
282  //---------------------------------------------------------------------------
283  // TODO: refactor Kokkos reductions to remove dependency on
284  // volatile member overloads since they are being deprecated in c++20
285  //---------------------------------------------------------------------------
286 
288  template <class RType,
289  typename std::enable_if<std::is_convertible<RType, RealType>::value,
290  int>::type = 0>
291  KOKKOS_INLINE_FUNCTION complex(const volatile complex<RType>& src) noexcept
292  // Intentionally do the conversions implicitly here so that users don't
293  // get any warnings about narrowing, etc., that they would expect to get
294  // otherwise.
295  : re_(src.re_), im_(src.im_) {}
296 
306  //
307  // Templated, so as not to be a copy assignment operator (Kokkos issue #2577)
308  // Intended to behave as
309  // void operator=(const complex&) volatile noexcept
310  //
311  // Use cases:
312  // complex r;
313  // const complex cr;
314  // volatile complex vl;
315  // vl = r;
316  // vl = cr;
317  template <class Complex,
318  typename std::enable_if<std::is_same<Complex, complex>::value,
319  int>::type = 0>
320  KOKKOS_INLINE_FUNCTION void operator=(const Complex& src) volatile noexcept {
321  re_ = src.re_;
322  im_ = src.im_;
323  // We deliberately do not return anything here. See explanation
324  // in public documentation above.
325  }
326 
328  // TODO Should this return void like the other volatile assignment operators?
329  //
330  // Templated, so as not to be a copy assignment operator (Kokkos issue #2577)
331  // Intended to behave as
332  // volatile complex& operator=(const volatile complex&) volatile noexcept
333  //
334  // Use cases:
335  // volatile complex vr;
336  // const volatile complex cvr;
337  // volatile complex vl;
338  // vl = vr;
339  // vl = cvr;
340  template <class Complex,
341  typename std::enable_if<std::is_same<Complex, complex>::value,
342  int>::type = 0>
343  KOKKOS_INLINE_FUNCTION volatile complex& operator=(
344  const volatile Complex& src) volatile noexcept {
345  re_ = src.re_;
346  im_ = src.im_;
347  return *this;
348  }
349 
351  //
352  // Templated, so as not to be a copy assignment operator (Kokkos issue #2577)
353  // Intended to behave as
354  // complex& operator=(const volatile complex&) noexcept
355  //
356  // Use cases:
357  // volatile complex vr;
358  // const volatile complex cvr;
359  // complex l;
360  // l = vr;
361  // l = cvr;
362  //
363  template <class Complex,
364  typename std::enable_if<std::is_same<Complex, complex>::value,
365  int>::type = 0>
366  KOKKOS_INLINE_FUNCTION complex& operator=(
367  const volatile Complex& src) noexcept {
368  re_ = src.re_;
369  im_ = src.im_;
370  return *this;
371  }
372 
373  // Mirroring the behavior of the assignment operators from complex RHS in the
374  // RealType RHS versions.
375 
377  KOKKOS_INLINE_FUNCTION void operator=(const volatile RealType& val) noexcept {
378  re_ = val;
379  im_ = RealType(0);
380  // We deliberately do not return anything here. See explanation
381  // in public documentation above.
382  }
383 
385  KOKKOS_INLINE_FUNCTION complex& operator=(
386  const RealType& val) volatile noexcept {
387  re_ = val;
388  im_ = RealType(0);
389  return *this;
390  }
391 
393  // TODO Should this return void like the other volatile assignment operators?
394  KOKKOS_INLINE_FUNCTION complex& operator=(
395  const volatile RealType& val) volatile noexcept {
396  re_ = val;
397  im_ = RealType(0);
398  return *this;
399  }
400 
402  KOKKOS_INLINE_FUNCTION
403  volatile RealType& imag() volatile noexcept { return im_; }
404 
406  KOKKOS_INLINE_FUNCTION
407  volatile RealType& real() volatile noexcept { return re_; }
408 
410  KOKKOS_INLINE_FUNCTION
411  RealType imag() const volatile noexcept { return im_; }
412 
414  KOKKOS_INLINE_FUNCTION
415  RealType real() const volatile noexcept { return re_; }
416 
417  KOKKOS_INLINE_FUNCTION void operator+=(
418  const volatile complex<RealType>& src) volatile noexcept {
419  re_ += src.re_;
420  im_ += src.im_;
421  }
422 
423  KOKKOS_INLINE_FUNCTION void operator+=(
424  const volatile RealType& src) volatile noexcept {
425  re_ += src;
426  }
427 
428  KOKKOS_INLINE_FUNCTION void operator*=(
429  const volatile complex<RealType>& src) volatile noexcept {
430  const RealType realPart = re_ * src.re_ - im_ * src.im_;
431  const RealType imagPart = re_ * src.im_ + im_ * src.re_;
432 
433  re_ = realPart;
434  im_ = imagPart;
435  }
436 
437  KOKKOS_INLINE_FUNCTION void operator*=(
438  const volatile RealType& src) volatile noexcept {
439  re_ *= src;
440  im_ *= src;
441  }
442 
443  // TODO DSH 2019-10-7 why are there no volatile /= and friends?
444 };
445 
446 //==============================================================================
447 // <editor-fold desc="Equality and inequality"> {{{1
448 
449 // Note that this is not the same behavior as std::complex, which doesn't allow
450 // implicit conversions, but since this is the way we had it before, we have
451 // to do it this way now.
452 
454 template <class RealType1, class RealType2>
455 KOKKOS_INLINE_FUNCTION bool operator==(complex<RealType1> const& x,
456  complex<RealType2> const& y) noexcept {
457  using common_type = typename std::common_type<RealType1, RealType2>::type;
458  return common_type(x.real()) == common_type(y.real()) &&
459  common_type(x.imag()) == common_type(y.imag());
460 }
461 
462 // TODO (here and elsewhere) decide if we should convert to a Kokkos::complex
463 // and do the comparison in a device-marked function
465 template <class RealType1, class RealType2>
466 inline bool operator==(std::complex<RealType1> const& x,
467  complex<RealType2> const& y) noexcept {
468  using common_type = typename std::common_type<RealType1, RealType2>::type;
469  return common_type(x.real()) == common_type(y.real()) &&
470  common_type(x.imag()) == common_type(y.imag());
471 }
472 
474 template <class RealType1, class RealType2>
475 inline bool operator==(complex<RealType1> const& x,
476  std::complex<RealType2> const& y) noexcept {
477  using common_type = typename std::common_type<RealType1, RealType2>::type;
478  return common_type(x.real()) == common_type(y.real()) &&
479  common_type(x.imag()) == common_type(y.imag());
480 }
481 
483 template <
484  class RealType1, class RealType2,
485  // Constraints to avoid participation in oparator==() for every possible RHS
486  typename std::enable_if<std::is_convertible<RealType2, RealType1>::value,
487  int>::type = 0>
488 KOKKOS_INLINE_FUNCTION bool operator==(complex<RealType1> const& x,
489  RealType2 const& y) noexcept {
490  using common_type = typename std::common_type<RealType1, RealType2>::type;
491  return common_type(x.real()) == common_type(y) &&
492  common_type(x.imag()) == common_type(0);
493 }
494 
496 template <
497  class RealType1, class RealType2,
498  // Constraints to avoid participation in oparator==() for every possible RHS
499  typename std::enable_if<std::is_convertible<RealType1, RealType2>::value,
500  int>::type = 0>
501 KOKKOS_INLINE_FUNCTION bool operator==(RealType1 const& x,
502  complex<RealType2> const& y) noexcept {
503  using common_type = typename std::common_type<RealType1, RealType2>::type;
504  return common_type(x) == common_type(y.real()) &&
505  common_type(0) == common_type(y.imag());
506 }
507 
509 template <class RealType1, class RealType2>
510 KOKKOS_INLINE_FUNCTION bool operator!=(complex<RealType1> const& x,
511  complex<RealType2> const& y) noexcept {
512  using common_type = typename std::common_type<RealType1, RealType2>::type;
513  return common_type(x.real()) != common_type(y.real()) ||
514  common_type(x.imag()) != common_type(y.imag());
515 }
516 
518 template <class RealType1, class RealType2>
519 inline bool operator!=(std::complex<RealType1> const& x,
520  complex<RealType2> const& y) noexcept {
521  using common_type = typename std::common_type<RealType1, RealType2>::type;
522  return common_type(x.real()) != common_type(y.real()) ||
523  common_type(x.imag()) != common_type(y.imag());
524 }
525 
527 template <class RealType1, class RealType2>
528 inline bool operator!=(complex<RealType1> const& x,
529  std::complex<RealType2> const& y) noexcept {
530  using common_type = typename std::common_type<RealType1, RealType2>::type;
531  return common_type(x.real()) != common_type(y.real()) ||
532  common_type(x.imag()) != common_type(y.imag());
533 }
534 
536 template <
537  class RealType1, class RealType2,
538  // Constraints to avoid participation in oparator==() for every possible RHS
539  typename std::enable_if<std::is_convertible<RealType2, RealType1>::value,
540  int>::type = 0>
541 KOKKOS_INLINE_FUNCTION bool operator!=(complex<RealType1> const& x,
542  RealType2 const& y) noexcept {
543  using common_type = typename std::common_type<RealType1, RealType2>::type;
544  return common_type(x.real()) != common_type(y) ||
545  common_type(x.imag()) != common_type(0);
546 }
547 
549 template <
550  class RealType1, class RealType2,
551  // Constraints to avoid participation in oparator==() for every possible RHS
552  typename std::enable_if<std::is_convertible<RealType1, RealType2>::value,
553  int>::type = 0>
554 KOKKOS_INLINE_FUNCTION bool operator!=(RealType1 const& x,
555  complex<RealType2> const& y) noexcept {
556  using common_type = typename std::common_type<RealType1, RealType2>::type;
557  return common_type(x) != common_type(y.real()) ||
558  common_type(0) != common_type(y.imag());
559 }
560 
561 // </editor-fold> end Equality and inequality }}}1
562 //==============================================================================
563 
565 template <class RealType1, class RealType2>
566 KOKKOS_INLINE_FUNCTION
567  complex<typename std::common_type<RealType1, RealType2>::type>
568  operator+(const complex<RealType1>& x,
569  const complex<RealType2>& y) noexcept {
570  return complex<typename std::common_type<RealType1, RealType2>::type>(
571  x.real() + y.real(), x.imag() + y.imag());
572 }
573 
575 template <class RealType1, class RealType2>
576 KOKKOS_INLINE_FUNCTION
577  complex<typename std::common_type<RealType1, RealType2>::type>
578  operator+(const complex<RealType1>& x, const RealType2& y) noexcept {
579  return complex<typename std::common_type<RealType1, RealType2>::type>(
580  x.real() + y, x.imag());
581 }
582 
584 template <class RealType1, class RealType2>
585 KOKKOS_INLINE_FUNCTION
586  complex<typename std::common_type<RealType1, RealType2>::type>
587  operator+(const RealType1& x, const complex<RealType2>& y) noexcept {
588  return complex<typename std::common_type<RealType1, RealType2>::type>(
589  x + y.real(), y.imag());
590 }
591 
593 template <class RealType>
594 KOKKOS_INLINE_FUNCTION complex<RealType> operator+(
595  const complex<RealType>& x) noexcept {
596  return complex<RealType>{+x.real(), +x.imag()};
597 }
598 
600 template <class RealType1, class RealType2>
601 KOKKOS_INLINE_FUNCTION
602  complex<typename std::common_type<RealType1, RealType2>::type>
603  operator-(const complex<RealType1>& x,
604  const complex<RealType2>& y) noexcept {
605  return complex<typename std::common_type<RealType1, RealType2>::type>(
606  x.real() - y.real(), x.imag() - y.imag());
607 }
608 
610 template <class RealType1, class RealType2>
611 KOKKOS_INLINE_FUNCTION
612  complex<typename std::common_type<RealType1, RealType2>::type>
613  operator-(const complex<RealType1>& x, const RealType2& y) noexcept {
614  return complex<typename std::common_type<RealType1, RealType2>::type>(
615  x.real() - y, x.imag());
616 }
617 
619 template <class RealType1, class RealType2>
620 KOKKOS_INLINE_FUNCTION
621  complex<typename std::common_type<RealType1, RealType2>::type>
622  operator-(const RealType1& x, const complex<RealType2>& y) noexcept {
623  return complex<typename std::common_type<RealType1, RealType2>::type>(
624  x - y.real(), -y.imag());
625 }
626 
628 template <class RealType>
629 KOKKOS_INLINE_FUNCTION complex<RealType> operator-(
630  const complex<RealType>& x) noexcept {
631  return complex<RealType>(-x.real(), -x.imag());
632 }
633 
635 template <class RealType1, class RealType2>
636 KOKKOS_INLINE_FUNCTION
637  complex<typename std::common_type<RealType1, RealType2>::type>
638  operator*(const complex<RealType1>& x,
639  const complex<RealType2>& y) noexcept {
640  return complex<typename std::common_type<RealType1, RealType2>::type>(
641  x.real() * y.real() - x.imag() * y.imag(),
642  x.real() * y.imag() + x.imag() * y.real());
643 }
644 
653 template <class RealType1, class RealType2>
654 inline complex<typename std::common_type<RealType1, RealType2>::type> operator*(
655  const std::complex<RealType1>& x, const complex<RealType2>& y) {
656  return complex<typename std::common_type<RealType1, RealType2>::type>(
657  x.real() * y.real() - x.imag() * y.imag(),
658  x.real() * y.imag() + x.imag() * y.real());
659 }
660 
665 template <class RealType1, class RealType2>
666 KOKKOS_INLINE_FUNCTION
667  complex<typename std::common_type<RealType1, RealType2>::type>
668  operator*(const RealType1& x, const complex<RealType2>& y) noexcept {
669  return complex<typename std::common_type<RealType1, RealType2>::type>(
670  x * y.real(), x * y.imag());
671 }
672 
677 template <class RealType1, class RealType2>
678 KOKKOS_INLINE_FUNCTION
679  complex<typename std::common_type<RealType1, RealType2>::type>
680  operator*(const complex<RealType1>& y, const RealType2& x) noexcept {
681  return complex<typename std::common_type<RealType1, RealType2>::type>(
682  x * y.real(), x * y.imag());
683 }
684 
686 template <class RealType>
687 KOKKOS_INLINE_FUNCTION RealType imag(const complex<RealType>& x) noexcept {
688  return x.imag();
689 }
690 
692 template <class RealType>
693 KOKKOS_INLINE_FUNCTION RealType real(const complex<RealType>& x) noexcept {
694  return x.real();
695 }
696 
698 template <class T>
699 KOKKOS_INLINE_FUNCTION complex<T> polar(const T& r, const T& theta = T()) {
700  using Kokkos::Experimental::cos;
701  using Kokkos::Experimental::sin;
702  KOKKOS_EXPECTS(r >= 0);
703  return complex<T>(r * cos(theta), r * sin(theta));
704 }
705 
707 template <class RealType>
708 KOKKOS_INLINE_FUNCTION RealType abs(const complex<RealType>& x) {
709  using Kokkos::Experimental::hypot;
710  return hypot(x.real(), x.imag());
711 }
712 
714 template <class T>
715 KOKKOS_INLINE_FUNCTION complex<T> pow(const complex<T>& x, const T& y) {
716  using Kokkos::Experimental::atan2;
717  using Kokkos::Experimental::pow;
718  T r = abs(x);
719  T theta = atan2(x.imag(), x.real());
720  return polar(pow(r, y), y * theta);
721 }
722 
723 template <class T>
724 KOKKOS_INLINE_FUNCTION complex<T> pow(const T& x, const complex<T>& y) {
725  return pow(complex<T>(x), y);
726 }
727 
728 template <class T>
729 KOKKOS_INLINE_FUNCTION complex<T> pow(const complex<T>& x,
730  const complex<T>& y) {
731  using Kokkos::Experimental::log;
732 
733  return x == T() ? T() : exp(y * log(x));
734 }
735 
736 namespace Impl {
737 // NOTE promote would also be useful for math functions
738 template <class T, bool = std::is_integral<T>::value>
739 struct promote {
740  using type = double;
741 };
742 template <class T>
743 struct promote<T, false> {};
744 template <>
745 struct promote<long double> {
746  using type = long double;
747 };
748 template <>
749 struct promote<double> {
750  using type = double;
751 };
752 template <>
753 struct promote<float> {
754  using type = float;
755 };
756 template <class T>
757 using promote_t = typename promote<T>::type;
758 template <class T, class U>
759 struct promote_2 {
760  using type = decltype(promote_t<T>() + promote_t<U>());
761 };
762 template <class T, class U>
763 using promote_2_t = typename promote_2<T, U>::type;
764 } // namespace Impl
765 
766 template <class T, class U,
767  class = std::enable_if_t<std::is_arithmetic<T>::value>>
768 KOKKOS_INLINE_FUNCTION complex<Impl::promote_2_t<T, U>> pow(
769  const T& x, const complex<U>& y) {
770  using type = Impl::promote_2_t<T, U>;
771  return pow(type(x), complex<type>(y));
772 }
773 
774 template <class T, class U,
775  class = std::enable_if_t<std::is_arithmetic<U>::value>>
776 KOKKOS_INLINE_FUNCTION complex<Impl::promote_2_t<T, U>> pow(const complex<T>& x,
777  const U& y) {
778  using type = Impl::promote_2_t<T, U>;
779  return pow(complex<type>(x), type(y));
780 }
781 
782 template <class T, class U>
783 KOKKOS_INLINE_FUNCTION complex<Impl::promote_2_t<T, U>> pow(
784  const complex<T>& x, const complex<U>& y) {
785  using type = Impl::promote_2_t<T, U>;
786  return pow(complex<type>(x), complex<type>(y));
787 }
788 
791 template <class RealType>
792 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> sqrt(
793  const complex<RealType>& x) {
794  using Kokkos::Experimental::fabs;
795  using Kokkos::Experimental::sqrt;
796 
797  RealType r = x.real();
798  RealType i = x.imag();
799 
800  if (r == RealType()) {
801  RealType t = sqrt(fabs(i) / 2);
802  return Kokkos::complex<RealType>(t, i < RealType() ? -t : t);
803  } else {
804  RealType t = sqrt(2 * (abs(x) + fabs(r)));
805  RealType u = t / 2;
806  return r > RealType() ? Kokkos::complex<RealType>(u, i / t)
807  : Kokkos::complex<RealType>(fabs(i) / t,
808  i < RealType() ? -u : u);
809  }
810 }
811 
813 template <class RealType>
814 KOKKOS_INLINE_FUNCTION complex<RealType> conj(
815  const complex<RealType>& x) noexcept {
816  return complex<RealType>(real(x), -imag(x));
817 }
818 
820 template <class RealType>
821 KOKKOS_INLINE_FUNCTION complex<RealType> exp(const complex<RealType>& x) {
822  using Kokkos::Experimental::cos;
823  using Kokkos::Experimental::exp;
824  using Kokkos::Experimental::sin;
825  return exp(x.real()) * complex<RealType>(cos(x.imag()), sin(x.imag()));
826 }
827 
829 template <class RealType>
830 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> log(
831  const complex<RealType>& x) {
832  using Kokkos::Experimental::atan2;
833  using Kokkos::Experimental::log;
834  RealType phi = atan2(x.imag(), x.real());
835  return Kokkos::complex<RealType>(log(abs(x)), phi);
836 }
837 
839 template <class RealType>
840 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> sin(
841  const complex<RealType>& x) {
842  using Kokkos::Experimental::cos;
843  using Kokkos::Experimental::cosh;
844  using Kokkos::Experimental::sin;
845  using Kokkos::Experimental::sinh;
846  return Kokkos::complex<RealType>(sin(x.real()) * cosh(x.imag()),
847  cos(x.real()) * sinh(x.imag()));
848 }
849 
851 template <class RealType>
852 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> cos(
853  const complex<RealType>& x) {
854  using Kokkos::Experimental::cos;
855  using Kokkos::Experimental::cosh;
856  using Kokkos::Experimental::sin;
857  using Kokkos::Experimental::sinh;
858  return Kokkos::complex<RealType>(cos(x.real()) * cosh(x.imag()),
859  -sin(x.real()) * sinh(x.imag()));
860 }
861 
863 template <class RealType>
864 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> tan(
865  const complex<RealType>& x) {
866  return sin(x) / cos(x);
867 }
868 
870 template <class RealType>
871 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> sinh(
872  const complex<RealType>& x) {
873  using Kokkos::Experimental::cos;
874  using Kokkos::Experimental::cosh;
875  using Kokkos::Experimental::sin;
876  using Kokkos::Experimental::sinh;
877  return Kokkos::complex<RealType>(sinh(x.real()) * cos(x.imag()),
878  cosh(x.real()) * sin(x.imag()));
879 }
880 
882 template <class RealType>
883 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> cosh(
884  const complex<RealType>& x) {
885  using Kokkos::Experimental::cos;
886  using Kokkos::Experimental::cosh;
887  using Kokkos::Experimental::sin;
888  using Kokkos::Experimental::sinh;
889  return Kokkos::complex<RealType>(cosh(x.real()) * cos(x.imag()),
890  sinh(x.real()) * sin(x.imag()));
891 }
892 
894 template <class RealType>
895 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> tanh(
896  const complex<RealType>& x) {
897  return sinh(x) / cosh(x);
898 }
899 
901 template <class RealType>
902 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> asinh(
903  const complex<RealType>& x) {
904  return log(x + sqrt(x * x + RealType(1.0)));
905 }
906 
908 template <class RealType>
909 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> acosh(
910  const complex<RealType>& x) {
911  return RealType(2.0) * log(sqrt(RealType(0.5) * (x + RealType(1.0))) +
912  sqrt(RealType(0.5) * (x - RealType(1.0))));
913 }
914 
916 template <class RealType>
917 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> atanh(
918  const complex<RealType>& x) {
919  using Kokkos::Experimental::atan2;
920  using Kokkos::Experimental::log;
921 
922  const RealType i2 = x.imag() * x.imag();
923  const RealType r = RealType(1.0) - i2 - x.real() * x.real();
924 
925  RealType p = RealType(1.0) + x.real();
926  RealType m = RealType(1.0) - x.real();
927 
928  p = i2 + p * p;
929  m = i2 + m * m;
930 
931  RealType phi = atan2(RealType(2.0) * x.imag(), r);
932  return Kokkos::complex<RealType>(RealType(0.25) * (log(p) - log(m)),
933  RealType(0.5) * phi);
934 }
935 
937 template <class RealType>
938 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> asin(
939  const complex<RealType>& x) {
941  asinh(Kokkos::complex<RealType>(-x.imag(), x.real()));
942  return Kokkos::complex<RealType>(t.imag(), -t.real());
943 }
944 
946 template <class RealType>
947 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> acos(
948  const complex<RealType>& x) {
949  using Kokkos::Experimental::acos;
950  Kokkos::complex<RealType> t = asin(x);
951  RealType pi_2 = acos(RealType(0.0));
952  return Kokkos::complex<RealType>(pi_2 - t.real(), -t.imag());
953 }
954 
956 template <class RealType>
957 KOKKOS_INLINE_FUNCTION Kokkos::complex<RealType> atan(
958  const complex<RealType>& x) {
959  using Kokkos::Experimental::atan2;
960  using Kokkos::Experimental::log;
961  const RealType r2 = x.real() * x.real();
962  const RealType i = RealType(1.0) - r2 - x.imag() * x.imag();
963 
964  RealType p = x.imag() + RealType(1.0);
965  RealType m = x.imag() - RealType(1.0);
966 
967  p = r2 + p * p;
968  m = r2 + m * m;
969 
971  RealType(0.5) * atan2(RealType(2.0) * x.real(), i),
972  RealType(0.25) * log(p / m));
973 }
974 
978 template <class RealType>
979 inline complex<RealType> exp(const std::complex<RealType>& c) {
980  return complex<RealType>(std::exp(c.real()) * std::cos(c.imag()),
981  std::exp(c.real()) * std::sin(c.imag()));
982 }
983 
985 template <class RealType1, class RealType2>
986 KOKKOS_INLINE_FUNCTION
987  complex<typename std::common_type<RealType1, RealType2>::type>
988  operator/(const complex<RealType1>& x,
989  const RealType2& y) noexcept(noexcept(RealType1{} /
990  RealType2{})) {
991  return complex<typename std::common_type<RealType1, RealType2>::type>(
992  real(x) / y, imag(x) / y);
993 }
994 
996 template <class RealType1, class RealType2>
997 KOKKOS_INLINE_FUNCTION
998  complex<typename std::common_type<RealType1, RealType2>::type>
999  operator/(const complex<RealType1>& x,
1000  const complex<RealType2>& y) noexcept(noexcept(RealType1{} /
1001  RealType2{})) {
1002  using Kokkos::Experimental::fabs;
1003  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
1004  // If the real part is +/-Inf and the imaginary part is -/+Inf,
1005  // this won't change the result.
1006  using common_real_type =
1007  typename std::common_type<RealType1, RealType2>::type;
1008  const common_real_type s = fabs(real(y)) + fabs(imag(y));
1009 
1010  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
1011  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
1012  // because y/s is NaN.
1013  if (s == 0.0) {
1014  return complex<common_real_type>(real(x) / s, imag(x) / s);
1015  } else {
1016  const complex<common_real_type> x_scaled(real(x) / s, imag(x) / s);
1017  const complex<common_real_type> y_conj_scaled(real(y) / s, -imag(y) / s);
1018  const RealType1 y_scaled_abs =
1019  real(y_conj_scaled) * real(y_conj_scaled) +
1020  imag(y_conj_scaled) * imag(y_conj_scaled); // abs(y) == abs(conj(y))
1021  complex<common_real_type> result = x_scaled * y_conj_scaled;
1022  result /= y_scaled_abs;
1023  return result;
1024  }
1025 }
1026 
1028 template <class RealType1, class RealType2>
1029 KOKKOS_INLINE_FUNCTION
1030  complex<typename std::common_type<RealType1, RealType2>::type>
1031  operator/(const RealType1& x,
1032  const complex<RealType2>& y) noexcept(noexcept(RealType1{} /
1033  RealType2{})) {
1034  return complex<typename std::common_type<RealType1, RealType2>::type>(x) / y;
1035 }
1036 
1037 template <class RealType>
1038 std::ostream& operator<<(std::ostream& os, const complex<RealType>& x) {
1039  const std::complex<RealType> x_std(Kokkos::real(x), Kokkos::imag(x));
1040  os << x_std;
1041  return os;
1042 }
1043 
1044 template <class RealType>
1045 std::istream& operator>>(std::istream& is, complex<RealType>& x) {
1046  std::complex<RealType> x_std;
1047  is >> x_std;
1048  x = x_std; // only assigns on success of above
1049  return is;
1050 }
1051 
1052 template <class T>
1053 struct reduction_identity<Kokkos::complex<T>> {
1054  using t_red_ident = reduction_identity<T>;
1055  KOKKOS_FORCEINLINE_FUNCTION constexpr static Kokkos::complex<T>
1056  sum() noexcept {
1057  return Kokkos::complex<T>(t_red_ident::sum(), t_red_ident::sum());
1058  }
1059  KOKKOS_FORCEINLINE_FUNCTION constexpr static Kokkos::complex<T>
1060  prod() noexcept {
1061  return Kokkos::complex<T>(t_red_ident::prod(), t_red_ident::sum());
1062  }
1063 };
1064 
1065 } // namespace Kokkos
1066 
1067 #endif // KOKKOS_COMPLEX_HPP
KOKKOS_INLINE_FUNCTION KOKKOS_CONSTEXPR_14 RealType & real() noexcept
The real part of this complex number.
KOKKOS_INLINE_FUNCTION KOKKOS_CONSTEXPR_14 RealType & imag() noexcept
The imaginary part of this complex number.
Partial reimplementation of std::complex that works as the result of a Kokkos::parallel_reduce.
KOKKOS_INLINE_FUNCTION constexpr RealType real() const noexcept
The real part of this complex number.
KOKKOS_INLINE_FUNCTION complex & operator=(const RealType &val) noexcept
Assignment operator (from a real number).
KOKKOS_INLINE_FUNCTION complex(const RealType &val) noexcept
Constructor that takes just the real part, and sets the imaginary part to zero.
complex & operator=(const std::complex< RealType > &src) noexcept
Assignment operator from std::complex.
KOKKOS_INLINE_FUNCTION volatile RealType & imag() volatile noexcept
The imaginary part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION complex(const RealType &re, const RealType &im) noexcept
Constructor that takes the real and imaginary parts.
KOKKOS_INLINE_FUNCTION constexpr RealType imag() const noexcept
The imaginary part of this complex number.
KOKKOS_INLINE_FUNCTION RealType imag() const volatile noexcept
The imaginary part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION void operator=(const Complex &src) volatile noexcept
Assignment operator, for volatile *this and nonvolatile input.
KOKKOS_INLINE_FUNCTION complex(const std::complex< RealType > &src) noexcept
Conversion constructor from std::complex.
KOKKOS_INLINE_FUNCTION KOKKOS_CONSTEXPR_14 void imag(RealType v) noexcept
Set the imaginary part of this complex number.
KOKKOS_INLINE_FUNCTION volatile complex & operator=(const volatile Complex &src) volatile noexcept
Assignment operator, volatile LHS and volatile RHS.
KOKKOS_INLINE_FUNCTION RealType real() const volatile noexcept
The real part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION complex & operator=(const volatile RealType &val) volatile noexcept
Assignment operator volatile LHS and volatile RHS.
Atomic functions.
KOKKOS_INLINE_FUNCTION complex(const volatile complex< RType > &src) noexcept
Copy constructor from volatile.
RealType value_type
The type of the real or imaginary parts of this complex number.
KOKKOS_INLINE_FUNCTION volatile RealType & real() volatile noexcept
The real part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION void operator=(const volatile RealType &val) noexcept
Assignment operator (from a volatile real number).
Definition: dummy.cpp:3
KOKKOS_INLINE_FUNCTION KOKKOS_CONSTEXPR_14 void real(RealType v) noexcept
Set the real part of this complex number.
KOKKOS_INLINE_FUNCTION complex & operator=(const volatile Complex &src) noexcept
Assignment operator, volatile RHS and non-volatile LHS.
KOKKOS_INLINE_FUNCTION complex & operator=(const RealType &val) volatile noexcept
Assignment operator volatile LHS and non-volatile RHS.