Tracker: Float16 support

This issue is meant to be an overview / track progress on supporting Float16 in CUDA.jl

Library APIs

We'll be able to do a lot already if we just dispatch to the correct library calls, e.g., cublasHgemm for half-precision, or cublasGemmEx with appropriate CUDA datatypes set. Note that this requires a matching math mode, https://docs.nvidia.com/cuda/cublas/index.html#cublassetmathmode, so https://github.com/JuliaGPU/CUDA.jl/issues/354 may be relevant.

Since we support many devices / toolkits, we'll need to check what the requirements for these calls are, and whether the APIs offer a fall-back (e.g., when tensor-cores are disabled or unavailable). If not, we may need to offer pure Julia fallbacks, convert to Float32 arrays, or require the caller (i.e. Knet, Flux) to check for device capabilities before using half-precision data. That may be the easiest solution.

Alternatively, we could look into the TF32 format which IIUC should make it possible to more easily use Float16 functions when available, or fall back to existing Float32 functionality (since the representation doesn't change). That way we also do not need kernel programming capabilities (see below) to implement non-library operations.

Kernel programming

Beyond API calls, we need to be able to compile and execute half-precision code. This requires several features:

Half-precision intrinsics

libdevice doesn't seem to provide half-precision intrinsics, but instead cuda_fp16.hpp contains inline assembly snippets:

__CUDA_FP16_DECL__ __half hsin(const __half a) {
    __half r = __hsin_internal(a);
    asm("{\n\t"
        "  .reg.b16 i,r,t;     \n\t"
        "  mov.b16 r, %0;      \n\t"
        "  mov.b16 i, %1;      \n\t"
        "  mov.b16 t, 0x8000;  \n\t"
        "  and.b16 t,r,t;      \n\t"
        __SPEC_CASE(i, r, 0X32B3, 0x0800)
        __SPEC_CASE(i, r, 0X5CB0, 0x1000)
        __SPEC_CASE(i, r, 0XB2B3, 0x8800)
        __SPEC_CASE(i, r, 0XDCB0, 0x9000)
        "  or.b16  r,r,t;      \n\t"
        "  mov.b16 %0, r;      \n"
        "}\n" : "+h"(__HALF_TO_US(r)) : "h"(__HALF_TO_CUS(a)));
    return r;
}

We can do inline assembly from Julia: https://github.com/JuliaGPU/CUDA.jl/blob/d9db2beeb6bb61f5d6b4c586b8de2e82dc43d78c/src/device/intrinsics/warp_vote.jl#L11-L21

half LLVM representation

Julia currently represents half-precision floating-point as i16 and implements essential operations itself. Ideally we'd switch to using half and let LLVM emit the appropriate GPU instructions, like https://github.com/JuliaLang/julia/pull/26381 does. Alternatively, we can use our own implementation, like https://github.com/vchuravy/IEEEFloat16.jl. @vchuravy, any status update here?

Note that we already can use LLVM's half with llvmcall after https://github.com/JuliaLang/julia/pull/33970; which makes it possible to, e.g., llvmcall the WMMA hardware intrinsics. We could use this to interface with other intrinsics, or even implement more essential operations, but it'd be better to update the Julia compiler and have it use half instead of reimplementing everything.

cc @denizyuret