/
common.hpp
198 lines (176 loc) · 6.18 KB
/
common.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
/** @file
* @copyright University of Warsaw
* @section LICENSE
* GPLv3+ (see the COPYING file or http://www.gnu.org/licenses/)
*/
#pragma once
#include <libmpdata++/formulae/idxperm.hpp>
#include <libmpdata++/opts.hpp>
namespace libmpdataxx
{
namespace formulae
{
// helper to cast floating literals to correct precision based on the blitz array underlaying type
template<class arr_t>
constexpr auto fconst(const double v)
{
return static_cast<typename arr_t::T_numtype>(v);
}
// overloads of abs/min/max/where that pick out the correct version based on ix_t
template<class ix_t, class arg_t>
forceinline_macro auto abs(const arg_t &a, typename std::enable_if<std::is_same<ix_t, int>::value>::type* = 0)
{
return std::abs(a);
}
template<class ix_t, class arg_t>
forceinline_macro auto abs(const arg_t &a, typename std::enable_if<std::is_same<ix_t, rng_t>::value>::type* = 0)
{
return blitz::abs(a);
}
template<class ix_t, class a_t, class b_t>
forceinline_macro auto min(const a_t &a, const b_t &b, typename std::enable_if<std::is_same<ix_t, int>::value>::type* = 0)
{
using cm_t = typename std::common_type<a_t, b_t>::type;
return std::min(static_cast<cm_t>(a), static_cast<cm_t>(b));
}
template<class ix_t, class a_t, class b_t>
forceinline_macro auto min(const a_t &a, const b_t &b, typename std::enable_if<std::is_same<ix_t, rng_t>::value>::type* = 0)
{
return blitz::min(a, b);
}
template<class ix_t, class a_t, class b_t>
forceinline_macro auto max(const a_t &a, const b_t &b, typename std::enable_if<std::is_same<ix_t, int>::value>::type* = 0)
{
using cm_t = typename std::common_type<a_t, b_t>::type;
return std::max(static_cast<cm_t>(a), static_cast<cm_t>(b));
}
template<class ix_t, class a_t, class b_t>
forceinline_macro auto max(const a_t &a, const b_t &b, typename std::enable_if<std::is_same<ix_t, rng_t>::value>::type* = 0)
{
return blitz::max(a, b);
}
// variadic max & min
template<class ix_t, class a_t, class... b_ts>
forceinline_macro auto max(const a_t &a, const b_ts & ... bs)
{
return max<ix_t>(a, max<ix_t>(bs...));
}
template<class ix_t, class a_t, class... b_ts>
forceinline_macro auto min(const a_t &a, const b_ts & ... bs)
{
return min<ix_t>(a, min<ix_t>(bs...));
}
template<class ix_t, class arg_t>
forceinline_macro auto where(bool c, const arg_t &a, const arg_t &b, typename std::enable_if<std::is_same<ix_t, int>::value>::type* = 0)
{
return c ? a : b;
}
template<class ix_t, class c_t, class a_t, class b_t>
forceinline_macro auto where(c_t c, const a_t &a, const b_t &b, typename std::enable_if<std::is_same<ix_t, rng_t>::value>::type* = 0)
{
return blitz::where(c, a, b);
}
// nprt: implemented using min
template<opts::opts_t opts, class ix_t, class arr_t>
forceinline_macro auto negpart(
const arr_t &x,
typename std::enable_if<!opts::isset(opts, opts::npa)>::type* = 0 // enabled if npa == false
)
{
return return_helper<ix_t>(min<ix_t>(0, x));
}
// nprt: implemented using abs
template<opts::opts_t opts, class ix_t, class arr_t>
forceinline_macro auto negpart(
const arr_t &x,
typename std::enable_if<opts::isset(opts, opts::npa)>::type* = 0 // enabled if npa == true
)
{
return return_helper<ix_t>((x - abs<ix_t>(x)) / 2);
}
// pprt: implemented using max
template<opts::opts_t opts, class ix_t, class arr_t>
forceinline_macro auto pospart(
const arr_t &x,
typename std::enable_if<!opts::isset(opts, opts::npa)>::type* = 0 // enabled if npa == false
)
{
return return_helper<ix_t>(max<ix_t>(0, x));
}
// pprt: implemented using abx
template<opts::opts_t opts, class ix_t, class arr_t>
forceinline_macro auto pospart(
const arr_t &x,
typename std::enable_if<opts::isset(opts, opts::npa)>::type* = 0 // enabled if npa == true
)
{
return return_helper<ix_t>((x + abs<ix_t>(x)) / 2);
}
// 1D or ND: G = const = 1
template<opts::opts_t opts, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &,
typename std::enable_if<!opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == false
) {
return 1;
}
// 2D: G = const = 1
template<opts::opts_t opts, int d, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &, const ix_t &,
typename std::enable_if<!opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == false
) {
return 1;
}
// 3D: G = const = 1
template<opts::opts_t opts, int d, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &, const ix_t &, const ix_t &,
typename std::enable_if<!opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == false
) {
return 1;
}
// 1D on ND: G != const
template<opts::opts_t opts, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &i,
typename std::enable_if<opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == true
)
{
return return_helper<ix_t>(
G(i) + 0 // return_macro includes a call to blitz::safeToReturn() which expects an expression as an arg
);
}
// 2D: G != const
template<opts::opts_t opts, int d, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &i,
const ix_t &j,
typename std::enable_if<opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == true
)
{
return return_helper<ix_t>(
G(idxperm::pi<d>(i, j)) + 0
);
}
// 3D: G != const
template<opts::opts_t opts, int d, class arr_t, class ix_t>
inline auto G(
const arr_t &G,
const ix_t &i,
const ix_t &j,
const ix_t &k,
typename std::enable_if<opts::isset(opts, opts::nug)>::type* = 0 // enabled if nug == true
)
{
return return_helper<ix_t>(
G(idxperm::pi<d>(i, j, k)) + 0
);
}
} // namespace formulae
} // namespace libmpdataxx