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) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com>
5// Copyright (C) 2021 C. Antonio Sanchez <cantonios@google.com>
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_COMPLEX_GPU_H
13#define EIGEN_COMPLEX_GPU_H
14
15// Many std::complex methods such as operator+, operator-, operator* and
16// operator/ are not constexpr. Due to this, GCC and older versions of clang do
17// not treat them as device functions and thus Eigen functors making use of
18// these operators fail to compile. Here, we manually specialize these
19// operators for complex types when building for CUDA to enable their use
20// on-device.
21//
22// Eigen/Core includes this header ahead of the Core headers (Meta.h's
23// equal_strict, MathFunctions.h, GenericPacketMath.h, the functors), which
24// apply these operators to a dependent Scalar: two-phase
25// lookup finds non-ADL candidates in the definition context only, and ADL for
26// std::complex<T> reaches namespace std alone, so an overload declared later
27// leaves those templates bound to the host-only std:: operators. nvcc before
28// CUDA 13.3 resolved them at the instantiation point, which masked the ordering.
29//
30// NOTES:
31// - Compound assignment operators +=,-=,*=,/=(Scalar) will not work on device,
32// since they are already specialized in the standard. Using them will result
33// in silent kernel failures.
34// - Compiling with MSVC and using +=,-=,*=,/=(std::complex<Scalar>) will lead
35// to duplicate definition errors, since these are already specialized in
36// Visual Studio's <complex> header (contrary to the standard). This is
37// preferable to removing such definitions, which will lead to silent kernel
38// failures.
39// - Compiling with ICC requires defining _USE_COMPLEX_SPECIALIZATION_ prior
40// to the first inclusion of <complex>.
41// - Device code outside namespace Eigen that applies these operators to a
42// dependent Scalar reaches them only through `using namespace Eigen;` or
43// using-declarations in scope (see test/gpu_basic.cu); ADL alone finds the
44// host-only std:: operators.
45
46#if defined(EIGEN_GPUCC) && defined(EIGEN_GPU_COMPILE_PHASE)
47
48// ICC already specializes std::complex<float> and std::complex<double>
49// operators, preventing us from making them device functions here.
50// This will lead to silent runtime errors if the operators are used on device.
51//
52// To allow std::complex operator use on device, define _OVERRIDE_COMPLEX_SPECIALIZATION_
53// prior to first inclusion of <complex>. This prevents ICC from adding
54// its own specializations, so our custom ones below can be used instead.
55#if !(EIGEN_COMP_ICC && defined(_USE_COMPLEX_SPECIALIZATION_))
56
57// Import Eigen's internal operator specializations.
58#define EIGEN_USING_STD_COMPLEX_OPERATORS \
59 using Eigen::complex_operator_detail::operator+; \
60 using Eigen::complex_operator_detail::operator-; \
61 using Eigen::complex_operator_detail::operator*; \
62 using Eigen::complex_operator_detail::operator/; \
63 using Eigen::complex_operator_detail::operator+=; \
64 using Eigen::complex_operator_detail::operator-=; \
65 using Eigen::complex_operator_detail::operator*=; \
66 using Eigen::complex_operator_detail::operator/=; \
67 using Eigen::complex_operator_detail::operator==; \
68 using Eigen::complex_operator_detail::operator!=;
69
70// IWYU pragma: private
71#include "../../InternalHeaderCheck.h"
72
73namespace Eigen {
74
75namespace internal {
76// Defined in MathFunctions.h, which follows this header.
77template <typename T>
78EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_multiply(const std::complex<T>& a,
79 const std::complex<T>& b);
80template <typename T>
81EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> complex_divide(const std::complex<T>& a,
82 const std::complex<T>& b);
83} // namespace internal
84
85// Specialized std::complex overloads.
86namespace complex_operator_detail {
87
88// NOTE: We cannot specialize compound assignment operators with Scalar T,
89// (i.e. operator@=(const T&), for @=+,-,*,/)
90// since they are already specialized for float/double/long double within
91// the standard <complex> header. We also do not specialize the stream
92// operators.
93// numext is not declared yet; the std::complex members used below are
94// constexpr, hence device-callable under EIGEN_CONSTEXPR_ARE_DEVICE_FUNC
95// (nvcc: --expt-relaxed-constexpr), as real_impl<std::complex<T>>
96// already assumes.
97#define EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(T) \
98 \
99 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator+(const std::complex<T>& a) { return a; } \
100 \
101 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator-(const std::complex<T>& a) { \
102 return std::complex<T>(-a.real(), -a.imag()); \
103 } \
104 \
105 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator+(const std::complex<T>& a, \
106 const std::complex<T>& b) { \
107 return std::complex<T>(a.real() + b.real(), a.imag() + b.imag()); \
108 } \
109 \
110 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator+(const std::complex<T>& a, const T& b) { \
111 return std::complex<T>(a.real() + b, a.imag()); \
112 } \
113 \
114 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator+(const T& a, const std::complex<T>& b) { \
115 return std::complex<T>(a + b.real(), b.imag()); \
116 } \
117 \
118 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator-(const std::complex<T>& a, \
119 const std::complex<T>& b) { \
120 return std::complex<T>(a.real() - b.real(), a.imag() - b.imag()); \
121 } \
122 \
123 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator-(const std::complex<T>& a, const T& b) { \
124 return std::complex<T>(a.real() - b, a.imag()); \
125 } \
126 \
127 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator-(const T& a, const std::complex<T>& b) { \
128 return std::complex<T>(a - b.real(), -b.imag()); \
129 } \
130 \
131 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator*(const std::complex<T>& a, \
132 const std::complex<T>& b) { \
133 return internal::complex_multiply(a, b); \
134 } \
135 \
136 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator*(const std::complex<T>& a, const T& b) { \
137 return std::complex<T>(a.real() * b, a.imag() * b); \
138 } \
139 \
140 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator*(const T& a, const std::complex<T>& b) { \
141 return std::complex<T>(a * b.real(), a * b.imag()); \
142 } \
143 \
144 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator/(const std::complex<T>& a, \
145 const std::complex<T>& b) { \
146 return internal::complex_divide(a, b); \
147 } \
148 \
149 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator/(const std::complex<T>& a, const T& b) { \
150 return std::complex<T>(a.real() / b, a.imag() / b); \
151 } \
152 \
153 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T> operator/(const T& a, const std::complex<T>& b) { \
154 return internal::complex_divide(std::complex<T>(a, 0), b); \
155 } \
156 \
157 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T>& operator+=(std::complex<T>& a, const std::complex<T>& b) { \
158 a = std::complex<T>(a.real() + b.real(), a.imag() + b.imag()); \
159 return a; \
160 } \
161 \
162 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T>& operator-=(std::complex<T>& a, const std::complex<T>& b) { \
163 a = std::complex<T>(a.real() - b.real(), a.imag() - b.imag()); \
164 return a; \
165 } \
166 \
167 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T>& operator*=(std::complex<T>& a, const std::complex<T>& b) { \
168 a = internal::complex_multiply(a, b); \
169 return a; \
170 } \
171 \
172 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE std::complex<T>& operator/=(std::complex<T>& a, const std::complex<T>& b) { \
173 a = internal::complex_divide(a, b); \
174 return a; \
175 } \
176 \
177 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator==(const std::complex<T>& a, const std::complex<T>& b) { \
178 return a.real() == b.real() && a.imag() == b.imag(); \
179 } \
180 \
181 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator==(const std::complex<T>& a, const T& b) { \
182 return a.real() == b && a.imag() == 0; \
183 } \
184 \
185 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator==(const T& a, const std::complex<T>& b) { \
186 return a == b.real() && 0 == b.imag(); \
187 } \
188 \
189 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator!=(const std::complex<T>& a, const std::complex<T>& b) { \
190 return !(a == b); \
191 } \
192 \
193 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator!=(const std::complex<T>& a, const T& b) { return !(a == b); } \
194 \
195 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator!=(const T& a, const std::complex<T>& b) { return !(a == b); }
196
197// Do not specialize for long double, since that reduces to double on device.
198EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(float)
199EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(double)
200
201#undef EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS
202
203} // namespace complex_operator_detail
204
205EIGEN_USING_STD_COMPLEX_OPERATORS
206
207namespace numext {
208EIGEN_USING_STD_COMPLEX_OPERATORS
209} // namespace numext
210
211namespace internal {
212EIGEN_USING_STD_COMPLEX_OPERATORS
213
214} // namespace internal
215} // namespace Eigen
216
217#undef EIGEN_USING_STD_COMPLEX_OPERATORS
218
219#endif // !(EIGEN_COMP_ICC && _USE_COMPLEX_SPECIALIZATION_)
220
221#endif // EIGEN_GPUCC && EIGEN_GPU_COMPILE_PHASE
222
223#endif // EIGEN_COMPLEX_GPU_H