[RfC] RISC-V port and universal binary

Performance on Spacemit K3, thanks @edolnx for testing

master:
Total time (ms) : 65200
Nodes searched  : 3493826
Nodes/second    : 53586

riscv-scalable-port:
Total time (ms) : 15834
Nodes searched  : 3493826
Nodes/second    : 220653

Also thanks to @camel-cdr for guidance on RVV programming, and https://cloud-v.co for supplying an RVV instance to test with

passed STC:
LLR: 2.81 (-2.94,2.94) <0.00,2.00>
Total: 1152 W: 527 L: 108 D: 517
Ptnml(0-2): 0, 17, 167, 348, 44
https://tests.stockfishchess.org/tests/view/6a39895b3036e45021aeb368

## Summary

We've had a `riscv64` target for a while, but haven't really optimized for it, in particular the vector extension (RVV).

RVV, like SVE, is based on a scalable vector system where the vector length ranges from 128 to 65536. In practice implementations are between 128 and 2048, and 256 bits is quite common (e.g. the Spacemit K3 system above). Unfortunately this doesn't fit well into the rest of our code which assumes a fixed vector length, so what I've done is bypass the `VECTOR` ifdef (which now basically means "FIXED_LENGTH_VECTOR") and just have RVV-specific paths.

The ability to explicitly control `vl` makes the code quite readable, in my opinion. We use LMUL>1 in most places to take advantage of multi-vector instructions. Generally the LMULs were chosen to best support a 256-bit vlen, which is very common, but by virtue of how the vlen control works, the code works with any vlen. In a couple places, i.e., `get_changed_pieces` and `AffineTransformSparseInput::propagate`, we have separate implementations depending on the vlen, because the optimal LMUL varies a lot between implementations.

One little wrinkle is that `load_as` is compiled to a sequence of byte loads, because although unaligned loads are legal in RVA23, the spec says that they *may* be extremely slow (even though they usually aren't, in actual hw), so compilers are conservative. Thus I aligned the relevant buffers and made the semantics of `load_as` that the operand is aligned, by adding a runtime assertion.

### Universal binary

Adding a universal binary is pretty easy and we can just cross-compile. There are two targets: baseline rv64gc and riscv64-rva23, which is actually a smaller subset of RVA23 that also works on some older processors that don't support the full thing. We use clang because GCC, until recently, has a nasty bug with LTO and RVV.

Like the universal ARM and x86 builds, we check all the builds in CI. In this case we run bench with multiple vlens, 128 through 1024.

In the meantime I deleted the existing broken and unused riscv64 tests.

### Follow-ups

- Optimizations
- zvdot4a8i path

closes https://github.com/official-stockfish/Stockfish/pull/6920

No functional change

src/nnue/layers/affine_transform.h

@@ -105,6 +105,26 @@ affine_transform_non_ssse3(i32* output, const i8* weights, const i32* biases, co
 
         #endif
     }
+    #elif defined(USE_RVV)
+    for (IndexType i = 0; i < OutputDimensions; ++i)
+    {
+        const i8*  row  = &weights[i * PaddedInputDimensions];
+        vint32m1_t vsum = __riscv_vmv_v_x_i32m1(0, __riscv_vsetvlmax_e32m1());
+
+        for (usize j = 0; j < InputDimensions;)
+        {
+            usize vl = __riscv_vsetvl_e8m4(InputDimensions - j);
+
+            vint8m4_t  w    = __riscv_vle8_v_i8m4(&row[j], vl);
+            vuint8m4_t x    = __riscv_vle8_v_u8m4(&input[j], vl);
+            vint16m8_t prod = __riscv_vwmulsu_vv_i16m8(w, x, vl);
+
+            vsum = __riscv_vwredsum_vs_i16m8_i32m1(prod, vsum, vl);
+            j += vl;
+        }
+
+        output[i] = biases[i] + __riscv_vmv_x_s_i32m1_i32(vsum);
+    }
     #else
     std::memcpy(output, biases, sizeof(i32) * OutputDimensions);
 

src/nnue/layers/affine_transform_sparse_input.h

@@ -59,7 +59,8 @@ class AffineTransformSparseInput {
     static constexpr IndexType PaddedOutputDimensions =
       ceil_to_multiple<IndexType>(OutputDimensions, MaxSimdWidth);
 
-#if (defined(USE_SSSE3) || defined(USE_LSX) || defined(USE_LASX) || (USE_NEON >= 8))
+#if (defined(USE_SSSE3) || defined(USE_LSX) || defined(USE_LASX) || (USE_NEON >= 8) \
+     || defined(USE_RVV))
     static constexpr IndexType ChunkSize = 4;
 #else
     static constexpr IndexType ChunkSize = 1;
@@ -82,7 +83,8 @@ class AffineTransformSparseInput {
     }
 
     static constexpr IndexType get_weight_index(IndexType i) {
-#if (defined(USE_SSSE3) || defined(USE_LSX) || defined(USE_LASX) || (USE_NEON >= 8))
+#if (defined(USE_SSSE3) || defined(USE_LSX) || defined(USE_LASX) || (USE_NEON >= 8) \
+     || defined(USE_RVV))
         return get_weight_index_scrambled(i);
 #else
         return i;
@@ -275,6 +277,51 @@ class AffineTransformSparseInput {
     #ifdef vec_add_32
         #undef vec_add_32
     #endif
+#elif defined(USE_RVV)
+        static_assert(InputDimensions % 256 == 0);
+
+        const i8* weights_cp = weights;
+
+    #define RVV_SPARSE_PROPAGATE(LMUL) \
+        do \
+        { \
+            const usize blk = __riscv_vsetvlmax_e32m##LMUL(); \
+            for (IndexType ob = 0; ob < OutputDimensions; ob += blk) \
+            { \
+                const usize       vl  = __riscv_vsetvl_e32m##LMUL(OutputDimensions - ob); \
+                vint32m##LMUL##_t acc = __riscv_vle32_v_i32m##LMUL(biases + ob, vl); \
+                for (IndexType k = 0; k < InputDimensions / 256; ++k) \
+                { \
+                    u64   bits         = load_as<u64>(nnzInfo.bitset + k * 8); \
+                    isize base         = k * 64; \
+                    auto* base_addr    = input + base * sizeof(i32); \
+                    auto* weights_base = &weights_cp[base * OutputDimensions * ChunkSize]; \
+                    while (bits) \
+                    { \
+                        isize             i = pop_lsb(bits); \
+                        vuint8m##LMUL##_t a = __riscv_vreinterpret_v_u32m##LMUL##_u8m##LMUL( \
+                          __riscv_vmv_v_x_u32m##LMUL(load_as<u32>(base_addr + i * sizeof(i32)), \
+                                                     vl)); \
+                        vint8m##LMUL##_t b = __riscv_vle8_v_i8m##LMUL( \
+                          &weights_base[i * OutputDimensions * ChunkSize + ob * ChunkSize], \
+                          vl * ChunkSize); \
+                        acc = \
+                          __riscv_vadd_vv_i32m##LMUL(acc, SIMD::rvv_dpbusd_m##LMUL(a, b, vl), vl); \
+                    } \
+                } \
+                __riscv_vse32_v_i32m##LMUL(output + ob, acc, vl); \
+            } \
+        } while (0)
+
+        // Select LMUL
+        if (__riscv_vsetvlmax_e32m1() >= OutputDimensions)
+            RVV_SPARSE_PROPAGATE(1);
+        else if (__riscv_vsetvlmax_e32m2() >= OutputDimensions)
+            RVV_SPARSE_PROPAGATE(2);
+        else
+            RVV_SPARSE_PROPAGATE(4);
+
+    #undef RVV_SPARSE_PROPAGATE
 #else
         // Use dense implementation for the other architectures.
         affine_transform_non_ssse3<InputDimensions, PaddedInputDimensions, OutputDimensions>(

src/nnue/layers/clipped_relu.h

@@ -147,6 +147,24 @@ class ClippedReLU {
         }
         constexpr IndexType Start = NumChunks * 16;
 
+#elif defined(USE_RVV)
+
+        for (usize j = 0; j < InputDimensions;)
+        {
+            usize vl = __riscv_vsetvl_e32m4(InputDimensions - j);
+
+            vint32m4_t in = __riscv_vle32_v_i32m4(&input[j], vl);
+            in            = __riscv_vmax_vx_i32m4(in, 0, vl);
+
+            vint16m2_t words =
+              __riscv_vnclip_wx_i16m2(in, WeightScaleBitsLocal, __RISCV_VXRM_RDN, vl);
+            vint8m1_t narrowed = __riscv_vnclip_wx_i8m1(words, 0, __RISCV_VXRM_RDN, vl);
+
+            __riscv_vse8_v_u8m1(&output[j], __riscv_vreinterpret_v_i8m1_u8m1(narrowed), vl);
+            j += vl;
+        }
+        constexpr IndexType Start = InputDimensions;
+
 #else
         constexpr IndexType Start = 0;
 #endif

src/nnue/layers/sqr_clipped_relu.h

@@ -143,6 +143,22 @@ class SqrClippedReLU {
         }
         constexpr IndexType Start = NumChunks * 16;
 
+#elif defined(USE_RVV)
+
+        for (usize j = 0; j < InputDimensions;)
+        {
+            usize      vl = __riscv_vsetvl_e32m4(InputDimensions - j);
+            vint32m4_t in = __riscv_vle32_v_i32m4(&input[j], vl);
+
+            vint16m2_t words   = __riscv_vnclip_wx_i16m2(in, 0, __RISCV_VXRM_RDN, vl);
+            vint16m2_t sqr     = __riscv_vmulh_vv_i16m2(words, words, vl);
+            vint8m1_t narrowed = __riscv_vnclip_wx_i8m1(sqr, SimdShiftAmount, __RISCV_VXRM_RDN, vl);
+
+            __riscv_vse8_v_u8m1(&output[j], __riscv_vreinterpret_v_i8m1_u8m1(narrowed), vl);
+            j += vl;
+        }
+        constexpr IndexType Start = InputDimensions;
+
 #else
         constexpr IndexType Start = 0;
 #endif

src/nnue/nnue_accumulator.cpp

@@ -288,6 +288,64 @@ void apply_combined(Color                              perspective,
             vec_store_psqt(&toTilePsqt[k], psqt[k]);
     }
 
+#elif defined(USE_RVV)
+
+    usize tileOffset = 0;
+
+    const auto* psqWeights        = &featureTransformer.weights[0];
+    const auto* threatWeights     = &featureTransformer.threatWeights[0];
+    const auto* psqtWeights       = &featureTransformer.psqtWeights[0];
+    const auto* threatPsqtWeights = &featureTransformer.threatPsqtWeights[0];
+
+    while (tileOffset < Dimensions)
+    {
+        usize vl = __riscv_vsetvl_e16m8(Dimensions - tileOffset);
+
+        vint16m8_t accum = __riscv_vle16_v_i16m8(&fromAcc[tileOffset], vl);
+        for (int i : psqRemoved)
+            accum = __riscv_vsub_vv_i16m8(
+              accum, __riscv_vle16_v_i16m8(&psqWeights[i * Dimensions + tileOffset], vl), vl);
+        for (int i : psqAdded)
+            accum = __riscv_vadd_vv_i16m8(
+              accum, __riscv_vle16_v_i16m8(&psqWeights[i * Dimensions + tileOffset], vl), vl);
+        for (int i : thrRemoved)
+            accum = __riscv_vwsub_wv_i16m8(
+              accum, __riscv_vle8_v_i8m4(&threatWeights[i * Dimensions + tileOffset], vl), vl);
+        for (int i : thrAdded)
+            accum = __riscv_vwadd_wv_i16m8(
+              accum, __riscv_vle8_v_i8m4(&threatWeights[i * Dimensions + tileOffset], vl), vl);
+        __riscv_vse16_v_i16m8(&toAcc[tileOffset], accum, vl);
+
+        tileOffset += vl;
+    }
+
+    tileOffset = 0;
+
+    while (tileOffset < PSQTBuckets)
+    {
+        usize vl = __riscv_vsetvl_e32m1(PSQTBuckets - tileOffset);
+
+        vint32m1_t accum = __riscv_vle32_v_i32m1(&fromPsqtAcc[tileOffset], vl);
+        for (int i : psqRemoved)
+            accum = __riscv_vsub_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&psqtWeights[i * PSQTBuckets + tileOffset], vl), vl);
+        for (int i : psqAdded)
+            accum = __riscv_vadd_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&psqtWeights[i * PSQTBuckets + tileOffset], vl), vl);
+        for (int i : thrRemoved)
+            accum = __riscv_vsub_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&threatPsqtWeights[i * PSQTBuckets + tileOffset], vl),
+              vl);
+        for (int i : thrAdded)
+            accum = __riscv_vadd_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&threatPsqtWeights[i * PSQTBuckets + tileOffset], vl),
+              vl);
+
+        __riscv_vse32_v_i32m1(&toPsqtAcc[tileOffset], accum, vl);
+
+        tileOffset += vl;
+    }
+
 #else
 
     toAcc     = fromAcc;
@@ -448,6 +506,25 @@ Bitboard get_changed_pieces(const std::array<Piece, SQUARE_NB>& oldPieces,
     }
 
     return ~sameBB;
+#elif defined(USE_RVV)
+
+    #define IMPL(mx, bx) \
+        return __riscv_vmv_x_s_u64m1_u64(__riscv_vreinterpret_v_u8m1_u64m1( \
+          __riscv_vreinterpret_v_b##bx##_u8m1(__riscv_vmsne_vv_i8m##mx##_b##bx( \
+            __riscv_vle8_v_i8m##mx(reinterpret_cast<const i8*>(oldPieces.data()), 64), \
+            __riscv_vle8_v_i8m##mx(reinterpret_cast<const i8*>(newPieces.data()), 64), 64))))
+
+
+    usize vl = __riscv_vsetvlmax_e8m1();
+    if (vl >= 64)
+        IMPL(1, 8);
+    else if (vl == 32)
+        IMPL(2, 4);
+    else
+        IMPL(4, 2);
+
+    #undef IMPL
+
 #else
     Bitboard changed = 0;
 
@@ -596,6 +673,64 @@ void update_accumulator_refresh_cache(Color                     perspective,
             vec_store_psqt(&accTilePsqt[k], psqt[k]);
     }
 
+#elif defined(USE_RVV)
+
+    const auto* weights           = &featureTransformer.weights[0];
+    const auto* threatWeights     = &featureTransformer.threatWeights[0];
+    const auto* psqtWeights       = &featureTransformer.psqtWeights[0];
+    const auto* threatPsqtWeights = &featureTransformer.threatPsqtWeights[0];
+
+    usize tileOffset = 0;
+
+    while (tileOffset < Dimensions)
+    {
+        usize vl = __riscv_vsetvl_e16m8(Dimensions - tileOffset);
+
+        vint16m8_t accum = __riscv_vle16_v_i16m8(&entry.accumulation[tileOffset], vl);
+        for (int i : removed)
+            accum = __riscv_vsub_vv_i16m8(
+              accum, __riscv_vle16_v_i16m8(&weights[i * Dimensions + tileOffset], vl), vl);
+        for (int i : added)
+            accum = __riscv_vadd_vv_i16m8(
+              accum, __riscv_vle16_v_i16m8(&weights[i * Dimensions + tileOffset], vl), vl);
+
+        __riscv_vse16_v_i16m8(&entry.accumulation[tileOffset], accum, vl);
+
+        for (int i : active)
+            accum = __riscv_vwadd_wv_i16m8(
+              accum, __riscv_vle8_v_i8m4(&threatWeights[i * Dimensions + tileOffset], vl), vl);
+
+        __riscv_vse16_v_i16m8(&accumulator.accumulation[perspective][tileOffset], accum, vl);
+
+        tileOffset += vl;
+    }
+
+    tileOffset = 0;
+
+    while (tileOffset < PSQTBuckets)
+    {
+        usize vl = __riscv_vsetvl_e32m1(PSQTBuckets - tileOffset);
+
+        vint32m1_t accum = __riscv_vle32_v_i32m1(&entry.psqtAccumulation[tileOffset], vl);
+        for (int i : removed)
+            accum = __riscv_vsub_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&psqtWeights[i * PSQTBuckets + tileOffset], vl), vl);
+        for (int i : added)
+            accum = __riscv_vadd_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&psqtWeights[i * PSQTBuckets + tileOffset], vl), vl);
+
+        __riscv_vse32_v_i32m1(&entry.psqtAccumulation[tileOffset], accum, vl);
+
+        for (int i : active)
+            accum = __riscv_vadd_vv_i32m1(
+              accum, __riscv_vle32_v_i32m1(&threatPsqtWeights[i * PSQTBuckets + tileOffset], vl),
+              vl);
+
+        __riscv_vse32_v_i32m1(&accumulator.psqtAccumulation[perspective][tileOffset], accum, vl);
+
+        tileOffset += vl;
+    }
+
 #else
 
     for (const auto index : removed)

src/nnue/nnue_feature_transformer.h

@@ -342,6 +342,40 @@ class FeatureTransformer {
                 cursor.record2(packed[0], packed[1]);
             }
 
+#elif defined(USE_RVV)
+
+            usize j = 0;
+
+            while (j < HalfDimensions / 2)
+            {
+                usize vl = __riscv_vsetvl_e16m8(HalfDimensions / 2 - j);
+
+                vint16m8_t acc0 = __riscv_vle16_v_i16m8(&accumulation[perspectives[p]][j], vl);
+                vint16m8_t acc1 =
+                  __riscv_vle16_v_i16m8(&accumulation[perspectives[p]][j + HalfDimensions / 2], vl);
+
+                acc0 = __riscv_vmax_vx_i16m8(acc0, 0, vl);
+                acc1 = __riscv_vmax_vx_i16m8(acc1, 0, vl);
+
+                vuint8m4_t pa = __riscv_vnclipu_wx_u8m4(__riscv_vreinterpret_v_i16m8_u16m8(acc0), 0,
+                                                        __RISCV_VXRM_RDN, vl);
+                vuint8m4_t pb = __riscv_vnclipu_wx_u8m4(__riscv_vreinterpret_v_i16m8_u16m8(acc1), 0,
+                                                        __RISCV_VXRM_RDN, vl);
+
+                vuint8m4_t hi     = __riscv_vmulhu_vv_u8m4(pa, pb, vl);
+                vuint8m4_t result = __riscv_vsrl_vx_u8m4(hi, 1, vl);
+
+                __riscv_vse8_v_u8m4(&output[offset + j], result, vl);
+                j += vl;
+
+                // Record NNZ
+                usize    vl32 = vl / 4;
+                vbool8_t nnzMask =
+                  __riscv_vmsne_vx_u32m4_b8(__riscv_vreinterpret_v_u8m4_u32m4(result), 0, vl32);
+                __riscv_vsm_v_b8(cursor.out, nnzMask, vl32);
+                cursor.out += vl32 / 8;
+            }
+
 #else
 
             for (IndexType j = 0; j < HalfDimensions / 2; ++j)

src/nnue/nnz_helper.h

@@ -105,7 +105,7 @@ struct NNZInfo {
     NNZCursor make_cursor(bool perspective) { return {*this, perspective, count}; }
 #else
     // Each 8-bit chunk
-    u8 bitset[(Dimensions + 31) / 32];
+    alignas(8) u8 bitset[(Dimensions + 31) / 32];
 
     struct NNZCursor {
         u8* out;

src/nnue/simd.h

@@ -40,6 +40,10 @@
 
 #elif defined(USE_LSX)
     #include <lsxintrin.h>
+
+#elif defined(USE_RVV)
+    #include <riscv_vector.h>
+
 #endif
 
 #include "../types.h"
@@ -635,6 +639,35 @@ class SIMDTiling {
     static_assert(PSQTBuckets % PsqtTileHeight == 0, "PsqtTileHeight must divide PSQTBuckets");
 #endif
 };
+
+#if defined(USE_RVV)
+
+    #define RVV_DEFINE_DPBUSD(A, W, H) \
+        inline vint32##A##_t rvv_dpbusd_##A(vuint8##A##_t a, vint8##A##_t b, usize n) { \
+            vint16##W##_t  prod   = __riscv_vwmulsu_vv_i16##W(b, a, 4 * n); \
+            vuint32##W##_t prod32 = __riscv_vreinterpret_v_u16##W##_u32##W( \
+              __riscv_vreinterpret_v_i16##W##_u16##W(prod)); \
+            vint16##A##_t even = \
+              __riscv_vreinterpret_v_u16##A##_i16##A(__riscv_vnsrl_wx_u16##A(prod32, 0, 2 * n)); \
+            vint16##A##_t odd = \
+              __riscv_vreinterpret_v_u16##A##_i16##A(__riscv_vnsrl_wx_u16##A(prod32, 16, 2 * n)); \
+            vuint32##A##_t pairs = __riscv_vreinterpret_v_u16##A##_u32##A( \
+              __riscv_vreinterpret_v_i16##A##_u16##A(__riscv_vadd_vv_i16##A(even, odd, 2 * n))); \
+            vint16##H##_t lo = \
+              __riscv_vreinterpret_v_u16##H##_i16##H(__riscv_vnsrl_wx_u16##H(pairs, 0, n)); \
+            vint16##H##_t hi = \
+              __riscv_vreinterpret_v_u16##H##_i16##H(__riscv_vnsrl_wx_u16##H(pairs, 16, n)); \
+            return __riscv_vwadd_vv_i32##A(lo, hi, n); \
+        }
+
+RVV_DEFINE_DPBUSD(m1, m2, mf2)
+RVV_DEFINE_DPBUSD(m2, m4, m1)
+RVV_DEFINE_DPBUSD(m4, m8, m2)
+
+    #undef RVV_DEFINE_DPBUSD
+
+#endif
+
 }
 
 #endif