mpbench — 多倍長計算ライブラリ性能評価のコード一式#

記事「多倍長数値計算の性能評価」で使ったベンチマークコードと測定データ。測定仕様は common/SPEC.md、共通乱数は common/xorshift.md。

benchmark_misc/
├── common/     共通仕様 (SPEC.md), 共通乱数 (xorshift.md, xorshift_ref.py)
├── cpp/        C/C++: FLINT (eig, fft), MPLAPACK (eig, gemm), Eigen+mpreal, MPFR+OpenMP FFT, スカラー演算, Makefile
├── julia/      Julia: GenericSchur, GenericLinearAlgebra, MultiFloats, DoubleFloats, Arblib, GenericFFT, FFTW
├── python/     Python: mpmath (+gmpy2), python-flint
├── matlab/     MATLAB + Advanpix: advanpix_bench.m (課題 A/B/C), advanpix_scaling.m (時間・メモリ vs N), advanpix_dense.m (密行列, MPLAPACK と同条件)
├── wolfram/    Mathematica (wolframscript)
├── scripts/    env_info.sh, run_all.sh, analyze.py, run_advanpix_scaling.sh, fit_advanpix.py
└── results/    results.csv (全測定値), env.txt (測定環境)

すべてのスクリプトは同一スキーマの CSV（timestamp,host,library,version,language,task,subtask,digits,n,threads,trial,time_s,check_value,notes）に追記し、scripts/analyze.py が表（summary.md）と図（fig_*.png）を生成する。

準備（Debian / Ubuntu 系の例）#

システムパッケージ#

apt install build-essential libgmp-dev libmpfr-dev libmpc-dev libflint-dev libeigen3-dev libboost-dev libfftw3-dev python3-pip

libflint-dev は 3.0 以上（arb 統合後）が必要。QD（double-double / quad-double）は MPLAPACK 同梱のものを使うので libqd-dev は不要（下記）。

mpreal.h（Eigen と scalar_bench で使用）#

Pavel Holoborodko 氏の mpreal.h（http://www.holoborodko.com/pavel/mpfr/ 、GPL）を cpp/mpreal/ に置く。

git clone https://github.com/advanpix/mpreal.git cpp/mpreal

MPLAPACK 3.0.0#

wget https://github.com/nakatamaho/mplapack/releases/download/v3.0.0/mplapack-3.0.0.tar.xz
tar xf mplapack-3.0.0.tar.xz && cd mplapack-3.0.0
./configure --prefix=$HOME/opt/mplapack --enable-gmp=yes --enable-mpfr=yes --enable-dd=yes --enable-qd=yes --enable-binary128=yes --enable-double=yes --enable-optimization=yes --enable-gmpfrxx-mkII --with-system-gmp=yes --with-system-mpfr=yes --with-system-mpc=yes --with-cxx-standard=17
make -j$(nproc) && make install

注意点:

  • configure は 1 行で実行する（\ 継続と行末コメントを混ぜると # が引数に渡って invalid host type: # になる）。

  • --with-system-qd は付けない。dd/qd の複素型（dd_complex, qd_complex）は MPLAPACK 同梱の QD にしか無く、システムの libqd では qd/qd_complex.h が見つからずビルドが失敗する。

  • --enable-optimization で OpenMP 版の libmplapack_<type>_opt.so ができる。3.0 では BLAS と LAPACK が一つのライブラリにまとまっており（libmpblas_* は無い）、cpp/Makefile はそれに合わせてある（MPL_OPT= を付けると参照版にリンク）。

  • gmpfrxx_mkII の GMP バックエンドは mpf_set_default_prec を無視する。精度は gmpxx::set_default_mpf_precision_bits(bits)（MPFR は mpfrxx::set_default_precision_bits）で設定する。cpp/mplapack_eig.cpp はそのようにしてある。

  • 型名は mpfr_class / mpc_class（MPFR）、mpf_class / mpfc_class（GMP）、mplapack_binary128_t（binary128、ヘッダは mpblas_binary128.h）。

Julia パッケージ#

julia -e 'using Pkg; Pkg.add(["GenericSchur","GenericLinearAlgebra","MultiFloats","DoubleFloats","Arblib","GenericFFT","FFTW"])'

MultiFloats は MultiFloat{Float64,2} と {Float64,4} を使う（N=8 は算術が未定義のバージョンがある）。超越関数は MultiFloats.use_bigfloat_transcendentals() で MPFR にフォールバックさせている。

Python パッケージ#

pip3 install mpmath gmpy2 python-flint pandas matplotlib tabulate

（mpmath は gmpy2 があると自動でバックエンドに使う。python-flint は conda より pip の方が確実。tabulate は analyze.py の Markdown 表出力に必要）

MATLAB / Mathematica#

Advanpix のパスは ~/Documents/MATLAB/startup.m で addpath しておく（ベンチマークの関数ファイルの先頭に書くと「Local function name must be different from the script name」になる）。matlab と wolframscript は PATH に無ければ MATLAB=...、WOLFRAMSCRIPT=... の環境変数で run_all.sh に渡す。

ビルド#

make -C cpp flint mpfrfft scalar eigen MPREAL_DIR=./mpreal HAVE_BOOST=1
make -C cpp mplapack mplapack_gemm MPLAPACK=$HOME/opt/mplapack

実行#

bash scripts/env_info.sh          # 環境記録 → results/env_host>.txt
bash scripts/run_all.sh smoke     # 全スクリプトの疎通確認 (数分)
bash scripts/run_all.sh scalar    # 課題 A
bash scripts/run_all.sh fft       # 課題 C
bash scripts/run_all.sh eig       # 課題 B (長時間。下記の環境変数で分割する)
python3 scripts/analyze.py        # results/summary.md と fig_*.png

run_all.sh の環境変数（既定値）:

変数

意味

THREADS=1,2,4,8,16,32,64,<nproc>

FLINT / Julia / Advanpix のスレッド数スイープ（Julia の eig と fft はスレッドを使わないので最初の段だけ）

MPL_THREADS=1,8,<nproc>

MPLAPACK のスレッド数（Rgeev/Cheev は並列 BLAS の恩恵をほぼ受けないので段数を絞ってある）

N_EIG=10,20,50,...,1000 / N_FFT=... / DIGITS=100,1000 / TRIALS=3

サイズ・桁数・試行回数

ONLY=mplapack,flint,julia,advanpix,seq

回すライブラリ群を絞る（seq = Eigen, GenericLinearAlgebra, mpmath, Mathematica）

VERIFIED=

FLINT の保証付き固有値（eig_multiple）を省く（1000 桁 N=300 で 3 時間かかる）

例: ONLY=advanpix THREADS=32,64 N_EIG=500,700,1000 bash scripts/run_all.sh eig

Advanpix の時間・メモリ vs N と見積り式#

MATLAB=... HERM=1 bash scripts/run_advanpix_scaling.sh      # 1 点 1 プロセスで N=100〜2000 (既定 64 スレッド)
python3 scripts/fit_advanpix.py results/advanpix_scaling_<host>.csv   # フィットと見積り表 → results/advanpix_fit.md, fig_advanpix_*.png

MPLAPACK の Rgemm / Rgehrd / Rgeev（密行列）と Advanpix の同条件比較#

for t in 1 8 64; do OMP_NUM_THREADS=$t cpp/mplapack_gemm --digits 100,1000 --n 300,1000 --threads $t --trials 1 --out results/results_<host>.csv; done
matlab -batch "addpath('matlab'); advanpix_dense('digits',[100 1000],'n',[300 1000],'threads',[1 8 64],'out','results/results_<host>.csv')"

個別実行の例#

OMP_NUM_THREADS=8 cpp/mplapack_eig_mpfr --digits 100 --n 100,200 --threads 8 --check --herm
cpp/flint_eig --digits 100,1000 --n 10,50,100 --threads 1 --verified
cpp/mpfr_fft_omp --digits 1000 --n 1048576 --threads 16
julia -t 1 julia/eig_bench.jl --libs genericschur,multifloats --digits 100 --n 100 --herm
python3 python/eig_bench.py --digits 100 --n 100 --verified --herm
matlab -batch "addpath('matlab'); advanpix_bench('task','eig','digits',100,'n',[100 200],'threads',8,'herm',true)"
wolframscript -file wolfram/bench.wls task=fft digits=100 n=65536

すべてのスクリプトは --smoke（MATLAB は 'smoke',true、Wolfram は smoke=1）で最小設定の動作確認ができる。

測定上の注意（記事本文の「落とし穴」も参照）#

  • スカラー演算と gemm の入力は sqrt(u) で仮数を目標精度いっぱいに埋めている。倍精度からの変換値をそのまま使うと arb や MPFR が短い仮数で速く見える。

  • Circus 行列は上 Hessenberg 形なので、対角化の時間に Hessenberg 化は含まれない。

  • Advanpix は物理コア数を超えるスレッド（SMT）で大きく遅くなる。maxNumCompThreads を物理コア数にする。

  • MPLAPACK の _opt 版は小さい行列でスレッド起動のオーバーヘッドが支配的になるので、OMP_NUM_THREADS を明示する（未設定だと全論理コアを使う）。

  • Julia の GenericSchur と GenericLinearAlgebra は同時に読み込むと衝突するので別プロセスで走らせる（run_all.sh は分けてある）。