# Matrix multiplication using SME2 intrinsics in C

## In this learning path

- [Introduction](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/)
- [Overview](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/overview/)
- [Set up your SME2 development environment](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/1-get-started/)
- [Test your SME2 development environment](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/2-check-your-environment/)
- [Streaming mode and ZA state in SME](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/3-streaming-mode/)
- [Vanilla matrix multiplication](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/4-vanilla-matmul/)
- [Outer product](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/5-outer-product/)
- [SME2 assembly matrix multiplication](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/6-sme2-matmul-asm/)
- [Matrix multiplication using SME2 intrinsics in C](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/7-sme2-matmul-intr/)
- [Benchmarking](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/8-benchmarking/)
- [Debugging](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/9-debugging/)
- [Going further](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/10-going-further/)
- [Next Steps](https://learn.arm.com/learning-paths/cross-platform/multiplying-matrices-with-sme2/_next-steps/)

In this section, you will write an SME2-optimized matrix multiplication routine in C using the intrinsics that the compiler provides.

## What are intrinsics?

*Intrinsics*, also known as *compiler intrinsics* or *intrinsic functions*, are the functions available to application developers that the compiler has intimate knowledge of. This enables the compiler to either translate the function to a specific instruction or to perform specific optimizations, or both.

You can learn more about intrinsics in this [Wikipedia Article on Intrinsic Function](https://en.wikipedia.org/wiki/Intrinsic_function).

Using intrinsics allows you to write performance-critical code in C while still using standard constructs like loops. This produces performance close to what can be reached with hand-written assembly whilst being significantly more maintainable and portable.

All Arm-specific intrinsics are specified in the [ACLE](https://github.com/ARM-software/acle), which is the Arm C Language Extension. ACLE is supported by the main compilers, most notably [GCC](https://gcc.gnu.org/) and [Clang](https://clang.llvm.org).

## Implementation

In this example, a top-level function named `matmul_intr`, defined in `matmul_intr.c`, brings together the preprocessing and matrix multiplication steps:

```c
__arm_new("za") __arm_locally_streaming void matmul_intr(
    uint64_t M, uint64_t K, uint64_t N,
    const float *restrict matLeft,
    const float *restrict matRight,
    float *restrict matLeft_mod,
    float *restrict matResult) {
    uint64_t SVL = svcntsw();
    preprocess_l_intr(M, K, SVL, matLeft, matLeft_mod);
    matmul_intr_impl(M, K, N, SVL, matLeft_mod, matRight, matResult);
}
```

Note the use of `__arm_new("za")` and `__arm_locally_streaming` at line 1. These attributes ensure that the compiler saves the ZA storage, allowing your function to use it safely without destroying its content if it was still in use by one of the callers.

`SVL`, the dimension of the ZA storage, is requested from the underlying hardware with the `svcntsw()` function call at line 5, and passed down to the `preprocess_l_intr` and `matmul_intr_impl` functions. `svcntsw()` is a function provided by the ACLE library.

### Matrix preprocessing

```c
void preprocess_l_intr(
    uint64_t M, uint64_t K, uint64_t SVL,
    const float *restrict a,
    float *restrict a_mod) __arm_streaming __arm_inout("za") {
    const uint64_t M_mod = SVL * (M / SVL + (M % SVL != 0 ? 1 : 0));

    // The outer loop, iterating over rows (M dimension)
    for (uint64_t row = 0; row < M; row += SVL) {
        svbool_t pMDim = svwhilelt_b32(row, M);

        // The inner loop, iterating on columns (K dimension).
        for (uint64_t col = 0; col < K; col += 2 * SVL) {
            svcount_t pKDim = svwhilelt_c32(col, K, 2);

            // Load-as-rows
            for (uint64_t trow = 0; trow < SVL; trow += 4) {
                svcount_t p0 = svpsel_lane_c32(pKDim, pMDim, trow + 0);
                svcount_t p1 = svpsel_lane_c32(pKDim, pMDim, trow + 1);
                svcount_t p2 = svpsel_lane_c32(pKDim, pMDim, trow + 2);
                svcount_t p3 = svpsel_lane_c32(pKDim, pMDim, trow + 3);

                const uint64_t tile_UL_corner = (row + trow) * K + col;
                svfloat32x2_t zp0 = svld1_x2(p0, &a[tile_UL_corner + 0 * K]);
                svfloat32x2_t zp1 = svld1_x2(p1, &a[tile_UL_corner + 1 * K]);
                svfloat32x2_t zp2 = svld1_x2(p2, &a[tile_UL_corner + 2 * K]);
                svfloat32x2_t zp3 = svld1_x2(p3, &a[tile_UL_corner + 3 * K]);

                svfloat32x4_t zq0 = svcreate4(svget2(zp0, 0), svget2(zp1, 0), svget2(zp2, 0), svget2(zp3, 0));
                svfloat32x4_t zq1 = svcreate4(svget2(zp0, 1), svget2(zp1, 1), svget2(zp2, 1), svget2(zp3, 1));
                svwrite_hor_za32_f32_vg4(/* tile: */ 0, /* slice: */ trow, zq0);
                svwrite_hor_za32_f32_vg4(/* tile: */ 1, /* slice: */ trow, zq1);
            }

            // Read-as-columns and store
            const uint64_t dest_0 = row * K + col * SVL;
            const uint64_t dest_1 = dest_0 + SVL * SVL;
            for (uint64_t tcol = 0; tcol < SVL; tcol += 4) {
                svcount_t p0 = svwhilelt_c32(dest_0 + tcol * SVL, K * M_mod, 4);
                svcount_t p1 = svwhilelt_c32(dest_1 + tcol * SVL, K * M_mod, 4);
                svfloat32x4_t zq0 = svread_ver_za32_f32_vg4(/* tile: */ 0, /* slice: */ tcol);
                svfloat32x4_t zq1 = svread_ver_za32_f32_vg4(/* tile: */ 1, /* slice: */ tcol);
                svst1(p0, &a_mod[dest_0 + tcol * SVL], zq0);
                svst1(p1, &a_mod[dest_1 + tcol * SVL], zq1);
            }
        }
    }
}
```

Note that `preprocess_l_intr` has been annotated at line 3 with:

- `__arm_streaming` - because this function is using streaming instructions
- `__arm_inout("za")` - because `preprocess_l_intr` reuses the ZA storage from its caller

The matrix preprocessing is performed in a double-nested loop, over the `M` (line 7) and `K` (line 12) dimensions of the input matrix `a`. Both loops have an `SVL` step increment, which corresponds to the horizontal and vertical dimensions of the ZA storage that will be used.

Using intrinsics simplifies function development, provided you have a good understanding of the SME2 instruction set. Predicates, which are fundamental to both SVE and SME, allow you to express algorithms cleanly while handling corner cases efficiently.

### Outer-product multiplication

```c
void matmul_intr_impl(
    uint64_t M, uint64_t K, uint64_t N, uint64_t SVL,
    const float *restrict matLeft_mod, const float *restrict matRight,
    float *restrict matResult) __arm_streaming __arm_inout("za") {
    
    // Build the result matrix tile by tile.
    for (uint64_t row = 0; row < M; row += SVL) {
        svbool_t pMDim = svwhilelt_b32(row, M);

        for (uint64_t col = 0; col < N; col += SVL) {
            svbool_t pNDim = svwhilelt_b32(col, N);

            // Outer product + accumulation
            svzero_za();
            const uint64_t matLeft_pos = row * K;
            const uint64_t matRight_UL_corner = col;
            for (uint64_t k = 0; k < K; k++) {
                svfloat32_t zL = svld1(pMDim, &matLeft_mod[matLeft_pos + k * SVL]);
                svfloat32_t zR = svld1(pNDim, &matRight[matRight_UL_corner + k * N]);
                svmopa_za32_m(0, pMDim, pNDim, zL, zR);
            }

            // Store ZA to matResult.
            const uint64_t result_tile_UL_corner = row * N + col;
            for (uint64_t trow = 0; trow < SVL && row + trow < M; trow += 4) {
                svbool_t p0 = svpsel_lane_b32(pNDim, pMDim, row + trow + 0);
                svbool_t p1 = svpsel_lane_b32(pNDim, pMDim, row + trow + 1);
                svbool_t p2 = svpsel_lane_b32(pNDim, pMDim, row + trow + 2);
                svbool_t p3 = svpsel_lane_b32(pNDim, pMDim, row + trow + 3);

                svst1_hor_za32(/* tile: */ 0, /* slice: */ trow + 0, p0, &matResult[result_tile_UL_corner + (trow + 0) * N]);
                svst1_hor_za32(/* tile: */ 0, /* slice: */ trow + 1, p1, &matResult[result_tile_UL_corner + (trow + 1) * N]);
                svst1_hor_za32(/* tile: */ 0, /* slice: */ trow + 2, p2, &matResult[result_tile_UL_corner + (trow + 2) * N]);
                svst1_hor_za32(/* tile: */ 0, /* slice: */ trow + 3, p3, &matResult[result_tile_UL_corner + (trow + 3) * N]);
            }
        }
    }
}
```

Note again that the `matmul_intr_impl` function has been annotated at line 4 with:

- `__arm_streaming`, because the function is using streaming instructions
- `__arm_inout("za")`, because the function reuses the ZA storage from its caller

The multiplication with the outer product is performed in a double-nested loop, over the `M` (line 7) and `N` (line 11) dimensions of the input matrices `matLeft_mod` and `matRight`. Both loops have an `SVL` step increment, which corresponds to the horizontal and vertical dimensions of the ZA storage.

Once again, intrinsics make it easy to fully leverage SME2, provided you have a solid understanding of its available instructions. 

### Compile and run

The main function is exactly the same that was used for the assembly version, with the `IMPL` macro defined to be `intr` in the `Makefile`.

First, make sure that the `sme2_matmul_intr` executable is up-to-date:

```bash
ninja -C build-native/ sme2_matmul_intr
```

Then execute `sme2_matmul_intr` either natively or on the FVP:

```bash
./build-native/sme2_matmul_intr
```

```bash
adb push build-android/sme2_matmul_intr /data/local/tmp
adb shell chmod 755 /data/local/tmp/sme2_matmul_intr
adb shell /data/local/tmp/sme2_matmul_intr
```

```bash
docker run --rm -v "$PWD:/work" armswdev/sme2-learning-path:sme2-environment-v3 ./run-fvp.sh build-baremetal/sme2_matmul_intr
```

As with the `sme2_matmul_asm` program, you can provide the `M`, `K` and `N` parameters on the command line to `sme2_matmul_intr`.