1747 lines
50 KiB
C++
1747 lines
50 KiB
C++
/*
|
|
* ******************************************************************************
|
|
* *
|
|
* *
|
|
* * This program and the accompanying materials are made available under the
|
|
* * terms of the Apache License, Version 2.0 which is available at
|
|
* * https://www.apache.org/licenses/LICENSE-2.0.
|
|
* *
|
|
* * See the NOTICE file distributed with this work for additional
|
|
* * information regarding copyright ownership.
|
|
* * Unless required by applicable law or agreed to in writing, software
|
|
* * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
|
|
* * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
|
|
* * License for the specific language governing permissions and limitations
|
|
* * under the License.
|
|
* *
|
|
* * SPDX-License-Identifier: Apache-2.0
|
|
* *****************************************************************************
|
|
*/
|
|
|
|
/*
|
|
* templatemath.h
|
|
*
|
|
* Created on: Jan 1, 2016
|
|
* Author: agibsonccc
|
|
*/
|
|
|
|
#ifndef TEMPLATEMATH_H_
|
|
#define TEMPLATEMATH_H_
|
|
|
|
#include <system/dll.h>
|
|
#include <system/pointercast.h>
|
|
#include <math/platformmath.h>
|
|
#include <array/DataTypeUtils.h>
|
|
|
|
#define BFLOAT16_MAX_VALUE 32737.
|
|
#define HALF_MAX_VALUE 65504.
|
|
#define FLOAT_MAX_VALUE 3.4028235E38
|
|
#define DOUBLE_MAX_VALUE 1.7976931348623157E308
|
|
#define FLOAT_MIN_NORMAL 1.17549435e-38
|
|
|
|
#ifndef M_E
|
|
#define M_E 2.718281828459
|
|
#endif
|
|
|
|
#ifndef M_PI
|
|
#define M_PI 3.14159265358979323846
|
|
#endif
|
|
|
|
namespace sd {
|
|
#ifdef __CUDACC__
|
|
|
|
#endif
|
|
|
|
namespace math {
|
|
template<typename T>
|
|
math_def inline T nd4j_abs(T value);
|
|
|
|
template<typename T>
|
|
math_def inline void nd4j_swap(T &val1, T &val2);
|
|
|
|
template<typename T>
|
|
math_def inline T nd4j_max(T val1, T val2);
|
|
|
|
template<typename T>
|
|
math_def inline T nd4j_min(T val1, T val2);
|
|
|
|
template <typename T>
|
|
math_def inline bool nd4j_eq(T val1, T val2, double eps);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_re(T val1, T val2);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_rint(T val1);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_copysign(T val1, T val2);
|
|
|
|
template <typename T, typename Z>
|
|
math_def inline Z nd4j_softplus(T val);
|
|
|
|
template <typename T>
|
|
math_def inline T nd4j_rotl(T val, T shift);
|
|
|
|
template <typename T>
|
|
math_def inline T nd4j_rotr(T val, T shift);
|
|
|
|
//#ifndef __CUDACC__
|
|
template<typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_dot(X *x, Y *y, int length);
|
|
//#endif
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_ceil(T val1);
|
|
|
|
template<typename T>
|
|
math_def inline bool nd4j_isnan(T val1);
|
|
|
|
template<typename T>
|
|
math_def inline bool nd4j_isinf(T val1);
|
|
|
|
template<typename T>
|
|
math_def inline bool nd4j_isfin(T val1);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_cos(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_cosh(T val);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_exp(X val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_floor(T val);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_log(X val);
|
|
|
|
template<typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_pow(X val, Y val2);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_round(T val);
|
|
|
|
template<typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_remainder(X num, Y denom);
|
|
|
|
template<typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_fmod(X num, Y denom);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_erf(T num);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_erfc(T num);
|
|
|
|
math_def inline int32_t floatToRawIntBits(float d) {
|
|
union {
|
|
float f;
|
|
int32_t i;
|
|
} tmp;
|
|
tmp.f = d;
|
|
return tmp.i;
|
|
}
|
|
|
|
math_def inline float intBitsToFloat(int32_t i) {
|
|
union {
|
|
float f;
|
|
int32_t i;
|
|
} tmp;
|
|
tmp.i = i;
|
|
return tmp.f;
|
|
}
|
|
|
|
math_def inline float mulsignf(float x, float y) {
|
|
return intBitsToFloat(floatToRawIntBits(x) ^ (floatToRawIntBits(y) & (1 << 31)));
|
|
}
|
|
|
|
math_def inline float copysignfk(float x, float y) {
|
|
return intBitsToFloat((floatToRawIntBits(x) & ~(1 << 31)) ^ (floatToRawIntBits(y) & (1 << 31)));
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_sigmoid(T val) {
|
|
return (Z) 1.0f / ((Z) 1.0f + nd4j_exp<T, Z>(-val));
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_elu(T val, T alpha) {
|
|
if (val >= (T) 0.f)
|
|
return val;
|
|
return static_cast<Z>(alpha) * (nd4j_exp<T, Z>(val) - static_cast<Z>(1.0f));
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_leakyrelu(T val,T alpha) {
|
|
if (val < (T) 0.0f)
|
|
return alpha * val;
|
|
else
|
|
return val;
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_eluderivative(T val, T alpha) {
|
|
if (val >= static_cast<T>(0.0f))
|
|
return static_cast<Z>(1.0f);
|
|
return static_cast<Z>(alpha) * nd4j_exp<T, Z>(val);
|
|
//return val >= 0.0 ? 1.0 : nd4j_exp(val);
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_sin(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_sinh(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_softplus(T val) {
|
|
return nd4j_log<T, Z>((Z) 1.0f + nd4j_exp<T, Z>(val));
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_softsign(T val) {
|
|
return val / ((T) 1.0f + sd::math::nd4j_abs<T>(val));
|
|
}
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_sqrt(X val);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_tanh(X val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_tan(T val);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_atan2(X val1, X val2);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_atan2(X val1, X val2) {
|
|
return p_atan2<Z>(static_cast<Z>(val1), static_cast<Z>(val2));
|
|
}
|
|
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_tan(T tval) {
|
|
return p_tan<Z>(static_cast<Z>(tval));
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_tanhderivative(T val) {
|
|
Z tanh = nd4j_tanh<T,Z>(val);
|
|
return (Z) 1.0f - tanh * tanh;
|
|
}
|
|
template <typename T, typename Z>
|
|
math_def inline T nd4j_sigmoidderivative(T val) {
|
|
Z sigmoid = nd4j_sigmoid<T,Z>(val);
|
|
return sigmoid * ((Z) 1.0f - sigmoid);
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline T nd4j_softsignderivative(T val) {
|
|
T y = (T) 1.0f + nd4j_abs(val);
|
|
return (Z) 1.0f / (y * y);
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline T nd4j_sgn(T val) {
|
|
return val < (T) 0.0f ? (Z) -1.0f : val > (T) 0.0f ? (Z) 1.0f : (Z) 0.0f;
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_sign(T val) {
|
|
return nd4j_sgn<T, Z>(val);
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_signum(T val) {
|
|
return nd4j_sgn<T, Z>(val);
|
|
}
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_gamma(X a);
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_lgamma(X x);
|
|
|
|
//#ifndef __CUDACC__
|
|
/*
|
|
template<>
|
|
math_def inline float16 nd4j_dot<float16>(float16 *x, float16 *y, int length) {
|
|
float16 dot = (float16) 0.0f;
|
|
|
|
// TODO: since we can't use simd on unions, we might use something else here.
|
|
for(int e = 0; e < length; e++) {
|
|
dot += x[e] * y[e];
|
|
}
|
|
|
|
return dot;
|
|
}
|
|
*/
|
|
|
|
template<typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_dot(X *x, Y *y, int length) {
|
|
Z dot = (Z)0.0f;
|
|
|
|
for(int e = 0; e < length; e++) {
|
|
dot += static_cast<Z>(x[e]) * static_cast<Z>(y[e]);
|
|
}
|
|
|
|
return dot;
|
|
}
|
|
//#endif
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_acos(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_sech(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_acosh(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_asin(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_asinh(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_asinh(T val) {
|
|
//Math.log(Math.sqrt(Math.pow(x, 2) + 1) + x)
|
|
return nd4j_log<Z, Z>(nd4j_sqrt<Z, Z>(nd4j_pow<T,T,Z>(val, (T) 2) + (Z) 1.f) + (Z) val);
|
|
}
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_atan(T val);
|
|
|
|
template<typename T, typename Z>
|
|
math_def inline Z nd4j_atanh(T val);
|
|
|
|
|
|
template<>
|
|
math_def inline float16 nd4j_abs<float16>(float16 value) {
|
|
#ifdef NATIVE_HALFS
|
|
if (value < (float16) 0.f) {
|
|
return float16(__hneg(value.data));
|
|
} else
|
|
return value;
|
|
#else
|
|
return (float16) fabsf((float) value);
|
|
#endif
|
|
}
|
|
template<>
|
|
math_def inline bfloat16 nd4j_abs<bfloat16>(bfloat16 value) {
|
|
return (bfloat16) fabsf((float) value);
|
|
}
|
|
template<>
|
|
math_def inline float nd4j_abs<float>(float value) {
|
|
return fabsf(value);
|
|
}
|
|
|
|
template<>
|
|
math_def inline double nd4j_abs<double>(double value) {
|
|
return fabs(value);
|
|
}
|
|
|
|
template<>
|
|
math_def inline int nd4j_abs<int>(int value) {
|
|
return abs(value);
|
|
}
|
|
|
|
template<>
|
|
math_def inline Nd4jLong nd4j_abs<Nd4jLong>(Nd4jLong value) {
|
|
return llabs(value);
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_abs<bool>(bool value) {
|
|
return value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline uint8_t nd4j_abs<uint8_t>(uint8_t value) {
|
|
return value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline uint16_t nd4j_abs<uint16_t>(uint16_t value) {
|
|
return value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline uint32_t nd4j_abs<uint32_t>(uint32_t value) {
|
|
return value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline Nd4jULong nd4j_abs<Nd4jULong>(Nd4jULong value) {
|
|
return value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline int8_t nd4j_abs<int8_t>(int8_t value) {
|
|
return value < 0 ? -value : value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline int16_t nd4j_abs<int16_t>(int16_t value) {
|
|
return value < 0 ? -value : value;
|
|
}
|
|
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<float16>(float16 value) {
|
|
return *(value.data.getXP()) == 0x7fffU;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<bfloat16>(bfloat16 value) {
|
|
return value == bfloat16::nan(); //0x7fffU;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<float>(float value) {
|
|
return value != value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<double>(double value) {
|
|
return value != value;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<int>(int value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<uint32_t>(uint32_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<uint16_t>(uint16_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<uint8_t>(uint8_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<int16_t>(int16_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<int8_t>(int8_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<bool>(bool value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<Nd4jLong>(Nd4jLong value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isnan<Nd4jULong>(Nd4jULong value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<float16>(float16 value) {
|
|
return value < (float16) -HALF_MAX_VALUE || value > (float16) HALF_MAX_VALUE;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<bfloat16>(bfloat16 value) {
|
|
return value < (bfloat16) -BFLOAT16_MAX_VALUE || value > (bfloat16) BFLOAT16_MAX_VALUE;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<float>(float value) {
|
|
#ifdef __CUDACC__
|
|
return isinf(value);
|
|
#else
|
|
return std::isinf(value);
|
|
#endif
|
|
//return value < -FLOAT_MAX_VALUE || value > FLOAT_MAX_VALUE;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<double>(double value) {
|
|
#ifdef __CUDACC__
|
|
return isinf(value);
|
|
#else
|
|
return std::isinf(value);
|
|
#endif
|
|
//return value < -DOUBLE_MAX_VALUE || value > DOUBLE_MAX_VALUE;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<int>(int value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<uint32_t>(uint32_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<uint16_t>(uint16_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<uint8_t>(uint8_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<int16_t>(int16_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<int8_t>(int8_t value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<bool>(bool value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<Nd4jLong>(Nd4jLong value) {
|
|
return false;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_isinf<Nd4jULong>(Nd4jULong value) {
|
|
return false;
|
|
}
|
|
|
|
template<typename T>
|
|
math_def inline bool nd4j_isfin(T value) {
|
|
return !nd4j_isnan<T>(value) && !nd4j_isinf<T>(value);
|
|
}
|
|
|
|
template<>
|
|
math_def inline float16 nd4j_copysign<float16>(float16 val1, float16 val2) {
|
|
return (float16) copysignf((float) val1, (float) val2);
|
|
}
|
|
|
|
template<>
|
|
math_def inline float nd4j_copysign<float>(float val1, float val2) {
|
|
return copysignf(val1, val2);
|
|
}
|
|
|
|
template<>
|
|
math_def inline double nd4j_copysign<double>(double val1, double val2) {
|
|
return copysign(val1, val2);
|
|
}
|
|
|
|
template<>
|
|
math_def inline int nd4j_copysign<int>(int val1, int val2) {
|
|
if (val2 < 0) return -(nd4j_abs<int>(val1));
|
|
else return nd4j_abs<int>(val1);
|
|
}
|
|
|
|
template<>
|
|
math_def inline Nd4jLong nd4j_copysign<Nd4jLong>(Nd4jLong val1, Nd4jLong val2) {
|
|
if (val2 < 0) return -(nd4j_abs<Nd4jLong>(val1));
|
|
else return nd4j_abs<Nd4jLong>(val1);
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_max(bool val1, bool val2) {
|
|
return (val1 || val2) ? true : false;
|
|
}
|
|
|
|
template<typename T>
|
|
math_def inline T nd4j_max(T val1, T val2) {
|
|
return val1 > val2 ? val1 : val2;
|
|
}
|
|
|
|
template<>
|
|
math_def inline bool nd4j_min(bool val1, bool val2) {
|
|
return (val1 && val2) ? true : false;
|
|
}
|
|
|
|
template<typename T>
|
|
math_def inline T nd4j_min(T val1, T val2) {
|
|
return val1 < val2 ? val1 : val2;
|
|
}
|
|
|
|
template <typename T>
|
|
math_def inline bool nd4j_eq(T d1, T d2, double eps) {
|
|
if (sd::math::nd4j_isinf<T>(d1) && sd::math::nd4j_isinf<T>(d2)) {
|
|
if (d1 > 0 && d2 > 0)
|
|
return true;
|
|
else if (d1 < 0 && d2 < 0)
|
|
return true;
|
|
else
|
|
return false;
|
|
}
|
|
|
|
auto diff = static_cast<double>(sd::math::nd4j_abs<T>(d1 - d2));
|
|
|
|
|
|
// works well except in the range of very large numbers
|
|
if (diff <= eps)
|
|
return true;
|
|
|
|
// Knuth approach
|
|
// works well except in the range of very small numbers
|
|
if (diff <= sd::math::nd4j_max<double>(sd::math::nd4j_abs<double>(static_cast<double>(d1)), sd::math::nd4j_abs<double>(static_cast<double>(d2))) * eps)
|
|
return true;
|
|
|
|
return false;
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_ceil(X val) {
|
|
return static_cast<Z>(p_ceil<X>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_round(X val) {
|
|
return static_cast<Z>(p_round<X>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_asin(X val) {
|
|
return p_asin<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_atan(X val) {
|
|
return p_atan<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_atanh(X val) {
|
|
return p_atanh<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_cosh(X val) {
|
|
return p_cosh<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_rint(X val) {
|
|
return p_rint<X>(val);
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_sinh(X val) {
|
|
return p_sinh<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_acos(X val) {
|
|
return p_acos<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_sech(X val) {
|
|
return static_cast<Z>(1) / nd4j_cosh<X,Z>(val);
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_acosh(X val) {
|
|
return p_acosh<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_cos(X val) {
|
|
return p_cos<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_exp(X val) {
|
|
return p_exp<X>(val);
|
|
}
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_floor(X val) {
|
|
return static_cast<Z>(p_floor<X>(val));
|
|
}
|
|
|
|
template<typename X, typename Z>
|
|
math_def inline Z nd4j_log(X val) {
|
|
return static_cast<Z>(p_log<X>(val));
|
|
}
|
|
|
|
/**
|
|
* This func is special case - it must return floating point value, and optionally Y arg can be floating point argument
|
|
* @tparam X
|
|
* @tparam Y
|
|
* @tparam Z
|
|
* @param val
|
|
* @param val2
|
|
* @return
|
|
*/
|
|
template <>
|
|
math_def inline float nd4j_pow(float val, float val2) {
|
|
return p_pow<float>(val, val2);
|
|
}
|
|
|
|
template <typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_pow(X val, Y val2) {
|
|
return p_pow<Z>(static_cast<Z>(val), static_cast<Z>(val2));
|
|
}
|
|
|
|
/**
|
|
* LogGamma(a) - float point extension of ln(n!)
|
|
**/
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_lgamma(X x) {
|
|
// if (x <= X(0.0))
|
|
// {
|
|
// std::stringstream os;
|
|
// os << "Logarithm of Gamma has sence only for positive values, but " << x << " was given.";
|
|
// throw std::invalid_argument( os.str() );
|
|
// }
|
|
|
|
if (x < X(12.0)) {
|
|
return nd4j_log<Z,Z>(nd4j_gamma<X,Z>(x));
|
|
}
|
|
|
|
// Abramowitz and Stegun 6.1.41
|
|
// Asymptotic series should be good to at least 11 or 12 figures
|
|
// For error analysis, see Whittiker and Watson
|
|
// A Course in Modern Analysis (1927), page 252
|
|
|
|
static const double c[8] = {
|
|
1.0/12.0,
|
|
-1.0/360.0,
|
|
1.0/1260.0,
|
|
-1.0/1680.0,
|
|
1.0/1188.0,
|
|
-691.0/360360.0,
|
|
1.0/156.0,
|
|
-3617.0/122400.0
|
|
};
|
|
|
|
double z = Z(1.0 / Z(x * x));
|
|
double sum = c[7];
|
|
|
|
for (int i = 6; i >= 0; i--) {
|
|
sum *= z;
|
|
sum += c[i];
|
|
}
|
|
|
|
double series = sum / Z(x);
|
|
|
|
static const double halfLogTwoPi = 0.91893853320467274178032973640562;
|
|
|
|
return Z((double(x) - 0.5) * nd4j_log<X,double>(x) - double(x) + halfLogTwoPi + series);
|
|
}
|
|
|
|
|
|
|
|
template<typename T>
|
|
math_def inline T nd4j_re(T val1, T val2) {
|
|
if (val1 == (T) 0.0f && val2 == (T) 0.0f)
|
|
return (T) 0.0f;
|
|
|
|
return nd4j_abs<T>(val1 - val2) / (nd4j_abs<T>(val1) + nd4j_abs<T>(val2));
|
|
}
|
|
|
|
|
|
template <typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_remainder(X val, Y val2) {
|
|
return p_remainder<Z>(static_cast<Z>(val), static_cast<Z>(val2));
|
|
}
|
|
|
|
template <typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_fmod(X val, Y val2) {
|
|
return p_fmod<Z>(static_cast<Z>(val), static_cast<Z>(val2));
|
|
}
|
|
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_sin(X val) {
|
|
return p_sin<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_sqrt(X val) {
|
|
return p_sqrt<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
|
|
template <typename X>
|
|
math_def inline X neg_tanh(X val) {
|
|
X o = static_cast<X>(1.0f);
|
|
X t = static_cast<X>(2.0f);
|
|
X e = static_cast<X>(M_E);
|
|
|
|
auto p = sd::math::nd4j_pow<X, X, X>(e, val * t);
|
|
return (p - o)/ (p + o);
|
|
}
|
|
|
|
template <typename X>
|
|
math_def inline X pos_tanh(X val) {
|
|
X o = static_cast<X>(1.0f);
|
|
X t = static_cast<X>(-2.0f);
|
|
X e = static_cast<X>(M_E);
|
|
|
|
auto p = sd::math::nd4j_pow<X, X, X>(e, val * t);
|
|
return (o - p) / (o + p);
|
|
}
|
|
|
|
|
|
math_def inline float neu_tanh(float val, float sign) {
|
|
float e(M_E);
|
|
float av = sign * val;
|
|
auto p = sd::math::nd4j_pow<float, float, float>(e, -av * 2.f);
|
|
return (1 - p) / (1 + p);
|
|
}
|
|
|
|
template <>
|
|
math_def inline float nd4j_tanh(float val) {
|
|
float sign = copysignfk(1.0f, val);
|
|
return sign * neu_tanh(val, sign);
|
|
}
|
|
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_tanh(X val) {
|
|
return val <= 0 ? neg_tanh(val) : pos_tanh(val);
|
|
}
|
|
|
|
template <typename T>
|
|
math_def inline T nd4j_rotl(T val, T shift) {
|
|
return p_rotl<T>(val, shift);
|
|
}
|
|
|
|
template <typename T>
|
|
math_def inline T nd4j_rotr(T val, T shift) {
|
|
return p_rotr<T>(val, shift);
|
|
}
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_erf(X val) {
|
|
return p_erf<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_erfc(X val) {
|
|
return p_erfc<Z>(static_cast<Z>(val));
|
|
}
|
|
|
|
template<typename T>
|
|
math_def inline void nd4j_swap(T &val1, T &val2) {
|
|
T temp = val1; val1=val2; val2=temp;
|
|
};
|
|
|
|
template <typename X, typename Z>
|
|
math_def inline Z nd4j_gamma(X a) {
|
|
// nd4j_lgamma<X,Z>(a);
|
|
// return (Z)std::tgamma(a);
|
|
// Split the function domain into three intervals:
|
|
// (0, 0.001), [0.001, 12), and (12, infinity)
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
// First interval: (0, 0.001)
|
|
//
|
|
// For small a, 1/Gamma(a) has power series a + gamma a^2 - ...
|
|
// So in this range, 1/Gamma(a) = a + gamma a^2 with error on the order of a^3.
|
|
// The relative error over this interval is less than 6e-7.
|
|
|
|
const double eulerGamma = 0.577215664901532860606512090; // Euler's gamma constant
|
|
|
|
if (a < X(0.001))
|
|
return Z(1.0 / ((double)a * (1.0 + eulerGamma * (double)a)));
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
// Second interval: [0.001, 12)
|
|
|
|
if (a < X(12.0)) {
|
|
// The algorithm directly approximates gamma over (1,2) and uses
|
|
// reduction identities to reduce other arguments to this interval.
|
|
|
|
double y = (double)a;
|
|
int n = 0;
|
|
bool argWasLessThanOne = y < 1.0;
|
|
|
|
// Add or subtract integers as necessary to bring y into (1,2)
|
|
// Will correct for this below
|
|
if (argWasLessThanOne) {
|
|
y += 1.0;
|
|
}
|
|
else {
|
|
n = static_cast<int>(floor(y)) - 1; // will use n later
|
|
y -= n;
|
|
}
|
|
|
|
// numerator coefficients for approximation over the interval (1,2)
|
|
static const double p[] = {
|
|
-1.71618513886549492533811E+0,
|
|
2.47656508055759199108314E+1,
|
|
-3.79804256470945635097577E+2,
|
|
6.29331155312818442661052E+2,
|
|
8.66966202790413211295064E+2,
|
|
-3.14512729688483675254357E+4,
|
|
-3.61444134186911729807069E+4,
|
|
6.64561438202405440627855E+4
|
|
};
|
|
|
|
// denominator coefficients for approximation over the interval (1,2)
|
|
static const double q[] = {
|
|
-3.08402300119738975254353E+1,
|
|
3.15350626979604161529144E+2,
|
|
-1.01515636749021914166146E+3,
|
|
-3.10777167157231109440444E+3,
|
|
2.25381184209801510330112E+4,
|
|
4.75584627752788110767815E+3,
|
|
-1.34659959864969306392456E+5,
|
|
-1.15132259675553483497211E+5
|
|
};
|
|
|
|
double num = 0.0;
|
|
double den = 1.0;
|
|
|
|
|
|
double z = y - 1;
|
|
for (auto i = 0; i < 8; i++) {
|
|
num = (num + p[i]) * z;
|
|
den = den * z + q[i];
|
|
}
|
|
double result = num / den + 1.0;
|
|
|
|
// Apply correction if argument was not initially in (1,2)
|
|
if (argWasLessThanOne) {
|
|
// Use identity gamma(z) = gamma(z+1)/z
|
|
// The variable "result" now holds gamma of the original y + 1
|
|
// Thus we use y-1 to get back the orginal y.
|
|
result /= (y - 1.0);
|
|
}
|
|
else {
|
|
// Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
|
|
for (auto i = 0; i < n; i++)
|
|
result *= y++;
|
|
}
|
|
|
|
return Z(result);
|
|
}
|
|
|
|
///////////////////////////////////////////////////////////////////////////
|
|
// Third interval: [12, infinity)
|
|
|
|
if (a > 171.624) {
|
|
// Correct answer too large to display. Force +infinity.
|
|
return Z(DOUBLE_MAX_VALUE);
|
|
// return DataTypeUtils::infOrMax<Z>();
|
|
}
|
|
|
|
return sd::math::nd4j_exp<Z,Z>(sd::math::nd4j_lgamma<X,Z>(a));
|
|
}
|
|
|
|
template <typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_igamma(X a, Y x) {
|
|
Z aim = nd4j_pow<X, X, Z>(x, a) / (nd4j_exp<X, Z>(x) * nd4j_gamma<Y, Z>(a));
|
|
auto sum = Z(0.);
|
|
auto denom = Z(1.);
|
|
if (a <= X(0.000001))
|
|
//throw std::runtime_error("Cannot calculate gamma for a zero val.");
|
|
return Z(0);
|
|
|
|
for (int i = 0; Z(1./denom) > Z(1.0e-12); i++) {
|
|
denom *= (a + i);
|
|
sum += nd4j_pow<X, int, Z>(x, i) / denom;
|
|
}
|
|
return aim * sum;
|
|
}
|
|
|
|
template <typename X, typename Y, typename Z>
|
|
math_def inline Z nd4j_igammac(X a, Y x) {
|
|
return Z(1.) - nd4j_igamma<X, Y, Z>(a, x);
|
|
}
|
|
|
|
#ifdef __CUDACC__
|
|
namespace atomics {
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicAdd(T* address, T val);
|
|
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicSub(T* address, T val);
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicMul(T* address, T val);
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicDiv(T* address, T val);
|
|
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicMin(T* address, T val);
|
|
template <typename T>
|
|
inline __device__ T nd4j_atomicMax(T* address, T val);
|
|
|
|
template <>
|
|
inline __device__ int32_t nd4j_atomicMin<int32_t>(int32_t* address, int32_t val) {
|
|
return atomicMin(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint32_t nd4j_atomicMin<uint32_t>(uint32_t* address, uint32_t val) {
|
|
return atomicMin(address, val);
|
|
}
|
|
template <>
|
|
inline __device__ float nd4j_atomicMin<float>(float* address, float val) {
|
|
int* address_as_ull = (int*)address;
|
|
int old = __float_as_int(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, __float_as_int(math::nd4j_min(val, __int_as_float(assumed))));
|
|
} while (assumed != old);
|
|
return __int_as_float(old);
|
|
}
|
|
template <>
|
|
inline __device__ double nd4j_atomicMin<double>(double* address, double val) {
|
|
unsigned long long int* address_as_ull = (unsigned long long int*)address;
|
|
unsigned long long int old = __double_as_longlong(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, __double_as_longlong(math::nd4j_min(val, __longlong_as_double(assumed))));
|
|
} while (assumed != old);
|
|
return __longlong_as_double(old);
|
|
}
|
|
template <>
|
|
inline __device__ uint64_t nd4j_atomicMin<uint64_t>(uint64_t* address, uint64_t val) {
|
|
#if __CUDA_ARCH__ >= 350
|
|
return atomicMin((unsigned long long*)address, (unsigned long long)val);
|
|
#else
|
|
unsigned long long int* address_as_ull = (unsigned long long int*)address;
|
|
unsigned long long int old = __double_as_longlong(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, math::nd4j_min((unsigned long long)val, assumed));
|
|
} while (assumed != old);
|
|
return old;
|
|
#endif
|
|
}
|
|
template <>
|
|
inline __device__ Nd4jLong nd4j_atomicMin<Nd4jLong>(Nd4jLong* address, Nd4jLong val) {
|
|
|
|
#if __CUDA_ARCH__ >= 350
|
|
return atomicMin((unsigned long long*)address, (unsigned long long)val);
|
|
#else
|
|
unsigned long long int* address_as_ull = (unsigned long long int*)address;
|
|
unsigned long long int old = (unsigned long long)val, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, math::nd4j_min(val, (Nd4jLong)assumed));
|
|
} while (assumed != old);
|
|
return old;
|
|
#endif
|
|
|
|
}
|
|
template <>
|
|
inline __device__ int16_t nd4j_atomicMin<int16_t>(int16_t* address, int16_t val) {
|
|
int32_t temp = *address;
|
|
*address = atomicMin(&temp, (int)val);
|
|
return *address;
|
|
}
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicMin<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
return bfloat16(nd4j_atomicMin<int16_t>(&address->_data, val._data));
|
|
}
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicMin<float16>(float16* address, float16 val) {
|
|
return float16(nd4j_atomicMin<int16_t>(reinterpret_cast<int16_t*>(&address->data), (int16_t)val.data));
|
|
}
|
|
template <>
|
|
inline __device__ int32_t nd4j_atomicMax<int32_t>(int32_t* address, int32_t val) {
|
|
return atomicMax(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint32_t nd4j_atomicMax<uint32_t>(uint32_t* address, uint32_t val) {
|
|
return atomicMax(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ double nd4j_atomicMax<double>(double* address, double val) {
|
|
unsigned long long int* address_as_ull = (unsigned long long int*)address;
|
|
unsigned long long int old = __double_as_longlong(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, __double_as_longlong(math::nd4j_max(val, __longlong_as_double(assumed))));
|
|
} while (assumed != old);
|
|
return __longlong_as_double(old);
|
|
}
|
|
template <>
|
|
inline __device__ float nd4j_atomicMax<float>(float* address, float val) {
|
|
int* address_as_ull = (int*)address;
|
|
int old = __float_as_int(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, __float_as_int(math::nd4j_max(val, __int_as_float(assumed))));
|
|
} while (assumed != old);
|
|
return __int_as_float(old);
|
|
}
|
|
template <>
|
|
inline __device__ uint8_t nd4j_atomicMin<uint8_t>(uint8_t* address, uint8_t val) {
|
|
uint32_t temp = *address;
|
|
*address = atomicMin(&temp, (uint32_t)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int8_t nd4j_atomicMin<int8_t>(int8_t* address, int8_t val) {
|
|
int32_t temp = *address;
|
|
*address = atomicMin(&temp, (int)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint16_t nd4j_atomicMin<uint16_t>(uint16_t* address, uint16_t val) {
|
|
uint32_t temp = *address;
|
|
*address = atomicMin(&temp, (uint32_t)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint8_t nd4j_atomicMax<uint8_t>(uint8_t* address, uint8_t val) {
|
|
uint32_t temp = *address;
|
|
*address = atomicMax(&temp, (uint32_t)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int8_t nd4j_atomicMax<int8_t>(int8_t* address, int8_t val) {
|
|
int32_t temp = *address;
|
|
*address = atomicMax(&temp, (int)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint16_t nd4j_atomicMax<uint16_t>(uint16_t* address, uint16_t val) {
|
|
uint32_t temp = *address;
|
|
*address = atomicMax(&temp, (uint32_t)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int16_t nd4j_atomicMax<int16_t>(int16_t* address, int16_t val) {
|
|
int32_t temp = *address;
|
|
*address = atomicMax(&temp, (int32_t)val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicMax<float16>(float16* address, float16 val) {
|
|
auto address_as_ull = (int*) address;
|
|
|
|
long addr = (long) address;
|
|
bool misaligned = addr & 0x3;
|
|
|
|
if (misaligned)
|
|
address_as_ull = (int *) (address - 1);
|
|
|
|
PAIR old, assumed, fresh;
|
|
|
|
old.W = *address_as_ull;
|
|
do {
|
|
|
|
if (!misaligned) {
|
|
float16 res = nd4j_max((float16) old.B.H, val);
|
|
fresh.B.H = res.data;
|
|
fresh.B.L = old.B.L;
|
|
} else {
|
|
float16 res = nd4j_max((float16) old.B.L, val);
|
|
fresh.B.L = res.data;
|
|
fresh.B.H = old.B.H;
|
|
}
|
|
|
|
assumed.W = old.W;
|
|
old.W = atomicCAS(address_as_ull, assumed.W, fresh.W);
|
|
} while (assumed.W != old.W);
|
|
|
|
if (!misaligned) return old.B.H;
|
|
else return old.B.L;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicMax<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
auto address_as_ull = (int*) address;
|
|
|
|
long addr = (long)(address);
|
|
bool misaligned = addr & 0x3;
|
|
|
|
if (misaligned)
|
|
address_as_ull = (int *) (address - 1);
|
|
|
|
BPAIR old, assumed, fresh;
|
|
|
|
old.W = *address_as_ull;
|
|
do {
|
|
|
|
if (!misaligned) {
|
|
bfloat16 res = nd4j_max(old.B.H, val);
|
|
fresh.B.H = res;
|
|
fresh.B.L = old.B.L;
|
|
} else {
|
|
bfloat16 res = nd4j_max(old.B.L, val);
|
|
fresh.B.L = res;
|
|
fresh.B.H = old.B.H;
|
|
}
|
|
|
|
assumed.W = old.W;
|
|
old.W = atomicCAS(address_as_ull, assumed.W, fresh.W);
|
|
} while (assumed.W != old.W);
|
|
|
|
if (!misaligned) return old.B.H;
|
|
else return old.B.L;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint64_t nd4j_atomicMax<uint64_t>(uint64_t* address, uint64_t val) {
|
|
#if __CUDA_ARCH__ >= 350
|
|
return atomicMax((unsigned long long*)address, (unsigned long long)val);
|
|
#else
|
|
unsigned long long int* address_as_ull = (unsigned long long int*)address;
|
|
unsigned long long int old = __double_as_longlong(val), assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, math::nd4j_max((unsigned long long)val, assumed));
|
|
} while (assumed != old);
|
|
return old;
|
|
#endif
|
|
}
|
|
|
|
template <>
|
|
inline __device__ Nd4jLong nd4j_atomicMax<Nd4jLong>(Nd4jLong* address, Nd4jLong val) {
|
|
unsigned long long int* address_as_ull = (unsigned long long int *) address;
|
|
|
|
//return (Nd4jLong) atomicAdd(address_as_ull, (unsigned long long int) val);
|
|
unsigned long long int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, (unsigned long long)nd4j_max(val, (Nd4jLong)assumed));
|
|
} while (assumed != old);
|
|
return old;
|
|
}
|
|
|
|
|
|
template <>
|
|
inline __device__ double nd4j_atomicAdd<double>(double* address, double val) {
|
|
unsigned long long int* address_as_ull =
|
|
(unsigned long long int *) address;
|
|
unsigned long long int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed,__double_as_longlong(val +
|
|
__longlong_as_double(assumed)));
|
|
} while (assumed != old);
|
|
return __longlong_as_double(old);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ Nd4jLong nd4j_atomicAdd<Nd4jLong>(Nd4jLong* address, Nd4jLong val) {
|
|
unsigned long long int* address_as_ull = (unsigned long long int *) address;
|
|
|
|
//return (Nd4jLong) atomicAdd(address_as_ull, (unsigned long long int) val);
|
|
unsigned long long int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, val + assumed);
|
|
} while (assumed != old);
|
|
return old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ long nd4j_atomicAdd<long>(long* address, long val) {
|
|
unsigned long long* address_as_ull = (unsigned long long int *) address;
|
|
|
|
// return atomicAdd(address, val);
|
|
unsigned long int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, val + assumed);
|
|
} while (assumed != old);
|
|
return old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint32_t nd4j_atomicAdd<uint32_t>(uint32_t* address, uint32_t val) {
|
|
return atomicAdd(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint64_t nd4j_atomicAdd<uint64_t>(uint64_t* address, uint64_t val) {
|
|
// unsigned long long* address_as_ull = (unsigned long long int *) address;
|
|
//
|
|
//// return atomicAdd(address, val);
|
|
// unsigned long int old = *address_as_ull, assumed;
|
|
// do {
|
|
// assumed = old;
|
|
// old = atomicCAS(address_as_ull, assumed, val + assumed);
|
|
// } while (assumed != old);
|
|
// return old;
|
|
return (uint64_t)atomicAdd((unsigned long long*)address, (unsigned long long)val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicAdd<float16>(float16* address, float16 val) {
|
|
#if __CUDA_ARCH__ >= 700 && CUDA_VERSION_MAJOR >=10
|
|
atomicAdd(reinterpret_cast<__half*>(address), val.data);
|
|
#else
|
|
auto address_as_ull = (int*) address;
|
|
|
|
long addr = (long) address;
|
|
bool misaligned = addr & 0x3;
|
|
|
|
if (misaligned)
|
|
address_as_ull = (int *) (address - 1);
|
|
|
|
PAIR old, assumed, fresh;
|
|
|
|
old.W = *address_as_ull;
|
|
do {
|
|
|
|
if (!misaligned) {
|
|
float16 res = ((float16) old.B.H) + val;
|
|
fresh.B.H = res.data;
|
|
fresh.B.L = old.B.L;
|
|
} else {
|
|
float16 res = ((float16) old.B.L) + val;
|
|
fresh.B.L = res.data;
|
|
fresh.B.H = old.B.H;
|
|
}
|
|
|
|
assumed.W = old.W;
|
|
old.W = atomicCAS(address_as_ull, assumed.W, fresh.W);
|
|
} while (assumed.W != old.W);
|
|
|
|
if (!misaligned) return old.B.H;
|
|
else return old.B.L;
|
|
#endif
|
|
}
|
|
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicAdd<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
auto address_as_ull = (int*) address;
|
|
|
|
auto addr = (long)(address);
|
|
bool misaligned = addr & 0x3;
|
|
|
|
if (misaligned)
|
|
address_as_ull = (int *) (address - 1);
|
|
|
|
BPAIR old, assumed, fresh;
|
|
|
|
old.W = *address_as_ull;
|
|
do {
|
|
|
|
if (!misaligned) {
|
|
bfloat16 res = old.B.H + val;
|
|
fresh.B.H = res;
|
|
fresh.B.L = old.B.L;
|
|
} else {
|
|
bfloat16 res = old.B.L + val;
|
|
fresh.B.L = res;
|
|
fresh.B.H = old.B.H;
|
|
}
|
|
|
|
assumed.W = old.W;
|
|
old.W = atomicCAS(address_as_ull, assumed.W, fresh.W);
|
|
} while (assumed.W != old.W);
|
|
|
|
if (!misaligned) return old.B.H;
|
|
else return old.B.L;
|
|
}
|
|
|
|
template <typename T>
|
|
static inline __device__ T internal_16bit_atomicAdd(T* address, T val) {
|
|
size_t shift = ((size_t)address & 2);
|
|
int *base_address = (int *)((char*)address - shift);
|
|
|
|
union I16PAIR {
|
|
struct {
|
|
T H;
|
|
T L;
|
|
} B;
|
|
int W;
|
|
|
|
__host__ __device__
|
|
I16PAIR() {};
|
|
|
|
__host__ __device__
|
|
~I16PAIR() {};
|
|
};
|
|
|
|
I16PAIR pairNew, pairOld, pairAssumed;
|
|
|
|
if (reinterpret_cast<int*>(address) == base_address) {
|
|
pairOld.B.L = val;
|
|
do {
|
|
|
|
pairNew.B.L = pairOld.B.L;
|
|
pairNew.B.H = pairOld.B.H + val;
|
|
pairAssumed.W = pairOld.W;
|
|
|
|
pairOld.W = atomicCAS(base_address, pairAssumed.W, pairNew.W);
|
|
} while (pairAssumed.W != pairOld.W);
|
|
|
|
return (T) pairOld.B.H;
|
|
} else {
|
|
pairOld.B.H = val;
|
|
do {
|
|
|
|
pairNew.B.H = pairOld.B.H;
|
|
pairNew.B.L = pairOld.B.L + val;
|
|
pairAssumed.W = pairOld.W;
|
|
pairOld.W = atomicCAS(base_address, pairAssumed.W, pairNew.W);
|
|
|
|
} while (pairAssumed.W != pairOld.W);
|
|
|
|
return (T) pairOld.B.L;
|
|
}
|
|
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int16_t nd4j_atomicAdd<int16_t>(int16_t* address, int16_t val) {
|
|
return internal_16bit_atomicAdd<int16_t>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint16_t nd4j_atomicAdd<uint16_t>(uint16_t* address, uint16_t val) {
|
|
return internal_16bit_atomicAdd<uint16_t>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int8_t nd4j_atomicAdd<int8_t>(int8_t* address, int8_t val) {
|
|
int res = *address;
|
|
atomicAdd(&res, (int)val);
|
|
*address = res;
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint8_t nd4j_atomicAdd<uint8_t>(uint8_t* address, uint8_t val) {
|
|
int res = *address;
|
|
atomicAdd(&res, (int)val);
|
|
*address = res;
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ bool nd4j_atomicAdd<bool>(bool* address, bool val) {
|
|
*address += (val);
|
|
return *address;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ double nd4j_atomicSub<double>(double* address, double val) {
|
|
return nd4j_atomicAdd<double>(address, -val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ double nd4j_atomicMul<double>(double* address, double val) {
|
|
unsigned long long int* address_as_ull =
|
|
(unsigned long long int*) address;
|
|
unsigned long long int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed,__double_as_longlong(val *
|
|
__longlong_as_double(assumed)));
|
|
} while (assumed != old);
|
|
return __longlong_as_double(old);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ double nd4j_atomicDiv<double>(double* address, double val) {
|
|
return nd4j_atomicMul<double>(address, 1./val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float nd4j_atomicAdd<float>(float* address, float val) {
|
|
return atomicAdd(address,val);
|
|
}
|
|
//template <>
|
|
//inline __device__ int nd4j_atomicAdd<int>(int* address, int val) {
|
|
// return atomicAdd(address, val);
|
|
//}
|
|
template <>
|
|
inline __device__ int32_t nd4j_atomicAdd<int32_t>(int32_t* address, int32_t val) {
|
|
return (int32_t)atomicAdd((int*)address, (int)val);
|
|
}
|
|
|
|
|
|
template <>
|
|
inline __device__ float nd4j_atomicSub<float>(float* address, float val) {
|
|
return nd4j_atomicAdd<float>(address, -val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicSub<float16>(float16* address, float16 val) {
|
|
return nd4j_atomicAdd<float16>(address, -val);
|
|
}
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicSub<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
return nd4j_atomicAdd<bfloat16>(address, -val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float nd4j_atomicMul<float>(float* address, float val) {
|
|
int* address_as_ull =
|
|
( int*)address;
|
|
int old = *address_as_ull, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(address_as_ull, assumed, __float_as_int(val *
|
|
__int_as_float(assumed)));
|
|
} while (assumed != old);
|
|
return __int_as_float(old);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int8_t nd4j_atomicMul<int8_t>(int8_t* address, int8_t val) {
|
|
unsigned int *base_address = (unsigned int *)((size_t)address & ~3);
|
|
unsigned int selectors[] = {0x3214, 0x3240, 0x3410, 0x4210};
|
|
unsigned int sel = selectors[(size_t)address & 3];
|
|
unsigned int old, assumed, mul, new_;
|
|
|
|
old = *base_address;
|
|
|
|
do {
|
|
|
|
assumed = old;
|
|
mul = val * (int8_t)__byte_perm(old, 0, ((size_t)address & 3) | 0x4440);
|
|
new_ = __byte_perm(old, mul, sel);
|
|
|
|
if (new_ == old)
|
|
break;
|
|
|
|
old = atomicCAS(base_address, assumed, new_);
|
|
} while (assumed != old);
|
|
return (int8_t)old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ unsigned char nd4j_atomicMul<unsigned char>(unsigned char* address, unsigned char val) {
|
|
unsigned int *base_address = (unsigned int *)((size_t)address & ~3);
|
|
unsigned int selectors[] = {0x3214, 0x3240, 0x3410, 0x4210};
|
|
unsigned int sel = selectors[(size_t)address & 3];
|
|
unsigned int old, assumed, mul, new_;
|
|
|
|
old = *base_address;
|
|
|
|
do {
|
|
|
|
assumed = old;
|
|
mul = val * (uint8_t)__byte_perm(old, 0, ((size_t)address & 3) | 0x4440);
|
|
new_ = __byte_perm(old, mul, sel);
|
|
|
|
if (new_ == old)
|
|
break;
|
|
|
|
old = atomicCAS(base_address, assumed, new_);
|
|
} while (assumed != old);
|
|
return (uint8_t)old;
|
|
}
|
|
|
|
template <typename T>
|
|
static inline __device__ T internal_16bit_atomicMul(T* address, T val) {
|
|
size_t shift = ((size_t)address & 2);
|
|
int *base_address = (int *)((char*)address - shift);
|
|
|
|
union I16PAIR {
|
|
struct {
|
|
T H;
|
|
T L;
|
|
} B;
|
|
int W;
|
|
|
|
__host__ __device__
|
|
I16PAIR() {};
|
|
|
|
__host__ __device__
|
|
~I16PAIR() {};
|
|
};
|
|
|
|
I16PAIR pairNew, pairOld, pairAssumed;
|
|
|
|
if (reinterpret_cast<int*>(address) == base_address) {
|
|
pairOld.B.L = val;
|
|
do {
|
|
pairNew.B.L = pairOld.B.L;
|
|
pairNew.B.H = pairOld.B.H * val;
|
|
pairAssumed.W = pairOld.W;
|
|
|
|
pairOld.W = atomicCAS(base_address, pairAssumed.W, pairNew.W);
|
|
} while (pairAssumed.W != pairOld.W);
|
|
|
|
return (T) pairOld.B.H;
|
|
} else {
|
|
pairOld.B.H = val;
|
|
do {
|
|
pairNew.B.H = pairOld.B.H;
|
|
pairNew.B.L = pairOld.B.L * val;
|
|
pairAssumed.W = pairOld.W;
|
|
pairOld.W = atomicCAS(base_address, pairAssumed.W, pairNew.W);
|
|
|
|
} while (pairAssumed.W != pairOld.W);
|
|
|
|
return (T) pairOld.B.L;
|
|
}
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int16_t nd4j_atomicMul<int16_t>(int16_t* address, int16_t val) {
|
|
return internal_16bit_atomicMul<int16_t>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint16_t nd4j_atomicMul<uint16_t>(uint16_t* address, uint16_t val) {
|
|
return internal_16bit_atomicMul<uint16_t>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int nd4j_atomicMul<int>(int* address, int val) {
|
|
int* res_address = address;
|
|
int old = *res_address, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(res_address, assumed, val * assumed);
|
|
} while (assumed != old);
|
|
return old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ unsigned int nd4j_atomicMul<unsigned int>(unsigned int* address, unsigned int val) {
|
|
unsigned int* res_address = address;
|
|
unsigned int old = *res_address, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(res_address, assumed, val * assumed);
|
|
} while (assumed != old);
|
|
return old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ int64_t nd4j_atomicMul<int64_t>(int64_t* address, int64_t val) {
|
|
unsigned long long int* res_address = (unsigned long long int*)address;
|
|
unsigned long long int old = *res_address, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(res_address, assumed, val * assumed);
|
|
} while (assumed != old);
|
|
return (int64_t)old;
|
|
}
|
|
|
|
template <>
|
|
inline __device__ uint64_t nd4j_atomicMul<uint64_t>(uint64_t* address, uint64_t val) {
|
|
unsigned long long int* res_address = (unsigned long long int*)address;
|
|
unsigned long long int old = *res_address, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(res_address, assumed, val * assumed);
|
|
} while (assumed != old);
|
|
return (uint64_t)old;
|
|
}
|
|
|
|
#if !defined(_WIN32) && !defined(_WIN64)
|
|
template <>
|
|
inline __device__ Nd4jLong nd4j_atomicMul<Nd4jLong>(Nd4jLong* address, Nd4jLong val) {
|
|
unsigned long long int* res_address = (unsigned long long*)address;
|
|
unsigned long long int old = *res_address, assumed;
|
|
do {
|
|
assumed = old;
|
|
old = atomicCAS(res_address, assumed, val * assumed);
|
|
} while (assumed != old);
|
|
return (Nd4jLong)old;
|
|
}
|
|
#endif
|
|
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicMul<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
return internal_16bit_atomicMul<bfloat16>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicMul<float16>(float16* address, float16 val) {
|
|
return internal_16bit_atomicMul<float16>(address, val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float nd4j_atomicDiv<float>(float* address, float val) {
|
|
return nd4j_atomicMul<float>(address, 1.f / val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ float16 nd4j_atomicDiv<float16>(float16* address, float16 val) {
|
|
return internal_16bit_atomicMul<float16>(address, (float16) 1.f / val);
|
|
}
|
|
|
|
template <>
|
|
inline __device__ bfloat16 nd4j_atomicDiv<bfloat16>(bfloat16* address, bfloat16 val) {
|
|
return internal_16bit_atomicMul<bfloat16>(address, (bfloat16) 1 / val);
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
}
|
|
|
|
#ifdef _OPENMP
|
|
|
|
#ifndef MAX_FLOAT
|
|
#define MAX_FLOAT 1e37
|
|
#endif
|
|
|
|
#pragma omp declare reduction(maxTF : float,double,float16,bfloat16 : \
|
|
omp_out = sd::math::nd4j_max(omp_in, omp_out) )\
|
|
initializer (omp_priv=-MAX_FLOAT)
|
|
|
|
#pragma omp declare reduction(minTF : float,double,float16,bfloat16 : \
|
|
omp_out = sd::math::nd4j_min(omp_in, omp_out) )\
|
|
initializer (omp_priv=MAX_FLOAT)
|
|
|
|
#pragma omp declare reduction(maxT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = sd::math::nd4j_max(omp_in, omp_out) )\
|
|
initializer (omp_priv=0)
|
|
|
|
#pragma omp declare reduction(minT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = sd::math::nd4j_min(omp_in, omp_out) )\
|
|
initializer (omp_priv=0)
|
|
|
|
#pragma omp declare reduction(amaxT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = sd::math::nd4j_max(sd::math::nd4j_abs(omp_in), sd::math::nd4j_abs(omp_out)) )
|
|
|
|
#pragma omp declare reduction(aminT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = sd::math::nd4j_min(sd::math::nd4j_abs(omp_in), sd::math::nd4j_abs(omp_out)) )
|
|
|
|
#pragma omp declare reduction(asumT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = sd::math::nd4j_abs(omp_in) + sd::math::nd4j_abs(omp_out))\
|
|
initializer (omp_priv=0)
|
|
|
|
#pragma omp declare reduction(sumT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = omp_in + omp_out)\
|
|
initializer (omp_priv=0)
|
|
|
|
#pragma omp declare reduction(prodT : float,double,float16,bfloat16,int,Nd4jLong,Nd4jULong,int8_t,uint8_t,bool,int16_t,uint16_t,uint32_t : \
|
|
omp_out = omp_in * omp_out)\
|
|
initializer (omp_priv=1)
|
|
|
|
#endif
|
|
|
|
#endif /* TEMPLATEMATH_H_ */
|