Kokkos Core Kernels Package  Version of the Day
Kokkos_Complex.hpp
1 /*
2 //@HEADER
3 // ************************************************************************
4 //
5 // Kokkos v. 2.0
6 // Copyright (2014) Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact H. Carter Edwards (hcedwar@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 */
43 #ifndef KOKKOS_COMPLEX_HPP
44 #define KOKKOS_COMPLEX_HPP
45 
46 #include <Kokkos_Atomic.hpp>
47 #include <complex>
48 #include <iostream>
49 
50 namespace Kokkos {
51 
59 template<class RealType>
60 class complex {
61 private:
62  RealType re_, im_;
63 
64 public:
66  typedef RealType value_type;
67 
69  KOKKOS_INLINE_FUNCTION complex () :
70  re_ (0.0), im_ (0.0)
71  {}
72 
74  KOKKOS_INLINE_FUNCTION complex (const complex<RealType>& src) :
75  re_ (src.re_), im_ (src.im_)
76  {}
77 
79  KOKKOS_INLINE_FUNCTION complex (const volatile complex<RealType>& src) :
80  re_ (src.re_), im_ (src.im_)
81  {}
82 
88  template<class InputRealType>
89  complex (const std::complex<InputRealType>& src) :
90  re_ (std::real (src)), im_ (std::imag (src))
91  {}
92 
98  operator std::complex<RealType> () const {
99  return std::complex<RealType> (re_, im_);
100  }
101 
104  template<class InputRealType>
105  KOKKOS_INLINE_FUNCTION complex (const InputRealType& val) :
106  re_ (val), im_ (0.0)
107  {}
108 
110  template<class RealType1, class RealType2>
111  KOKKOS_INLINE_FUNCTION complex (const RealType1& re, const RealType2& im) :
112  re_ (re), im_ (im)
113  {}
114 
116  template<class InputRealType>
117  KOKKOS_INLINE_FUNCTION
119  re_ = src.re_;
120  im_ = src.im_;
121  return *this;
122  }
123 
125  template<class InputRealType>
126  KOKKOS_INLINE_FUNCTION
127  volatile complex<RealType>& operator= (const complex<InputRealType>& src) volatile {
128  re_ = src.re_;
129  im_ = src.im_;
130  return *this;
131  }
132 
134  template<class InputRealType>
135  KOKKOS_INLINE_FUNCTION
136  volatile complex<RealType>& operator= (const volatile complex<InputRealType>& src) volatile {
137  re_ = src.re_;
138  im_ = src.im_;
139  return *this;
140  }
141 
143  template<class InputRealType>
144  KOKKOS_INLINE_FUNCTION
146  re_ = src.re_;
147  im_ = src.im_;
148  return *this;
149  }
150 
152  template<class InputRealType>
153  KOKKOS_INLINE_FUNCTION
154  complex<RealType>& operator= (const InputRealType& val) {
155  re_ = val;
156  im_ = static_cast<RealType> (0.0);
157  return *this;
158  }
159 
161  template<class InputRealType>
162  KOKKOS_INLINE_FUNCTION
163  void operator= (const InputRealType& val) volatile {
164  re_ = val;
165  im_ = static_cast<RealType> (0.0);
166  }
167 
173  template<class InputRealType>
174  complex<RealType>& operator= (const std::complex<InputRealType>& src) {
175  re_ = std::real (src);
176  im_ = std::imag (src);
177  return *this;
178  }
179 
181  KOKKOS_INLINE_FUNCTION RealType& imag () {
182  return im_;
183  }
184 
186  KOKKOS_INLINE_FUNCTION RealType& real () {
187  return re_;
188  }
189 
191  KOKKOS_INLINE_FUNCTION const RealType imag () const {
192  return im_;
193  }
194 
196  KOKKOS_INLINE_FUNCTION const RealType real () const {
197  return re_;
198  }
199 
201  KOKKOS_INLINE_FUNCTION volatile RealType& imag () volatile {
202  return im_;
203  }
204 
206  KOKKOS_INLINE_FUNCTION volatile RealType& real () volatile {
207  return re_;
208  }
209 
211  KOKKOS_INLINE_FUNCTION const RealType imag () const volatile {
212  return im_;
213  }
214 
216  KOKKOS_INLINE_FUNCTION const RealType real () const volatile {
217  return re_;
218  }
219 
220  KOKKOS_INLINE_FUNCTION
221  complex<RealType>& operator += (const complex<RealType>& src) {
222  re_ += src.re_;
223  im_ += src.im_;
224  return *this;
225  }
226 
227  KOKKOS_INLINE_FUNCTION
228  void operator += (const volatile complex<RealType>& src) volatile {
229  re_ += src.re_;
230  im_ += src.im_;
231  }
232 
233  KOKKOS_INLINE_FUNCTION
234  complex<RealType>& operator += (const RealType& src) {
235  re_ += src;
236  return *this;
237  }
238 
239  KOKKOS_INLINE_FUNCTION
240  void operator += (const volatile RealType& src) volatile {
241  re_ += src;
242  }
243 
244  KOKKOS_INLINE_FUNCTION
245  complex<RealType>& operator -= (const complex<RealType>& src) {
246  re_ -= src.re_;
247  im_ -= src.im_;
248  return *this;
249  }
250 
251  KOKKOS_INLINE_FUNCTION
252  complex<RealType>& operator -= (const RealType& src) {
253  re_ -= src;
254  return *this;
255  }
256 
257  KOKKOS_INLINE_FUNCTION
258  complex<RealType>& operator *= (const complex<RealType>& src) {
259  const RealType realPart = re_ * src.re_ - im_ * src.im_;
260  const RealType imagPart = re_ * src.im_ + im_ * src.re_;
261  re_ = realPart;
262  im_ = imagPart;
263  return *this;
264  }
265 
266  KOKKOS_INLINE_FUNCTION
267  void operator *= (const volatile complex<RealType>& src) volatile {
268  const RealType realPart = re_ * src.re_ - im_ * src.im_;
269  const RealType imagPart = re_ * src.im_ + im_ * src.re_;
270  re_ = realPart;
271  im_ = imagPart;
272  }
273 
274  KOKKOS_INLINE_FUNCTION
275  complex<RealType>& operator *= (const RealType& src) {
276  re_ *= src;
277  im_ *= src;
278  return *this;
279  }
280 
281  KOKKOS_INLINE_FUNCTION
282  void operator *= (const volatile RealType& src) volatile {
283  re_ *= src;
284  im_ *= src;
285  }
286 
287  KOKKOS_INLINE_FUNCTION
288  complex<RealType>& operator /= (const complex<RealType>& y) {
289  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
290  // If the real part is +/-Inf and the imaginary part is -/+Inf,
291  // this won't change the result.
292  const RealType s = ::fabs (y.real ()) + ::fabs (y.imag ());
293 
294  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
295  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
296  // because y/s is NaN.
297  if (s == 0.0) {
298  this->re_ /= s;
299  this->im_ /= s;
300  }
301  else {
302  const complex<RealType> x_scaled (this->re_ / s, this->im_ / s);
303  const complex<RealType> y_conj_scaled (y.re_ / s, -(y.im_) / s);
304  const RealType y_scaled_abs = y_conj_scaled.re_ * y_conj_scaled.re_ +
305  y_conj_scaled.im_ * y_conj_scaled.im_; // abs(y) == abs(conj(y))
306  *this = x_scaled * y_conj_scaled;
307  *this /= y_scaled_abs;
308  }
309  return *this;
310  }
311 
312  KOKKOS_INLINE_FUNCTION
313  complex<RealType>& operator /= (const RealType& src) {
314  re_ /= src;
315  im_ /= src;
316  return *this;
317  }
318 };
319 
321 template<class RealType>
322 KOKKOS_INLINE_FUNCTION
325  return complex<RealType> (x.real () + y.real (), x.imag () + y.imag ());
326 }
327 
329 template<class RealType>
330 KOKKOS_INLINE_FUNCTION
333  return x;
334 }
335 
337 template<class RealType>
338 KOKKOS_INLINE_FUNCTION
341  return complex<RealType> (x.real () - y.real (), x.imag () - y.imag ());
342 }
343 
345 template<class RealType>
346 KOKKOS_INLINE_FUNCTION
349  return complex<RealType> (-x.real (), -x.imag ());
350 }
351 
353 template<class RealType>
354 KOKKOS_INLINE_FUNCTION
357  return complex<RealType> (x.real () * y.real () - x.imag () * y.imag (),
358  x.real () * y.imag () + x.imag () * y.real ());
359 }
360 
371 template<class RealType>
373 operator * (const std::complex<RealType>& x, const complex<RealType>& y) {
374  return complex<RealType> (x.real () * y.real () - x.imag () * y.imag (),
375  x.real () * y.imag () + x.imag () * y.real ());
376 }
377 
382 template<class RealType>
383 KOKKOS_INLINE_FUNCTION
385 operator * (const RealType& x, const complex<RealType>& y) {
386  return complex<RealType> (x * y.real (), x * y.imag ());
387 }
388 
389 
391 template<class RealType>
392 KOKKOS_INLINE_FUNCTION
393 RealType imag (const complex<RealType>& x) {
394  return x.imag ();
395 }
396 
398 template<class RealType>
399 KOKKOS_INLINE_FUNCTION
400 RealType real (const complex<RealType>& x) {
401  return x.real ();
402 }
403 
405 template<class RealType>
406 KOKKOS_INLINE_FUNCTION
407 RealType abs (const complex<RealType>& x) {
408  // FIXME (mfh 31 Oct 2014) Scale to avoid unwarranted overflow.
409  return ::sqrt (real (x) * real (x) + imag (x) * imag (x));
410 }
411 
413 template<class RealType>
414 KOKKOS_INLINE_FUNCTION
416  return complex<RealType> (real (x), -imag (x));
417 }
418 
419 
421 template<class RealType1, class RealType2>
422 KOKKOS_INLINE_FUNCTION
424 operator / (const complex<RealType1>& x, const RealType2& y) {
425  return complex<RealType1> (real (x) / y, imag (x) / y);
426 }
427 
429 template<class RealType>
430 KOKKOS_INLINE_FUNCTION
433  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
434  // If the real part is +/-Inf and the imaginary part is -/+Inf,
435  // this won't change the result.
436  const RealType s = ::fabs (real (y)) + ::fabs (imag (y));
437 
438  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
439  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
440  // because y/s is NaN.
441  if (s == 0.0) {
442  return complex<RealType> (real (x) / s, imag (x) / s);
443  }
444  else {
445  const complex<RealType> x_scaled (real (x) / s, imag (x) / s);
446  const complex<RealType> y_conj_scaled (real (y) / s, -imag (y) / s);
447  const RealType y_scaled_abs = real (y_conj_scaled) * real (y_conj_scaled) +
448  imag (y_conj_scaled) * imag (y_conj_scaled); // abs(y) == abs(conj(y))
449  complex<RealType> result = x_scaled * y_conj_scaled;
450  result /= y_scaled_abs;
451  return result;
452  }
453 }
454 
456 template<class RealType>
457 KOKKOS_INLINE_FUNCTION
459  return real (x) == real (y) && imag (x) == imag (y);
460 }
461 
463 template<class RealType>
464 KOKKOS_INLINE_FUNCTION
465 bool operator == (const std::complex<RealType>& x, const complex<RealType>& y) {
466  return std::real (x) == real (y) && std::imag (x) == imag (y);
467 }
468 
470 template<class RealType1, class RealType2>
471 KOKKOS_INLINE_FUNCTION
472 bool operator == (const complex<RealType1>& x, const RealType2& y) {
473  return real (x) == y && imag (x) == static_cast<RealType1> (0.0);
474 }
475 
477 template<class RealType>
478 KOKKOS_INLINE_FUNCTION
479 bool operator == (const RealType& x, const complex<RealType>& y) {
480  return y == x;
481 }
482 
484 template<class RealType>
485 KOKKOS_INLINE_FUNCTION
487  return real (x) != real (y) || imag (x) != imag (y);
488 }
489 
491 template<class RealType>
492 KOKKOS_INLINE_FUNCTION
493 bool operator != (const std::complex<RealType>& x, const complex<RealType>& y) {
494  return std::real (x) != real (y) || std::imag (x) != imag (y);
495 }
496 
498 template<class RealType1, class RealType2>
499 KOKKOS_INLINE_FUNCTION
500 bool operator != (const complex<RealType1>& x, const RealType2& y) {
501  return real (x) != y || imag (x) != static_cast<RealType1> (0.0);
502 }
503 
505 template<class RealType>
506 KOKKOS_INLINE_FUNCTION
507 bool operator != (const RealType& x, const complex<RealType>& y) {
508  return y != x;
509 }
510 
511 template<class RealType>
512 std::ostream& operator << (std::ostream& os, const complex<RealType>& x) {
513  const std::complex<RealType> x_std (Kokkos::real (x), Kokkos::imag (x));
514  os << x_std;
515  return os;
516 }
517 
518 template<class RealType>
519 std::ostream& operator >> (std::ostream& os, complex<RealType>& x) {
520  std::complex<RealType> x_std;
521  os >> x_std;
522  x = x_std; // only assigns on success of above
523  return os;
524 }
525 
526 
527 } // namespace Kokkos
528 
529 #endif // KOKKOS_COMPLEX_HPP
KOKKOS_INLINE_FUNCTION RealType & imag()
The imaginary part of this complex number.
KOKKOS_INLINE_FUNCTION complex(const RealType1 &re, const RealType2 &im)
Constructor that takes the real and imaginary parts.
KOKKOS_INLINE_FUNCTION complex< RealType > operator*(const complex< RealType > &x, const complex< RealType > &y)
Binary * operator for complex.
Partial reimplementation of std::complex that works as the result of a Kokkos::parallel_reduce.
KOKKOS_INLINE_FUNCTION complex(const InputRealType &val)
Constructor that takes just the real part, and sets the imaginary part to zero.
KOKKOS_INLINE_FUNCTION bool operator!=(const complex< RealType > &x, const complex< RealType > &y)
Inequality operator for two complex numbers.
KOKKOS_INLINE_FUNCTION complex(const volatile complex< RealType > &src)
Copy constructor from volatile.
KOKKOS_INLINE_FUNCTION const RealType real() const
The real part of this complex number.
KOKKOS_INLINE_FUNCTION complex(const complex< RealType > &src)
Copy constructor.
KOKKOS_INLINE_FUNCTION complex()
Default constructor (initializes both real and imaginary parts to zero).
KOKKOS_INLINE_FUNCTION RealType abs(const complex< RealType > &x)
Absolute value (magnitude) of a complex number.
KOKKOS_INLINE_FUNCTION RealType & real()
The real part of this complex number.
complex(const std::complex< InputRealType > &src)
Conversion constructor from std::complex.
KOKKOS_INLINE_FUNCTION volatile RealType & imag() volatile
The imaginary part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION const RealType imag() const
The imaginary part of this complex number.
KOKKOS_INLINE_FUNCTION complex< RealType > & operator=(const complex< InputRealType > &src)
Assignment operator.
KOKKOS_INLINE_FUNCTION complex< RealType > operator+(const complex< RealType > &x, const complex< RealType > &y)
Binary + operator for complex.
KOKKOS_INLINE_FUNCTION complex< RealType > operator-(const complex< RealType > &x, const complex< RealType > &y)
Binary - operator for complex.
RealType value_type
The type of the real or imaginary parts of this complex number.
KOKKOS_INLINE_FUNCTION RealType real(const complex< RealType > &x)
Real part of a complex number.
KOKKOS_INLINE_FUNCTION volatile RealType & real() volatile
The real part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION const RealType imag() const volatile
The imaginary part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION const RealType real() const volatile
The real part of this complex number (volatile overload).
KOKKOS_INLINE_FUNCTION complex< RealType > conj(const complex< RealType > &x)
Conjugate of a complex number.
Atomic functions.
KOKKOS_INLINE_FUNCTION bool operator==(const complex< RealType > &x, const complex< RealType > &y)
Equality operator for two complex numbers.
KOKKOS_INLINE_FUNCTION complex< RealType1 > operator/(const complex< RealType1 > &x, const RealType2 &y)
Binary operator / for complex and real numbers.
KOKKOS_INLINE_FUNCTION RealType imag(const complex< RealType > &x)
Imaginary part of a complex number.