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 
133  template<class InputRealType>
134  KOKKOS_INLINE_FUNCTION
135  void operator= (const complex<InputRealType>& src) volatile {
136  re_ = src.re_;
137  im_ = src.im_;
138  // We deliberately do not return anything here. See explanation
139  // in public documentation above.
140  }
141 
143  template<class InputRealType>
144  KOKKOS_INLINE_FUNCTION
145  volatile complex<RealType>& operator= (const volatile complex<InputRealType>& src) volatile {
146  re_ = src.re_;
147  im_ = src.im_;
148  return *this;
149  }
150 
152  template<class InputRealType>
153  KOKKOS_INLINE_FUNCTION
155  re_ = src.re_;
156  im_ = src.im_;
157  return *this;
158  }
159 
161  template<class InputRealType>
162  KOKKOS_INLINE_FUNCTION
163  complex<RealType>& operator= (const InputRealType& val) {
164  re_ = val;
165  im_ = static_cast<RealType> (0.0);
166  return *this;
167  }
168 
170  template<class InputRealType>
171  KOKKOS_INLINE_FUNCTION
172  void operator= (const InputRealType& val) volatile {
173  re_ = val;
174  im_ = static_cast<RealType> (0.0);
175  }
176 
182  template<class InputRealType>
183  complex<RealType>& operator= (const std::complex<InputRealType>& src) {
184  re_ = std::real (src);
185  im_ = std::imag (src);
186  return *this;
187  }
188 
190  KOKKOS_INLINE_FUNCTION RealType& imag () {
191  return im_;
192  }
193 
195  KOKKOS_INLINE_FUNCTION RealType& real () {
196  return re_;
197  }
198 
200  KOKKOS_INLINE_FUNCTION const RealType imag () const {
201  return im_;
202  }
203 
205  KOKKOS_INLINE_FUNCTION const RealType real () const {
206  return re_;
207  }
208 
210  KOKKOS_INLINE_FUNCTION volatile RealType& imag () volatile {
211  return im_;
212  }
213 
215  KOKKOS_INLINE_FUNCTION volatile RealType& real () volatile {
216  return re_;
217  }
218 
220  KOKKOS_INLINE_FUNCTION const RealType imag () const volatile {
221  return im_;
222  }
223 
225  KOKKOS_INLINE_FUNCTION const RealType real () const volatile {
226  return re_;
227  }
228 
229  KOKKOS_INLINE_FUNCTION
230  complex<RealType>& operator += (const complex<RealType>& src) {
231  re_ += src.re_;
232  im_ += src.im_;
233  return *this;
234  }
235 
236  KOKKOS_INLINE_FUNCTION
237  void operator += (const volatile complex<RealType>& src) volatile {
238  re_ += src.re_;
239  im_ += src.im_;
240  }
241 
242  KOKKOS_INLINE_FUNCTION
243  complex<RealType>& operator += (const RealType& src) {
244  re_ += src;
245  return *this;
246  }
247 
248  KOKKOS_INLINE_FUNCTION
249  void operator += (const volatile RealType& src) volatile {
250  re_ += src;
251  }
252 
253  KOKKOS_INLINE_FUNCTION
254  complex<RealType>& operator -= (const complex<RealType>& src) {
255  re_ -= src.re_;
256  im_ -= src.im_;
257  return *this;
258  }
259 
260  KOKKOS_INLINE_FUNCTION
261  complex<RealType>& operator -= (const RealType& src) {
262  re_ -= src;
263  return *this;
264  }
265 
266  KOKKOS_INLINE_FUNCTION
267  complex<RealType>& operator *= (const complex<RealType>& src) {
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  return *this;
273  }
274 
275  KOKKOS_INLINE_FUNCTION
276  void operator *= (const volatile complex<RealType>& src) volatile {
277  const RealType realPart = re_ * src.re_ - im_ * src.im_;
278  const RealType imagPart = re_ * src.im_ + im_ * src.re_;
279  re_ = realPart;
280  im_ = imagPart;
281  }
282 
283  KOKKOS_INLINE_FUNCTION
284  complex<RealType>& operator *= (const RealType& src) {
285  re_ *= src;
286  im_ *= src;
287  return *this;
288  }
289 
290  KOKKOS_INLINE_FUNCTION
291  void operator *= (const volatile RealType& src) volatile {
292  re_ *= src;
293  im_ *= src;
294  }
295 
296  KOKKOS_INLINE_FUNCTION
297  complex<RealType>& operator /= (const complex<RealType>& y) {
298  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
299  // If the real part is +/-Inf and the imaginary part is -/+Inf,
300  // this won't change the result.
301  const RealType s = ::fabs (y.real ()) + ::fabs (y.imag ());
302 
303  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
304  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
305  // because y/s is NaN.
306  if (s == 0.0) {
307  this->re_ /= s;
308  this->im_ /= s;
309  }
310  else {
311  const complex<RealType> x_scaled (this->re_ / s, this->im_ / s);
312  const complex<RealType> y_conj_scaled (y.re_ / s, -(y.im_) / s);
313  const RealType y_scaled_abs = y_conj_scaled.re_ * y_conj_scaled.re_ +
314  y_conj_scaled.im_ * y_conj_scaled.im_; // abs(y) == abs(conj(y))
315  *this = x_scaled * y_conj_scaled;
316  *this /= y_scaled_abs;
317  }
318  return *this;
319  }
320 
321  KOKKOS_INLINE_FUNCTION
322  complex<RealType>& operator /= (const RealType& src) {
323  re_ /= src;
324  im_ /= src;
325  return *this;
326  }
327 };
328 
330 template<class RealType>
331 KOKKOS_INLINE_FUNCTION
332 complex<RealType>
334  return complex<RealType> (x.real () + y.real (), x.imag () + y.imag ());
335 }
336 
338 template<class RealType>
339 KOKKOS_INLINE_FUNCTION
340 complex<RealType>
342  return x;
343 }
344 
346 template<class RealType>
347 KOKKOS_INLINE_FUNCTION
348 complex<RealType>
350  return complex<RealType> (x.real () - y.real (), x.imag () - y.imag ());
351 }
352 
354 template<class RealType>
355 KOKKOS_INLINE_FUNCTION
356 complex<RealType>
358  return complex<RealType> (-x.real (), -x.imag ());
359 }
360 
362 template<class RealType>
363 KOKKOS_INLINE_FUNCTION
364 complex<RealType>
366  return complex<RealType> (x.real () * y.real () - x.imag () * y.imag (),
367  x.real () * y.imag () + x.imag () * y.real ());
368 }
369 
380 template<class RealType>
381 complex<RealType>
382 operator * (const std::complex<RealType>& x, const complex<RealType>& y) {
383  return complex<RealType> (x.real () * y.real () - x.imag () * y.imag (),
384  x.real () * y.imag () + x.imag () * y.real ());
385 }
386 
391 template<class RealType>
392 KOKKOS_INLINE_FUNCTION
393 complex<RealType>
394 operator * (const RealType& x, const complex<RealType>& y) {
395  return complex<RealType> (x * y.real (), x * y.imag ());
396 }
397 
398 
400 template<class RealType>
401 KOKKOS_INLINE_FUNCTION
402 RealType imag (const complex<RealType>& x) {
403  return x.imag ();
404 }
405 
407 template<class RealType>
408 KOKKOS_INLINE_FUNCTION
409 RealType real (const complex<RealType>& x) {
410  return x.real ();
411 }
412 
414 template<class RealType>
415 KOKKOS_INLINE_FUNCTION
416 RealType abs (const complex<RealType>& x) {
417  // FIXME (mfh 31 Oct 2014) Scale to avoid unwarranted overflow.
418  return ::sqrt (real (x) * real (x) + imag (x) * imag (x));
419 }
420 
422 template<class RealType>
423 KOKKOS_INLINE_FUNCTION
425  return complex<RealType> (real (x), -imag (x));
426 }
427 
428 
430 template<class RealType1, class RealType2>
431 KOKKOS_INLINE_FUNCTION
432 complex<RealType1>
433 operator / (const complex<RealType1>& x, const RealType2& y) {
434  return complex<RealType1> (real (x) / y, imag (x) / y);
435 }
436 
438 template<class RealType>
439 KOKKOS_INLINE_FUNCTION
440 complex<RealType>
442  // Scale (by the "1-norm" of y) to avoid unwarranted overflow.
443  // If the real part is +/-Inf and the imaginary part is -/+Inf,
444  // this won't change the result.
445  const RealType s = ::fabs (real (y)) + ::fabs (imag (y));
446 
447  // If s is 0, then y is zero, so x/y == real(x)/0 + i*imag(x)/0.
448  // In that case, the relation x/y == (x/s) / (y/s) doesn't hold,
449  // because y/s is NaN.
450  if (s == 0.0) {
451  return complex<RealType> (real (x) / s, imag (x) / s);
452  }
453  else {
454  const complex<RealType> x_scaled (real (x) / s, imag (x) / s);
455  const complex<RealType> y_conj_scaled (real (y) / s, -imag (y) / s);
456  const RealType y_scaled_abs = real (y_conj_scaled) * real (y_conj_scaled) +
457  imag (y_conj_scaled) * imag (y_conj_scaled); // abs(y) == abs(conj(y))
458  complex<RealType> result = x_scaled * y_conj_scaled;
459  result /= y_scaled_abs;
460  return result;
461  }
462 }
463 
465 template<class RealType>
466 KOKKOS_INLINE_FUNCTION
468  return real (x) == real (y) && imag (x) == imag (y);
469 }
470 
472 template<class RealType>
473 KOKKOS_INLINE_FUNCTION
474 bool operator == (const std::complex<RealType>& x, const complex<RealType>& y) {
475  return std::real (x) == real (y) && std::imag (x) == imag (y);
476 }
477 
479 template<class RealType1, class RealType2>
480 KOKKOS_INLINE_FUNCTION
481 bool operator == (const complex<RealType1>& x, const RealType2& y) {
482  return real (x) == y && imag (x) == static_cast<RealType1> (0.0);
483 }
484 
486 template<class RealType>
487 KOKKOS_INLINE_FUNCTION
488 bool operator == (const RealType& x, const complex<RealType>& y) {
489  return y == x;
490 }
491 
493 template<class RealType>
494 KOKKOS_INLINE_FUNCTION
496  return real (x) != real (y) || imag (x) != imag (y);
497 }
498 
500 template<class RealType>
501 KOKKOS_INLINE_FUNCTION
502 bool operator != (const std::complex<RealType>& x, const complex<RealType>& y) {
503  return std::real (x) != real (y) || std::imag (x) != imag (y);
504 }
505 
507 template<class RealType1, class RealType2>
508 KOKKOS_INLINE_FUNCTION
509 bool operator != (const complex<RealType1>& x, const RealType2& y) {
510  return real (x) != y || imag (x) != static_cast<RealType1> (0.0);
511 }
512 
514 template<class RealType>
515 KOKKOS_INLINE_FUNCTION
516 bool operator != (const RealType& x, const complex<RealType>& y) {
517  return y != x;
518 }
519 
520 template<class RealType>
521 std::ostream& operator << (std::ostream& os, const complex<RealType>& x) {
522  const std::complex<RealType> x_std (Kokkos::real (x), Kokkos::imag (x));
523  os << x_std;
524  return os;
525 }
526 
527 template<class RealType>
528 std::ostream& operator >> (std::ostream& os, complex<RealType>& x) {
529  std::complex<RealType> x_std;
530  os >> x_std;
531  x = x_std; // only assigns on success of above
532  return os;
533 }
534 
535 
536 } // namespace Kokkos
537 
538 #endif // KOKKOS_COMPLEX_HPP
KOKKOS_INLINE_FUNCTION const RealType imag() const
The imaginary part of this complex number.
KOKKOS_INLINE_FUNCTION const RealType real() const
The real part of this complex number.
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 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 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.