ryujin 2.1.1 revision ee5cbcbf2346c1299c942d0e1f13b46449973c18
Loading...
Searching...
No Matches
simd.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include "simd.h"
9#include "simd_fast_pow.template.h"
10
11#include <cmath>
12
13#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
14// Import the vectorlib library prefixed in the vcl namespace
15#define VCL_NAMESPACE vcl
16DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
17#include "../simd-math/vectorclass.h"
18#include "../simd-math/vectormath_exp.h"
19DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
20#undef VCL_NAMESPACE
21#else
22// Make ryujin::pow known as vcl::ryujin
23namespace vcl = ryujin;
24#endif
25
26namespace ryujin
27{
28 /*****************************************************************************
29 * Helper typetraits for dealing with vcl: *
30 ****************************************************************************/
31
32 namespace
33 {
34 /*
35 * A type trait to select the correct VCL type:
36 */
37 template <typename T, std::size_t width>
38 struct VectorClassType {
39 };
40
41#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)
42 template <>
43 struct VectorClassType<float, 16> {
44 using value_type = vcl::Vec16f;
45 };
46
47 template <>
48 struct VectorClassType<double, 8> {
49 using value_type = vcl::Vec8d;
50 };
51#endif
52
53#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)
54 template <>
55 struct VectorClassType<float, 8> {
56 using value_type = vcl::Vec8f;
57 };
58
59 template <>
60 struct VectorClassType<double, 4> {
61 using value_type = vcl::Vec4d;
62 };
63#endif
64
65#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
66 template <>
67 struct VectorClassType<float, 4> {
68 using value_type = vcl::Vec4f;
69 };
70
71 template <>
72 struct VectorClassType<double, 2> {
73 using value_type = vcl::Vec2d;
74 };
75
76 template <>
77 struct VectorClassType<float, 1> {
78 using value_type = vcl::Vec4f;
79 };
80
81 template <>
82 struct VectorClassType<double, 1> {
83 using value_type = vcl::Vec2d;
84 };
85
86#else
87 template <>
88 struct VectorClassType<float, 1> {
89 using value_type = float;
90 };
91
92 template <>
93 struct VectorClassType<double, 1> {
94 using value_type = double;
95 };
96#endif
97
98 /*
99 * Convert a dealii::VectorizedArray to a VCL container type:
100 */
101 template <typename T, std::size_t width>
102 DEAL_II_ALWAYS_INLINE inline typename VectorClassType<T, width>::value_type
103 to_vcl(const dealii::VectorizedArray<T, width> x)
104 {
105 return typename VectorClassType<T, width>::value_type(x.data);
106 }
107
108
109 /*
110 * Convert a VCL container type to a dealii::VectorizedArray:
111 */
112 template <typename T, std::size_t width>
113 DEAL_II_ALWAYS_INLINE inline dealii::VectorizedArray<T, width>
114 from_vcl(typename VectorClassType<T, width>::value_type x)
115 {
116 dealii::VectorizedArray<T, width> result;
117#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
118 if constexpr (width == 1)
119 result.data = x.extract(0);
120 else
121#endif
122 result.data = x;
123 return result;
124 }
125
126
127 /*
128 * Helper functions to convert to float arrays and back.
129 */
130 template <typename T, std::size_t width>
131 struct FC {
132 };
133
134 template <std::size_t width>
135 struct FC<double, width> {
136 // There is no Vec2f, so make sure to use Vec4f instead
137#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
138 static constexpr std::size_t float_width = (width <= 2 ? 4 : width);
139#else
140 static_assert(width == 1, "internal error");
141 static constexpr std::size_t float_width = width;
142#endif
143
144 static DEAL_II_ALWAYS_INLINE inline
145 typename VectorClassType<float, float_width>::value_type
146 to_float(typename VectorClassType<double, width>::value_type x)
147 {
148#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
149 return vcl::to_float(x);
150#else
151 return x;
152#endif
153 }
154
155 static DEAL_II_ALWAYS_INLINE inline
156 typename VectorClassType<double, width>::value_type
157 to_double(typename VectorClassType<float, float_width>::value_type x)
158 {
159#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
160 if constexpr (width == 1) {
161 return static_cast<double>(x.extract(0));
162 } else if constexpr (width == 2) {
163 const vcl::Vec4d temp = vcl::to_double(x);
164 return vcl::Vec2d(temp.extract(0), temp.extract(1));
165 } else {
166 return vcl::to_double(x);
167 }
168#else
169 return x;
170#endif
171 }
172 };
173
174 template <std::size_t width>
175 struct FC<float, width> {
176 static DEAL_II_ALWAYS_INLINE inline
177 typename VectorClassType<float, width>::value_type
178 to_float(typename VectorClassType<float, width>::value_type x)
179 {
180 return x;
181 }
182 static DEAL_II_ALWAYS_INLINE inline
183 typename VectorClassType<float, width>::value_type
184 to_double(typename VectorClassType<float, width>::value_type x)
185 {
186 return x;
187 }
188 };
189 } // namespace
190
191
192 /*****************************************************************************
193 * pow() implementation: *
194 ****************************************************************************/
195
196#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
197 template <typename T, std::size_t width>
198 // DEAL_II_ALWAYS_INLINE inline
199 dealii::VectorizedArray<T, width>
200 pow(const dealii::VectorizedArray<T, width> x, const T b)
201 {
202 return from_vcl<T, width>(vcl::pow(to_vcl(x), b));
203 }
204
205
206 template <typename T, std::size_t width>
207 // DEAL_II_ALWAYS_INLINE inline
208 dealii::VectorizedArray<T, width>
209 pow(const dealii::VectorizedArray<T, width> x,
210 const dealii::VectorizedArray<T, width> b)
211 {
212 return from_vcl<T, width>(vcl::pow(to_vcl(x), to_vcl(b)));
213 }
214
215#else
216
217 template <typename T, std::size_t width>
218 // DEAL_II_ALWAYS_INLINE inline
219 dealii::VectorizedArray<T, width>
220 pow(const dealii::VectorizedArray<T, width> x, const T b)
221 {
222 // Call generic deal.II implementation
223 return std::pow(x, b);
224 }
225
226
227 template <typename T, std::size_t width>
228 // DEAL_II_ALWAYS_INLINE inline
229 dealii::VectorizedArray<T, width>
230 pow(const dealii::VectorizedArray<T, width> x,
231 const dealii::VectorizedArray<T, width> b)
232 {
233 // Call generic deal.II implementation
234 return std::pow(x, b);
235 }
236#endif
237
238
239 /*****************************************************************************
240 * Fast pow() implementation: *
241 ****************************************************************************/
242
243#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
244 template <typename T, std::size_t width>
245 // DEAL_II_ALWAYS_INLINE inline
246 dealii::VectorizedArray<T, width> fast_pow(
247 const dealii::VectorizedArray<T, width> x, const T b, const Bias bias)
248 {
249 using vcl_type = decltype(FC<T, width>::to_float(to_vcl(x)));
250 return from_vcl<T, width>(FC<T, width>::to_double(
251 fast_pow_impl(FC<T, width>::to_float(to_vcl(x)), vcl_type(b), bias)));
252 }
253
254
255 template <typename T, std::size_t width>
256 // DEAL_II_ALWAYS_INLINE inline
257 dealii::VectorizedArray<T, width>
258 fast_pow(const dealii::VectorizedArray<T, width> x,
259 const dealii::VectorizedArray<T, width> b,
260 const Bias bias)
261 {
262 return from_vcl<T, width>(
263 FC<T, width>::to_double(fast_pow_impl(FC<T, width>::to_float(to_vcl(x)),
264 FC<T, width>::to_float(to_vcl(b)),
265 bias)));
266 }
267
268#else
269
270 template <typename T, std::size_t width>
271 // DEAL_II_ALWAYS_INLINE inline
272 dealii::VectorizedArray<T, width>
273 fast_pow(const dealii::VectorizedArray<T, width> x, const T b, const Bias)
274 {
275 // Call generic deal.II implementation
276 return std::pow(x, b);
277 }
278
279
280 template <typename T, std::size_t width>
281 // DEAL_II_ALWAYS_INLINE inline
282 dealii::VectorizedArray<T, width>
283 fast_pow(const dealii::VectorizedArray<T, width> x,
284 const dealii::VectorizedArray<T, width> b,
285 const Bias)
286 {
287 // Call generic deal.II implementation
288 return std::pow(x, b);
289 }
290#endif
291
292
293} // namespace ryujin
ryujin::pow
DEAL_II_HOST_DEVICE T pow(const T x, const T b)
ryujin::fast_pow
DEAL_II_HOST_DEVICE T fast_pow(const T x, const T b, const Bias bias=Bias::none)
ryujin::fast_pow_impl
T fast_pow_impl(const T x, const T b, const Bias)
Definition simd_fast_pow.template.h:25
ryujin::Bias
Bias
Definition simd.h:300
simd_fast_pow.template.h
Generated by doxygen 1.9.8