Bfloat16

.github/workflows/windows-mingw.yml

@@ -50,7 +50,7 @@ jobs:
         mkdir build
         cd build
         cmake -G "MinGW Makefiles" -DBUILD_SHARED_LIBS=${{ matrix.config.shared }} -DCMAKE_BUILD_TYPE=${{ matrix.config.build_type }} -DGINKGO_COMPILER_FLAGS=${{ matrix.config.cflags }} ..
-        cmake --build . -j4
+        cmake --build . -j2
       shell: cmd
 
     - name: install

.github/workflows/windows-msvc-ref.yml

@@ -27,8 +27,9 @@ jobs:
       fail-fast: false
       matrix:
         config:
-        - {shared: "ON", build_type: "Debug", name: "reference/debug/shared"}
-        - {shared: "OFF", build_type: "Release", name: "reference/release/static"}
+        # Debug with half precision has the issue "library limit of 65535 objects exceeded"
+        - {shared: "ON", build_type: "Debug", name: "reference/debug/shared", half: "OFF"}
+        - {shared: "OFF", build_type: "Release", name: "reference/release/static", half: "ON"}
         # Debug static needs too much storage
         # - {shared: "OFF", build_type: "Debug", name: "reference/debug/static"}
     name: msvc/${{ matrix.config.name }}
@@ -47,7 +48,7 @@ jobs:
       run: |
         mkdir build
         cd build
-        cmake  -DBUILD_SHARED_LIBS=${{ matrix.config.shared }} -DCMAKE_CXX_FLAGS_DEBUG='/MDd /Zi /Ob1 /O1 /Od /RTC1' -DGINKGO_BUILD_CUDA=OFF -DGINKGO_BUILD_OMP=OFF ..
+        cmake  -DBUILD_SHARED_LIBS=${{ matrix.config.shared }} -DCMAKE_CXX_FLAGS_DEBUG='/MDd /Zi /Ob1 /O1 /Od /RTC1' -DGINKGO_BUILD_CUDA=OFF -DGINKGO_BUILD_OMP=OFF -DGINKGO_ENABLE_HALF=${{ matrix.config.half }}..
         cmake --build . -j4 --config ${{ matrix.config.build_type }}
         ctest . -C ${{ matrix.config.build_type }} --output-on-failure
 

CMakeLists.txt

@@ -29,6 +29,7 @@ option(GINKGO_BUILD_DOC "Generate documentation" OFF)
 option(GINKGO_FAST_TESTS "Reduces the input size for a few tests known to be time-intensive" OFF)
 option(GINKGO_TEST_NONDEFAULT_STREAM "Uses non-default streams in CUDA and HIP tests" OFF)
 option(GINKGO_MIXED_PRECISION "Instantiate true mixed-precision kernels (otherwise they will be conversion-based using implicit temporary storage)" OFF)
+option(GINKGO_ENABLE_HALF "Enable the half operation" ON)
 option(GINKGO_SKIP_DEPENDENCY_UPDATE
     "Do not update dependencies each time the project is rebuilt" ON)
 option(GINKGO_EXPORT_BUILD_DIR

accessor/cuda_helper.hpp

@@ -47,6 +47,9 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "utils.hpp"
 
 
+struct __half;
+
+
 namespace gko {
 namespace acc {
 namespace detail {
@@ -57,6 +60,16 @@ struct cuda_type {
     using type = T;
 };
 
+template <>
+struct cuda_type<gko::half> {
+    using type = __half;
+};
+
+template <>
+struct cuda_type<gko::bfloat16> {
+    using type = __nv_bfloat16;
+};
+
 // Unpack cv and reference / pointer qualifiers
 template <typename T>
 struct cuda_type<const T> {
@@ -87,7 +100,7 @@ struct cuda_type<T&&> {
 // Transform std::complex to thrust::complex
 template <typename T>
 struct cuda_type<std::complex<T>> {
-    using type = thrust::complex<T>;
+    using type = thrust::complex<typename cuda_type<T>::type>;
 };
 
 

accessor/hip_helper.hpp

@@ -47,6 +47,9 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include "utils.hpp"
 
 
+struct __half;
+
+
 namespace gko {
 namespace acc {
 namespace detail {
@@ -83,11 +86,20 @@ struct hip_type<T&&> {
     using type = typename hip_type<T>::type&&;
 };
 
+template <>
+struct hip_type<gko::half> {
+    using type = __half;
+};
+
+template <>
+struct hip_type<gko::bfloat16> {
+    using type = hip_bfloat16;
+};
 
 // Transform std::complex to thrust::complex
 template <typename T>
 struct hip_type<std::complex<T>> {
-    using type = thrust::complex<T>;
+    using type = thrust::complex<typename hip_type<T>::type>;
 };
 
 

accessor/reduced_row_major_reference.hpp

@@ -103,7 +103,7 @@ class reduced_storage
     operator=(arithmetic_type val) &&
     {
         storage_type* const GKO_ACC_RESTRICT r_ptr = ptr_;
-        *r_ptr = val;
+        *r_ptr = detail::implicit_explicit_conversion<storage_type>(val);
         return val;
     }
 
@@ -115,7 +115,8 @@ class reduced_storage
     }
 
     constexpr GKO_ACC_ATTRIBUTES arithmetic_type
-    operator=(reduced_storage&& ref) && noexcept
+        operator=(reduced_storage&& ref) &&
+        noexcept
     {
         std::move(*this) = ref.implicit_conversion();
         return *this;

accessor/reference_helper.hpp

@@ -43,8 +43,10 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 // CUDA TOOLKIT < 11 does not support constexpr in combination with
 // thrust::complex, which is why constexpr is only present in later versions
-#if defined(__CUDA_ARCH__) && defined(__CUDACC_VER_MAJOR__) && \
-    (__CUDACC_VER_MAJOR__ < 11)
+// TODO: NVC++ constexpr
+#if (defined(__CUDA_ARCH__) && defined(__CUDACC_VER_MAJOR__) && \
+     (__CUDACC_VER_MAJOR__ < 11)) ||                            \
+    (defined(__NVCOMPILER) && GINKGO_ENABLE_HALF)
 
 #define GKO_ACC_ENABLE_REFERENCE_CONSTEXPR
 

benchmark/CMakeLists.txt

@@ -24,15 +24,12 @@ function(ginkgo_benchmark_cusparse_linops type def)
 endfunction()
 
 function(ginkgo_benchmark_hipsparse_linops type def)
-    add_library(hipsparse_linops_${type} utils/hip_linops.hip.cpp)
+    set_source_files_properties(utils/hip_linops.hip.cpp PROPERTIES HIP_SOURCE_PROPERTY_FORMAT TRUE)
+    hip_add_library(hipsparse_linops_${type} utils/hip_linops.hip.cpp
+    HIPCC_OPTIONS ${GINKGO_HIPCC_OPTIONS} -D${def}
+    CLANG_OPTIONS ${GINKGO_HIP_CLANG_OPTIONS}
+    NVCC_OPTIONS ${GINKGO_HIP_NVCC_OPTIONS})
     target_compile_definitions(hipsparse_linops_${type} PUBLIC ${def})
-    EXECUTE_PROCESS(COMMAND ${HIP_PATH}/bin/hipconfig --cpp_config OUTPUT_VARIABLE HIP_CXX_FLAGS)
-    set_target_properties(hipsparse_linops_${type} PROPERTIES COMPILE_FLAGS ${HIP_CXX_FLAGS})
-    # use Thrust C++ device just for compilation, we don't use thrust::complex in the benchmarks
-    target_compile_definitions(hipsparse_linops_${type} PUBLIC -DTHRUST_DEVICE_SYSTEM=THRUST_DEVICE_SYSTEM_CPP)
-    target_include_directories(hipsparse_linops_${type} SYSTEM PRIVATE
-        ${HSA_HEADER} ${HIP_INCLUDE_DIRS}
-        ${HIPBLAS_INCLUDE_DIRS} ${HIPSPARSE_INCLUDE_DIRS})
     target_link_libraries(hipsparse_linops_${type} Ginkgo::ginkgo ${HIPSPARSE_LIBRARIES})
 endfunction()
 
@@ -79,17 +76,25 @@ function(ginkgo_add_single_benchmark_executable name use_lib_linops macro_def ty
     target_compile_options("${name}" PRIVATE ${GINKGO_COMPILER_FLAGS})
     ginkgo_benchmark_add_tuning_maybe("${name}")
     if("${use_lib_linops}")
-        if (GINKGO_BUILD_CUDA)
-            target_compile_definitions("${name}" PRIVATE HAS_CUDA=1)
-            target_link_libraries("${name}" cusparse_linops_${type})
-        endif()
-        if (GINKGO_BUILD_HIP)
-            target_compile_definitions("${name}" PRIVATE HAS_HIP=1)
-            target_link_libraries("${name}" hipsparse_linops_${type})
-        endif()
-        if (GINKGO_BUILD_DPCPP)
-            target_compile_definitions("${name}" PRIVATE HAS_DPCPP=1)
-            target_link_libraries("${name}" onemkl_linops_${type})
+        if ("${type}" STREQUAL "h")
+            # only cuda supports half currently
+            if (GINKGO_BUILD_CUDA)
+                target_compile_definitions("${name}" PRIVATE HAS_CUDA=1)
+                target_link_libraries("${name}" cusparse_linops_${type})
+            endif()
+        else()
+            if (GINKGO_BUILD_CUDA)
+                target_compile_definitions("${name}" PRIVATE HAS_CUDA=1)
+                target_link_libraries("${name}" cusparse_linops_${type})
+            endif()
+            if (GINKGO_BUILD_HIP)
+                target_compile_definitions("${name}" PRIVATE HAS_HIP=1)
+                target_link_libraries("${name}" hipsparse_linops_${type})
+            endif()
+            if (GINKGO_BUILD_DPCPP)
+                target_compile_definitions("${name}" PRIVATE HAS_DPCPP=1)
+                target_link_libraries("${name}" onemkl_linops_${type})
+            endif()
         endif()
     endif()
 endfunction(ginkgo_add_single_benchmark_executable)
@@ -119,6 +124,9 @@ if (GINKGO_BUILD_CUDA)
     ginkgo_benchmark_cusparse_linops(s GKO_BENCHMARK_USE_SINGLE_PRECISION)
     ginkgo_benchmark_cusparse_linops(z GKO_BENCHMARK_USE_DOUBLE_COMPLEX_PRECISION)
     ginkgo_benchmark_cusparse_linops(c GKO_BENCHMARK_USE_SINGLE_COMPLEX_PRECISION)
+    if (GINKGO_ENABLE_HALF)
+        ginkgo_benchmark_cusparse_linops(h GKO_BENCHMARK_USE_HALF_PRECISION)
+    endif()
     add_library(cuda_timer utils/cuda_timer.cpp)
     target_link_libraries(cuda_timer ginkgo CUDA::cudart)
 endif()

benchmark/run_all_benchmarks.sh

@@ -110,6 +110,8 @@ elif [ "${BENCHMARK_PRECISION}" == "dcomplex" ]; then
     BENCH_SUFFIX="_dcomplex"
 elif [ "${BENCHMARK_PRECISION}" == "scomplex" ]; then
     BENCH_SUFFIX="_scomplex"
+elif [ "${BENCHMARK_PRECISION}" == "half" ]; then
+    BENCH_SUFFIX="_half"
 else
     echo "BENCHMARK_PRECISION is set to the not supported \"${BENCHMARK_PRECISION}\"." 1>&2
     echo "Currently supported values: \"double\", \"single\", \"dcomplex\" and \"scomplex\"" 1>&2
@@ -216,9 +218,16 @@ keep_latest() {
 compute_matrix_statistics() {
     [ "${DRY_RUN}" == "true" ] && return
     cp "$1" "$1.imd" # make sure we're not loosing the original input
-    ./matrix_statistics/matrix_statistics${BENCH_SUFFIX} \
-        --backup="$1.bkp" --double_buffer="$1.bkp2" \
-        <"$1.imd" 2>&1 >"$1"
+    if [ "${BENCH_SUFFIX}" == "_half" ]; then
+        # half precision benchmark still uses single for statistics
+        ./matrix_statistics/matrix_statistics_single \
+            --backup="$1.bkp" --double_buffer="$1.bkp2" \
+            <"$1.imd" 2>&1 >"$1"
+    else
+        ./matrix_statistics/matrix_statistics${BENCH_SUFFIX} \
+            --backup="$1.bkp" --double_buffer="$1.bkp2" \
+            <"$1.imd" 2>&1 >"$1"
+    fi
     keep_latest "$1" "$1.bkp" "$1.bkp2" "$1.imd"
 }
 

benchmark/spmv/CMakeLists.txt

@@ -1,4 +1,9 @@
 ginkgo_add_typed_benchmark_executables(spmv "YES" spmv.cpp)
+# TODO: move to all benchmark
+if (GINKGO_ENABLE_HALF)
+    ginkgo_add_single_benchmark_executable(
+        "spmv_half" "YES" "GKO_BENCHMARK_USE_HALF_PRECISION" "h" spmv.cpp)
+endif()
 if(GINKGO_BUILD_MPI)
     add_subdirectory(distributed)
 endif()

benchmark/spmv/spmv_common.hpp

@@ -144,7 +144,9 @@ struct SpmvBenchmark : Benchmark<spmv_benchmark_state<Generator>> {
             exec->synchronize();
             auto max_relative_norm2 =
                 compute_max_relative_norm2(x_clone.get(), state.answer.get());
-            format_case["max_relative_norm2"] = max_relative_norm2;
+            format_case["max_relative_norm2"] =
+                static_cast<typename gko::detail::arth_type<rc_etype>::type>(
+                    max_relative_norm2);
         }
 
         IterationControl ic{timer};

benchmark/utils/cuda_linops.cpp

@@ -558,14 +558,19 @@ class CusparseHybrid
     ((CUDA_VERSION >= 10020) && !(defined(_WIN32) || defined(__CYGWIN__)))
 
 
+// cuSPARSE does not support 16 bit compute for full 16 bit floating point
+// input. Also, the scalar must be the compute type, i.e. float.
 template <typename ValueType>
-void cusparse_generic_spmv(std::shared_ptr<const gko::CudaExecutor> gpu_exec,
-                           const cusparseSpMatDescr_t mat,
-                           const gko::array<ValueType>& scalars,
-                           const gko::LinOp* b, gko::LinOp* x,
-                           cusparseOperation_t trans, cusparseSpMVAlg_t alg)
+void cusparse_generic_spmv(
+    std::shared_ptr<const gko::CudaExecutor> gpu_exec,
+    const cusparseSpMatDescr_t mat,
+    const gko::array<typename gko::detail::arth_type<ValueType>::type>& scalars,
+    const gko::LinOp* b, gko::LinOp* x, cusparseOperation_t trans,
+    cusparseSpMVAlg_t alg)
 {
     cudaDataType_t cu_value = gko::kernels::cuda::cuda_data_type<ValueType>();
+    cudaDataType_t compute_value = gko::kernels::cuda::cuda_data_type<
+        typename gko::detail::arth_type<ValueType>::type>();
     using gko::kernels::cuda::as_culibs_type;
     auto dense_b = gko::as<gko::matrix::Dense<ValueType>>(b);
     auto dense_x = gko::as<gko::matrix::Dense<ValueType>>(x);
@@ -584,13 +589,14 @@ void cusparse_generic_spmv(std::shared_ptr<const gko::CudaExecutor> gpu_exec,
     gko::size_type buffer_size = 0;
     GKO_ASSERT_NO_CUSPARSE_ERRORS(cusparseSpMV_bufferSize(
         gpu_exec->get_cusparse_handle(), trans, &scalars.get_const_data()[0],
-        mat, vecb, &scalars.get_const_data()[1], vecx, cu_value, alg,
+        mat, vecb, &scalars.get_const_data()[1], vecx, compute_value, alg,
         &buffer_size));
     gko::array<char> buffer_array(gpu_exec, buffer_size);
     auto dbuffer = buffer_array.get_data();
     GKO_ASSERT_NO_CUSPARSE_ERRORS(cusparseSpMV(
         gpu_exec->get_cusparse_handle(), trans, &scalars.get_const_data()[0],
-        mat, vecb, &scalars.get_const_data()[1], vecx, cu_value, alg, dbuffer));
+        mat, vecb, &scalars.get_const_data()[1], vecx, compute_value, alg,
+        dbuffer));
     GKO_ASSERT_NO_CUSPARSE_ERRORS(cusparseDestroyDnVec(vecx));
     GKO_ASSERT_NO_CUSPARSE_ERRORS(cusparseDestroyDnVec(vecb));
 }
@@ -669,8 +675,8 @@ class CusparseGenericCsr
 protected:
     void apply_impl(const gko::LinOp* b, gko::LinOp* x) const override
     {
-        cusparse_generic_spmv(this->get_gpu_exec(), mat_, scalars, b, x, trans_,
-                              Alg);
+        cusparse_generic_spmv<ValueType>(this->get_gpu_exec(), mat_, scalars, b,
+                                         x, trans_, Alg);
     }
 
     void apply_impl(const gko::LinOp* alpha, const gko::LinOp* b,
@@ -686,9 +692,11 @@ class CusparseGenericCsr
     {}
 
 private:
+    using compute_type = typename gko::detail::arth_type<ValueType>::type;
     // Contains {alpha, beta}
-    gko::array<ValueType> scalars{
-        this->get_executor(), {gko::one<ValueType>(), gko::zero<ValueType>()}};
+    gko::array<compute_type> scalars{
+        this->get_executor(),
+        {gko::one<compute_type>(), gko::zero<compute_type>()}};
     std::shared_ptr<csr> csr_;
     cusparseOperation_t trans_;
     cusparseSpMatDescr_t mat_;
@@ -761,8 +769,8 @@ class CusparseGenericCoo
 protected:
     void apply_impl(const gko::LinOp* b, gko::LinOp* x) const override
     {
-        cusparse_generic_spmv(this->get_gpu_exec(), mat_, scalars, b, x, trans_,
-                              default_csr_alg);
+        cusparse_generic_spmv<ValueType>(this->get_gpu_exec(), mat_, scalars, b,
+                                         x, trans_, default_csr_alg);
     }
 
     void apply_impl(const gko::LinOp* alpha, const gko::LinOp* b,
@@ -777,9 +785,11 @@ class CusparseGenericCoo
     {}
 
 private:
+    using compute_type = typename gko::detail::arth_type<ValueType>::type;
     // Contains {alpha, beta}
-    gko::array<ValueType> scalars{
-        this->get_executor(), {gko::one<ValueType>(), gko::zero<ValueType>()}};
+    gko::array<compute_type> scalars{
+        this->get_executor(),
+        {gko::one<compute_type>(), gko::zero<compute_type>()}};
     std::shared_ptr<coo> coo_;
     cusparseOperation_t trans_;
     cusparseSpMatDescr_t mat_;

benchmark/utils/generator.hpp

@@ -158,10 +158,7 @@ struct DefaultSystemGenerator {
     {
         auto res = Vec::create(exec);
         res->read(gko::matrix_data<ValueType, itype>(
-            size,
-            std::uniform_real_distribution<gko::remove_complex<ValueType>>(-1.0,
-                                                                           1.0),
-            get_engine()));
+            size, std::uniform_real_distribution<>(-1.0, 1.0), get_engine()));
         return res;
     }