前回の続き。

msyksphinz.hatenablog.com

2. Indirect Memory Access(ギャザー)を RVV 化してみる

次は少し難しいケース。配列 x をインデックスとして、y[x[i]] を足し込む。これは 連続アクセスではない ので、通常の vle ではなく gather(indexed load) が必要になる。

floatima(size_t n,size_t x[],float y[]) {float sum =0;for (size_t i =0; i < n; i++) {
        sum += y[x[i]];
    }return sum;
}

とりあえず同様にAIに依頼すると、最初は次のような実装が出てきた(抜粋)。

vy = __riscv_vluxei64_v_f32m1(y, vx, vl);
vsum = __riscv_vfredusum_vs_f32m1_f32m1(vy, vsum, vl);

しかしこれは そのままだと不完全。理由は vluxei* の「index が何を表すか」を取り違えやすいからだ。

vluxei* の index は「要素番号」ではなく「バイトオフセット」

RVVの vluxei* 系(indexed load)は、(少なくとも一般的には)base アドレス + index(バイト単位のオフセット) でアドレス計算する。

つまり x[i] が「要素番号」なら、sizeof(float) を掛けてバイトオフセットへ変換してから渡す必要がある。

この点を指摘すると、AIは以下のように修正してきた。

vx = __riscv_vle64_v_u64m1(x, vl);// 要素番号 -> バイトオフセットへ変換(float は 4 bytes)vuint64m1_t vx_bytes = __riscv_vmul_vx_u64m1(vx,sizeof(float), vl);

vy = __riscv_vluxei64_v_f32m1(y, vx_bytes, vl);

方向性は正しい。さらに sizeof(float)=4=2^2 なので、定数倍は乗算より左シフトのほうが軽い(場合が多い)。そこで次の形にするのが自然。

vuint64m1_t vx_bytes = __riscv_vsll_vi_u64m1(vx,2, vl);

---

2.2 Reduction(総和)を「どこで」やるか：正しさと速さのトレードオフ

次はここ。AIの初期案はループ内で vfredusum していた。

• ループごとに gather した vy を ベクトルreduce して vsum へ畳み込む
• 最後に vfmv_f_s でスカラに取り出して返す

これは機能的には分かりやすいが、性能面で気になる点がある：

• vfredusum は(実装にもよるが)レイテンシがそれなりに重い 操作になりやすい
• ループ内で毎回 reduce すると、依存鎖(vsum が次の反復に依存) が強くなってパイプラインが詰まりやすい

一方で、浮動小数点の足し算は結合法則が成り立たないため、reduce のやり方を変えると 最終結果が僅かに変わる 可能性がある。

今回の用途では「スカラ版とビット一致が必要」ではなく、多少の差は許容する前提にする。

制約：加算の順序はスカラ版と厳密に同一である必要はない(リダクション順序で数値が多少変わり得る)

---

2.3 追加最適化案：複数アキュムレータで依存を緩和

上の制約を明示したうえで、AIに「さらに最適化」を依頼すると、複数アキュムレータ(vsum0〜vsum3) を使う案が出てきた。狙いは：

• 依存鎖を短くして ILP を稼ぐ
• reduce を 1本に集中させず、反復ごとに別アキュムレータへ分散

考え方としては理解できる(もちろん、実機の reduce 実装やスケジューリング次第で効く/効かないはある)。

最終的に生成されたコードは以下のようになった：

#if __riscv_v_intrinsic >= 1000000
#include <riscv_vector.h>

float ima_rvv(size_t n, const size_t *x, const float *y) {
    size_t vl;
    vfloat32m1_t vsum;
    float sum;

    // reduction 初期値(ベクトル長1でOK)
    size_t vl_init = __riscv_vsetvl_e32m1(1);
    vsum = __riscv_vfmv_v_f_f32m1(0.0f, vl_init);

    for (size_t avl = n; avl > 0; avl -= vl) {
        vl = __riscv_vsetvl_e64m1(avl);

        // index(要素番号)をロード
        vuint64m1_t vidx = __riscv_vle64_v_u64m1(x, vl);

        // 要素番号 -> バイトオフセット(float=4 bytes)
        vuint64m1_t vidx_bytes = __riscv_vsll_vi_u64m1(vidx, 2, vl);

        // gather load: y[ x[i] ]
        vfloat32m1_t vy = __riscv_vluxei64_v_f32m1(y, vidx_bytes, vl);

        // reduce(順序は変わり得る)
        vsum = __riscv_vfredusum_vs_f32m1_f32m1(vy, vsum, vl);

        x += vl;
    }

    sum = __riscv_vfmv_f_s_f32m1_f32(vsum);
    return sum;
}
#endif