diff --git a/blas/CMakeLists.txt b/blas/CMakeLists.txt new file mode 100644 index 000000000..bf3554d1e --- /dev/null +++ b/blas/CMakeLists.txt @@ -0,0 +1,22 @@ +#------------------------------------------------------------------------------ +# Copyright (c) 2026 Ainekko, Co. +# SPDX-License-Identifier: Apache-2.0 +#------------------------------------------------------------------------------ + +function(add_etsoc_blas_kernel TARGET_NAME) + set(options) + set(oneValueArgs INSTALL_DESTINATION) + set(multiValueArgs) + cmake_parse_arguments(ARGS "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) + + if(NOT DEFINED ARGS_INSTALL_DESTINATION) + set(ARGS_INSTALL_DESTINATION kernels/blas) + endif() + + add_etsoc_riscv_executable(${TARGET_NAME} ${ARGS_UNPARSED_ARGUMENTS}) + target_link_libraries(${TARGET_NAME} etsoc_crt0) + install(TARGETS ${TARGET_NAME} ${TARGET_NAME}_dbg DESTINATION ${ARGS_INSTALL_DESTINATION}) +endfunction() + +add_subdirectory(reference) +add_subdirectory(optimized) diff --git a/blas/README.md b/blas/README.md new file mode 100644 index 000000000..c8bfa8f5b --- /dev/null +++ b/blas/README.md @@ -0,0 +1,34 @@ +# BLAS Kernels + +This directory holds custom GP-SDK device kernels for the ET-SOC1 BLAS work. + +It is intended to be consumed through the existing `gp-sdk/device` custom-kernel +hook by setting: + +- `CUSTOM_KERNELS_SRC_DIR=/blas` +- `CUSTOM_KERNELS_BIN_DIR=/blas` + +Layout conventions: + +- `blas/reference/` holds numerically conservative scalar baselines +- `blas/optimized/` holds vector and tensor implementations +- datatype directories sit below those roots +- BLAS level and operation directories sit below the datatype + +Datatype rollout order: + +- `fp32` first +- `fp16` +- `bf16` +- `int16` +- `int8` + +Implementation policy: + +- every optimized kernel must have a matching reference kernel +- reference and optimized kernels must keep the same external argument contract +- numerical comparisons are always made against the reference implementation + +The intent is to always keep a trustworthy reference kernel beside each +optimized family so we can measure error growth as operation ordering and +reduction strategies change. diff --git a/blas/optimized/CMakeLists.txt b/blas/optimized/CMakeLists.txt new file mode 100644 index 000000000..f18746e78 --- /dev/null +++ b/blas/optimized/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(fp32) diff --git a/blas/optimized/README.md b/blas/optimized/README.md new file mode 100644 index 000000000..3385f9e28 --- /dev/null +++ b/blas/optimized/README.md @@ -0,0 +1,59 @@ +# Optimized Kernels + +Optimized kernels hold vector and tensor implementations that are validated +against the matching reference kernels. + +Allowed here: + +- vector and tensor implementations +- tiling, blocking, fusion, and staging for performance +- controlled operation reordering when justified by performance + +Required here: + +- keep the same external kernel contract as the matching reference kernel +- document any intentional numerical tradeoff in code comments or commit history +- validate against the matching reference implementation + +Not allowed here: + +- silently changing kernel semantics +- changing datatype interpretation relative to the reference path +- landing an optimized kernel before a reference kernel exists + +Current kernels: + +- `fp32/level1/axpy`: scalar inner loop unrolled by 4 with the same minion + partitioning and external contract as the reference kernel +- `fp32/level1/dot`: ET packed-SIMD loads plus packed multiply/add in the main + loop, with a masked packed-SIMD tail in the default artifact and a scalar-tail + comparison artifact kept alongside it +- `fp32/level2/gemv`: ET packed-SIMD vectorization across contiguous output-row + blocks for the non-transpose, unit-stride path, with scalar fallback for + transpose and non-unit-stride cases +- `fp32/level3/gemm`: + - `blas_gemm_optimized_fp32_vector`: single-minion 4x4 blocked `NN` + micro-kernel with K strip-mining for lower-latency execution + - `blas_gemm_optimized_fp32_tensor`: tensor-engine tiled `NN` path for + higher-throughput execution, with scalar fallback for unsupported tails and + non-`NN` cases + - current tensor bring-up keeps the tensor engine on the bulk of each tile + and applies a post-store scalar correction for the dropped final inner term + observed on the current hardware path + +Current validation policy: + +- prefer `*_dbg` kernel artifacts when validating optimized kernels +- treat release-ELF bring-up as a separate runtime/toolchain issue +- do not treat a release-ELF failure by itself as evidence that the optimized + kernel math is wrong when the matching `_dbg` artifact passes + +Current ET-SIMD bring-up notes: + +- `fp32/level1/dot` now uses ET packed-SIMD loads and packed multiply/add in + the main loop +- the default optimized `dot` artifact uses a masked packed-SIMD tail to avoid + switching back into scalar FP arithmetic while packed state is still live in + the overlaid `f` register file +- `blas_dot_optimized_fp32_scalar_tail.elf{,_dbg}` is retained as a comparison + artifact for validating the scalar-tail alternative diff --git a/blas/optimized/fp32/CMakeLists.txt b/blas/optimized/fp32/CMakeLists.txt new file mode 100644 index 000000000..8fd9c533b --- /dev/null +++ b/blas/optimized/fp32/CMakeLists.txt @@ -0,0 +1,3 @@ +add_subdirectory(level1) +add_subdirectory(level2) +add_subdirectory(level3) diff --git a/blas/optimized/fp32/level1/CMakeLists.txt b/blas/optimized/fp32/level1/CMakeLists.txt new file mode 100644 index 000000000..4a11fbf8d --- /dev/null +++ b/blas/optimized/fp32/level1/CMakeLists.txt @@ -0,0 +1,2 @@ +add_subdirectory(axpy) +add_subdirectory(dot) diff --git a/blas/optimized/fp32/level1/axpy/CMakeLists.txt b/blas/optimized/fp32/level1/axpy/CMakeLists.txt new file mode 100644 index 000000000..952d4fe64 --- /dev/null +++ b/blas/optimized/fp32/level1/axpy/CMakeLists.txt @@ -0,0 +1,6 @@ +add_etsoc_blas_kernel( + blas_axpy_optimized_fp32.elf + axpy.cpp + INSTALL_DESTINATION kernels/blas/optimized/fp32/level1/axpy +) +target_include_directories(blas_axpy_optimized_fp32.elf PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/blas/optimized/fp32/level1/axpy/axpy.cpp b/blas/optimized/fp32/level1/axpy/axpy.cpp new file mode 100644 index 000000000..7ae22bb80 --- /dev/null +++ b/blas/optimized/fp32/level1/axpy/axpy.cpp @@ -0,0 +1,51 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include +#include + +#include + +#include "axpy_kernel_arguments.h" +#include "entryPoint.h" + +int entryPoint_0(KernelArguments* args); +DECLARE_KERNEL_ENTRY_POINTS(entryPoint_0, nullptr); + +int entryPoint_0(KernelArguments* args) { + const auto minionId = get_relative_thread_id(); + const size_t numWorkers = SOC_MINIONS_PER_SHIRE; + + if (args->numElements == 0) { + return 0; + } + + size_t elemsPerWorker = (args->numElements + numWorkers - 1) / numWorkers; + if (elemsPerWorker % 16) { + elemsPerWorker += 16 - (elemsPerWorker % 16); + } + + const size_t begin = elemsPerWorker * minionId; + const size_t end = std::min(elemsPerWorker * (minionId + 1), static_cast(args->numElements)); + if (begin > (args->numElements - 1)) { + return 0; + } + + size_t i = begin; + const size_t loopEnd = begin + ((end - begin) / 4) * 4; + for (; i < loopEnd; i += 4) { + args->y[i + 0] = args->alpha * args->x[i + 0] + args->y[i + 0]; + args->y[i + 1] = args->alpha * args->x[i + 1] + args->y[i + 1]; + args->y[i + 2] = args->alpha * args->x[i + 2] + args->y[i + 2]; + args->y[i + 3] = args->alpha * args->x[i + 3] + args->y[i + 3]; + } + + for (; i < end; ++i) { + args->y[i] = args->alpha * args->x[i] + args->y[i]; + } + + return 0; +} diff --git a/blas/optimized/fp32/level1/axpy/axpy_kernel_arguments.h b/blas/optimized/fp32/level1/axpy/axpy_kernel_arguments.h new file mode 100644 index 000000000..42078b93b --- /dev/null +++ b/blas/optimized/fp32/level1/axpy/axpy_kernel_arguments.h @@ -0,0 +1,13 @@ +#ifndef ET_BLAS_OPTIMIZED_FP32_AXPY_KERNEL_ARGUMENTS_H +#define ET_BLAS_OPTIMIZED_FP32_AXPY_KERNEL_ARGUMENTS_H + +#include + +struct KernelArguments { + uint64_t numElements; + const float* x; + float* y; + float alpha; +} __attribute__((packed)); + +#endif diff --git a/blas/optimized/fp32/level1/dot/CMakeLists.txt b/blas/optimized/fp32/level1/dot/CMakeLists.txt new file mode 100644 index 000000000..addc23795 --- /dev/null +++ b/blas/optimized/fp32/level1/dot/CMakeLists.txt @@ -0,0 +1,12 @@ +function(add_dot_variant target_name tail_define) + add_etsoc_blas_kernel( + ${target_name} + dot.cpp + INSTALL_DESTINATION kernels/blas/optimized/fp32/level1/dot + ) + target_include_directories(${target_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + target_compile_definitions(${target_name} PRIVATE ${tail_define}=1) +endfunction() + +add_dot_variant(blas_dot_optimized_fp32.elf ET_BLAS_DOT_TAIL_MASKED) +add_dot_variant(blas_dot_optimized_fp32_scalar_tail.elf ET_BLAS_DOT_TAIL_SCALAR) diff --git a/blas/optimized/fp32/level1/dot/dot.cpp b/blas/optimized/fp32/level1/dot/dot.cpp new file mode 100644 index 000000000..ee9aa4252 --- /dev/null +++ b/blas/optimized/fp32/level1/dot/dot.cpp @@ -0,0 +1,184 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include +#include +#include + +#include +#include +#include + +#include "CommonCode.h" +#include "dot_kernel_arguments.h" +#include "entryPoint.h" +#include "sync.h" + +#if !defined(ET_BLAS_DOT_TAIL_SCALAR) && !defined(ET_BLAS_DOT_TAIL_MASKED) +#define ET_BLAS_DOT_TAIL_MASKED 1 +#endif + +namespace { + +constexpr size_t kVectorLanes = 8; + +inline __attribute__((always_inline)) void setMaskForLaneCount(size_t activeLanes) { + switch (activeLanes) { + case 0: + mask_set(0, 0x00); + break; + case 1: + mask_set(0, 0x01); + break; + case 2: + mask_set(0, 0x03); + break; + case 3: + mask_set(0, 0x07); + break; + case 4: + mask_set(0, 0x0f); + break; + case 5: + mask_set(0, 0x1f); + break; + case 6: + mask_set(0, 0x3f); + break; + case 7: + mask_set(0, 0x7f); + break; + default: + mask_set(0, 0xff); + break; + } +} + +inline __attribute__((always_inline)) float packedHorizontalSum(float packedValue) { + alignas(32) float lanes[kVectorLanes]; + __asm__ __volatile__("fsw.ps %[packedValue], 0(%[dst])\n" + : + : [dst] "r"(lanes), [packedValue] "f"(packedValue) + : "memory"); + + float result = 0.0f; + for (size_t lane = 0; lane < kVectorLanes; ++lane) { + result += lanes[lane]; + } + return result; +} + +inline __attribute__((always_inline)) float sdotVectorCore( + const float* x, const float* y, size_t begin, size_t end) { + float packedAcc; + uint32_t zeroWord = 0; + __asm__ __volatile__("fbcx.ps %[packedAcc], %[zeroWord]\n" + : [packedAcc] "=&f"(packedAcc) + : [zeroWord] "r"(zeroWord)); + + size_t i = begin; + const size_t loopEnd = begin + ((end - begin) / kVectorLanes) * kVectorLanes; + for (; i < loopEnd; i += kVectorLanes) { + float xVec; + float yVec; + float prodVec; + const float* xPtr = x + i; + const float* yPtr = y + i; + + __asm__ __volatile__("flw.ps %[xVec], 0(%[xPtr])\n" + "flw.ps %[yVec], 0(%[yPtr])\n" + : [xVec] "=&f"(xVec), [yVec] "=&f"(yVec) + : [xPtr] "r"(xPtr), [yPtr] "r"(yPtr) + : "memory"); + + __asm__ __volatile__("fmul.ps %[prodVec], %[xVec], %[yVec]\n" + "fadd.ps %[packedAcc], %[packedAcc], %[prodVec]\n" + : [packedAcc] "+&f"(packedAcc), [prodVec] "=&f"(prodVec) + : [xVec] "f"(xVec), [yVec] "f"(yVec)); + } + +#if defined(ET_BLAS_DOT_TAIL_MASKED) + const size_t remaining = end - i; + if (remaining > 0) { + alignas(32) float tailX[kVectorLanes] = {}; + alignas(32) float tailY[kVectorLanes] = {}; + for (size_t lane = 0; lane < remaining; ++lane) { + tailX[lane] = x[i + lane]; + tailY[lane] = y[i + lane]; + } + + setMaskForLaneCount(remaining); + + float xVec; + float yVec; + float prodVec; + __asm__ __volatile__("flw.ps %[xVec], 0(%[xPtr])\n" + "flw.ps %[yVec], 0(%[yPtr])\n" + : [xVec] "=&f"(xVec), [yVec] "=&f"(yVec) + : [xPtr] "r"(tailX), [yPtr] "r"(tailY) + : "memory"); + + __asm__ __volatile__("fmul.ps %[prodVec], %[xVec], %[yVec]\n" + "fadd.ps %[packedAcc], %[packedAcc], %[prodVec]\n" + : [packedAcc] "+&f"(packedAcc), [prodVec] "=&f"(prodVec) + : [xVec] "f"(xVec), [yVec] "f"(yVec)); + + mask_set(0, 0xff); + } +#endif + + float localSum = packedHorizontalSum(packedAcc); + +#if defined(ET_BLAS_DOT_TAIL_SCALAR) + for (; i < end; ++i) { + localSum += x[i] * y[i]; + } +#endif + + return localSum; +} + +} // namespace + +int entryPoint_0(KernelArguments* args); +DECLARE_KERNEL_ENTRY_POINTS(entryPoint_0, nullptr); + +int entryPoint_0(KernelArguments* args) { + const auto minionId = get_relative_thread_id(); + const size_t numWorkers = SOC_MINIONS_PER_SHIRE; + + if (args->numElements == 0) { + if (minionId == 0) { + *(args->res) = 0.0f; + } + return 0; + } + + size_t elemsPerWorker = (args->numElements + numWorkers - 1) / numWorkers; + if (elemsPerWorker % 16) { + elemsPerWorker += 16 - (elemsPerWorker % 16); + } + + const size_t begin = elemsPerWorker * minionId; + const size_t end = std::min(elemsPerWorker * (minionId + 1), static_cast(args->numElements)); + + if (begin <= (args->numElements - 1)) { + const float localSum = sdotVectorCore(args->x, args->y, begin, end); + args->partials[begin] = localSum; + evictCacheLine(0x1ULL, reinterpret_cast(&args->partials[begin])); + } + hart::barrier(); + + if (minionId == 0) { + float result = 0.0f; + for (size_t i = 0; i < args->numElements; i += elemsPerWorker) { + result += args->partials[i]; + } + *(args->res) = result; + } + + return 0; +} diff --git a/blas/optimized/fp32/level1/dot/dot_kernel_arguments.h b/blas/optimized/fp32/level1/dot/dot_kernel_arguments.h new file mode 100644 index 000000000..f9770fb2b --- /dev/null +++ b/blas/optimized/fp32/level1/dot/dot_kernel_arguments.h @@ -0,0 +1,15 @@ +#ifndef ET_BLAS_OPTIMIZED_FP32_DOT_KERNEL_ARGUMENTS_H +#define ET_BLAS_OPTIMIZED_FP32_DOT_KERNEL_ARGUMENTS_H + +#include + +// Keep this ABI aligned with the reference dot kernel. +struct KernelArguments { + uint64_t numElements; + const float* x; + const float* y; + float* partials; + float* res; +} __attribute__((packed)); + +#endif diff --git a/blas/optimized/fp32/level2/CMakeLists.txt b/blas/optimized/fp32/level2/CMakeLists.txt new file mode 100644 index 000000000..5f1a80a85 --- /dev/null +++ b/blas/optimized/fp32/level2/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(gemv) diff --git a/blas/optimized/fp32/level2/gemv/CMakeLists.txt b/blas/optimized/fp32/level2/gemv/CMakeLists.txt new file mode 100644 index 000000000..d17cdf75d --- /dev/null +++ b/blas/optimized/fp32/level2/gemv/CMakeLists.txt @@ -0,0 +1,6 @@ +add_etsoc_blas_kernel( + blas_gemv_optimized_fp32.elf + gemv.cpp + INSTALL_DESTINATION kernels/blas/optimized/fp32/level2/gemv +) +target_include_directories(blas_gemv_optimized_fp32.elf PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/blas/optimized/fp32/level2/gemv/gemv.cpp b/blas/optimized/fp32/level2/gemv/gemv.cpp new file mode 100644 index 000000000..82f22a4c6 --- /dev/null +++ b/blas/optimized/fp32/level2/gemv/gemv.cpp @@ -0,0 +1,204 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include +#include +#include + +#include +#include + +#include "entryPoint.h" +#include "gemv_kernel_arguments.h" + +namespace { + +constexpr size_t kVectorLanes = 8; + +bool isTranspose(char trans) { + return trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'; +} + +int startIndex(int length, int inc) { + return inc >= 0 ? 0 : (length - 1) * (-inc); +} + +inline __attribute__((always_inline)) void setMaskForLaneCount(size_t activeLanes) { + switch (activeLanes) { + case 0: + mask_set(0, 0x00); + break; + case 1: + mask_set(0, 0x01); + break; + case 2: + mask_set(0, 0x03); + break; + case 3: + mask_set(0, 0x07); + break; + case 4: + mask_set(0, 0x0f); + break; + case 5: + mask_set(0, 0x1f); + break; + case 6: + mask_set(0, 0x3f); + break; + case 7: + mask_set(0, 0x7f); + break; + default: + mask_set(0, 0xff); + break; + } +} + +inline __attribute__((always_inline)) float broadcastPackedFloat(float value) { + uint32_t bits = 0; + __builtin_memcpy(&bits, &value, sizeof(bits)); + + float packedValue; + __asm__ __volatile__("fbcx.ps %[packedValue], %[bits]\n" + : [packedValue] "=&f"(packedValue) + : [bits] "r"(bits)); + return packedValue; +} + +} // namespace + +int entryPoint_0(KernelArguments* args); +DECLARE_KERNEL_ENTRY_POINTS(entryPoint_0, nullptr); + +int entryPoint_0(KernelArguments* args) { + if (args->m <= 0 || args->n <= 0) { + return 0; + } + + const bool transpose = isTranspose(args->trans); + const int outputCount = transpose ? args->n : args->m; + const int reductionCount = transpose ? args->m : args->n; + const int xStart = startIndex(reductionCount, args->incx); + const int yStart = startIndex(outputCount, args->incy); + + const auto minionId = get_relative_thread_id(); + const size_t numWorkers = SOC_MINIONS_PER_SHIRE; + size_t elemsPerWorker = (static_cast(outputCount) + numWorkers - 1) / numWorkers; + if (elemsPerWorker % 16) { + elemsPerWorker += 16 - (elemsPerWorker % 16); + } + + const size_t begin = elemsPerWorker * minionId; + const size_t end = std::min(elemsPerWorker * (minionId + 1), static_cast(outputCount)); + if (begin > static_cast(outputCount - 1)) { + return 0; + } + + const bool usePackedVectorPath = !transpose && args->incx == 1 && args->incy == 1; + + if (usePackedVectorPath) { + const float packedAlpha = broadcastPackedFloat(args->alpha); + const float packedBeta = broadcastPackedFloat(args->beta); + + for (size_t rowBase = begin; rowBase < end; rowBase += kVectorLanes) { + const size_t rowsInBlock = std::min(kVectorLanes, end - rowBase); + setMaskForLaneCount(rowsInBlock); + + float accVec = broadcastPackedFloat(0.0f); + for (int column = 0; column < args->n; ++column) { + const float* aPtr = args->a + rowBase + static_cast(column) * args->lda; + const float* xPtr = args->x + xStart + column; + + float aVec; + float xVec; + float prodVec; + + uint32_t xBits = 0; + __builtin_memcpy(&xBits, xPtr, sizeof(xBits)); + + __asm__ __volatile__("flw.ps %[aVec], 0(%[aPtr])\n" + "fbcx.ps %[xVec], %[xBits]\n" + "fmul.ps %[prodVec], %[aVec], %[xVec]\n" + "fadd.ps %[accVec], %[accVec], %[prodVec]\n" + : [accVec] "+&f"(accVec), [aVec] "=&f"(aVec), [xVec] "=&f"(xVec), + [prodVec] "=&f"(prodVec) + : [aPtr] "r"(aPtr), [xBits] "r"(xBits) + : "memory"); + } + + float yVec; + float scaledAccVec; + float scaledYVec; + float outVec; + float* yPtr = args->y + yStart + static_cast(rowBase); + + __asm__ __volatile__("flw.ps %[yVec], 0(%[yPtr])\n" + "fmul.ps %[scaledAccVec], %[accVec], %[packedAlpha]\n" + "fmul.ps %[scaledYVec], %[yVec], %[packedBeta]\n" + "fadd.ps %[outVec], %[scaledAccVec], %[scaledYVec]\n" + "fsw.ps %[outVec], 0(%[yPtr])\n" + : [yVec] "=&f"(yVec), [scaledAccVec] "=&f"(scaledAccVec), + [scaledYVec] "=&f"(scaledYVec), [outVec] "=&f"(outVec) + : [accVec] "f"(accVec), [packedAlpha] "f"(packedAlpha), + [packedBeta] "f"(packedBeta), [yPtr] "r"(yPtr) + : "memory"); + } + + mask_set(0, 0xff); + return 0; + } + + for (size_t outputIndex = begin; outputIndex < end; ++outputIndex) { + float sum0 = 0.0f; + float sum1 = 0.0f; + float sum2 = 0.0f; + float sum3 = 0.0f; + int inner = 0; + + if (!transpose) { + const float* matrixRow = args->a + outputIndex; + const int loopEnd = (args->n / 4) * 4; + for (; inner < loopEnd; inner += 4) { + sum0 += matrixRow[static_cast(inner + 0) * args->lda] * + args->x[xStart + (inner + 0) * args->incx]; + sum1 += matrixRow[static_cast(inner + 1) * args->lda] * + args->x[xStart + (inner + 1) * args->incx]; + sum2 += matrixRow[static_cast(inner + 2) * args->lda] * + args->x[xStart + (inner + 2) * args->incx]; + sum3 += matrixRow[static_cast(inner + 3) * args->lda] * + args->x[xStart + (inner + 3) * args->incx]; + } + float sum = sum0 + sum1 + sum2 + sum3; + for (; inner < args->n; ++inner) { + sum += matrixRow[static_cast(inner) * args->lda] * args->x[xStart + inner * args->incx]; + } + float& yValue = args->y[yStart + static_cast(outputIndex) * args->incy]; + yValue = args->alpha * sum + args->beta * yValue; + } else { + const int loopEnd = (args->m / 4) * 4; + for (; inner < loopEnd; inner += 4) { + sum0 += args->a[static_cast(inner + 0) + outputIndex * args->lda] * + args->x[xStart + (inner + 0) * args->incx]; + sum1 += args->a[static_cast(inner + 1) + outputIndex * args->lda] * + args->x[xStart + (inner + 1) * args->incx]; + sum2 += args->a[static_cast(inner + 2) + outputIndex * args->lda] * + args->x[xStart + (inner + 2) * args->incx]; + sum3 += args->a[static_cast(inner + 3) + outputIndex * args->lda] * + args->x[xStart + (inner + 3) * args->incx]; + } + float sum = sum0 + sum1 + sum2 + sum3; + for (; inner < args->m; ++inner) { + sum += args->a[static_cast(inner) + outputIndex * args->lda] * + args->x[xStart + inner * args->incx]; + } + float& yValue = args->y[yStart + static_cast(outputIndex) * args->incy]; + yValue = args->alpha * sum + args->beta * yValue; + } + } + + return 0; +} diff --git a/blas/optimized/fp32/level2/gemv/gemv_kernel_arguments.h b/blas/optimized/fp32/level2/gemv/gemv_kernel_arguments.h new file mode 100644 index 000000000..9f3844c65 --- /dev/null +++ b/blas/optimized/fp32/level2/gemv/gemv_kernel_arguments.h @@ -0,0 +1,20 @@ +#ifndef ET_BLAS_OPTIMIZED_FP32_GEMV_KERNEL_ARGUMENTS_H +#define ET_BLAS_OPTIMIZED_FP32_GEMV_KERNEL_ARGUMENTS_H + +#include + +struct KernelArguments { + char trans; + int32_t m; + int32_t n; + float alpha; + const float* a; + int32_t lda; + const float* x; + int32_t incx; + float beta; + float* y; + int32_t incy; +}; + +#endif diff --git a/blas/optimized/fp32/level3/CMakeLists.txt b/blas/optimized/fp32/level3/CMakeLists.txt new file mode 100644 index 000000000..0cf2c16da --- /dev/null +++ b/blas/optimized/fp32/level3/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(gemm) diff --git a/blas/optimized/fp32/level3/gemm/CMakeLists.txt b/blas/optimized/fp32/level3/gemm/CMakeLists.txt new file mode 100644 index 000000000..1c166dcb3 --- /dev/null +++ b/blas/optimized/fp32/level3/gemm/CMakeLists.txt @@ -0,0 +1,11 @@ +function(add_gemm_variant target_name source_file) + add_etsoc_blas_kernel( + ${target_name} + ${source_file} + INSTALL_DESTINATION kernels/blas/optimized/fp32/level3/gemm + ) + target_include_directories(${target_name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) +endfunction() + +add_gemm_variant(blas_gemm_optimized_fp32_vector.elf gemm.cpp) +add_gemm_variant(blas_gemm_optimized_fp32_tensor.elf gemm_tensor.cpp) diff --git a/blas/optimized/fp32/level3/gemm/gemm.cpp b/blas/optimized/fp32/level3/gemm/gemm.cpp new file mode 100644 index 000000000..e3228deb7 --- /dev/null +++ b/blas/optimized/fp32/level3/gemm/gemm.cpp @@ -0,0 +1,99 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include +#include +#include + +#include + +#include "entryPoint.h" +#include "gemm_kernel_arguments.h" + +namespace { + +bool isTranspose(char trans) { + return trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'; +} + +float generalElement(const float* matrix, int ld, int row, int col, bool transpose) { + return transpose ? matrix[col + static_cast(row) * ld] : matrix[row + static_cast(col) * ld]; +} + +} // namespace + +int entryPoint_0(KernelArguments* args); +DECLARE_KERNEL_ENTRY_POINTS(entryPoint_0, nullptr); + +int entryPoint_0(KernelArguments* args) { + if (args->m <= 0 || args->n <= 0 || args->k < 0) { + return 0; + } + + if (get_relative_thread_id() != 0) { + return 0; + } + + const bool transA = isTranspose(args->transa); + const bool transB = isTranspose(args->transb); + + if (!(transA || transB)) { + // Baseline optimized path: + // - 4x4 register-resident micro-kernel over the M/N dimensions + // - strip-mine K so the active A/B working set stays bounded + constexpr int kBlock = 32; + for (int colBlock = 0; colBlock < args->n; colBlock += 4) { + const int colsInBlock = std::min(4, args->n - colBlock); + for (int rowBlock = 0; rowBlock < args->m; rowBlock += 4) { + const int rowsInBlock = std::min(4, args->m - rowBlock); + float acc[4][4] = {}; + + for (int kBase = 0; kBase < args->k; kBase += kBlock) { + const int kEnd = std::min(args->k, kBase + kBlock); + for (int inner = kBase; inner < kEnd; ++inner) { + float aVals[4] = {0.0f, 0.0f, 0.0f, 0.0f}; + float bVals[4] = {0.0f, 0.0f, 0.0f, 0.0f}; + + for (int r = 0; r < rowsInBlock; ++r) { + aVals[r] = args->a[(rowBlock + r) + static_cast(inner) * args->lda]; + } + for (int c = 0; c < colsInBlock; ++c) { + bVals[c] = args->b[inner + static_cast(colBlock + c) * args->ldb]; + } + + for (int r = 0; r < rowsInBlock; ++r) { + for (int c = 0; c < colsInBlock; ++c) { + acc[r][c] += aVals[r] * bVals[c]; + } + } + } + } + + for (int c = 0; c < colsInBlock; ++c) { + for (int r = 0; r < rowsInBlock; ++r) { + float& cValue = args->c[(rowBlock + r) + static_cast(colBlock + c) * args->ldc]; + cValue = args->alpha * acc[r][c] + args->beta * cValue; + } + } + } + } + return 0; + } + + for (int col = 0; col < args->n; ++col) { + for (int row = 0; row < args->m; ++row) { + float sum = 0.0f; + for (int inner = 0; inner < args->k; ++inner) { + sum += generalElement(args->a, args->lda, row, inner, transA) * + generalElement(args->b, args->ldb, inner, col, transB); + } + float& cValue = args->c[row + static_cast(col) * args->ldc]; + cValue = args->alpha * sum + args->beta * cValue; + } + } + + return 0; +} diff --git a/blas/optimized/fp32/level3/gemm/gemm_kernel_arguments.h b/blas/optimized/fp32/level3/gemm/gemm_kernel_arguments.h new file mode 100644 index 000000000..e6f4baefd --- /dev/null +++ b/blas/optimized/fp32/level3/gemm/gemm_kernel_arguments.h @@ -0,0 +1,22 @@ +#ifndef ET_BLAS_OPTIMIZED_FP32_GEMM_KERNEL_ARGUMENTS_H +#define ET_BLAS_OPTIMIZED_FP32_GEMM_KERNEL_ARGUMENTS_H + +#include + +struct KernelArguments { + char transa; + char transb; + int32_t m; + int32_t n; + int32_t k; + float alpha; + const float* a; + int32_t lda; + const float* b; + int32_t ldb; + float beta; + float* c; + int32_t ldc; +}; + +#endif diff --git a/blas/optimized/fp32/level3/gemm/gemm_tensor.cpp b/blas/optimized/fp32/level3/gemm/gemm_tensor.cpp new file mode 100644 index 000000000..905bd11b5 --- /dev/null +++ b/blas/optimized/fp32/level3/gemm/gemm_tensor.cpp @@ -0,0 +1,206 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include +#include +#include + +#include +#include +#include + +#include "entryPoint.h" +#include "gemm_kernel_arguments.h" + +namespace { + +constexpr int kTensorTileRows = 16; +constexpr int kTensorTileCols = 16; +constexpr int kTensorTileK = 16; +constexpr int kTensorColGranularity = 4; +constexpr uint64_t kTensorStrideBytes = 64; +constexpr uint64_t kTensorOpcodeFp32 = 0; + +alignas(64) float gATile[kTensorTileRows * kTensorTileCols]; +alignas(64) float gBTile[kTensorTileRows * kTensorTileCols]; +alignas(64) float gCTile[kTensorTileRows * kTensorTileCols]; + +bool isTranspose(char trans) { + return trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'; +} + +float generalElement(const float* matrix, int ld, int row, int col, bool transpose) { + return transpose ? matrix[col + static_cast(row) * ld] : matrix[row + static_cast(col) * ld]; +} + +void drainCoalescingBuffers() { + constexpr uint32_t l2CacheBanks = 4; + + for (uint32_t i = 0; i < l2CacheBanks; ++i) { + volatile uint64_t* control = reinterpret_cast(ESR_CACHE(SHIRE_OWN, i, SC_IDX_COP_SM_CTL_USER)); + uint64_t state = 0; + do { + state = (*control >> 24) & 0xff; + } while (state != 4); + } + + for (uint32_t i = 0; i < l2CacheBanks; ++i) { + volatile uint64_t* control = reinterpret_cast(ESR_CACHE(SHIRE_OWN, i, SC_IDX_COP_SM_CTL_USER)); + *control = (1ULL << 0) | (10ULL << 8); + } + + for (uint32_t i = 0; i < l2CacheBanks; ++i) { + volatile uint64_t* control = reinterpret_cast(ESR_CACHE(SHIRE_OWN, i, SC_IDX_COP_SM_CTL_USER)); + uint64_t state = 0; + do { + state = (*control >> 24) & 0xff; + } while (state != 4); + } +} + +void runGenericGemm(KernelArguments* args, bool transA, bool transB) { + for (int col = 0; col < args->n; ++col) { + for (int row = 0; row < args->m; ++row) { + float sum = 0.0f; + for (int inner = 0; inner < args->k; ++inner) { + sum += generalElement(args->a, args->lda, row, inner, transA) * + generalElement(args->b, args->ldb, inner, col, transB); + } + float& cValue = args->c[row + static_cast(col) * args->ldc]; + cValue = args->alpha * sum + args->beta * cValue; + } + } +} + +void runScalarNNTile(KernelArguments* args, int rowBegin, int rowEnd, int colBegin, int colEnd) { + for (int col = colBegin; col < colEnd; ++col) { + for (int row = rowBegin; row < rowEnd; ++row) { + float sum = 0.0f; + for (int inner = 0; inner < args->k; ++inner) { + sum += args->a[row + static_cast(inner) * args->lda] * + args->b[inner + static_cast(col) * args->ldb]; + } + float& cValue = args->c[row + static_cast(col) * args->ldc]; + cValue = args->alpha * sum + args->beta * cValue; + } + } +} + +void clearTile(float* tile) { + for (int i = 0; i < kTensorTileRows * kTensorTileCols; ++i) { + tile[i] = 0.0f; + } +} + +void packATile(float* tile, const KernelArguments* args, int rowBlock, int rowsInTile, int kBlock, int kInTile) { + clearTile(tile); + for (int row = 0; row < rowsInTile; ++row) { + for (int inner = 0; inner < kInTile; ++inner) { + tile[row * kTensorTileCols + inner] = + args->a[(rowBlock + row) + static_cast(kBlock + inner) * args->lda]; + } + } +} + +void packBTile(float* tile, const KernelArguments* args, int kBlock, int kInTile, int colBlock, int colsInTile) { + clearTile(tile); + for (int inner = 0; inner < kInTile; ++inner) { + for (int col = 0; col < colsInTile; ++col) { + tile[inner * kTensorTileCols + col] = + args->b[(kBlock + inner) + static_cast(colBlock + col) * args->ldb]; + } + } +} + +void runTensorNNTile(KernelArguments* args, int rowBlock, int rowsInTile, int colBlock, int colsInTile) { + clearTile(gCTile); + + bool clearRf = true; + for (int kBlock = 0; kBlock < args->k; kBlock += kTensorTileK) { + const int kInTile = std::min(kTensorTileK, args->k - kBlock); + packATile(gATile, args, rowBlock, rowsInTile, kBlock, kInTile); + packBTile(gBTile, args, kBlock, kInTile, colBlock, colsInTile); + + tensor_load(false, false, 0, 0, 0, reinterpret_cast(gATile), 0, static_cast(rowsInTile - 1), + kTensorStrideBytes, 0); + tensor_load(false, false, 32, 0, 1, reinterpret_cast(gBTile), 0, static_cast(kInTile - 1), + kTensorStrideBytes, 1); + tensor_wait(TENSOR_LOAD_WAIT_0); + + tensor_fma(false, static_cast(colsInTile / kTensorColGranularity - 1), + static_cast(rowsInTile - 1), static_cast(kInTile - 1), 0, + false, false, false, true, 32, 0, kTensorOpcodeFp32, clearRf ? 1 : 0); + tensor_wait(TENSOR_FMA_WAIT); + clearRf = false; + } + + tensor_store(0, 0, static_cast(colsInTile - 1), static_cast(rowsInTile - 1), + reinterpret_cast(gCTile), 0, kTensorStrideBytes); + tensor_wait(TENSOR_STORE_WAIT); + drainCoalescingBuffers(); + + for (int col = 0; col < colsInTile; ++col) { + for (int row = 0; row < rowsInTile; ++row) { + float tileValue = gCTile[row * kTensorTileCols + col]; + + // The current tensor-FMA hardware path drops the last inner term of each + // tensor pass. Apply the scalar correction only after tensor_store so the + // overlaid F-register file is no longer live. + for (int kBlock = 0; kBlock < args->k; kBlock += kTensorTileK) { + const int kInTile = std::min(kTensorTileK, args->k - kBlock); + const int correctionInner = kBlock + kInTile - 1; + tileValue += args->a[(rowBlock + row) + static_cast(correctionInner) * args->lda] * + args->b[correctionInner + static_cast(colBlock + col) * args->ldb]; + } + + float& cValue = args->c[(rowBlock + row) + static_cast(colBlock + col) * args->ldc]; + cValue = args->alpha * tileValue + args->beta * cValue; + } + } +} + +} // namespace + +int entryPoint_0(KernelArguments* args); +DECLARE_KERNEL_ENTRY_POINTS(entryPoint_0, nullptr); + +int entryPoint_0(KernelArguments* args) { + if (args->m <= 0 || args->n <= 0 || args->k < 0) { + return 0; + } + + if (get_relative_thread_id() != 0) { + return 0; + } + + const bool transA = isTranspose(args->transa); + const bool transB = isTranspose(args->transb); + + if (transA || transB || args->k == 0) { + runGenericGemm(args, transA, transB); + return 0; + } + + for (int rowBlock = 0; rowBlock < args->m; rowBlock += kTensorTileRows) { + const int rowsInTile = std::min(kTensorTileRows, args->m - rowBlock); + int colBlock = 0; + for (; colBlock < args->n;) { + int colsInTensor = std::min(kTensorTileCols, args->n - colBlock); + colsInTensor -= colsInTensor % kTensorColGranularity; + + if (colsInTensor >= kTensorColGranularity) { + runTensorNNTile(args, rowBlock, rowsInTile, colBlock, colsInTensor); + colBlock += colsInTensor; + continue; + } + + runScalarNNTile(args, rowBlock, rowBlock + rowsInTile, colBlock, args->n); + break; + } + } + + return 0; +} diff --git a/blas/optimized/host/CMakeLists.txt b/blas/optimized/host/CMakeLists.txt new file mode 100644 index 000000000..ba2c19b82 --- /dev/null +++ b/blas/optimized/host/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 3.20) + +project(blas-optimized-host LANGUAGES CXX) + +include(GNUInstallDirs) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +set(ET_PLATFORM_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/../../.." CACHE PATH "Path to the et-platform source tree") +set(BUILD_TESTS OFF CACHE BOOL "Build GP-SDK example launchers" FORCE) + +list(APPEND CMAKE_MODULE_PATH /opt/et/lib/cmake /opt/et/lib/cmake/cmake-modules) +list(APPEND CMAKE_PREFIX_PATH /opt/et/lib/cmake /opt/et/lib/cmake/cmake-modules) + +add_subdirectory(${ET_PLATFORM_SOURCE_DIR}/gp-sdk/host gp-sdk-host) + +add_executable(blas_optimized_sdot_verifier sdot_verifier.cpp) +target_link_libraries(blas_optimized_sdot_verifier PRIVATE etsoc_gpsdk) + +add_executable(blas_optimized_proof_of_life_verifier proof_of_life_verifier.cpp) +target_link_libraries(blas_optimized_proof_of_life_verifier PRIVATE etsoc_gpsdk) + +install(TARGETS + blas_optimized_proof_of_life_verifier + blas_optimized_sdot_verifier + DESTINATION ${CMAKE_INSTALL_BINDIR} +) diff --git a/blas/optimized/host/proof_of_life_verifier.cpp b/blas/optimized/host/proof_of_life_verifier.cpp new file mode 100644 index 000000000..d26133a5d --- /dev/null +++ b/blas/optimized/host/proof_of_life_verifier.cpp @@ -0,0 +1,624 @@ +/*------------------------------------------------------------------------- + * Copyright (c) 2026 Ainekko, Co. + * SPDX-License-Identifier: Apache-2.0 + *------------------------------------------------------------------------- + */ + +#include "GenericLauncher.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +struct AxpyArgs { + uint64_t numElements; + const float* x; + float* y; + float alpha; +} __attribute__((packed)); + +struct GemvArgs { + char trans; + int32_t m; + int32_t n; + float alpha; + const float* a; + int32_t lda; + const float* x; + int32_t incx; + float beta; + float* y; + int32_t incy; +}; + +struct GemmArgs { + char transa; + char transb; + int32_t m; + int32_t n; + int32_t k; + float alpha; + const float* a; + int32_t lda; + const float* b; + int32_t ldb; + float beta; + float* c; + int32_t ldc; +}; + +struct Options { + fs::path kernel_root = ""; + fs::path host_results_path = "blas/optimized/proof_of_life_host.md"; + fs::path device_results_path = "blas/optimized/proof_of_life_device.md"; + int kernel_launch_timeout = 30; + std::string device_type = "silicon"; + std::string selected_case = "all"; + int saxpy_num_elements = 259; + int gemv_m = 17; + int gemv_n = 19; + int gemm_m = 17; + int gemm_n = 15; + int gemm_k = 13; + double epsilon = 1.0e-5; +}; + +struct VectorRecord { + std::string label; + std::vector values; +}; + +struct MatrixRecord { + std::string label; + int rows = 0; + int cols = 0; + std::vector values; +}; + +struct CaseResult { + std::string name; + std::string params; + fs::path kernel_artifact; + bool ok = false; + std::vector vectors; + std::vector matrices; +}; + +struct VerificationPair { + CaseResult host; + CaseResult device; +}; + +Options parseArgs(int argc, char* const* argv, std::vector& nextlevel) { + static constexpr const char* helpMsg = + "Usage: [options]\n\n" + "Proof-of-life verifier for optimized SAXPY/SGEMV/SGEMM kernels.\n\n" + "Optional switches:\n" + " -k, --kernel_root root containing optimized fp32 kernels\n" + " --host_results_path markdown file for host reference outputs\n" + " --device_results_path markdown file for device outputs\n" + " --case case to run (all, saxpy, sgemv_n, sgemm_vector_nn, sgemm_tensor_nn)\n" + " -t, --kernel_launch_timeout timeout (in seconds) to wait for kernel completion\n" + " -d, --device_type device type to use (sysemu, fake, silicon)\n" + " --saxpy_num_elements vector length for SAXPY\n" + " --gemv_m matrix rows for GEMV\n" + " --gemv_n matrix columns for GEMV\n" + " --gemm_m output rows for GEMM\n" + " --gemm_n output columns for GEMM\n" + " --gemm_k reduction dimension for GEMM\n" + " -e, --epsilon comparison tolerance for float results\n"; + + static constexpr const char* shortOpts = "k:t:d:e:h"; + static const std::vector