Medial Code Documentation
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
11#ifndef EIGEN_COMPLEX_CUDA_H
12#define EIGEN_COMPLEX_CUDA_H
13
14// Many std::complex methods such as operator+, operator-, operator* and
15// operator/ are not constexpr. Due to this, GCC and older versions of clang do
16// not treat them as device functions and thus Eigen functors making use of
17// these operators fail to compile. Here, we manually specialize these
18// operators and functors for complex types when building for CUDA to enable
19// their use on-device.
20//
21// NOTES:
22// - Compound assignment operators +=,-=,*=,/=(Scalar) will not work on device,
23// since they are already specialized in the standard. Using them will result
24// in silent kernel failures.
25// - Compiling with MSVC and using +=,-=,*=,/=(std::complex<Scalar>) will lead
26// to duplicate definition errors, since these are already specialized in
27// Visual Studio's <complex> header (contrary to the standard). This is
28// preferable to removing such definitions, which will lead to silent kernel
29// failures.
30// - Compiling with ICC requires defining _USE_COMPLEX_SPECIALIZATION_ prior
31// to the first inclusion of <complex>.
32
33#if defined(EIGEN_CUDACC) && defined(EIGEN_GPU_COMPILE_PHASE)
34
35// ICC already specializes std::complex<float> and std::complex<double>
36// operators, preventing us from making them device functions here.
37// This will lead to silent runtime errors if the operators are used on device.
38//
39// To allow std::complex operator use on device, define _OVERRIDE_COMPLEX_SPECIALIZATION_
40// prior to first inclusion of <complex>. This prevents ICC from adding
41// its own specializations, so our custom ones below can be used instead.
42#if !(defined(EIGEN_COMP_ICC) && defined(_USE_COMPLEX_SPECIALIZATION_))
43
44// Import Eigen's internal operator specializations.
45#define EIGEN_USING_STD_COMPLEX_OPERATORS \
46 using Eigen::complex_operator_detail::operator+; \
47 using Eigen::complex_operator_detail::operator-; \
48 using Eigen::complex_operator_detail::operator*; \
49 using Eigen::complex_operator_detail::operator/; \
50 using Eigen::complex_operator_detail::operator+=; \
51 using Eigen::complex_operator_detail::operator-=; \
52 using Eigen::complex_operator_detail::operator*=; \
53 using Eigen::complex_operator_detail::operator/=; \
54 using Eigen::complex_operator_detail::operator==; \
55 using Eigen::complex_operator_detail::operator!=;
56
57namespace Eigen {
58
59// Specialized std::complex overloads.
60namespace complex_operator_detail {
61
62template<typename T>
63EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
64std::complex<T> complex_multiply(const std::complex<T>& a, const std::complex<T>& b) {
65 const T a_real = numext::real(a);
66 const T a_imag = numext::imag(a);
67 const T b_real = numext::real(b);
68 const T b_imag = numext::imag(b);
69 return std::complex<T>(
70 a_real * b_real - a_imag * b_imag,
71 a_imag * b_real + a_real * b_imag);
72}
73
74template<typename T>
75EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
76std::complex<T> complex_divide_fast(const std::complex<T>& a, const std::complex<T>& b) {
77 const T a_real = numext::real(a);
78 const T a_imag = numext::imag(a);
79 const T b_real = numext::real(b);
80 const T b_imag = numext::imag(b);
81 const T norm = (b_real * b_real + b_imag * b_imag);
82 return std::complex<T>((a_real * b_real + a_imag * b_imag) / norm,
83 (a_imag * b_real - a_real * b_imag) / norm);
84}
85
86template<typename T>
87EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
88std::complex<T> complex_divide_stable(const std::complex<T>& a, const std::complex<T>& b) {
89 const T a_real = numext::real(a);
90 const T a_imag = numext::imag(a);
91 const T b_real = numext::real(b);
92 const T b_imag = numext::imag(b);
93 // Smith's complex division (https://arxiv.org/pdf/1210.4539.pdf),
94 // guards against over/under-flow.
95 const bool scale_imag = numext::abs(b_imag) <= numext::abs(b_real);
96 const T rscale = scale_imag ? T(1) : b_real / b_imag;
97 const T iscale = scale_imag ? b_imag / b_real : T(1);
98 const T denominator = b_real * rscale + b_imag * iscale;
99 return std::complex<T>((a_real * rscale + a_imag * iscale) / denominator,
100 (a_imag * rscale - a_real * iscale) / denominator);
101}
102
103template<typename T>
104EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
105std::complex<T> complex_divide(const std::complex<T>& a, const std::complex<T>& b) {
106#if EIGEN_FAST_MATH
107 return complex_divide_fast(a, b);
108#else
109 return complex_divide_stable(a, b);
110#endif
111}
112
113// NOTE: We cannot specialize compound assignment operators with Scalar T,
114// (i.e. operator@=(const T&), for @=+,-,*,/)
115// since they are already specialized for float/double/long double within
116// the standard <complex> header. We also do not specialize the stream
117// operators.
118#define EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(T) \
119 \
120EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
121std::complex<T> operator+(const std::complex<T>& a) { return a; } \
122 \
123EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
124std::complex<T> operator-(const std::complex<T>& a) { \
125 return std::complex<T>(-numext::real(a), -numext::imag(a)); \
126} \
127 \
128EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
129std::complex<T> operator+(const std::complex<T>& a, const std::complex<T>& b) { \
130 return std::complex<T>(numext::real(a) + numext::real(b), numext::imag(a) + numext::imag(b)); \
131} \
132 \
133EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
134std::complex<T> operator+(const std::complex<T>& a, const T& b) { \
135 return std::complex<T>(numext::real(a) + b, numext::imag(a)); \
136} \
137 \
138EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
139std::complex<T> operator+(const T& a, const std::complex<T>& b) { \
140 return std::complex<T>(a + numext::real(b), numext::imag(b)); \
141} \
142 \
143EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
144std::complex<T> operator-(const std::complex<T>& a, const std::complex<T>& b) { \
145 return std::complex<T>(numext::real(a) - numext::real(b), numext::imag(a) - numext::imag(b)); \
146} \
147 \
148EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
149std::complex<T> operator-(const std::complex<T>& a, const T& b) { \
150 return std::complex<T>(numext::real(a) - b, numext::imag(a)); \
151} \
152 \
153EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
154std::complex<T> operator-(const T& a, const std::complex<T>& b) { \
155 return std::complex<T>(a - numext::real(b), -numext::imag(b)); \
156} \
157 \
158EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
159std::complex<T> operator*(const std::complex<T>& a, const std::complex<T>& b) { \
160 return complex_multiply(a, b); \
161} \
162 \
163EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
164std::complex<T> operator*(const std::complex<T>& a, const T& b) { \
165 return std::complex<T>(numext::real(a) * b, numext::imag(a) * b); \
166} \
167 \
168EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
169std::complex<T> operator*(const T& a, const std::complex<T>& b) { \
170 return std::complex<T>(a * numext::real(b), a * numext::imag(b)); \
171} \
172 \
173EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
174std::complex<T> operator/(const std::complex<T>& a, const std::complex<T>& b) { \
175 return complex_divide(a, b); \
176} \
177 \
178EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
179std::complex<T> operator/(const std::complex<T>& a, const T& b) { \
180 return std::complex<T>(numext::real(a) / b, numext::imag(a) / b); \
181} \
182 \
183EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
184std::complex<T> operator/(const T& a, const std::complex<T>& b) { \
185 return complex_divide(std::complex<T>(a, 0), b); \
186} \
187 \
188EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
189std::complex<T>& operator+=(std::complex<T>& a, const std::complex<T>& b) { \
190 numext::real_ref(a) += numext::real(b); \
191 numext::imag_ref(a) += numext::imag(b); \
192 return a; \
193} \
194 \
195EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
196std::complex<T>& operator-=(std::complex<T>& a, const std::complex<T>& b) { \
197 numext::real_ref(a) -= numext::real(b); \
198 numext::imag_ref(a) -= numext::imag(b); \
199 return a; \
200} \
201 \
202EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
203std::complex<T>& operator*=(std::complex<T>& a, const std::complex<T>& b) { \
204 a = complex_multiply(a, b); \
205 return a; \
206} \
207 \
208EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
209std::complex<T>& operator/=(std::complex<T>& a, const std::complex<T>& b) { \
210 a = complex_divide(a, b); \
211 return a; \
212} \
213 \
214EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
215bool operator==(const std::complex<T>& a, const std::complex<T>& b) { \
216 return numext::real(a) == numext::real(b) && numext::imag(a) == numext::imag(b); \
217} \
218 \
219EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
220bool operator==(const std::complex<T>& a, const T& b) { \
221 return numext::real(a) == b && numext::imag(a) == 0; \
222} \
223 \
224EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
225bool operator==(const T& a, const std::complex<T>& b) { \
226 return a == numext::real(b) && 0 == numext::imag(b); \
227} \
228 \
229EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
230bool operator!=(const std::complex<T>& a, const std::complex<T>& b) { \
231 return !(a == b); \
232} \
233 \
234EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
235bool operator!=(const std::complex<T>& a, const T& b) { \
236 return !(a == b); \
237} \
238 \
239EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE \
240bool operator!=(const T& a, const std::complex<T>& b) { \
241 return !(a == b); \
242}
243
244// Do not specialize for long double, since that reduces to double on device.
245EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(float)
246EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS(double)
247
248#undef EIGEN_CREATE_STD_COMPLEX_OPERATOR_SPECIALIZATIONS
249
250
251} // namespace complex_operator_detail
252
253EIGEN_USING_STD_COMPLEX_OPERATORS
254
255namespace numext {
256EIGEN_USING_STD_COMPLEX_OPERATORS
257} // namespace numext
258
259namespace internal {
260EIGEN_USING_STD_COMPLEX_OPERATORS
261
262} // namespace internal
263} // namespace Eigen
264
265#endif // !(EIGEN_COMP_ICC && _USE_COMPLEX_SPECIALIZATION_)
266
267#endif // EIGEN_CUDACC && EIGEN_GPU_COMPILE_PHASE
268
269#endif // EIGEN_COMPLEX_CUDA_H
Eigen::MatrixBase::norm
EIGEN_DEVICE_FUNC RealScalar norm() const
Definition Dot.h:103
Eigen
Namespace containing all symbols from the Eigen library.
Definition LDLT.h:16
Generated on Thu Dec 18 2025 18:56:06 for Medial Code Documentation by doxygen 1.9.8