RISC-V Vectorized FFTW3を試す (2. RISC-V向に、GCC13を使ってビルド)

前回、OfficialのリリースパッケージでRISC-V GCCを作成したのは、最新版に対応したRVV 1.0 Intrinsicを使いたかったからだった。

msyksphinz.hatenablog.com

これを使って、FFTW3をビルドし直すことにする。

新しいGCCを使用して前回のブログで使用したFFTW3をビルドすると、とりあえず何も修正なしでFFTW3をビルド可能だ。

git clone https://github.com/sh-zheng/fftw3.git fftw3-zheng
cd fftw3-zheng
./bootstrap.sh --host=riscv64-unknown-elf --disable-threads
./configure --enable-maintainer-mode --enable-rvv --host=riscv64-unknown-elf "CC=riscv64-unknown-elf-gcc -march=rv64gcv -O0 -g" "CXX=riscv64-unknown-elf-g++ -march=rv64gcv -O0 -g"
make -j$(nproc)

これを使ってサンプルアプリケーションを作ってみる。以下の記事を参考にして、FFTW3をライブラリとして呼び出してみる。

auewe.hatenablog.com

• fourier.c

/***************************/
/* FourierTransformation.c */
/***************************/

// 要素数Nの数列 a[m] をフーリエ変換で b[n] にする。添字は 0,...,N-1
// b[n] = sum[m=0,N-1]  a[m]*exp(-2*PI*I* (m/N) *n)

#include <stdio.h>
#include <complex.h> // complex.h は fftw3.h より先に include する
#include <fftw3.h>   // windows環境では #include "C:/path/to/fftw3.h"
                     // あるいは        #include "./相対パス/fftw3.h"
#include "count_utils.h"


int main( void )
{
  
  int N=32;

  // a,b は double _Complex 型のC99標準複素配列と実質的に同じ
  // double _Complex a[4] としても動くけど計算速度が低下する可能性あり
  fftw_complex *a, *b;
  a = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  b = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

  // プランの生成
  // フーリエ逆変換つまり位相因子を exp(-k)ではなく exp(+k)とする場合は
  // FFTW_FORWARD の代わりに FFTW_BACKWARD とする
  fftw_plan plan;
  plan = fftw_plan_dft_1d( N, a, b, FFTW_FORWARD, FFTW_ESTIMATE );
  
  // フーリエ変換前の数列値を設定
  for (int i = 0; i < N; i++) {
    a[i] = (float)i + 0.0*I;
  }
 
  // a[0] = 1.0 + 0.0*I;
  // a[1] = 2.0 + 0.0*I;
  // a[2] = 5.0 + 0.0*I;
  // a[3] = 3.0 + 0.0*I;

  // フーリエ変換実行   b[n]に計算結果が入る
  long long start_cycle   = get_cycle();
  // long long start_vecinst = get_vecinst();
  long long start_inst    = get_instret();

  fftw_execute(plan);

  long long end_cycle   = get_cycle();
  // long long end_vecinst = get_vecinst();
  long long end_inst    = get_instret();

  // b[n]の値を表示
  int n;
  for( n=0; n<N; n++ ){
    printf("b_%d = %+lf %+lf*i\n", n, creal(b[n]), cimag(b[n]) );
  }

  // ここで a[m] の値を変えて再度 fftw_execute(plan) を実行すれば、
  // b[n] が再計算される。

  // 計算終了時、メモリ開放を忘れないように
  if(plan) fftw_destroy_plan(plan);
  fftw_free(a); fftw_free(b);

  printf("allinst = %lld\n", end_inst - start_inst);
  // printf("vecinst = %lld\n", end_vecinst - start_vecinst);
  printf("cycles = %lld\n", end_cycle - start_cycle);

  return 0;
}

riscv64-unknown-elf-gcc -I. -Wall -Werror -O2 \
    -I/home/msyksphinz/work/riscv/fftw/fftw3-zheng/api \
    -o fourier fourier.o \
    /home/msyksphinz/work/riscv/fftw/fftw3-zheng/.libs/libfftw3.a \
    -L/home/msyksphinz/work/riscv/fftw/fftw3-zheng/.libs \
    -lfftw3 -lm

Spikeでの実行方法は以下。

spike --isa=rv64gcv --varch="vlen:512,elen:64" pk fourier

bbl loader
b_0 = +496.000000 +0.000000*i
b_1 = -16.000000 +162.450726*i
b_2 = -16.000000 +80.437432*i
b_3 = -16.000000 +52.744931*i
b_4 = -16.000000 +38.627417*i
b_5 = -16.000000 +29.933895*i
b_6 = -16.000000 +23.945692*i
b_7 = -16.000000 +19.496056*i
b_8 = -16.000000 +16.000000*i
b_9 = -16.000000 +13.130861*i
b_10 = -16.000000 +10.690858*i
b_11 = -16.000000 +8.552178*i
b_12 = -16.000000 +6.627417*i
b_13 = -16.000000 +4.853547*i
b_14 = -16.000000 +3.182598*i
b_15 = -16.000000 +1.575862*i
b_16 = -16.000000 +0.000000*i
b_17 = -16.000000 -1.575862*i
b_18 = -16.000000 -3.182598*i
b_19 = -16.000000 -4.853547*i
b_20 = -16.000000 -6.627417*i
b_21 = -16.000000 -8.552178*i
b_22 = -16.000000 -10.690858*i
b_23 = -16.000000 -13.130861*i
b_24 = -16.000000 -16.000000*i
b_25 = -16.000000 -19.496056*i
b_26 = -16.000000 -23.945692*i
b_27 = -16.000000 -29.933895*i
b_28 = -16.000000 -38.627417*i
b_29 = -16.000000 -52.744931*i
b_30 = -16.000000 -80.437432*i
b_31 = -16.000000 -162.450726*i
allinst = 369
cycles = 369