[libclc] Move fma to the CLC library

libclc/CMakeLists.txt

@@ -28,6 +28,8 @@ set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS
   spirv/lib/SOURCES;
   # CLC internal libraries
   clc/lib/generic/SOURCES;
+  clc/lib/clspv/SOURCES;
+  clc/lib/spirv/SOURCES;
 )
 
 set( LIBCLC_MIN_LLVM 3.9.0 )

libclc/generic/include/math/clc_fma.h → libclc/clc/include/clc/internal/math/clc_sw_fma.h

@@ -1,11 +1,12 @@
-#define __CLC_FUNCTION __clc_fma
-#define __CLC_INTRINSIC "llvm.fma"
-#include "math/ternary_intrin.inc"
+#ifndef __CLC_INTERNAL_MATH_CLC_SW_FMA_H__
+#define __CLC_INTERNAL_MATH_CLC_SW_FMA_H__
 
-#define __FLOAT_ONLY
 #define __CLC_FUNCTION __clc_sw_fma
 #define __CLC_BODY <clc/shared/ternary_decl.inc>
+
 #include <clc/math/gentype.inc>
+
 #undef __CLC_BODY
 #undef __CLC_FUNCTION
-#undef __FLOAT_ONLY
+
+#endif // __CLC_INTERNAL_MATH_CLC_SW_FMA_H__

libclc/clc/include/clc/math/clc_fma.h

@@ -0,0 +1,12 @@
+#ifndef __CLC_MATH_CLC_FMA_H__
+#define __CLC_MATH_CLC_FMA_H__
+
+#define __CLC_FUNCTION __clc_fma
+#define __CLC_BODY <clc/shared/ternary_decl.inc>
+
+#include <clc/math/gentype.inc>
+
+#undef __CLC_BODY
+#undef __CLC_FUNCTION
+
+#endif // __CLC_MATH_CLC_FMA_H__

libclc/clc/include/clc/math/math.h

@@ -39,12 +39,12 @@
 #define PINF 0x200
 
 #if (defined __AMDGCN__ || defined __R600__) && !defined __HAS_FMAF__
-#define HAVE_HW_FMA32() (0)
+#define __CLC_HAVE_HW_FMA32() (0)
 #elif defined(CLC_SPIRV)
 bool __attribute__((noinline)) __clc_runtime_has_hw_fma32(void);
-#define HAVE_HW_FMA32() __clc_runtime_has_hw_fma32()
+#define __CLC_HAVE_HW_FMA32() __clc_runtime_has_hw_fma32()
 #else
-#define HAVE_HW_FMA32() (1)
+#define __CLC_HAVE_HW_FMA32() (1)
 #endif
 
 #define HAVE_BITALIGN() (0)

libclc/clc/lib/clspv/SOURCES

@@ -0,0 +1 @@
+math/clc_sw_fma.cl

libclc/clc/lib/clspv/math/clc_sw_fma.cl

@@ -0,0 +1,287 @@
+/*
+ * Copyright (c) 2014 Advanced Micro Devices, Inc.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ */
+
+// This version is derived from the generic fma software implementation
+// (__clc_sw_fma), but avoids the use of ulong in favor of uint2. The logic has
+// been updated as appropriate.
+
+#include <clc/clc_as_type.h>
+#include <clc/clcmacro.h>
+#include <clc/float/definitions.h>
+#include <clc/integer/clc_abs.h>
+#include <clc/integer/clc_clz.h>
+#include <clc/integer/clc_hadd.h>
+#include <clc/integer/clc_mul_hi.h>
+#include <clc/integer/definitions.h>
+#include <clc/math/clc_mad.h>
+#include <clc/math/math.h>
+#include <clc/relational/clc_isinf.h>
+#include <clc/relational/clc_isnan.h>
+#include <clc/shared/clc_max.h>
+
+struct fp {
+  uint2 mantissa;
+  int exponent;
+  uint sign;
+};
+
+static uint2 u2_set(uint hi, uint lo) {
+  uint2 res;
+  res.lo = lo;
+  res.hi = hi;
+  return res;
+}
+
+static uint2 u2_set_u(uint val) { return u2_set(0, val); }
+
+static uint2 u2_mul(uint a, uint b) {
+  uint2 res;
+  res.hi = __clc_mul_hi(a, b);
+  res.lo = a * b;
+  return res;
+}
+
+static uint2 u2_sll(uint2 val, uint shift) {
+  if (shift == 0)
+    return val;
+  if (shift < 32) {
+    val.hi <<= shift;
+    val.hi |= val.lo >> (32 - shift);
+    val.lo <<= shift;
+  } else {
+    val.hi = val.lo << (shift - 32);
+    val.lo = 0;
+  }
+  return val;
+}
+
+static uint2 u2_srl(uint2 val, uint shift) {
+  if (shift == 0)
+    return val;
+  if (shift < 32) {
+    val.lo >>= shift;
+    val.lo |= val.hi << (32 - shift);
+    val.hi >>= shift;
+  } else {
+    val.lo = val.hi >> (shift - 32);
+    val.hi = 0;
+  }
+  return val;
+}
+
+static uint2 u2_or(uint2 a, uint b) {
+  a.lo |= b;
+  return a;
+}
+
+static uint2 u2_and(uint2 a, uint2 b) {
+  a.lo &= b.lo;
+  a.hi &= b.hi;
+  return a;
+}
+
+static uint2 u2_add(uint2 a, uint2 b) {
+  uint carry = (__clc_hadd(a.lo, b.lo) >> 31) & 0x1;
+  a.lo += b.lo;
+  a.hi += b.hi + carry;
+  return a;
+}
+
+static uint2 u2_add_u(uint2 a, uint b) { return u2_add(a, u2_set_u(b)); }
+
+static uint2 u2_inv(uint2 a) {
+  a.lo = ~a.lo;
+  a.hi = ~a.hi;
+  return u2_add_u(a, 1);
+}
+
+static uint u2_clz(uint2 a) {
+  uint leading_zeroes = __clc_clz(a.hi);
+  if (leading_zeroes == 32) {
+    leading_zeroes += __clc_clz(a.lo);
+  }
+  return leading_zeroes;
+}
+
+static bool u2_eq(uint2 a, uint2 b) { return a.lo == b.lo && a.hi == b.hi; }
+
+static bool u2_zero(uint2 a) { return u2_eq(a, u2_set_u(0)); }
+
+static bool u2_gt(uint2 a, uint2 b) {
+  return a.hi > b.hi || (a.hi == b.hi && a.lo > b.lo);
+}
+
+_CLC_DEF _CLC_OVERLOAD float __clc_sw_fma(float a, float b, float c) {
+  /* special cases */
+  if (__clc_isnan(a) || __clc_isnan(b) || __clc_isnan(c) || __clc_isinf(a) ||
+      __clc_isinf(b)) {
+    return __clc_mad(a, b, c);
+  }
+
+  /* If only c is inf, and both a,b are regular numbers, the result is c*/
+  if (__clc_isinf(c)) {
+    return c;
+  }
+
+  a = __clc_flush_denormal_if_not_supported(a);
+  b = __clc_flush_denormal_if_not_supported(b);
+  c = __clc_flush_denormal_if_not_supported(c);
+
+  if (a == 0.0f || b == 0.0f) {
+    return c;
+  }
+
+  if (c == 0) {
+    return a * b;
+  }
+
+  struct fp st_a, st_b, st_c;
+
+  st_a.exponent = a == .0f ? 0 : ((__clc_as_uint(a) & 0x7f800000) >> 23) - 127;
+  st_b.exponent = b == .0f ? 0 : ((__clc_as_uint(b) & 0x7f800000) >> 23) - 127;
+  st_c.exponent = c == .0f ? 0 : ((__clc_as_uint(c) & 0x7f800000) >> 23) - 127;
+
+  st_a.mantissa =
+      u2_set_u(a == .0f ? 0 : (__clc_as_uint(a) & 0x7fffff) | 0x800000);
+  st_b.mantissa =
+      u2_set_u(b == .0f ? 0 : (__clc_as_uint(b) & 0x7fffff) | 0x800000);
+  st_c.mantissa =
+      u2_set_u(c == .0f ? 0 : (__clc_as_uint(c) & 0x7fffff) | 0x800000);
+
+  st_a.sign = __clc_as_uint(a) & 0x80000000;
+  st_b.sign = __clc_as_uint(b) & 0x80000000;
+  st_c.sign = __clc_as_uint(c) & 0x80000000;
+
+  // Multiplication.
+  // Move the product to the highest bits to maximize precision
+  // mantissa is 24 bits => product is 48 bits, 2bits non-fraction.
+  // Add one bit for future addition overflow,
+  // add another bit to detect subtraction underflow
+  struct fp st_mul;
+  st_mul.sign = st_a.sign ^ st_b.sign;
+  st_mul.mantissa = u2_sll(u2_mul(st_a.mantissa.lo, st_b.mantissa.lo), 14);
+  st_mul.exponent =
+      !u2_zero(st_mul.mantissa) ? st_a.exponent + st_b.exponent : 0;
+
+  // FIXME: Detecting a == 0 || b == 0 above crashed GCN isel
+  if (st_mul.exponent == 0 && u2_zero(st_mul.mantissa))
+    return c;
+
+// Mantissa is 23 fractional bits, shift it the same way as product mantissa
+#define C_ADJUST 37ul
+
+  // both exponents are bias adjusted
+  int exp_diff = st_mul.exponent - st_c.exponent;
+
+  st_c.mantissa = u2_sll(st_c.mantissa, C_ADJUST);
+  uint2 cutoff_bits = u2_set_u(0);
+  uint2 cutoff_mask = u2_add(u2_sll(u2_set_u(1), __clc_abs(exp_diff)),
+                             u2_set(0xffffffff, 0xffffffff));
+  if (exp_diff > 0) {
+    cutoff_bits =
+        exp_diff >= 64 ? st_c.mantissa : u2_and(st_c.mantissa, cutoff_mask);
+    st_c.mantissa =
+        exp_diff >= 64 ? u2_set_u(0) : u2_srl(st_c.mantissa, exp_diff);
+  } else {
+    cutoff_bits = -exp_diff >= 64 ? st_mul.mantissa
+                                  : u2_and(st_mul.mantissa, cutoff_mask);
+    st_mul.mantissa =
+        -exp_diff >= 64 ? u2_set_u(0) : u2_srl(st_mul.mantissa, -exp_diff);
+  }
+
+  struct fp st_fma;
+  st_fma.sign = st_mul.sign;
+  st_fma.exponent = __clc_max(st_mul.exponent, st_c.exponent);
+  if (st_c.sign == st_mul.sign) {
+    st_fma.mantissa = u2_add(st_mul.mantissa, st_c.mantissa);
+  } else {
+    // cutoff bits borrow one
+    st_fma.mantissa =
+        u2_add(u2_add(st_mul.mantissa, u2_inv(st_c.mantissa)),
+               (!u2_zero(cutoff_bits) && (st_mul.exponent > st_c.exponent)
+                    ? u2_set(0xffffffff, 0xffffffff)
+                    : u2_set_u(0)));
+  }
+
+  // underflow: st_c.sign != st_mul.sign, and magnitude switches the sign
+  if (u2_gt(st_fma.mantissa, u2_set(0x7fffffff, 0xffffffff))) {
+    st_fma.mantissa = u2_inv(st_fma.mantissa);
+    st_fma.sign = st_mul.sign ^ 0x80000000;
+  }
+
+  // detect overflow/underflow
+  int overflow_bits = 3 - u2_clz(st_fma.mantissa);
+
+  // adjust exponent
+  st_fma.exponent += overflow_bits;
+
+  // handle underflow
+  if (overflow_bits < 0) {
+    st_fma.mantissa = u2_sll(st_fma.mantissa, -overflow_bits);
+    overflow_bits = 0;
+  }
+
+  // rounding
+  uint2 trunc_mask = u2_add(u2_sll(u2_set_u(1), C_ADJUST + overflow_bits),
+                            u2_set(0xffffffff, 0xffffffff));
+  uint2 trunc_bits =
+      u2_or(u2_and(st_fma.mantissa, trunc_mask), !u2_zero(cutoff_bits));
+  uint2 last_bit =
+      u2_and(st_fma.mantissa, u2_sll(u2_set_u(1), C_ADJUST + overflow_bits));
+  uint2 grs_bits = u2_sll(u2_set_u(4), C_ADJUST - 3 + overflow_bits);
+
+  // round to nearest even
+  if (u2_gt(trunc_bits, grs_bits) ||
+      (u2_eq(trunc_bits, grs_bits) && !u2_zero(last_bit))) {
+    st_fma.mantissa =
+        u2_add(st_fma.mantissa, u2_sll(u2_set_u(1), C_ADJUST + overflow_bits));
+  }
+
+  // Shift mantissa back to bit 23
+  st_fma.mantissa = u2_srl(st_fma.mantissa, C_ADJUST + overflow_bits);
+
+  // Detect rounding overflow
+  if (u2_gt(st_fma.mantissa, u2_set_u(0xffffff))) {
+    ++st_fma.exponent;
+    st_fma.mantissa = u2_srl(st_fma.mantissa, 1);
+  }
+
+  if (u2_zero(st_fma.mantissa)) {
+    return 0.0f;
+  }
+
+  // Flating point range limit
+  if (st_fma.exponent > 127) {
+    return __clc_as_float(__clc_as_uint(INFINITY) | st_fma.sign);
+  }
+
+  // Flush denormals
+  if (st_fma.exponent <= -127) {
+    return __clc_as_float(st_fma.sign);
+  }
+
+  return __clc_as_float(st_fma.sign | ((st_fma.exponent + 127) << 23) |
+                        ((uint)st_fma.mantissa.lo & 0x7fffff));
+}
+
+_CLC_TERNARY_VECTORIZE(_CLC_DEF _CLC_OVERLOAD, float, __clc_sw_fma, float,
+                       float, float)

libclc/clc/lib/generic/SOURCES

@@ -20,12 +20,14 @@ integer/clc_upsample.cl
 math/clc_ceil.cl
 math/clc_copysign.cl
 math/clc_fabs.cl
+math/clc_fma.cl
 math/clc_floor.cl
 math/clc_frexp.cl
 math/clc_mad.cl
 math/clc_modf.cl
 math/clc_nextafter.cl
 math/clc_rint.cl
+math/clc_sw_fma.cl
 math/clc_trunc.cl
 relational/clc_all.cl
 relational/clc_any.cl

libclc/clc/lib/generic/math/clc_fma.cl

@@ -0,0 +1,6 @@
+#include <clc/internal/clc.h>
+#include <clc/internal/math/clc_sw_fma.h>
+#include <clc/math/math.h>
+
+#define __CLC_BODY <clc_fma.inc>
+#include <clc/math/gentype.inc>

libclc/clc/lib/generic/math/clc_fma.inc

@@ -0,0 +1,8 @@
+_CLC_DEF _CLC_OVERLOAD __CLC_GENTYPE __clc_fma(__CLC_GENTYPE a, __CLC_GENTYPE b,
+                                               __CLC_GENTYPE c) {
+#if __CLC_FPSIZE == 32
+  if (!__CLC_HAVE_HW_FMA32())
+    return __clc_sw_fma(a, b, c);
+#endif
+  return __builtin_elementwise_fma(a, b, c);
+}