Skip to content

opencl: keep half-precision matrix-vector products on the GEMM path #3727

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
melonakos wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

opencl: keep half-precision matrix-vector products on the GEMM path #3727

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8243caa
Keep half-precision matrix-vector products on the GEMM path in OpenCL
melonakos Sep 11, 2026
fa5a75c
State the observed GEMV failure rather than an unproven cause in the …
melonakos Sep 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
State the observed GEMV failure rather than an unproven cause in the … 
…comment
  • Loading branch information
@melonakos
melonakos committed Sep 11, 2026
commit fa5a75c328bd89149879f2b25ff14d4ee26239f0
436 changes: 218 additions & 218 deletions README.md
View file
Open in desktop

Large diffs are not rendered by default.

367 changes: 184 additions & 183 deletions src/backend/opencl/blas.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,183 +1,184 @@
/*******************************************************
* Copyright (c) 2014, ArrayFire
* All rights reserved.
*
* This file is distributed under 3-clause BSD license.
* The complete license agreement can be obtained at:
* http://arrayfire.com/licenses/BSD-3-Clause
********************************************************/

#include <blas.hpp>

#include <Array.hpp>
#include <arith.hpp>
#include <common/half.hpp>
#include <common/traits.hpp>
#include <complex.hpp>
#include <err_opencl.hpp>
#include <math.hpp>
#include <reduce.hpp>
#include <transpose.hpp>

#include <complex>
#include <vector>

// Includes one of the supported OpenCL BLAS back-ends (e.g. clBLAS, CLBlast)
#include <cpu/cpu_blas.hpp>
#include <magma/magma_blas.h>
#include <type_traits>

using arrayfire::common::half;

namespace arrayfire {
namespace opencl {

void initBlas() { gpu_blas_init(); }

void deInitBlas() { gpu_blas_deinit(); }

// Converts an af_mat_prop options to a transpose type for one of the OpenCL
// BLAS back-ends
OPENCL_BLAS_TRANS_T
toBlasTranspose(af_mat_prop opt) {
switch (opt) {
case AF_MAT_NONE: return OPENCL_BLAS_NO_TRANS;
case AF_MAT_TRANS: return OPENCL_BLAS_TRANS;
case AF_MAT_CTRANS: return OPENCL_BLAS_CONJ_TRANS;
default: AF_ERROR("INVALID af_mat_prop", AF_ERR_ARG);
}
}

template<typename T>
void gemm_fallback(Array<T> &out, af_mat_prop optLhs, af_mat_prop optRhs,
const T *alpha, const Array<T> &lhs, const Array<T> &rhs,
const T *beta) {
cpu::gemm(out, optLhs, optRhs, alpha, lhs, rhs, beta);
}

template<>
void gemm_fallback<half>(Array<half> & /*out*/, af_mat_prop /*optLhs*/,
af_mat_prop /*optRhs*/, const half * /*alpha*/,
const Array<half> & /*lhs*/,
const Array<half> & /*rhs*/, const half * /*beta*/) {
assert(false && "CPU fallback not implemented for f16");
}

template<typename Ti, typename To>
void gemm(Array<To> &out, af_mat_prop optLhs, af_mat_prop optRhs,
const To *alpha, const Array<Ti> &lhs, const Array<Ti> &rhs,
const To *beta) {
#if defined(WITH_LINEAR_ALGEBRA)
// Do not force offload gemm on OSX Intel devices
if (OpenCLCPUOffload(false) &&
static_cast<af_dtype>(dtype_traits<Ti>::af_type) != f16) {
gemm_fallback(out, optLhs, optRhs, alpha, lhs, rhs, beta);
return;
}
#endif
const auto lOpts = toBlasTranspose(optLhs);
const auto rOpts = toBlasTranspose(optRhs);

const auto aRowDim = (lOpts == OPENCL_BLAS_NO_TRANS) ? 0 : 1;
const auto aColDim = (lOpts == OPENCL_BLAS_NO_TRANS) ? 1 : 0;
const auto bColDim = (rOpts == OPENCL_BLAS_NO_TRANS) ? 1 : 0;

const dim4 &lDims = lhs.dims();
const dim4 &rDims = rhs.dims();
const int M = lDims[aRowDim];
const int N = rDims[bColDim];
const int K = lDims[aColDim];
const dim4 oDims = out.dims();

const dim4 &lStrides = lhs.strides();
const dim4 &rStrides = rhs.strides();
const dim4 oStrides = out.strides();

int batchSize = static_cast<int>(oDims[2] * oDims[3]);

bool is_l_d2_batched = oDims[2] == lDims[2];
bool is_l_d3_batched = oDims[3] == lDims[3];
bool is_r_d2_batched = oDims[2] == rDims[2];
bool is_r_d3_batched = oDims[3] == rDims[3];

for (int n = 0; n < batchSize; n++) {
int w = static_cast<int>(n / oDims[2]);
int z = static_cast<int>(n - w * oDims[2]);

int loff = z * (is_l_d2_batched * lStrides[2]) +
w * (is_l_d3_batched * lStrides[3]);
int roff = z * (is_r_d2_batched * rStrides[2]) +
w * (is_r_d3_batched * rStrides[3]);

dim_t lOffset = lhs.getOffset() + loff;
dim_t rOffset = rhs.getOffset() + roff;
dim_t oOffset = out.getOffset() + z * oStrides[2] + w * oStrides[3];

cl::Event event;
// CLBlast's half-precision GEMV returns wrong values (CNugteren/CLBlast
// issue 561); its GEMM does not, so keep f16 on the GEMM path even
// when the right-hand side is a single column.
const bool useGemv =
rDims[bColDim] == 1 && !std::is_same<Ti, half>::value;
if (useGemv) {
dim_t incr = (optRhs == AF_MAT_NONE) ? rStrides[0] : rStrides[1];
gpu_blas_gemv_func<Ti> gemv;
OPENCL_BLAS_CHECK(gemv(lOpts, lDims[0], lDims[1], *alpha,
(*lhs.get())(), lOffset, lStrides[1],
(*rhs.get())(), rOffset, incr, *beta,
(*out.get())(), oOffset, oStrides[0], 1,
&getQueue()(), 0, nullptr, &event()));
} else {
gpu_blas_gemm_func<Ti> gemm;
OPENCL_BLAS_CHECK(gemm(lOpts, rOpts, M, N, K, *alpha,
(*lhs.get())(), lOffset, lStrides[1],
(*rhs.get())(), rOffset, rStrides[1], *beta,
(*out.get())(), oOffset, oStrides[1], 1,
&getQueue()(), 0, nullptr, &event()));
}
}
}

template<>
void gemm<schar, float>(Array<float> &out, af_mat_prop optLhs,
af_mat_prop optRhs, const float *alpha,
const Array<schar> &lhs, const Array<schar> &rhs,
const float *beta) {
TYPE_ERROR(3, af_dtype::s8);
}

template<typename T>
Array<T> dot(const Array<T> &lhs, const Array<T> &rhs, af_mat_prop optLhs,
af_mat_prop optRhs) {
const Array<T> lhs_ = (optLhs == AF_MAT_NONE ? lhs : conj<T>(lhs));
const Array<T> rhs_ = (optRhs == AF_MAT_NONE ? rhs : conj<T>(rhs));

const Array<T> temp = arithOp<T, af_mul_t>(lhs_, rhs_, lhs_.dims());
return reduce<af_add_t, T, T>(temp, 0, false, 0);
}

#define INSTANTIATE_GEMM(TYPE) \
template void gemm<TYPE>(Array<TYPE> & out, af_mat_prop optLhs, \
af_mat_prop optRhs, const TYPE *alpha, \
const Array<TYPE> &lhs, const Array<TYPE> &rhs, \
const TYPE *beta);

INSTANTIATE_GEMM(float)
INSTANTIATE_GEMM(cfloat)
INSTANTIATE_GEMM(double)
INSTANTIATE_GEMM(cdouble)
INSTANTIATE_GEMM(half)

#define INSTANTIATE_DOT(TYPE) \
template Array<TYPE> dot<TYPE>(const Array<TYPE> &lhs, \
const Array<TYPE> &rhs, af_mat_prop optLhs, \
af_mat_prop optRhs);

INSTANTIATE_DOT(float)
INSTANTIATE_DOT(double)
INSTANTIATE_DOT(cfloat)
INSTANTIATE_DOT(cdouble)
INSTANTIATE_DOT(half)

} // namespace opencl
} // namespace arrayfire
/*******************************************************
* Copyright (c) 2014, ArrayFire
* All rights reserved.
*
* This file is distributed under 3-clause BSD license.
* The complete license agreement can be obtained at:
* http://arrayfire.com/licenses/BSD-3-Clause
********************************************************/

#include <blas.hpp>

#include <Array.hpp>
#include <arith.hpp>
#include <common/half.hpp>
#include <common/traits.hpp>
#include <complex.hpp>
#include <err_opencl.hpp>
#include <math.hpp>
#include <reduce.hpp>
#include <transpose.hpp>

#include <complex>
#include <vector>

// Includes one of the supported OpenCL BLAS back-ends (e.g. clBLAS, CLBlast)
#include <cpu/cpu_blas.hpp>
#include <magma/magma_blas.h>
#include <type_traits>

using arrayfire::common::half;

namespace arrayfire {
namespace opencl {

void initBlas() { gpu_blas_init(); }

void deInitBlas() { gpu_blas_deinit(); }

// Converts an af_mat_prop options to a transpose type for one of the OpenCL
// BLAS back-ends
OPENCL_BLAS_TRANS_T
toBlasTranspose(af_mat_prop opt) {
switch (opt) {
case AF_MAT_NONE: return OPENCL_BLAS_NO_TRANS;
case AF_MAT_TRANS: return OPENCL_BLAS_TRANS;
case AF_MAT_CTRANS: return OPENCL_BLAS_CONJ_TRANS;
default: AF_ERROR("INVALID af_mat_prop", AF_ERR_ARG);
}
}

template<typename T>
void gemm_fallback(Array<T> &out, af_mat_prop optLhs, af_mat_prop optRhs,
const T *alpha, const Array<T> &lhs, const Array<T> &rhs,
const T *beta) {
cpu::gemm(out, optLhs, optRhs, alpha, lhs, rhs, beta);
}

template<>
void gemm_fallback<half>(Array<half> & /*out*/, af_mat_prop /*optLhs*/,
af_mat_prop /*optRhs*/, const half * /*alpha*/,
const Array<half> & /*lhs*/,
const Array<half> & /*rhs*/, const half * /*beta*/) {
assert(false && "CPU fallback not implemented for f16");
}

template<typename Ti, typename To>
void gemm(Array<To> &out, af_mat_prop optLhs, af_mat_prop optRhs,
const To *alpha, const Array<Ti> &lhs, const Array<Ti> &rhs,
const To *beta) {
#if defined(WITH_LINEAR_ALGEBRA)
// Do not force offload gemm on OSX Intel devices
if (OpenCLCPUOffload(false) &&
static_cast<af_dtype>(dtype_traits<Ti>::af_type) != f16) {
gemm_fallback(out, optLhs, optRhs, alpha, lhs, rhs, beta);
return;
}
#endif
const auto lOpts = toBlasTranspose(optLhs);
const auto rOpts = toBlasTranspose(optRhs);

const auto aRowDim = (lOpts == OPENCL_BLAS_NO_TRANS) ? 0 : 1;
const auto aColDim = (lOpts == OPENCL_BLAS_NO_TRANS) ? 1 : 0;
const auto bColDim = (rOpts == OPENCL_BLAS_NO_TRANS) ? 1 : 0;

const dim4 &lDims = lhs.dims();
const dim4 &rDims = rhs.dims();
const int M = lDims[aRowDim];
const int N = rDims[bColDim];
const int K = lDims[aColDim];
const dim4 oDims = out.dims();

const dim4 &lStrides = lhs.strides();
const dim4 &rStrides = rhs.strides();
const dim4 oStrides = out.strides();

int batchSize = static_cast<int>(oDims[2] * oDims[3]);

bool is_l_d2_batched = oDims[2] == lDims[2];
bool is_l_d3_batched = oDims[3] == lDims[3];
bool is_r_d2_batched = oDims[2] == rDims[2];
bool is_r_d3_batched = oDims[3] == rDims[3];

for (int n = 0; n < batchSize; n++) {
int w = static_cast<int>(n / oDims[2]);
int z = static_cast<int>(n - w * oDims[2]);

int loff = z * (is_l_d2_batched * lStrides[2]) +
w * (is_l_d3_batched * lStrides[3]);
int roff = z * (is_r_d2_batched * rStrides[2]) +
w * (is_r_d3_batched * rStrides[3]);

dim_t lOffset = lhs.getOffset() + loff;
dim_t rOffset = rhs.getOffset() + roff;
dim_t oOffset = out.getOffset() + z * oStrides[2] + w * oStrides[3];

cl::Event event;
// With fp16 data CLBlast's GEMV path returns wrong values on Intel GPUs
// while its GEMM path is correct (#3674; possibly related to
// CNugteren/CLBlast#561), so keep half on GEMM even for a single
// column.
const bool useGemv =
rDims[bColDim] == 1 && !std::is_same<Ti, half>::value;
if (useGemv) {
dim_t incr = (optRhs == AF_MAT_NONE) ? rStrides[0] : rStrides[1];
gpu_blas_gemv_func<Ti> gemv;
OPENCL_BLAS_CHECK(gemv(lOpts, lDims[0], lDims[1], *alpha,
(*lhs.get())(), lOffset, lStrides[1],
(*rhs.get())(), rOffset, incr, *beta,
(*out.get())(), oOffset, oStrides[0], 1,
&getQueue()(), 0, nullptr, &event()));
} else {
gpu_blas_gemm_func<Ti> gemm;
OPENCL_BLAS_CHECK(gemm(lOpts, rOpts, M, N, K, *alpha,
(*lhs.get())(), lOffset, lStrides[1],
(*rhs.get())(), rOffset, rStrides[1], *beta,
(*out.get())(), oOffset, oStrides[1], 1,
&getQueue()(), 0, nullptr, &event()));
}
}
}

template<>
void gemm<schar, float>(Array<float> &out, af_mat_prop optLhs,
af_mat_prop optRhs, const float *alpha,
const Array<schar> &lhs, const Array<schar> &rhs,
const float *beta) {
TYPE_ERROR(3, af_dtype::s8);
}

template<typename T>
Array<T> dot(const Array<T> &lhs, const Array<T> &rhs, af_mat_prop optLhs,
af_mat_prop optRhs) {
const Array<T> lhs_ = (optLhs == AF_MAT_NONE ? lhs : conj<T>(lhs));
const Array<T> rhs_ = (optRhs == AF_MAT_NONE ? rhs : conj<T>(rhs));

const Array<T> temp = arithOp<T, af_mul_t>(lhs_, rhs_, lhs_.dims());
return reduce<af_add_t, T, T>(temp, 0, false, 0);
}

#define INSTANTIATE_GEMM(TYPE) \
template void gemm<TYPE>(Array<TYPE> & out, af_mat_prop optLhs, \
af_mat_prop optRhs, const TYPE *alpha, \
const Array<TYPE> &lhs, const Array<TYPE> &rhs, \
const TYPE *beta);

INSTANTIATE_GEMM(float)
INSTANTIATE_GEMM(cfloat)
INSTANTIATE_GEMM(double)
INSTANTIATE_GEMM(cdouble)
INSTANTIATE_GEMM(half)

#define INSTANTIATE_DOT(TYPE) \
template Array<TYPE> dot<TYPE>(const Array<TYPE> &lhs, \
const Array<TYPE> &rhs, af_mat_prop optLhs, \
af_mat_prop optRhs);

INSTANTIATE_DOT(float)
INSTANTIATE_DOT(double)
INSTANTIATE_DOT(cfloat)
INSTANTIATE_DOT(cdouble)
INSTANTIATE_DOT(half)

} // namespace opencl
} // namespace arrayfire
Loading
Loading