1 // Copyright (c) 2018-2019, NVIDIA CORPORATION. All rights reserved.
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // http://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 // This file defines host runtime functions that can be used for folding
16 // intrinsic functions.
17 // The default HostIntrinsicProceduresLibrary is built with <cmath> and
18 // <complex> functions that are guaranteed to exist from the C++ standard.
19
20 #include "intrinsics-library-templates.h"
21 #include <cmath>
22 #include <complex>
23
24 namespace Fortran::evaluate {
25
26 // Note: argument passing is ignored in equivalence
HasEquivalentProcedure(const IntrinsicProcedureRuntimeDescription & sym) const27 bool HostIntrinsicProceduresLibrary::HasEquivalentProcedure(
28 const IntrinsicProcedureRuntimeDescription &sym) const {
29 const auto rteProcRange{procedures_.equal_range(sym.name)};
30 const size_t nargs{sym.argumentsType.size()};
31 for (auto iter{rteProcRange.first}; iter != rteProcRange.second; ++iter) {
32 if (nargs == iter->second.argumentsType.size() &&
33 sym.returnType == iter->second.returnType &&
34 (sym.isElemental || iter->second.isElemental)) {
35 bool match{true};
36 int pos{0};
37 for (const auto &type : sym.argumentsType) {
38 if (type != iter->second.argumentsType[pos++]) {
39 match = false;
40 break;
41 }
42 }
43 if (match) {
44 return true;
45 }
46 }
47 }
48 return false;
49 }
50
51 // Map numerical intrinsic to <cmath>/<complex> functions
52
53 // Define which host runtime functions will be used for folding
54
55 template<typename HostT>
AddLibmRealHostProcedures(HostIntrinsicProceduresLibrary & hostIntrinsicLibrary)56 static void AddLibmRealHostProcedures(
57 HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
58 using F = FuncPointer<HostT, HostT>;
59 using F2 = FuncPointer<HostT, HostT, HostT>;
60 HostRuntimeIntrinsicProcedure libmSymbols[]{
61 {"acos", F{std::acos}, true},
62 {"acosh", F{std::acosh}, true},
63 {"asin", F{std::asin}, true},
64 {"asinh", F{std::asinh}, true},
65 {"atan", F{std::atan}, true},
66 {"atan", F2{std::atan2}, true},
67 {"atanh", F{std::atanh}, true},
68 {"cos", F{std::cos}, true},
69 {"cosh", F{std::cosh}, true},
70 {"erf", F{std::erf}, true},
71 {"erfc", F{std::erfc}, true},
72 {"exp", F{std::exp}, true},
73 {"gamma", F{std::tgamma}, true},
74 {"hypot", F2{std::hypot}, true},
75 {"log", F{std::log}, true},
76 {"log10", F{std::log10}, true},
77 {"log_gamma", F{std::lgamma}, true},
78 {"mod", F2{std::fmod}, true},
79 {"pow", F2{std::pow}, true},
80 {"sin", F{std::sin}, true},
81 {"sinh", F{std::sinh}, true},
82 {"sqrt", F{std::sqrt}, true},
83 {"tan", F{std::tan}, true},
84 {"tanh", F{std::tanh}, true},
85 };
86 // Note: cmath does not have modulo and erfc_scaled equivalent
87
88 // Note regarding lack of bessel function support:
89 // C++17 defined standard Bessel math functions std::cyl_bessel_j
90 // and std::cyl_neumann that can be used for Fortran j and y
91 // bessel functions. However, they are not yet implemented in
92 // clang libc++ (ok in GNU libstdc++). C maths functions j0...
93 // are not C standard but a GNU extension so they are not used
94 // to avoid introducing incompatibilities.
95 // Use libpgmath to get bessel function folding support.
96 // TODO: Add Bessel functions when possible.
97
98 for (auto sym : libmSymbols) {
99 if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
100 hostIntrinsicLibrary.AddProcedure(std::move(sym));
101 }
102 }
103 }
104
105 template<typename HostT>
AddLibmComplexHostProcedures(HostIntrinsicProceduresLibrary & hostIntrinsicLibrary)106 static void AddLibmComplexHostProcedures(
107 HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
108 using F = FuncPointer<std::complex<HostT>, const std::complex<HostT> &>;
109 using F2 = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
110 const std::complex<HostT> &>;
111 using F2a = FuncPointer<std::complex<HostT>, const HostT &,
112 const std::complex<HostT> &>;
113 using F2b = FuncPointer<std::complex<HostT>, const std::complex<HostT> &,
114 const HostT &>;
115 HostRuntimeIntrinsicProcedure libmSymbols[]{
116 {"abs", FuncPointer<HostT, const std::complex<HostT> &>{std::abs}, true},
117 {"acos", F{std::acos}, true},
118 {"acosh", F{std::acosh}, true},
119 {"asin", F{std::asin}, true},
120 {"asinh", F{std::asinh}, true},
121 {"atan", F{std::atan}, true},
122 {"atanh", F{std::atanh}, true},
123 {"cos", F{std::cos}, true},
124 {"cosh", F{std::cosh}, true},
125 {"exp", F{std::exp}, true},
126 {"log", F{std::log}, true},
127 {"pow", F2{std::pow}, true},
128 {"pow", F2a{std::pow}, true},
129 {"pow", F2b{std::pow}, true},
130 {"sin", F{std::sin}, true},
131 {"sinh", F{std::sinh}, true},
132 {"sqrt", F{std::sqrt}, true},
133 {"tan", F{std::tan}, true},
134 {"tanh", F{std::tanh}, true},
135 };
136
137 for (auto sym : libmSymbols) {
138 if (!hostIntrinsicLibrary.HasEquivalentProcedure(sym)) {
139 hostIntrinsicLibrary.AddProcedure(std::move(sym));
140 }
141 }
142 }
143
InitHostIntrinsicLibraryWithLibm(HostIntrinsicProceduresLibrary & lib)144 void InitHostIntrinsicLibraryWithLibm(HostIntrinsicProceduresLibrary &lib) {
145 if constexpr (host::FortranTypeExists<float>()) {
146 AddLibmRealHostProcedures<float>(lib);
147 }
148 if constexpr (host::FortranTypeExists<double>()) {
149 AddLibmRealHostProcedures<double>(lib);
150 }
151 if constexpr (host::FortranTypeExists<long double>()) {
152 AddLibmRealHostProcedures<long double>(lib);
153 }
154
155 if constexpr (host::FortranTypeExists<std::complex<float>>()) {
156 AddLibmComplexHostProcedures<float>(lib);
157 }
158 if constexpr (host::FortranTypeExists<std::complex<double>>()) {
159 AddLibmComplexHostProcedures<double>(lib);
160 }
161 if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
162 AddLibmComplexHostProcedures<long double>(lib);
163 }
164 }
165
166 #if LINK_WITH_LIBPGMATH
167 namespace pgmath {
168 // Define mapping between numerical intrinsics and libpgmath symbols
169 // namespace is used to have shorter names on repeated patterns.
170 // A class would be better to hold all these defs, but GCC does not
171 // support specialization of template variables inside class even
172 // if it is C++14 standard compliant here because there are only full
173 // specializations.
174
175 // List of intrinsics that have libpgmath implementations that can be used for
176 // constant folding The tag names must match the name used inside libpgmath name
177 // so that the macro below work.
178 enum class I {
179 acos,
180 acosh,
181 asin,
182 asinh,
183 atan,
184 atan2,
185 atanh,
186 bessel_j0,
187 bessel_j1,
188 bessel_jn,
189 bessel_y0,
190 bessel_y1,
191 bessel_yn,
192 cos,
193 cosh,
194 erf,
195 erfc,
196 erfc_scaled,
197 exp,
198 gamma,
199 hypot,
200 log,
201 log10,
202 log_gamma,
203 mod,
204 pow,
205 sin,
206 sinh,
207 sqrt,
208 tan,
209 tanh
210 };
211
212 // Library versions: P for Precise, R for Relaxed, F for Fast
213 enum class L { F, R, P };
214
215 struct NoSuchRuntimeSymbol {};
216 template<L, I, typename> constexpr auto Sym{NoSuchRuntimeSymbol{}};
217
218 // Macros to declare fast/relaxed/precise libpgmath variants.
219 #define DECLARE_PGMATH_FAST_REAL(func) \
220 extern "C" float __fs_##func##_1(float); \
221 extern "C" double __fd_##func##_1(double); \
222 template<> constexpr auto Sym<L::F, I::func, float>{__fs_##func##_1}; \
223 template<> constexpr auto Sym<L::F, I::func, double>{__fd_##func##_1};
224
225 #define DECLARE_PGMATH_FAST_COMPLEX(func) \
226 extern "C" float _Complex __fc_##func##_1(float _Complex); \
227 extern "C" double _Complex __fz_##func##_1(double _Complex); \
228 template<> \
229 constexpr auto Sym<L::F, I::func, std::complex<float>>{__fc_##func##_1}; \
230 template<> \
231 constexpr auto Sym<L::F, I::func, std::complex<double>>{__fz_##func##_1};
232
233 #define DECLARE_PGMATH_FAST_ALL_FP(func) \
234 DECLARE_PGMATH_FAST_REAL(func) \
235 DECLARE_PGMATH_FAST_COMPLEX(func)
236
237 #define DECLARE_PGMATH_PRECISE_REAL(func) \
238 extern "C" float __ps_##func##_1(float); \
239 extern "C" double __pd_##func##_1(double); \
240 template<> constexpr auto Sym<L::P, I::func, float>{__ps_##func##_1}; \
241 template<> constexpr auto Sym<L::P, I::func, double>{__pd_##func##_1};
242
243 #define DECLARE_PGMATH_PRECISE_COMPLEX(func) \
244 extern "C" float _Complex __pc_##func##_1(float _Complex); \
245 extern "C" double _Complex __pz_##func##_1(double _Complex); \
246 template<> \
247 constexpr auto Sym<L::P, I::func, std::complex<float>>{__pc_##func##_1}; \
248 template<> \
249 constexpr auto Sym<L::P, I::func, std::complex<double>>{__pz_##func##_1};
250
251 #define DECLARE_PGMATH_PRECISE_ALL_FP(func) \
252 DECLARE_PGMATH_PRECISE_REAL(func) \
253 DECLARE_PGMATH_PRECISE_COMPLEX(func)
254
255 #define DECLARE_PGMATH_RELAXED_REAL(func) \
256 extern "C" float __rs_##func##_1(float); \
257 extern "C" double __rd_##func##_1(double); \
258 template<> constexpr auto Sym<L::R, I::func, float>{__rs_##func##_1}; \
259 template<> constexpr auto Sym<L::R, I::func, double>{__rd_##func##_1};
260
261 #define DECLARE_PGMATH_RELAXED_COMPLEX(func) \
262 extern "C" float _Complex __rc_##func##_1(float _Complex); \
263 extern "C" double _Complex __rz_##func##_1(double _Complex); \
264 template<> \
265 constexpr auto Sym<L::R, I::func, std::complex<float>>{__rc_##func##_1}; \
266 template<> \
267 constexpr auto Sym<L::R, I::func, std::complex<double>>{__rz_##func##_1};
268
269 #define DECLARE_PGMATH_RELAXED_ALL_FP(func) \
270 DECLARE_PGMATH_RELAXED_REAL(func) \
271 DECLARE_PGMATH_RELAXED_COMPLEX(func)
272
273 #define DECLARE_PGMATH_REAL(func) \
274 DECLARE_PGMATH_FAST_REAL(func) \
275 DECLARE_PGMATH_PRECISE_REAL(func) \
276 DECLARE_PGMATH_RELAXED_REAL(func)
277
278 #define DECLARE_PGMATH_COMPLEX(func) \
279 DECLARE_PGMATH_FAST_COMPLEX(func) \
280 DECLARE_PGMATH_PRECISE_COMPLEX(func) \
281 DECLARE_PGMATH_RELAXED_COMPLEX(func)
282
283 #define DECLARE_PGMATH_ALL(func) \
284 DECLARE_PGMATH_REAL(func) \
285 DECLARE_PGMATH_COMPLEX(func)
286
287 // Macros to declare fast/relaxed/precise libpgmath variants with two arguments.
288 #define DECLARE_PGMATH_FAST_REAL2(func) \
289 extern "C" float __fs_##func##_1(float, float); \
290 extern "C" double __fd_##func##_1(double, double); \
291 template<> constexpr auto Sym<L::F, I::func, float>{__fs_##func##_1}; \
292 template<> constexpr auto Sym<L::F, I::func, double>{__fd_##func##_1};
293
294 #define DECLARE_PGMATH_FAST_COMPLEX2(func) \
295 extern "C" float _Complex __fc_##func##_1(float _Complex, float _Complex); \
296 extern "C" double _Complex __fz_##func##_1( \
297 double _Complex, double _Complex); \
298 template<> \
299 constexpr auto Sym<L::F, I::func, std::complex<float>>{__fc_##func##_1}; \
300 template<> \
301 constexpr auto Sym<L::F, I::func, std::complex<double>>{__fz_##func##_1};
302
303 #define DECLARE_PGMATH_FAST_ALL_FP2(func) \
304 DECLARE_PGMATH_FAST_REAL2(func) \
305 DECLARE_PGMATH_FAST_COMPLEX2(func)
306
307 #define DECLARE_PGMATH_PRECISE_REAL2(func) \
308 extern "C" float __ps_##func##_1(float, float); \
309 extern "C" double __pd_##func##_1(double, double); \
310 template<> constexpr auto Sym<L::P, I::func, float>{__ps_##func##_1}; \
311 template<> constexpr auto Sym<L::P, I::func, double>{__pd_##func##_1};
312
313 #define DECLARE_PGMATH_PRECISE_COMPLEX2(func) \
314 extern "C" float _Complex __pc_##func##_1(float _Complex, float _Complex); \
315 extern "C" double _Complex __pz_##func##_1( \
316 double _Complex, double _Complex); \
317 template<> \
318 constexpr auto Sym<L::P, I::func, std::complex<float>>{__pc_##func##_1}; \
319 template<> \
320 constexpr auto Sym<L::P, I::func, std::complex<double>>{__pz_##func##_1};
321
322 #define DECLARE_PGMATH_PRECISE_ALL_FP2(func) \
323 DECLARE_PGMATH_PRECISE_REAL2(func) \
324 DECLARE_PGMATH_PRECISE_COMPLEX2(func)
325
326 #define DECLARE_PGMATH_RELAXED_REAL2(func) \
327 extern "C" float __rs_##func##_1(float, float); \
328 extern "C" double __rd_##func##_1(double, double); \
329 template<> constexpr auto Sym<L::R, I::func, float>{__rs_##func##_1}; \
330 template<> constexpr auto Sym<L::R, I::func, double>{__rd_##func##_1};
331
332 #define DECLARE_PGMATH_RELAXED_COMPLEX2(func) \
333 extern "C" float _Complex __rc_##func##_1(float _Complex, float _Complex); \
334 extern "C" double _Complex __rz_##func##_1( \
335 double _Complex, double _Complex); \
336 template<> \
337 constexpr auto Sym<L::R, I::func, std::complex<float>>{__rc_##func##_1}; \
338 template<> \
339 constexpr auto Sym<L::R, I::func, std::complex<double>>{__rz_##func##_1};
340
341 #define DECLARE_PGMATH_RELAXED_ALL_FP2(func) \
342 DECLARE_PGMATH_RELAXED_REAL2(func) \
343 DECLARE_PGMATH_RELAXED_COMPLEX2(func)
344
345 #define DECLARE_PGMATH_REAL2(func) \
346 DECLARE_PGMATH_FAST_REAL2(func) \
347 DECLARE_PGMATH_PRECISE_REAL2(func) \
348 DECLARE_PGMATH_RELAXED_REAL2(func)
349
350 #define DECLARE_PGMATH_COMPLEX2(func) \
351 DECLARE_PGMATH_FAST_COMPLEX2(func) \
352 DECLARE_PGMATH_PRECISE_COMPLEX2(func) \
353 DECLARE_PGMATH_RELAXED_COMPLEX2(func)
354
355 #define DECLARE_PGMATH_ALL2(func) \
356 DECLARE_PGMATH_REAL2(func) \
357 DECLARE_PGMATH_COMPLEX2(func)
358
359 // Marcos to declare __mth_i libpgmath variants
360 #define DECLARE_PGMATH_MTH_VERSION_REAL(func) \
361 extern "C" float __mth_i_##func(float); \
362 extern "C" double __mth_i_d##func(double); \
363 template<> constexpr auto Sym<L::F, I::func, float>{__mth_i_##func}; \
364 template<> constexpr auto Sym<L::F, I::func, double>{__mth_i_d##func}; \
365 template<> constexpr auto Sym<L::P, I::func, float>{__mth_i_##func}; \
366 template<> constexpr auto Sym<L::P, I::func, double>{__mth_i_d##func}; \
367 template<> constexpr auto Sym<L::R, I::func, float>{__mth_i_##func}; \
368 template<> constexpr auto Sym<L::R, I::func, double>{__mth_i_d##func};
369
370 // Actual libpgmath declarations
371 DECLARE_PGMATH_ALL(acos)
372 DECLARE_PGMATH_MTH_VERSION_REAL(acosh)
373 DECLARE_PGMATH_ALL(asin)
374 DECLARE_PGMATH_MTH_VERSION_REAL(asinh)
375 DECLARE_PGMATH_ALL(atan)
376 DECLARE_PGMATH_REAL2(atan2)
377 DECLARE_PGMATH_MTH_VERSION_REAL(atanh)
378 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_j0)
379 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_j1)
380 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_y0)
381 DECLARE_PGMATH_MTH_VERSION_REAL(bessel_y1)
382 // bessel_jn and bessel_yn takes an int as first arg
383 extern "C" float __mth_i_bessel_jn(int, float);
384 extern "C" double __mth_i_dbessel_jn(int, double);
385 template<> constexpr auto Sym<L::F, I::bessel_jn, float>{__mth_i_bessel_jn};
386 template<> constexpr auto Sym<L::F, I::bessel_jn, double>{__mth_i_dbessel_jn};
387 template<> constexpr auto Sym<L::P, I::bessel_jn, float>{__mth_i_bessel_jn};
388 template<> constexpr auto Sym<L::P, I::bessel_jn, double>{__mth_i_dbessel_jn};
389 template<> constexpr auto Sym<L::R, I::bessel_jn, float>{__mth_i_bessel_jn};
390 template<> constexpr auto Sym<L::R, I::bessel_jn, double>{__mth_i_dbessel_jn};
391 extern "C" float __mth_i_bessel_yn(int, float);
392 extern "C" double __mth_i_dbessel_yn(int, double);
393 template<> constexpr auto Sym<L::F, I::bessel_yn, float>{__mth_i_bessel_yn};
394 template<> constexpr auto Sym<L::F, I::bessel_yn, double>{__mth_i_dbessel_yn};
395 template<> constexpr auto Sym<L::P, I::bessel_yn, float>{__mth_i_bessel_yn};
396 template<> constexpr auto Sym<L::P, I::bessel_yn, double>{__mth_i_dbessel_yn};
397 template<> constexpr auto Sym<L::R, I::bessel_yn, float>{__mth_i_bessel_yn};
398 template<> constexpr auto Sym<L::R, I::bessel_yn, double>{__mth_i_dbessel_yn};
399 DECLARE_PGMATH_ALL(cos)
400 DECLARE_PGMATH_ALL(cosh)
401 DECLARE_PGMATH_MTH_VERSION_REAL(erf)
402 DECLARE_PGMATH_MTH_VERSION_REAL(erfc)
403 DECLARE_PGMATH_MTH_VERSION_REAL(erfc_scaled)
404 DECLARE_PGMATH_ALL(exp)
405 DECLARE_PGMATH_MTH_VERSION_REAL(gamma)
406 extern "C" float __mth_i_hypot(float, float);
407 extern "C" double __mth_i_dhypot(double, double);
408 template<> constexpr auto Sym<L::F, I::hypot, float>{__mth_i_hypot};
409 template<> constexpr auto Sym<L::F, I::hypot, double>{__mth_i_dhypot};
410 template<> constexpr auto Sym<L::P, I::hypot, float>{__mth_i_hypot};
411 template<> constexpr auto Sym<L::P, I::hypot, double>{__mth_i_dhypot};
412 template<> constexpr auto Sym<L::R, I::hypot, float>{__mth_i_hypot};
413 template<> constexpr auto Sym<L::R, I::hypot, double>{__mth_i_dhypot};
414 DECLARE_PGMATH_ALL(log)
415 DECLARE_PGMATH_REAL(log10)
416 DECLARE_PGMATH_MTH_VERSION_REAL(log_gamma)
417 // no function for modulo in libpgmath
418 extern "C" float __fs_mod_1(float, float);
419 extern "C" double __fd_mod_1(double, double);
420 template<> constexpr auto Sym<L::F, I::mod, float>{__fs_mod_1};
421 template<> constexpr auto Sym<L::F, I::mod, double>{__fd_mod_1};
422 template<> constexpr auto Sym<L::P, I::mod, float>{__fs_mod_1};
423 template<> constexpr auto Sym<L::P, I::mod, double>{__fd_mod_1};
424 template<> constexpr auto Sym<L::R, I::mod, float>{__fs_mod_1};
425 template<> constexpr auto Sym<L::R, I::mod, double>{__fd_mod_1};
426 DECLARE_PGMATH_ALL2(pow)
DECLARE_PGMATH_ALL(sin)427 DECLARE_PGMATH_ALL(sin)
428 DECLARE_PGMATH_ALL(sinh)
429 DECLARE_PGMATH_MTH_VERSION_REAL(sqrt)
430 DECLARE_PGMATH_COMPLEX(sqrt) // real versions are __mth_i...
431 DECLARE_PGMATH_ALL(tan)
432 DECLARE_PGMATH_ALL(tanh)
433
434 // Fill the function map used for folding with libpgmath symbols
435 template<L Lib, typename HostT>
436 static void AddLibpgmathRealHostProcedures(
437 HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
438 static_assert(std::is_same_v<HostT, float> || std::is_same_v<HostT, double>);
439 HostRuntimeIntrinsicProcedure pgmathSymbols[]{
440 {"acos", Sym<Lib, I::acos, HostT>, true},
441 {"acosh", Sym<Lib, I::acosh, HostT>, true},
442 {"asin", Sym<Lib, I::asin, HostT>, true},
443 {"asinh", Sym<Lib, I::asinh, HostT>, true},
444 {"atan", Sym<Lib, I::atan, HostT>, true},
445 {"atan", Sym<Lib, I::atan2, HostT>,
446 true}, // atan is also the generic name for atan2
447 {"atanh", Sym<Lib, I::atanh, HostT>, true},
448 {"bessel_j0", Sym<Lib, I::bessel_j0, HostT>, true},
449 {"bessel_j1", Sym<Lib, I::bessel_j1, HostT>, true},
450 {"bessel_jn", Sym<Lib, I::bessel_jn, HostT>, true},
451 {"bessel_y0", Sym<Lib, I::bessel_y0, HostT>, true},
452 {"bessel_y1", Sym<Lib, I::bessel_y1, HostT>, true},
453 {"bessel_yn", Sym<Lib, I::bessel_yn, HostT>, true},
454 {"cos", Sym<Lib, I::cos, HostT>, true},
455 {"cosh", Sym<Lib, I::cosh, HostT>, true},
456 {"erf", Sym<Lib, I::erf, HostT>, true},
457 {"erfc", Sym<Lib, I::erfc, HostT>, true},
458 {"erfc_scaled", Sym<Lib, I::erfc_scaled, HostT>, true},
459 {"exp", Sym<Lib, I::exp, HostT>, true},
460 {"gamma", Sym<Lib, I::gamma, HostT>, true},
461 {"hypot", Sym<Lib, I::hypot, HostT>, true},
462 {"log", Sym<Lib, I::log, HostT>, true},
463 {"log10", Sym<Lib, I::log10, HostT>, true},
464 {"log_gamma", Sym<Lib, I::log_gamma, HostT>, true},
465 {"mod", Sym<Lib, I::mod, HostT>, true},
466 {"pow", Sym<Lib, I::pow, HostT>, true},
467 {"sin", Sym<Lib, I::sin, HostT>, true},
468 {"sinh", Sym<Lib, I::sinh, HostT>, true},
469 {"sqrt", Sym<Lib, I::sqrt, HostT>, true},
470 {"tan", Sym<Lib, I::tan, HostT>, true},
471 {"tanh", Sym<Lib, I::tanh, HostT>, true},
472 };
473
474 for (auto sym : pgmathSymbols) {
475 hostIntrinsicLibrary.AddProcedure(std::move(sym));
476 }
477 }
478
479 // Note: std::complex and _complex are layout compatible but are not guaranteed
480 // to be linkage compatible. For instance, on i386, float _Complex is returned
481 // by a pair of register but std::complex<float> is returned by structure
482 // address. To fix the issue, wrapper around C _Complex functions are defined
483 // below.
484 template<FuncPointer<float _Complex, float _Complex> func>
ComplexCFuncWrapper(std::complex<float> & arg)485 static std::complex<float> ComplexCFuncWrapper(std::complex<float> &arg) {
486 float _Complex res{func(*reinterpret_cast<float _Complex *>(&arg))};
487 return *reinterpret_cast<std::complex<float> *>(&res);
488 }
489
490 template<FuncPointer<double _Complex, double _Complex> func>
ComplexCFuncWrapper(std::complex<double> & arg)491 static std::complex<double> ComplexCFuncWrapper(std::complex<double> &arg) {
492 double _Complex res{func(*reinterpret_cast<double _Complex *>(&arg))};
493 return *reinterpret_cast<std::complex<double> *>(&res);
494 }
495
496 template<FuncPointer<float _Complex, float _Complex, float _Complex> func>
ComplexCFuncWrapper(std::complex<float> & arg1,std::complex<float> & arg2)497 static std::complex<float> ComplexCFuncWrapper(
498 std::complex<float> &arg1, std::complex<float> &arg2) {
499 float _Complex res{func(*reinterpret_cast<float _Complex *>(&arg1),
500 *reinterpret_cast<float _Complex *>(&arg2))};
501 return *reinterpret_cast<std::complex<float> *>(&res);
502 }
503
504 template<FuncPointer<double _Complex, double _Complex, double _Complex> func>
ComplexCFuncWrapper(std::complex<double> & arg1,std::complex<double> & arg2)505 static std::complex<double> ComplexCFuncWrapper(
506 std::complex<double> &arg1, std::complex<double> &arg2) {
507 double _Complex res{func(*reinterpret_cast<double _Complex *>(&arg1),
508 *reinterpret_cast<double _Complex *>(&arg2))};
509 return *reinterpret_cast<std::complex<double> *>(&res);
510 }
511
512 template<L Lib, typename HostT>
AddLibpgmathComplexHostProcedures(HostIntrinsicProceduresLibrary & hostIntrinsicLibrary)513 static void AddLibpgmathComplexHostProcedures(
514 HostIntrinsicProceduresLibrary &hostIntrinsicLibrary) {
515 static_assert(std::is_same_v<HostT, float> || std::is_same_v<HostT, double>);
516 using CHostT = std::complex<HostT>;
517 // cmath is used to complement pgmath when symbols are not available
518 using CmathF = FuncPointer<CHostT, const CHostT &>;
519 HostRuntimeIntrinsicProcedure pgmathSymbols[]{
520 {"abs", FuncPointer<HostT, const CHostT &>{std::abs}, true},
521 {"acos", ComplexCFuncWrapper<Sym<Lib, I::acos, CHostT>>, true},
522 {"acosh", CmathF{std::acosh}, true},
523 {"asin", ComplexCFuncWrapper<Sym<Lib, I::asin, CHostT>>, true},
524 {"asinh", CmathF{std::asinh}, true},
525 {"atan", ComplexCFuncWrapper<Sym<Lib, I::atan, CHostT>>, true},
526 {"atanh", CmathF{std::atanh}, true},
527 {"cos", ComplexCFuncWrapper<Sym<Lib, I::cos, CHostT>>, true},
528 {"cosh", ComplexCFuncWrapper<Sym<Lib, I::cosh, CHostT>>, true},
529 {"exp", ComplexCFuncWrapper<Sym<Lib, I::exp, CHostT>>, true},
530 {"log", ComplexCFuncWrapper<Sym<Lib, I::log, CHostT>>, true},
531 {"pow", ComplexCFuncWrapper<Sym<Lib, I::pow, CHostT>>, true},
532 {"sin", ComplexCFuncWrapper<Sym<Lib, I::sin, CHostT>>, true},
533 {"sinh", ComplexCFuncWrapper<Sym<Lib, I::sinh, CHostT>>, true},
534 {"sqrt", ComplexCFuncWrapper<Sym<Lib, I::sqrt, CHostT>>, true},
535 {"tan", ComplexCFuncWrapper<Sym<Lib, I::tan, CHostT>>, true},
536 {"tanh", ComplexCFuncWrapper<Sym<Lib, I::tanh, CHostT>>, true},
537 };
538
539 for (auto sym : pgmathSymbols) {
540 hostIntrinsicLibrary.AddProcedure(std::move(sym));
541 }
542 }
543
544 template<L Lib>
InitHostIntrinsicLibraryWithLibpgmath(HostIntrinsicProceduresLibrary & lib)545 static void InitHostIntrinsicLibraryWithLibpgmath(
546 HostIntrinsicProceduresLibrary &lib) {
547 if constexpr (host::FortranTypeExists<float>()) {
548 AddLibpgmathRealHostProcedures<Lib, float>(lib);
549 }
550 if constexpr (host::FortranTypeExists<double>()) {
551 AddLibpgmathRealHostProcedures<Lib, double>(lib);
552 }
553 if constexpr (host::FortranTypeExists<std::complex<float>>()) {
554 AddLibpgmathComplexHostProcedures<Lib, float>(lib);
555 }
556 if constexpr (host::FortranTypeExists<std::complex<double>>()) {
557 AddLibpgmathComplexHostProcedures<Lib, double>(lib);
558 }
559 // No long double functions in libpgmath
560 if constexpr (host::FortranTypeExists<long double>()) {
561 AddLibmRealHostProcedures<long double>(lib);
562 }
563 if constexpr (host::FortranTypeExists<std::complex<long double>>()) {
564 AddLibmComplexHostProcedures<long double>(lib);
565 }
566 }
567 }
568 #endif // LINK_WITH_LIBPGMATH
569
570 // Define which host runtime functions will be used for folding
HostIntrinsicProceduresLibrary()571 HostIntrinsicProceduresLibrary::HostIntrinsicProceduresLibrary() {
572 // TODO: When command line options regarding targeted numerical library is
573 // available, this needs to be revisited to take it into account. So far,
574 // default to libpgmath if F18 is built with it.
575 #if LINK_WITH_LIBPGMATH
576 // This looks and is stupid for now (until TODO above), but it is needed
577 // to silence clang warnings on unused symbols if all declared pgmath
578 // symbols are not used somewhere.
579 if (true) {
580 pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::P>(*this);
581 } else if (false) {
582 pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::F>(*this);
583 } else {
584 pgmath::InitHostIntrinsicLibraryWithLibpgmath<pgmath::L::R>(*this);
585 }
586 #else
587 InitHostIntrinsicLibraryWithLibm(*this);
588 #endif
589 }
590 }
591