Skip to content

Commit 9b558d9

Browse files
authored
Optimize the custom SDPA kernel for BF16 (pytorch#21789)
Differential Revision: D116040710 Pull Request resolved: pytorch#21789
1 parent 7a6c1c4 commit 9b558d9

6 files changed

Lines changed: 1115 additions & 23 deletions

File tree

‎extension/llm/custom_ops/op_sdpa_impl.h‎

Lines changed: 118 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,8 @@ void _q_at_k_gemm(
7171
const int64_t q_stride_m,
7272
const MaybeQuantizedMatrixData& k_data,
7373
const int64_t k_stride_n,
74-
accum_t* qk_data) {
74+
accum_t* qk_data,
75+
accum_t* widen_scratch) {
7576
ET_CHECK_MSG(q_data.dtype == k_data.dtype, "q and k must have same dtype");
7677
ET_CHECK_MSG(
7778
q_data.dtype == ScalarType::Char || q_data.dtype == ScalarType::Float ||
@@ -109,6 +110,45 @@ void _q_at_k_gemm(
109110
q_data.dtype == ScalarType::BFloat16 ||
110111
q_data.dtype == ScalarType::Half) {
111112
if constexpr (std::is_same<accum_t, float>::value) {
113+
if (widen_scratch != nullptr && q_data.dtype == ScalarType::BFloat16) {
114+
// Reduction is headSize here, short enough that a per-output dot pays
115+
// a cross-lane reduction it cannot amortize, while a packed BLAS can.
116+
// BLAS has no bf16 entry point, so widen and use the fp32 one; the
117+
// conversion is O(mk + nk) against O(mnk) of multiply.
118+
accum_t* k_f32 = widen_scratch;
119+
accum_t* q_f32 = widen_scratch + k_n * qk_k;
120+
const auto* k_src =
121+
static_cast<const ::executorch::aten::BFloat16*>(k_data.data);
122+
const auto* q_src =
123+
static_cast<const ::executorch::aten::BFloat16*>(q_data.data);
124+
for (int64_t i = 0; i < k_n; ++i) {
125+
for (int64_t l = 0; l < qk_k; ++l) {
126+
k_f32[i * qk_k + l] =
127+
static_cast<accum_t>(k_src[i * k_stride_n + l]);
128+
}
129+
}
130+
for (int64_t j = 0; j < q_m; ++j) {
131+
for (int64_t l = 0; l < qk_k; ++l) {
132+
q_f32[j * qk_k + l] =
133+
static_cast<accum_t>(q_src[j * q_stride_m + l]);
134+
}
135+
}
136+
::executorch::cpublas::gemm(
137+
::executorch::cpublas::TransposeType::Transpose,
138+
::executorch::cpublas::TransposeType::NoTranspose,
139+
k_n,
140+
q_m,
141+
qk_k,
142+
static_cast<accum_t>(1),
143+
k_f32,
144+
qk_k,
145+
q_f32,
146+
qk_k,
147+
static_cast<accum_t>(0),
148+
qk_data,
149+
k_n);
150+
return;
151+
}
112152
auto do_gemm = [&](auto rt_tag) {
113153
using rt = decltype(rt_tag);
114154
::executorch::cpublas::gemm(
@@ -289,7 +329,8 @@ void _qk_at_v_gemm(
289329
accum_t* o_data,
290330
const int64_t o_stride_m,
291331
const accum_t beta,
292-
accum_t* buf_qdq_ptr) {
332+
accum_t* buf_qdq_ptr,
333+
const bool widen_v) {
293334
if (v_data.dtype == ScalarType::Char) {
294335
if constexpr (std::is_same<accum_t, float>::value) {
295336
const float* qk = static_cast<const float*>(qk_data);
@@ -339,6 +380,37 @@ void _qk_at_v_gemm(
339380
// qk has been cast to the activation dtype (see qk_reduced_data); both
340381
// operands are reduced precision and accumulate into the float output.
341382
if constexpr (std::is_same<accum_t, float>::value) {
383+
if (widen_v) {
384+
ET_CHECK_MSG(
385+
v_data.dtype == ScalarType::BFloat16,
386+
"widen_v requires BFloat16 V data");
387+
// Weights are still accum_t, so widen V to match and let BLAS do the
388+
// multiply. buf_qdq_ptr is free here: per-thread, ctx-allocated, the
389+
// right size, and only the quantized branch above uses it.
390+
const auto* v_src =
391+
static_cast<const ::executorch::aten::BFloat16*>(v_data.data);
392+
for (int64_t kk = 0; kk < k; ++kk) {
393+
for (int64_t nn = 0; nn < n; ++nn) {
394+
buf_qdq_ptr[kk * n + nn] =
395+
static_cast<accum_t>(v_src[kk * v_stride_n + nn]);
396+
}
397+
}
398+
::executorch::cpublas::gemm(
399+
::executorch::cpublas::TransposeType::NoTranspose,
400+
::executorch::cpublas::TransposeType::NoTranspose,
401+
n,
402+
m,
403+
k,
404+
static_cast<accum_t>(1),
405+
buf_qdq_ptr,
406+
n,
407+
static_cast<const accum_t*>(qk_data),
408+
qk_stride_m,
409+
beta,
410+
o_data,
411+
o_stride_m);
412+
return;
413+
}
342414
auto do_gemm = [&](auto rt_tag) {
343415
using rt = decltype(rt_tag);
344416
::executorch::cpublas::gemm(
@@ -848,6 +920,34 @@ void cpu_flash_attention(
848920
buf_reduced = scratch_reduced.get();
849921
}
850922
}
923+
// Widening only pays for a short reduction feeding many output columns; both
924+
// bounds are empirical.
925+
constexpr int64_t kMaxHeadSizeForWidenedQK = 128;
926+
constexpr int64_t kMinQBlockForWidenedQK = 32;
927+
// Scratch for widening q@K.T to fp32 (see _q_at_k_gemm): one K block plus one
928+
// q block. qBlockSize cannot exceed qSplitSize, so include the runtime bounds
929+
// that determine whether any block can use the widened path.
930+
const bool widen_qk =
931+
std::is_same<scalar_t, ::executorch::aten::BFloat16>::value &&
932+
::executorch::cpublas::gemm_uses_blas() &&
933+
headSize <= kMaxHeadSizeForWidenedQK &&
934+
qSplitSize >= kMinQBlockForWidenedQK;
935+
int64_t size_per_thread_widen =
936+
widen_qk ? (kvSplitSize + qSplitSize) * headSize : 0;
937+
std::unique_ptr<char[]> allocated_buf_widen;
938+
accum_t* widen_buf = nullptr;
939+
if (widen_qk) {
940+
int64_t size_widen_bytes =
941+
size_per_thread_widen * num_thread * sizeof(accum_t);
942+
Result<void*> scratch_widen = ctx.allocate_temp(size_widen_bytes, 64);
943+
if (!scratch_widen.ok()) {
944+
allocated_buf_widen = std::make_unique<char[]>(size_widen_bytes);
945+
widen_buf = reinterpret_cast<accum_t*>(allocated_buf_widen.get());
946+
} else {
947+
widen_buf = reinterpret_cast<accum_t*>(scratch_widen.get());
948+
}
949+
}
950+
851951
int64_t size_per_thread_qdq_vec = kvSplitSize * headSize;
852952
// Lets align size_per_thread_qdq_vec to 64 bytes, for coalesced cache reads,
853953
// by padding with right number of per thread elements
@@ -891,6 +991,8 @@ void cpu_flash_attention(
891991
: nullptr;
892992
accum_t* buf_qdq_ptr =
893993
scratch_for_quant_dequant + ompIdx * size_per_thread_qdq_vec;
994+
accum_t* widen_ptr =
995+
widen_buf ? widen_buf + ompIdx * size_per_thread_widen : nullptr;
894996

895997
for (int64_t z = begin; z < end; z++) {
896998
int64_t m = k * qSplitSize;
@@ -987,7 +1089,9 @@ void cpu_flash_attention(
9871089
qStrideM,
9881090
k_sub_matrix_data,
9891091
kStrideN,
990-
qk_data);
1092+
qk_data,
1093+
(widen_qk && qBlockSize >= kMinQBlockForWidenedQK) ? widen_ptr
1094+
: nullptr);
9911095

9921096
// There are 4 cases that is_causal has to cover to fill
9931097
// not-attendable-position with -inf
@@ -1139,9 +1243,17 @@ void cpu_flash_attention(
11391243
// For reduced-precision activations the attention weights are cast to
11401244
// the activation dtype so that Softmax(q @ k.T) @ v runs as a
11411245
// reduced-in/float-out gemm matching the value matrix.
1246+
// Casting the weights down only pays when the bf16 attn@V kernel is
1247+
// the one that runs. Above kMinQBlockForWidenedAV it is faster to keep
1248+
// them in accum_t, widen V and let BLAS multiply -- also one rounding
1249+
// step fewer. Below it the widening stops amortizing.
1250+
constexpr int64_t kMinQBlockForWidenedAV = 64;
1251+
const bool widen_v = is_reduced_type && !is_quantized_sdpa &&
1252+
std::is_same<scalar_t, ::executorch::aten::BFloat16>::value &&
1253+
qBlockSize >= kMinQBlockForWidenedAV;
11421254
const void* qk_gemm_data = qk_data;
11431255
if constexpr (is_reduced_type) {
1144-
if (!is_quantized_sdpa) {
1256+
if (!is_quantized_sdpa && !widen_v) {
11451257
vec::convert<accum_t, scalar_t>(
11461258
qk_data, qk_reduced_data, qBlockSize * kvBlockSize);
11471259
qk_gemm_data = qk_reduced_data;
@@ -1159,7 +1271,8 @@ void cpu_flash_attention(
11591271
dst_data,
11601272
headSize,
11611273
n == 0 ? static_cast<accum_t>(0) : static_cast<accum_t>(1),
1162-
buf_qdq_ptr);
1274+
buf_qdq_ptr,
1275+
widen_v);
11631276
}
11641277
// dst <- dst / sum[row]
11651278
// reorder MHA output with strides

0 commit comments

Comments
 (0)