Fortran 2018でGPUを高速化する Part 3 ― GPUカーネルのメモリアクセスとFP32・FP64を比較する
はじめに
Part 1では、Fortran 2018のDO CONCURRENTとOpenACCを使って、2次元の格子計算をNVIDIA GPUへ実行させる仕組みを整理しました。Part 2では、各時間ステップでCPUとGPUの間を往復していた配列を、GPUメモリに置いたまま更新する方法へ変更しました。
そこでPart 3では、GPUメモリに置いた後、どう計算するかを調べます。
転送を減らしても、GPUで実行するループの書き方によって速度が変わることがあります。とくにFortranは配列の第1添字がメモリ上で連続する言語です。この特徴がGPU上のスレッド配置とどう関係するか、同じ差分式を3種類のループで計算して比較します。
もう一つのテーマは、単精度real32（FP32）と倍精度real64（FP64）の使い分けです。FP32では配列サイズを小さくできますが、計算精度も変わります。水深や河川水位のような数値を扱うときには、速くなる可能性だけで判断できません。
今回も前回と同じ2次元拡散方程式の5点差分を使い、計算モデルを複雑にしないまま、GPUカーネルの配置と数値誤差に焦点を絞ります。拡散方程式は洪水氾濫の支配方程式そのものではありませんが、格子上で近傍セルを参照し、次の時刻へ進める考え方は共通しています。
検証状況（2026年11月3日）：GNU Fortran 14.2.0でCPU基準版とOpenACCホスト版をコンパイル・実行しました。FP32/FP64の2種類、7/8/31/32ステップの4条件、OpenACCの3方式で、CPU基準版との全格子照合24件が成功しました。Pythonの回帰テスト10件も成功しています。一方、
nvfortranでのGPUコンパイル、Docker Composeでの実機起動、NVIDIA GPU上の性能測定は、この作成環境では行っていません。GPU性能に関する数値は掲載せず、実機で測定するためのコードと手順を用意しました。今後検証する予定です。
1. Part 1・Part 2・Part 3のつながり
Part 2までの実装は、GPU上にuとvの2配列を確保して、奇数ステップと偶数ステップで読み書きする役割を交換する方式でした。
- Part 1：Fortranの並列化構文とGPUオフロードを理解する
- Part 2：CPU↔GPUの転送を減らし、配列をGPUに常駐させる
- Part 3：GPU上のループ配置を比較し、FP32とFP64の精度を測る
図1 Part 1〜3の流れ
今回の目的は「この書き方が必ず最速です」と断定することではありません。コンパイラがGPU上のループをどう割り当てたか、実際の計算時間と誤差がどうなったかを、測れる形にすることです。
2. 今回比較するプログラム
配布ZIPには、CPU基準版とOpenACC版を収録しています。OpenACC版は、実行時の引数で3種類のループに切り替えます。
| 種類 | 実行モード | ループの特徴 | 主な確認対象 |
|---|---|---|---|
| CPU基準版 | cpu |
jの内側でiを更新、全セルで境界チェック |
数値結果の基準 |
| OpenACCベースライン | baseline |
collapse(2)で全セル、境界を分岐処理 |
Part 2の方式に近い構成 |
| OpenACC連続方向候補 | contiguous |
外側jをgang、内側iをvector
|
連続方向へのスレッド配置 |
| OpenACC比較候補 | strided |
外側iをgang、内側jをvector
|
不連続方向への配置 |
FP32とFP64は共通のFortranソースからビルドします。精度の違い以外で、初期値や時間ステップ、境界条件を変更しません。
なお、後の2方式は最適化の候補です。contiguousがベースラインより高速になるとは限りませんし、stridedが必ず遅いとも限りません。GPUのスレッド配置、キャッシュ、コンパイラの判断を含めて、実機で確認する必要があります。
3. 前回と同じ2次元5点差分を使用する
計算対象は次の拡散方程式です。
\frac{\partial u}{\partial t}
=D\left(\frac{\partial^2u}{\partial x^2}
+\frac{\partial^2u}{\partial y^2}\right)
5点差分では、ひとつ前の時刻の中心セルと上下左右の4セルを使って次の値を求めます。
u^{n+1}_{i,j}=u^n_{i,j}+\alpha
\bigl(u^n_{i-1,j}+u^n_{i+1,j}+u^n_{i,j-1}+u^n_{i,j+1}-4u^n_{i,j}\bigr)
今回もalpha=0.1とし、四辺は0固定、最初の中央1セルだけ1とします。等間隔格子の陽的差分として、alpha <= 1/4の安定条件を満たす設定です。
! 各セルの更新式（FP32/FP64共通）
v(i,j) = u(i,j) + alpha * ( &
u(i-1,j) + u(i+1,j) + &
u(i,j-1) + u(i,j+1) - 4.0_rk*u(i,j))
uを読み取りながらvへ書きます。次のステップではvを読んでuを更新します。途中で同じ配列を上書きしてしまうと、計算順序の依存が生じるため、今回も2配列を使います。
前回から変更した点として、Part 3の3方式では、初回転送条件をそろえるために両配列をcopyin(u,v)でGPUへ送ります。vも境界値0で初期化されているため、内部セルだけを計算する方式でも境界を保持できます。Part 2のcopyin(u) create(v)と転送量が異なるので、Part 2とPart 3のTOTAL_Sをそのまま比較しないでください。
4. Fortranの配列順序を確認する
Fortranの2次元配列u(nx,ny)では、第1添字iが先に連続します。C言語の行優先を前提とした書き方と混同しないことが大切です。
例えばu(4,3)のメモリ配置は、概念的には次の順番です。
u(1,1), u(2,1), u(3,1), u(4,1),
u(1,2), u(2,2), u(3,2), u(4,2),
u(1,3), u(2,3), u(3,3), u(4,3)
したがってCPUで配列を連続方向に走査するときには、次の順番が自然です。
! 第1添字iを内側で回す例
do j = 2, ny-1
do i = 2, nx-1
! u(i,j)へのアクセス
end do
end do
図2 Fortranの列優先メモリ配置
GPUの場合、隣接するスレッドが連続したアドレスをアクセスできると、メモリ帯域を使いやすくなることがあります。ただしOpenACCが実際にどのループをGPUのgangやvectorに割り当てるかは、コンパイル時の診断を見て確認します。単にFortranソースのループ順だけからGPU上の実行順を断定しません。
5. GPUカーネルを3方式で比較する
5.1 baseline：collapse(2)で全格子を処理する
前回のGPU常駐版に近い処理です。
!$acc parallel loop collapse(2) present(u,v)
do j = 1, ny
do i = 1, nx
if (i==1 .or. i==nx .or. j==1 .or. j==ny) then
v(i,j) = 0.0_rk
else
v(i,j) = u(i,j) + alpha*( &
u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) &
-4.0_rk*u(i,j))
end if
end do
end do
!$acc end parallel loop
2重ループをまとめて並列化し、境界を各セルで判定します。単純で読みやすい実装です。コンパイラがどのようなスレッド割付を選んだかは-Minfo=accelで確認します。
5.2 contiguous：第1添字をvectorにする候補
!$acc parallel loop gang present(u,v)
do j = 2, ny-1
!$acc loop vector
do i = 2, nx-1
v(i,j) = u(i,j) + alpha*( &
u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) &
-4.0_rk*u(i,j))
end do
end do
!$acc end parallel loop
ここでは境界セルを計算しません。uとvの四辺は初期値0のまま保持されるため、各ステップで境界の分岐を行う必要がありません。ループ配置に加えて、境界分岐もなくしているので、baselineとの差をメモリアクセスだけの効果と解釈してはいけません。
5.3 strided：第2添字をvectorにする比較用コード
!$acc parallel loop gang present(u,v)
do i = 2, nx-1
!$acc loop vector
do j = 2, ny-1
v(i,j) = u(i,j) + alpha*( &
u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1) &
-4.0_rk*u(i,j))
end do
end do
!$acc end parallel loop
こちらも内部セルだけを更新しますが、vectorの対象をjに変えています。contiguousとstridedは演算と境界処理が同じなので、ループの向きによる差を確認しやすい組になっています。
図3 GPUカーネルの3方式
5.4 時間ステップごとの配列切り替えは継続する
Part 2と同じように、奇数ステップはu → v、偶数ステップはv → uです。
初期状態 u = 入力場, v = 0
Step 1 uを読んでvへ書く
Step 2 vを読んでuへ書く
Step 3 uを読んでvへ書く
...
最終値 奇数ステップ後はv、偶数ステップ後はu
最終値が入っている方だけ!$acc update self(...)でCPUへ戻します。毎ステップ全面コピーする処理はありません。
また、この実装では基本的に時間ステップごとにOpenACC計算領域を呼ぶ設計です。1回の時間ステップに対するGPUカーネルの実際の起動数は、後述する実行時ログで確認します。時間ステップをまたいで複数回の差分更新を1つのカーネルにまとめる「時間方向の融合」は、今回は行いません。近傍セルの依存関係を考慮しない融合は、正しい計算にならないためです。
6. FP32とFP64の違い
Fortranでは次の2つを使います。
use, intrinsic :: iso_fortran_env, only: real32, real64
real(real32) :: value_fp32
real(real64) :: value_fp64
| 項目 | FP32（単精度） | FP64（倍精度） |
|---|---|---|
| 1要素あたりの記憶容量 | 4 bytes | 8 bytes |
| 有効桁数の目安 | 約7桁 | 約15～16桁 |
| 2配列のメモリ量 |
2 × nx × ny × 4 bytes |
2 × nx × ny × 8 bytes |
| 想定される用途 | 一定の誤差が許容される処理 | 厳密な数値検証や精度を優先する処理 |
FP32のデータ量はFP64のおおむね半分です。メモリ帯域が支配的な計算では有利になる場合があります。ただしFP32が常に2倍高速になるわけではありません。GPUの演算器、キャッシュ、カーネル起動の負担、格子の大きさによって結果は変わります。
特に民生用NVIDIA GPUでは、FP64の演算能力がFP32と比べて制約される製品があります。それでも、5点差分のようにメモリアクセスの比重が大きい処理では、単純に理論演算性能だけで速度差を説明できません。
6.1 同じソースを精度別にビルドする
src/diffusion_common.F90の精度指定は次のようにしています。
#ifdef SINGLE_PRECISION
integer, parameter, public :: rk = real32
character(len=*), parameter, public :: precision_name = 'FP32'
#else
integer, parameter, public :: rk = real64
character(len=*), parameter, public :: precision_name = 'FP64'
#endif
real(rk), parameter, public :: alpha = 0.1_rk
FP32をビルドするときは-DSINGLE_PRECISIONを渡します。FP64をビルドするときは指定しません。このため、alphaと配列・演算の精度をまとめて変更できます。
重要なのは、コンパイルフラグだけで暗黙に型を変えないことです。-r8のように既定型を一括変更すると、別のモジュールや外部ライブラリとの型の整合性を崩すことがあります。今回はreal(rk)を明示しています。
7. GPUの転送量とカーネル内部のアクセス量
7.1 GPUに置く配列の大きさ
単一の配列u(nx,ny)が必要とするバイト数は、次のとおりです。
S=n_x n_y p
ここでp=4（FP32）またはp=8（FP64）です。GPU上にはuとvを保持するので、最低限の配列本体の容量は2Sです。
今回の3方式は、初めにcopyin(u,v)で2配列を転送し、最後に最終配列だけをCPUへ戻します。そのため、明示的な配列転送の概算は2S+S=3Sになります。例えば2048×2048格子では、次のようになります。
| 項目 | FP32 | FP64 |
|---|---|---|
| 配列1つの大きさ | 16 MiB | 32 MiB |
| GPU上の2配列 | 32 MiB | 64 MiB |
| 最初の転送（u・v） | 32 MiB | 64 MiB |
| 最後の転送（結果のみ） | 16 MiB | 32 MiB |
| 明示的な転送合計 | 48 MiB | 96 MiB |
これはデータ本体だけの概算です。配列の管理情報、CUDAの初期化、転送速度、PCIe上での実効効率などは別です。
7.2 各セルが参照する値は6要素相当
5点差分では、中心と上下左右の5点を参照し、新しい値を1点へ書きます。単純に考えると1セルあたり5回の読み込みと1回の書き込みで、FP32なら24 bytes、FP64なら48 bytesのデータアクセスに相当します。
ただし隣接するセルどうしで入力値を再利用でき、GPUのキャッシュに載ることもあります。実際のGPU DRAM転送量が必ず24/48 bytesになるわけではありません。 キャッシュヒット率やメモリ帯域の観察には、Nsight Computeなどのプロファイラを使います。
8. 数値結果はどこまで一致させるか
最初に精度を確認し、その後で速度を比較するという順番は、Part 2から変えません。
同精度のCPUとGPUを照合するときでも、コンパイラによる浮動小数点演算の再配置、FMA（積和演算）の使用、丸め方の違いがあり得ます。そのため、「全セルでビット単位の一致」を必須にはせず、許容差を明示します。
|u_{\mathrm{ref}}-u_{\mathrm{test}}|
\leq \mathrm{atol}+\mathrm{rtol}\,|u_{\mathrm{ref}}|
atolは絶対許容差、rtolは相対許容差です。参照値が0に近い場所では、相対誤差だけで判定すると過剰に厳しくなります。
| 比較 | 配布スクリプトの既定値 |
|---|---|
| 同じFP64でCPUとOpenACCを照合 |
atol=1e-12, rtol=1e-10
|
| 同じFP32でCPUとOpenACCを照合 |
atol=1e-6, rtol=1e-5
|
| FP64 CPUを参照にしてFP32 CPUを照合 |
atol=1e-6, rtol=1e-4
|
これらは、今回の小さな拡散計算用の初期検証値です。水位・水深・流量の計算へそのまま転用してよい基準ではありません。実務モデルでは、許容水位誤差、体積保存、最大浸水深、ピーク時刻など、目的に合わせた条件が必要です。
8.1 最大誤差だけではなく誤差分布も見る
compare_fields.pyでは、全格子について次を計算します。
- 最大絶対誤差：局所的に最も差が大きい場所を確認するため
- RMSE：各セルの誤差を二乗平均でまとめるため
- 相対L2誤差：全体の誤差の大きさを参照値に対して確認するため
- 質量差：領域内の値の合計がどれほど違うか確認するため
- 許容差に違反したセル数：局所的な不一致を見逃さないため
図4 FP32・FP64の数値検証
8.2 この環境で確認できた値
31×27格子を32ステップ進めたCPU基準版で、FP64とFP32を比較しました。
| 指標 | CPUホストでの検証結果 |
|---|---|
| 最大絶対誤差 | 3.949045×10⁻⁹ |
| 相対L2誤差 | 1.138846×10⁻⁷ |
| FP32−FP64の質量差 | 1.036993×10⁻⁷ |
| 検証したセル数 | 837 |
| 許容差違反セル数 | 0 |
このケースでは指定した許容差を満たしました。ただし、この結果をもって、FP32が長時間の水理・流体計算に十分だとは判断できません。 計算の種類、積分時間、境界条件、最小水深付近の扱いなどが変われば、誤差の影響も変わります。
9. 開発環境：Ubuntu 24.04 / WSL2 / Docker Compose
Part 2と同じ系統の環境です。GPU向けコンパイラにはNVIDIA HPC SDK 26.9のnvfortranを使います。2026年11月1日時点では、26.9が公式リリース一覧の現行版です。
NVIDIA公式のHPC SDK 26.9開発用コンテナは大容量です。以前のイメージが手元にある場合は再利用できます。十分な空き容量を確認してください。
必要なものは次のとおりです。
- Windows 11 + WSL2 Ubuntu 24.04、またはUbuntu 24.04
- NVIDIA GPUと対応するNVIDIAドライバ
- Docker Engine・Compose、NVIDIA Container Toolkit
-
nvfortranが入ったNVIDIA HPC SDKコンテナ - 比較用の
gfortranとPython 3
前回と同様に、GPUドライバはホスト側で管理します。 WSL2内部へNVIDIAのLinux用ディスプレイドライバを別途インストールする方法は採りません。
9.1 ZIPの構成
fortran2018_gpu_part3_r003/
├── Dockerfile
├── compose.yaml
├── README.md
├── TEST_RESULTS.md
├── src/
│ ├── diffusion_common.F90
│ ├── cpu_reference.f90
│ └── acc_stencil.f90
├── scripts/
│ ├── build_cpu.sh
│ ├── build_gpu.sh
│ ├── check_cpu.sh
│ ├── check_gpu.sh
│ ├── compare_fields.py
│ └── benchmark.py
├── tests/
│ └── test_compare_fields.py
├── images/
│ └── README.md
└── results/
└── README.md
resultsフォルダーには、実行後に格子の出力ファイル、比較ログ、測定CSVが作成されます。GPUの数値を記入したサンプルCSVは同梱していません。
9.2 コンパイラとGPU認識の確認
# compose.yamlのある場所で実行します。
docker compose build nvhpc
docker compose run --rm nvhpc nvfortran --version
docker compose run --rm nvhpc nvidia-smi
nvfortran --versionはコンパイラが利用できるかの確認であり、nvidia-smiはコンテナからGPUが見えるかの確認です。両方が成功しても、実際のGPUカーネルが正しい数値を返したことにはなりません。 この後にコンパイル・実行・全格子照合を行います。
9.3 GPU向けのコンパイルフラグ
nvfortran -O3 -Mpreprocess -acc=gpu \
-gpu=cc75,cc89,mem:separate -Minfo=accel \
-module build/gpu/mod_fp64 \
src/diffusion_common.F90 src/acc_stencil.f90 \
-o build/gpu/fp64_acc
このコマンドは、RTX 2080系向けのcc75とRTX 4060 Ti向けのcc89を例にしています。-gpu=mem:separateはGPU側のデータ転送を明示して測るための指定、-Minfo=accelはカーネル化したループのコンパイラ診断です。
実際のビルドスクリプトではFP32/FP64を両方コンパイルし、診断結果をbuild/gpu/へ保存します。
10. 計算の正しさを確認する
最初に、ローカルのGNU Fortranでホスト実行を確認します。
bash scripts/check_cpu.sh
gfortran -fopenaccは、GPUのない環境ではOpenACCコードのホスト代替実行になります。これはループ構造が正しく動くかを調べるための試験であり、NVIDIA GPUでの実行確認ではありません。
GPUが利用できる場合は、次を実行します。
docker compose run --rm nvhpc bash scripts/check_gpu.sh
check_gpu.shでは、nvfortranが作ったOpenACC版を本物のGPUで実行し、gfortranのCPU版と全格子を照合します。合格したときだけresults/check_gpu/PASS.markerが作られます。失敗時には作りません。
例えば手動で1ケースだけ照合するには、次のようにします。
# GPU版は「格子数 格子数 ステップ数 比較方式 保存先」の順です。
build/gpu/fp64_acc 31 27 32 contiguous results/fp64_gpu.dat
# CPU版は「格子数 格子数 ステップ数 保存先」の順です。
build/host/fp64_cpu 31 27 32 results/fp64_cpu.dat
python3 scripts/compare_fields.py \
results/fp64_cpu.dat results/fp64_gpu.dat \
--atol 1e-12 --rtol 1e-10
注意：上のコマンドは、ビルドを済ませ、results/が存在する状態で実行します。
11. GPUカーネルの動きを診断する
11.1 コンパイル診断：-Minfo=accel
build/gpu/fp32_compile.logとbuild/gpu/fp64_compile.logには、カーネル生成やループ配置に関する情報が出力されます。
cat build/gpu/fp64_compile.log
ここで、gangやvectorがどのループへ割り当てられたかを確認します。「contiguousを使ったので絶対に連続アクセスになった」と決めつけないための確認です。
11.2 実行時：カーネル起動数と転送を調べる
# カーネル起動の通知（1）を使う。確認用の短い計算です。
NVCOMPILER_ACC_NOTIFY=1 \
build/gpu/fp64_acc 128 128 7 contiguous
# 起動（1）＋データ転送（2）を見たい場合。
NVCOMPILER_ACC_NOTIFY=3 \
build/gpu/fp64_acc 128 128 7 baseline
通知の出力を見れば、実際にGPUカーネルが起動したかや、予想外の転送が入っていないかを確認できます。通知を有効にした状態は計時条件が変わるため、通常のベンチマークとは分けて実施します。
11.3 実行時：OpenACCの区間計測を調べる
NVCOMPILER_ACC_TIME=1 \
build/gpu/fp64_acc 1024 1024 100 contiguous
この情報とプログラムのLOOP_Sを混同しないでください。LOOP_SはCPU側の壁時計で測った時間で、カーネル起動や同期も含みます。GPUのカーネルだけにかかった時間を見る場合には、OpenACCランタイムの情報やNsight Systems / Nsight Computeを組み合わせます。
12. 性能比較は3つの時間を分ける
今回のプログラムは、少なくとも次の区間を出力します。
| 項目 | 意味 | 注意点 |
|---|---|---|
TOTAL_S |
ソルバー内部の対象区間の壁時計 | GPU版は入力転送・ループ・結果回収・解放を含む |
LOOP_S |
GPU上の時間ループ全体の壁時計 | GPU専用のカーネル時間ではない |
H2D_S |
GPUデータ領域の開始時刻まで | 初回のGPU初期化を含む可能性がある |
D2H_S |
最終結果を回収する区間 | 片方の配列だけを戻す |
EXIT_S |
GPUデータ領域を抜ける区間 | 解放処理などを含み得る |
process_s |
Pythonからのプロセス起動〜終了 | CUDA初期化や起動の負担も見える |
CPU版は計算ループの時間をTOTAL_Sとして表示し、GPU版はGPUデータ領域の開始から終了までを計測します。したがって、両方のTOTAL_Sは完全に同じ内訳ではありません。 ただしGPU側で転送まで含む「実務上のコスト」を見るには意味のある値です。CPU側が逐次版である点も含めて、性能比較の前提を書き残します。
12.1 ウォームアップと中央値
短い計算はばらつきが大きくなるため、測定前に1回実行し、以降を複数回測ります。中央値を代表値として用い、試行ごとの値は必ず保存します。
# GPU実機でcheck_gpu.shが成功した後に実行します。
docker compose run --rm nvhpc \
python3 scripts/benchmark.py \
--target gpu --nx 1024 --ny 1024 --steps 400 \
--warmups 1 --repeats 5
ベンチマーク用Pythonはそれぞれの実行ファイルを毎回新しいプロセスとして起動します。ウォームアップの別プロセスを使っても、同じCUDAコンテキストが使い回されるわけではありません。今回のTOTAL_Sはこうした初期化の影響を含み得る計測として扱います。純粋なカーネル時間だけを厳密に測りたい場合は、同一プロセスで繰り返し実行しGPUイベント計測を行う構成に拡張する必要があります。
12.2 出力されるファイル
results/
bench_part3_gpu_1024x1024_400steps.csv
bench_part3_gpu_1024x1024_400steps_median.csv
前者が各試行の生データ、後者が方法別・精度別の中央値です。CPUホスト代替実行の場合はファイル名のgpu部分がhostになります。
12.3 記事へ追記する測定表
GPU上で検証・測定してから埋めるため、今は空欄にしています。
| GPU実機・1024×1024・400 step | FP32 TOTAL_S中央値 |
FP64 TOTAL_S中央値 |
FP32 LOOP_S中央値 |
FP64 LOOP_S中央値 |
|---|---|---|---|---|
GPU baseline
|
― | ― | ― | ― |
GPU contiguous
|
― | ― | ― | ― |
GPU strided
|
― | ― | ― | ― |
| CPU基準（逐次） | ― | ― | ― | ― |
baselineに対する高速化率は、例えば次で計算します。
\mathrm{Speedup}=\frac{T_{\mathrm{baseline}}}{T_{\mathrm{candidate}}}
またFP64に対してFP32が何倍速いかを比較する場合は、同じ方式・同じ格子・同じステップ数のTOTAL_SまたはLOOP_Sを使います。別GPUや別コンパイラの結果を同じ表で混ぜないようにします。
12.4 小さい格子だけで速度を判断しない
小さい格子ではカーネル起動の負担が目立ち、格子が大きくなるとGPUのメモリ帯域や演算能力の影響が増えることがあります。128×128・512×512・1024×1024・2048×2048など、複数の格子サイズを測ると変化が分かります。
2048×2048以上では、GPUメモリだけでなく、全格子をテキスト出力する際のディスク容量にも注意が必要です。ベンチマーク本体はファイル出力を計測区間から外し、必要な検証時だけ出力するようにしています。
13. 洪水・河川解析に応用するときの注意
今回の格子計算は、単純な拡散方程式です。これだけでは、洪水氾濫解析の数値精度を保証できません。それでも、GPU化の順番を考える材料になります。
水深や流量の計算で先に確認することは、次のような項目です。
- CPUとGPUで水深・流量・水位の時系列がどの程度一致するか
- 全格子・全時刻の誤差だけでなく、ピーク水位と発生時刻の誤差はどの程度か
- 流入・流出・貯留量の水収支残差が十分小さいか
- 濡れ・乾き境界、微小水深、ゼロ割回避で精度を落としていないか
- FP32へ変更すると、CFL条件やタイムステップ制御に悪影響が出ないか
さらに1D–2D連成では、格子ごとの更新とは別に、河道・水路との交換流量や境界処理があります。GPU化する前後で同じ収支になるよう、更新手順と同期の位置を決める必要があります。
FP32を試す場合でも、全量をFP32へ一度に置き換えず、精度の必要な集計・保存量をFP64に残すといった混合精度の検討ができます。ただし、水量保存が改善するかどうかはアルゴリズム次第であり、精度を混ぜれば必ずよくなるわけではありません。
14. 今回のまとめ
今回は、Part 2で作ったGPU常駐方式を土台に、3種類のGPUカーネルを比較できる構成へ進めました。
- Fortranの
u(i,j)はi方向が連続することを確認しました。 -
collapse(2)、jを外側にしたgang/vector、iを外側にしたgang/vectorを用意しました。 - FP32とFP64を同じFortranソースからビルドできるようにしました。
- CPU基準版・OpenACCホスト版では全格子照合に成功しました。
- GPU実行時には、数値確認後に転送込み時間・ループ時間・プロセス時間を分けて計測します。
CPU実行で3方式が一致したことは確認できましたが、それはGPUでの高速化が完了したという意味ではありません。次はGPU実機でのコンパイルと数値照合、その後のプロファイリングです。
次のPart 4では、計測でカーネル起動の負担が大きかった場合の対策として、OpenACCのasyncとwait、GPUカーネル実行時間の可視化、時間ループの起動負担を調べます。精度と水収支を先に固定してから改善していきます。
ソースコードの配布
ソース一式はfortran2018_gpu_part3_sources_r003.zipにまとめています。Docker Compose、Fortran 2018、Python標準ライブラリによる全格子比較、ベンチマーク、回帰テストが含まれています。
参考のために、ZIPファイルを置いておきますので、ダウンロードして活用してください。現時点では実機GPUの測定結果やGPUの合格ログは配布していません。
参考資料
- NVIDIA HPC SDK 26.9のリリース一覧
- NVIDIA HPC Compilers User's Guide 26.9
- NVIDIA OpenACC Getting Started Guide 26.9
- NVIDIA HPC SDKコンテナ（NGC）
付録A：配布ソースの全文
実際のコードを確認しやすいように、主要ソースを折りたたみ形式で掲載します。Qiita上で長く表示される場合は閉じたまま読み進められます。記事の説明とソースが食い違わないよう、配布ファイルから内容をそのまま取り込んでいます。
src/diffusion_common.F90 ― 共通の精度・初期値・検証・入出力
!> Part 3で使用する精度切替・初期値・数値検査・入出力モジュールです。
!!
!! SINGLE_PRECISIONを定義した場合だけreal32を使います。
!! 未定義時はreal64を使います。CPU版とOpenACC版のどちらも
!! 必ず同じモジュールを使うため、初期値や判定基準を統一できます。
module diffusion_common
use, intrinsic :: iso_fortran_env, only : real32, real64, int64, error_unit
use, intrinsic :: ieee_arithmetic, only : ieee_is_finite
implicit none
private
#ifdef SINGLE_PRECISION
integer, parameter, public :: rk = real32
character(len=*), parameter, public :: precision_name = 'FP32'
#else
integer, parameter, public :: rk = real64
character(len=*), parameter, public :: precision_name = 'FP64'
#endif
real(rk), parameter, public :: alpha = 0.1_rk
public :: parse_grid, get_path, set_initial_condition, seconds_between
public :: print_summary, save_field, validate_counts
contains
!> 共通の3つの整数引数(nx, ny, steps)を読みます。
!! 引数を省略した場合は256x256格子・100ステップです。
subroutine parse_grid(nx, ny, steps)
integer, intent(out) :: nx, ny, steps
nx = 256
ny = 256
steps = 100
if (command_argument_count() >= 1) call parse_posint(1, nx)
if (command_argument_count() >= 2) call parse_posint(2, ny)
if (command_argument_count() >= 3) call parse_posint(3, steps)
if (nx < 5 .or. ny < 5 .or. steps < 1) &
error stop 'nx/ny must be >=5; steps must be >=1.'
end subroutine parse_grid
!> 不正な引数で静かに進まないよう、整数と範囲を検査します。
subroutine parse_posint(index, number)
integer, intent(in) :: index
integer, intent(out) :: number
integer :: ios
character(len=80) :: value
call get_command_argument(index, value, status=ios)
if (ios /= 0) error stop 'Failed to read integer argument.'
read(value, *, iostat=ios) number
if (ios /= 0 .or. number < 1) error stop 'Expected a positive integer.'
end subroutine parse_posint
!> 指定した引数位置の保存先を読みます。省略ならファイル出力しません。
subroutine get_path(index, path)
integer, intent(in) :: index
character(len=*), intent(out) :: path
integer :: ios
path = ''
if (command_argument_count() < index) return
call get_command_argument(index, path, status=ios)
if (ios /= 0 .or. len_trim(path) == 0) error stop 'Invalid output path.'
end subroutine get_path
!> コマンド引数の数が期待する範囲かを確認します。
subroutine validate_counts(min_args, max_args)
integer, intent(in) :: min_args, max_args
integer :: nargs
nargs = command_argument_count()
if (nargs > max_args .or. nargs < min_args) then
write(error_unit,'(a,i0,a,i0,a)') &
'Expected ', min_args, ' to ', max_args, ' arguments.'
error stop 2
end if
end subroutine validate_counts
!> 初期場は中央1セル=1、ほかは0。境界は全期間0で固定します。
subroutine set_initial_condition(u, v)
real(rk), intent(out) :: u(:,:), v(:,:)
integer :: nx, ny
if (any(shape(u) /= shape(v))) error stop 'Array shape mismatch.'
nx = size(u, 1)
ny = size(u, 2)
u = 0.0_rk
v = 0.0_rk
u((nx+1)/2, (ny+1)/2) = 1.0_rk
end subroutine set_initial_condition
!> SYSTEM_CLOCK(int64)の差から壁時計秒を計算します。
pure function seconds_between(t0, t1, rate) result(value)
integer(int64), intent(in) :: t0, t1, rate
real(real64) :: value
if (rate <= 0_int64 .or. t1 < t0) then
value = -1.0_real64
else
value = real(t1-t0, real64)/real(rate, real64)
end if
end function seconds_between
!> 実行結果を検査し、Pythonで読める形式で出力します。
!!
!! 境界から質量が流出し得るためMASS=1の一致判定はしません。
!! GPU速度とは独立した安定性検査であり、CPUとの一致は
!! compare_fields.pyで全格子を確認します。
subroutine print_summary(method, field, total_s, loop_s, h2d_s, d2h_s, exit_s, steps)
character(len=*), intent(in) :: method
real(rk), intent(in) :: field(:,:)
real(real64), intent(in) :: total_s, loop_s, h2d_s, d2h_s, exit_s
integer, intent(in) :: steps
real(real64) :: mass, boundary_tolerance
real(rk) :: min_u, max_u
integer :: nx, ny
nx = size(field,1)
ny = size(field,2)
if (.not. all(ieee_is_finite(field))) error stop 'Nonfinite solution.'
if (total_s < 0 .or. loop_s < 0 .or. h2d_s < 0 .or. d2h_s < 0 .or. exit_s < 0) &
error stop 'Invalid timing.'
mass = sum(real(field, real64))
min_u = minval(field)
max_u = maxval(field)
#ifdef SINGLE_PRECISION
boundary_tolerance = 1.0e-4_real64
#else
boundary_tolerance = 1.0e-10_real64
#endif
if (real(min_u, real64) < -boundary_tolerance .or. &
real(max_u, real64) > 1.0_real64 + boundary_tolerance) &
error stop 'Field outside physical bounds.'
if (mass < -boundary_tolerance .or. mass > 1.0_real64+boundary_tolerance) &
error stop 'Unstable mass.'
print '(a)', 'CHECK=PASS'
print '(a)', 'METHOD='//trim(method)
print '(a)', 'PRECISION='//precision_name
print '(a,i0)', 'NX=', nx
print '(a,i0)', 'NY=', ny
print '(a,i0)', 'STEPS=', steps
print '(a,es24.16)', 'MASS=', mass
print '(a,es24.16)', 'CENTRE=', real(field((nx+1)/2,(ny+1)/2), real64)
print '(a,es24.16)', 'MIN=', real(min_u,real64)
print '(a,es24.16)', 'MAX=', real(max_u,real64)
print '(a,es24.16)', 'TOTAL_S=', total_s
print '(a,es24.16)', 'LOOP_S=', loop_s
print '(a,es24.16)', 'H2D_S=', h2d_s
print '(a,es24.16)', 'D2H_S=', d2h_s
print '(a,es24.16)', 'EXIT_S=', exit_s
end subroutine print_summary
!> Fortran配列順(iが先に増加)の全格子値をテキスト保存します。
!! 最初の行はnx ny、その後は1セル1行です。計測終了後のみ呼びます。
subroutine save_field(path, field)
character(len=*), intent(in) :: path
real(rk), intent(in) :: field(:,:)
integer :: unit, ios, i, j
if (len_trim(path) == 0) return
open(newunit=unit,file=trim(path),status='replace',action='write',iostat=ios)
if (ios /= 0) error stop 'Failed to open field file.'
write(unit,'(i0,1x,i0)') size(field,1), size(field,2)
do j = 1, size(field,2)
do i = 1, size(field,1)
write(unit,'(es26.17e3)') real(field(i,j),real64)
end do
end do
close(unit)
print '(a)', 'FIELD_FILE='//trim(path)
end subroutine save_field
end module diffusion_common
src/cpu_reference.f90 ― CPU基準版
!> Part 2と同じ5点差分をCPU上で逐次計算する基準実装です。
!! CPU最速版ではありません。GPU版の全格子結果との比較が目的です。
program cpu_reference
use, intrinsic :: iso_fortran_env, only : int64, real64
use diffusion_common, only : rk, alpha, parse_grid, get_path, &
set_initial_condition, seconds_between, print_summary, save_field, validate_counts
implicit none
real(rk), allocatable :: u(:,:), v(:,:)
integer :: nx, ny, steps, step, i, j
integer(int64) :: t0, t1, rate
real(real64) :: total_s
character(len=512) :: field_path
call validate_counts(0, 4)
call parse_grid(nx, ny, steps)
call get_path(4, field_path)
allocate(u(nx,ny),v(nx,ny))
call set_initial_condition(u,v)
call system_clock(t0,rate)
do step = 1, steps
if (mod(step,2) == 1) then
do j = 1, ny
do i = 1, nx
if (i==1 .or. i==nx .or. j==1 .or. j==ny) then
v(i,j) = 0.0_rk
else
v(i,j) = u(i,j) + alpha*(u(i-1,j)+u(i+1,j)+ &
u(i,j-1)+u(i,j+1)-4.0_rk*u(i,j))
end if
end do
end do
else
do j = 1, ny
do i = 1, nx
if (i==1 .or. i==nx .or. j==1 .or. j==ny) then
u(i,j) = 0.0_rk
else
u(i,j) = v(i,j) + alpha*(v(i-1,j)+v(i+1,j)+ &
v(i,j-1)+v(i,j+1)-4.0_rk*v(i,j))
end if
end do
end do
end if
end do
call system_clock(t1)
total_s = seconds_between(t0,t1,rate)
if (mod(steps,2) == 1) then
call print_summary('CPU',v,total_s,total_s,0.0_real64,0.0_real64,0.0_real64,steps)
call save_field(field_path,v)
else
call print_summary('CPU',u,total_s,total_s,0.0_real64,0.0_real64,0.0_real64,steps)
call save_field(field_path,u)
end if
end program cpu_reference
src/acc_stencil.f90 ― OpenACCの3方式
!> GPU常駐の5点差分。配列と初期値はPart 2と同じです。
!!
!! 実験する3つのループ構成:
!! - baseline : 全セルをcollapse(2)し境界を条件分岐で更新
!! - contiguous : j外側(gang), i内側(vector)、内部セルだけ更新
!! - strided : i外側(gang), j内側(vector)、内部セルだけ更新
!!
!! Fortranのu(i,j)ではi方向がメモリ上で連続しています。
!! ただし、実際のGPUのスレッド割付・メモリ処理はコンパイラ次第です。
!! stridedが必ず遅い、contiguousが必ず速いとは決めつけません。
!!
!! CPUとGPUを分離する-gpu=mem:separateを用い、
!! コピーとループを別の区間として計測します。
program acc_stencil
use, intrinsic :: iso_fortran_env, only : int64, real64
use diffusion_common, only : rk, alpha, parse_grid, get_path, &
set_initial_condition, seconds_between, print_summary, save_field, validate_counts
implicit none
real(rk), allocatable :: u(:,:), v(:,:)
integer :: nx, ny, steps, step, i, j, ios
integer(int64) :: t0, t_enter, t_loop, t_back, t_end, rate
real(real64) :: total_s, loop_s, h2d_s, d2h_s, exit_s
character(len=512) :: field_path
character(len=32) :: variant
call validate_counts(4, 5)
call parse_grid(nx,ny,steps)
call get_command_argument(4,variant,status=ios)
if (ios /= 0) error stop 'Invalid variant.'
select case(trim(variant))
case ('baseline','contiguous','strided')
! 有効な比較方法です。
case default
error stop 'Variant must be baseline, contiguous or strided.'
end select
call get_path(5,field_path)
allocate(u(nx,ny),v(nx,ny))
call set_initial_condition(u,v)
call system_clock(t0,rate)
! 両配列の境界ゼロを保証するため、ここではu/vの両方を転送します。
! 3方式の初期転送量を揃え、比較条件を固定します。
!$acc data copyin(u,v)
call system_clock(t_enter)
do step=1,steps
if (mod(step,2) == 1) then
! ---- 奇数ステップ: u -> v ----
select case (trim(variant))
case ('baseline')
! 全セルを対象とし、境界もゼロを書き込みます。
!$acc parallel loop collapse(2) present(u,v)
do j=1,ny
do i=1,nx
if (i==1 .or. i==nx .or. j==1 .or. j==ny) then
v(i,j)=0.0_rk
else
v(i,j)=u(i,j)+alpha*(u(i-1,j)+u(i+1,j)+ &
u(i,j-1)+u(i,j+1)-4.0_rk*u(i,j))
end if
end do
end do
!$acc end parallel loop
case ('contiguous')
! i方向の隣り合う要素をvectorに割り当てる候補。
! 境界は初期化時からゼロで、内部の更新で変更しません。
!$acc parallel loop gang present(u,v)
do j=2,ny-1
!$acc loop vector
do i=2,nx-1
v(i,j)=u(i,j)+alpha*(u(i-1,j)+u(i+1,j)+ &
u(i,j-1)+u(i,j+1)-4.0_rk*u(i,j))
end do
end do
!$acc end parallel loop
case ('strided')
! 比較用にjを内側へ移動。iが連続しないため
! 実GPUでのメモリアクセスの違いを調べます。
!$acc parallel loop gang present(u,v)
do i=2,nx-1
!$acc loop vector
do j=2,ny-1
v(i,j)=u(i,j)+alpha*(u(i-1,j)+u(i+1,j)+ &
u(i,j-1)+u(i,j+1)-4.0_rk*u(i,j))
end do
end do
!$acc end parallel loop
end select
else
! ---- 偶数ステップ: v -> u ----
select case (trim(variant))
case ('baseline')
!$acc parallel loop collapse(2) present(u,v)
do j=1,ny
do i=1,nx
if (i==1 .or. i==nx .or. j==1 .or. j==ny) then
u(i,j)=0.0_rk
else
u(i,j)=v(i,j)+alpha*(v(i-1,j)+v(i+1,j)+ &
v(i,j-1)+v(i,j+1)-4.0_rk*v(i,j))
end if
end do
end do
!$acc end parallel loop
case ('contiguous')
!$acc parallel loop gang present(u,v)
do j=2,ny-1
!$acc loop vector
do i=2,nx-1
u(i,j)=v(i,j)+alpha*(v(i-1,j)+v(i+1,j)+ &
v(i,j-1)+v(i,j+1)-4.0_rk*v(i,j))
end do
end do
!$acc end parallel loop
case ('strided')
!$acc parallel loop gang present(u,v)
do i=2,nx-1
!$acc loop vector
do j=2,ny-1
u(i,j)=v(i,j)+alpha*(v(i-1,j)+v(i+1,j)+ &
v(i,j-1)+v(i,j+1)-4.0_rk*v(i,j))
end do
end do
!$acc end parallel loop
end select
end if
end do
! ループの非同期設定を追加しても、計時終了前に完了を保証します。
!$acc wait
call system_clock(t_loop)
! 使われていない中間配列はCPUへ戻しません。
if (mod(steps,2) == 1) then
!$acc update self(v)
else
!$acc update self(u)
end if
!$acc wait
call system_clock(t_back)
!$acc end data
call system_clock(t_end)
total_s = seconds_between(t0,t_end,rate)
h2d_s = seconds_between(t0,t_enter,rate)
loop_s = seconds_between(t_enter,t_loop,rate)
d2h_s = seconds_between(t_loop,t_back,rate)
exit_s = seconds_between(t_back,t_end,rate)
if (mod(steps,2) == 1) then
call print_summary('ACC_'//trim(variant),v,total_s,loop_s,h2d_s,d2h_s,exit_s,steps)
call save_field(field_path,v)
else
call print_summary('ACC_'//trim(variant),u,total_s,loop_s,h2d_s,d2h_s,exit_s,steps)
call save_field(field_path,u)
end if
end program acc_stencil
scripts/build_cpu.sh ― CPUビルド
#!/usr/bin/env bash
# gfortranによる検証用ホストビルド。-fopenaccはCPU代替実行です。
set -euo pipefail
cd "$(dirname "$0")/.."
FC="${FC:-gfortran}"
FCFLAGS="${FCFLAGS:--std=f2018 -O3 -Wall -Wextra}"
mkdir -p build/host/mod_fp32 build/host/mod_fp64
for precision in fp32 fp64; do
define=()
if [[ "$precision" == fp32 ]]; then define=(-DSINGLE_PRECISION); fi
moddir="build/host/mod_${precision}"
# CPU版とOpenACC版で.modファイルを混用しないよう、同一精度で揃えます。
"$FC" -cpp "${define[@]}" $FCFLAGS -J"$moddir" -I"$moddir" \
src/diffusion_common.F90 src/cpu_reference.f90 -o "build/host/${precision}_cpu"
"$FC" -cpp "${define[@]}" $FCFLAGS -fopenacc -J"$moddir" -I"$moddir" \
src/diffusion_common.F90 src/acc_stencil.f90 -o "build/host/${precision}_acc"
done
printf 'HOST_BUILD=PASS (gfortran; GPU performance NOT measured)\n'
scripts/build_gpu.sh ― GPUビルド
#!/usr/bin/env bash
# NVIDIA HPC SDK 26.9のnvfortranでOpenACC版だけをGPU向けビルド。
# NVIDIA RTX 2080 (cc75)とRTX 4060 Ti (cc89)を両方ターゲットに含めます。
# GPU_TARGETSを指定すれば環境に合わせて変更できます。
set -euo pipefail
cd "$(dirname "$0")/.."
command -v nvfortran >/dev/null || { echo 'nvfortran not found' >&2; exit 1; }
GPU_TARGETS="${GPU_TARGETS:-cc75,cc89}"
mkdir -p build/gpu/mod_fp32 build/gpu/mod_fp64
for precision in fp32 fp64; do
define=()
if [[ "$precision" == fp32 ]]; then define=(-DSINGLE_PRECISION); fi
moddir="build/gpu/mod_${precision}"
# モジュールは精度ごとに分け、OpenACCのGPUオフロードを有効化します。
nvfortran -O3 -Mpreprocess "${define[@]}" -acc=gpu \
-gpu="${GPU_TARGETS},mem:separate" -Minfo=accel \
-module "$moddir" \
src/diffusion_common.F90 src/acc_stencil.f90 \
-o "build/gpu/${precision}_acc" \
2>"build/gpu/${precision}_compile.log"
printf 'Compiled %s; diagnostic: %s\n' "$precision" "build/gpu/${precision}_compile.log"
done
printf 'GPU_BUILD=PASS (only compilation; GPU execution requires check_gpu.sh)\n'
scripts/check_cpu.sh ― CPUホストでの検証
#!/usr/bin/env bash
# 精度・奇偶・境界を含めたCPU側の全格子検証です。
# gfortran -fopenaccはホスト代替実行でありGPU動作確認ではありません。
set -euo pipefail
cd "$(dirname "$0")/.."
rm -f results/check_cpu/PASS.marker
bash scripts/build_cpu.sh
mkdir -p results/check_cpu
for precision in fp32 fp64; do
for steps in 7 8 31 32; do
tag="${precision}_31x27_${steps}"
"build/host/${precision}_cpu" 31 27 "$steps" "results/check_cpu/${tag}_cpu.dat" \
>"results/check_cpu/${tag}_cpu.log"
for method in baseline contiguous strided; do
"build/host/${precision}_acc" 31 27 "$steps" "$method" \
"results/check_cpu/${tag}_${method}.dat" \
>"results/check_cpu/${tag}_${method}.log"
done
# 全格子をCPU基準版と照合します。精度別に許容値を変えます。
if [[ "$precision" == fp32 ]]; then
atol=1e-6; rtol=1e-5
else
atol=1e-12; rtol=1e-10
fi
python3 scripts/compare_fields.py "results/check_cpu/${tag}_cpu.dat" \
"results/check_cpu/${tag}_baseline.dat" \
"results/check_cpu/${tag}_contiguous.dat" \
"results/check_cpu/${tag}_strided.dat" --atol "$atol" --rtol "$rtol"
done
done
# FP32とFP64の差は異精度比較なので大きめの許容差を指定します。
python3 scripts/compare_fields.py results/check_cpu/fp64_31x27_32_cpu.dat \
results/check_cpu/fp32_31x27_32_cpu.dat --atol 1e-6 --rtol 1e-4 \
--csv results/check_cpu/precision_comparison.csv
python3 -m unittest discover -s tests -v
touch results/check_cpu/PASS.marker
printf 'CPU_VALIDATION=PASS (2 precisions x 4 steps x 3 OpenACC methods)\n'
scripts/check_gpu.sh ― GPU実機での検証
#!/usr/bin/env bash
# gfortranで作ったCPU結果とnvfortran実GPU結果を全セルで比較します。
# 失敗時はPASS.markerを生成しません。
set -euo pipefail
cd "$(dirname "$0")/.."
rm -f results/check_gpu/PASS.marker
command -v nvidia-smi >/dev/null || { echo 'nvidia-smi not available' >&2; exit 1; }
nvidia-smi -L
bash scripts/build_cpu.sh
bash scripts/build_gpu.sh
mkdir -p results/check_gpu
for precision in fp32 fp64; do
for steps in 7 8 31 32; do
tag="${precision}_31x27_${steps}"
"build/host/${precision}_cpu" 31 27 "$steps" "results/check_gpu/${tag}_cpu.dat" \
>"results/check_gpu/${tag}_cpu.log"
for method in baseline contiguous strided; do
"build/gpu/${precision}_acc" 31 27 "$steps" "$method" \
"results/check_gpu/${tag}_${method}.dat" \
>"results/check_gpu/${tag}_${method}.log"
done
if [[ "$precision" == fp32 ]]; then
atol=1e-6; rtol=1e-5
else
atol=1e-12; rtol=1e-10
fi
python3 scripts/compare_fields.py "results/check_gpu/${tag}_cpu.dat" \
"results/check_gpu/${tag}_baseline.dat" \
"results/check_gpu/${tag}_contiguous.dat" \
"results/check_gpu/${tag}_strided.dat" --atol "$atol" --rtol "$rtol"
done
done
python3 scripts/compare_fields.py results/check_gpu/fp64_31x27_32_cpu.dat \
results/check_gpu/fp32_31x27_32_cpu.dat --atol 1e-6 --rtol 1e-4 \
--csv results/check_gpu/precision_comparison.csv
python3 -m unittest discover -s tests -v
touch results/check_gpu/PASS.marker
printf 'GPU_VALIDATION=PASS (real NVIDIA GPU + all grid comparisons)\n'
scripts/compare_fields.py ― 全格子の誤差計算
#!/usr/bin/env python3
"""Compare complete grids written by the Part 3 Fortran programs.
The first line of each file must contain nx and ny. Remaining lines are
column-major cell values, one number per line. Compare all cells using
|reference-candidate| <= atol + rtol*|reference|.
Pass through --csv to retain numeric errors for FP32/FP64 experiments.
A small tolerance is necessary when FP32 and FP64 or different compilers
are compared. A zero tolerance can be used for strict host comparisons.
"""
from __future__ import annotations
import argparse
import csv
import math
from pathlib import Path
def load_field(path: Path) -> tuple[int, int, list[float]]:
"""Read field, rejecting a wrong header, missing elements or NaN/Inf."""
with path.open(encoding="utf-8") as handle:
dims = handle.readline().split()
if len(dims) != 2:
raise ValueError(f"{path}: bad grid header")
try:
nx, ny = (int(item) for item in dims)
except ValueError as exc:
raise ValueError(f"{path}: noninteger header") from exc
if nx < 5 or ny < 5:
raise ValueError(f"{path}: nx and ny must be >=5")
data: list[float] = []
for lineno, row in enumerate(handle, start=2):
if not row.strip():
raise ValueError(f"{path}:{lineno}: blank row")
try:
val = float(row.strip().replace("D", "E"))
except ValueError as exc:
raise ValueError(f"{path}:{lineno}: malformed number") from exc
if not math.isfinite(val):
raise ValueError(f"{path}:{lineno}: nonfinite number")
data.append(val)
if len(data) != nx * ny:
raise ValueError(f"{path}: {nx*ny} values expected, got {len(data)}")
return nx, ny, data
def metrics(
reference: tuple[int, int, list[float]],
candidate: tuple[int, int, list[float]],
*,
atol: float,
rtol: float,
) -> dict[str, float | int]:
"""Calculate max abs error, RMSE, relative L2, mass error and failures."""
if reference[:2] != candidate[:2]:
raise ValueError("Field dimensions do not match")
ref, cand = reference[2], candidate[2]
if len(ref) != len(cand):
raise ValueError("Field lengths do not match")
differences = [abs(a - b) for a, b in zip(ref, cand)]
failures = sum(
delta > atol + rtol * abs(a)
for a, delta in zip(ref, differences)
)
sq_error = math.fsum(delta * delta for delta in differences)
sq_ref = math.fsum(a * a for a in ref)
mass_ref = math.fsum(ref)
mass_cand = math.fsum(cand)
max_difference = max(differences)
worst = differences.index(max_difference)
return {
"max_abs": max_difference,
"rmse": math.sqrt(sq_error / len(ref)),
"relative_l2": math.sqrt(sq_error / sq_ref) if sq_ref else math.sqrt(sq_error),
"mass_difference": mass_cand - mass_ref,
"worst_index_0based": worst,
"failures": failures,
"cells": len(ref),
}
def main() -> None:
"""Perform CLI comparisons and optionally save row-based CSV diagnostics."""
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("reference", type=Path)
parser.add_argument("candidates", type=Path, nargs="+")
parser.add_argument("--atol", type=float, default=1e-12)
parser.add_argument("--rtol", type=float, default=1e-10)
parser.add_argument("--csv", type=Path)
args = parser.parse_args()
if args.atol < 0 or args.rtol < 0 or not all(
math.isfinite(x) for x in (args.atol, args.rtol)
):
parser.error("atol and rtol must be finite and nonnegative")
expected = load_field(args.reference)
records: list[dict[str, str | float | int]] = []
failed = False
for filename in args.candidates:
measures = metrics(expected, load_field(filename), atol=args.atol, rtol=args.rtol)
rec = {"reference": str(args.reference), "candidate": str(filename), **measures}
records.append(rec)
failed = failed or measures["failures"] > 0
print(
f"COMPARE_FILE={filename} MAX_ABS={measures['max_abs']:.6e} "
f"REL_L2={measures['relative_l2']:.6e} "
f"MASS_DIFF={measures['mass_difference']:.6e} "
f"FAILURES={measures['failures']}"
)
if args.csv:
args.csv.parent.mkdir(parents=True, exist_ok=True)
with args.csv.open("w", encoding="utf-8", newline="") as handle:
writer = csv.DictWriter(handle, fieldnames=list(records[0]))
writer.writeheader()
writer.writerows(records)
if failed:
raise SystemExit("COMPARE=FAIL")
print("COMPARE=PASS")
if __name__ == "__main__":
main()
scripts/benchmark.py ― 中央値を出力する計測プログラム
#!/usr/bin/env python3
"""Run reproducible CPU and OpenACC benchmarks after full-grid verification.
For GPU target, CPU binaries are compiled by gfortran and OpenACC binaries
by nvfortran. The CSV records TOTAL_S (allocations + transfers + loop),
LOOP_S (host-wall including launches and synchronization), process wall,
H2D/D2H times, precision, variant, host name, and optional GPU name.
Each executable is a fresh process. CUDA initialization time can therefore
be included in per-process measurements. No GPU-only kernel time is claimed.
"""
from __future__ import annotations
import argparse
import csv
import math
import platform
import statistics
import subprocess
import time
from collections import defaultdict
from pathlib import Path
NUMERIC_KEYS = ("TOTAL_S", "LOOP_S", "H2D_S", "D2H_S", "EXIT_S", "MASS")
PRECISIONS = ("fp32", "fp64")
VARIANTS = ("cpu", "baseline", "contiguous", "strided")
def parse_metrics(stdout: str) -> dict[str, float | str]:
"""Parse success status and typed metrics without accepting invalid values."""
parsed: dict[str, float | str] = {}
for line in stdout.splitlines():
if "=" not in line:
continue
key, raw = line.strip().split("=", 1)
if key in NUMERIC_KEYS:
try:
value = float(raw.replace("D", "E"))
except ValueError as exc:
raise ValueError(f"Invalid {key}={raw}") from exc
if not math.isfinite(value):
raise ValueError(f"Nonfinite {key}")
parsed[key] = value
elif key in ("CHECK", "METHOD", "PRECISION"):
parsed[key] = raw.strip()
if parsed.get("CHECK") != "PASS":
raise ValueError("CHECK=PASS not found")
if not all(x in parsed for x in NUMERIC_KEYS):
raise ValueError("Missing timing field or MASS")
if any(float(parsed[x]) < 0 for x in NUMERIC_KEYS[:-1]):
raise ValueError("Negative elapsed time")
return parsed
def run_program(binary: Path, params: list[str]) -> tuple[dict[str, float | str], float]:
"""Execute one case, returning solver measurements and total process wall."""
start = time.perf_counter()
result = subprocess.run(
[str(binary.resolve()), *params], capture_output=True, text=True, check=False
)
duration = time.perf_counter() - start
if result.returncode:
raise RuntimeError(
f"{binary} exited {result.returncode}\n{result.stdout}\n{result.stderr}"
)
return parse_metrics(result.stdout), duration
def gpu_name() -> str:
"""Read a connected GPU's name if nvidia-smi is accessible."""
try:
result = subprocess.run(
["nvidia-smi", "--query-gpu=name", "--format=csv,noheader"],
capture_output=True, text=True, check=False, timeout=5,
)
if result.returncode == 0 and result.stdout.strip():
return result.stdout.splitlines()[0].strip()
except (OSError, subprocess.TimeoutExpired):
pass
return "unavailable"
def select_executables(base: Path, target: str) -> dict[tuple[str, str], Path]:
"""Map precision and method to built binaries for the chosen execution target."""
mapping = {}
for precision in PRECISIONS:
for method in VARIANTS:
subdir = "host" if method == "cpu" or target == "host" else "gpu"
suffix = "cpu" if method == "cpu" else "acc"
binary = base / subdir / f"{precision}_{suffix}"
if not binary.is_file():
raise FileNotFoundError(f"Missing binary: {binary}")
mapping[(precision, method)] = binary
return mapping
def main() -> None:
"""Run warmups, rotate variant order and record raw and median timing data."""
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--target", choices=("host", "gpu"), required=True)
parser.add_argument("--nx", type=int, default=1024)
parser.add_argument("--ny", type=int, default=1024)
parser.add_argument("--steps", type=int, default=400)
parser.add_argument("--warmups", type=int, default=1)
parser.add_argument("--repeats", type=int, default=5)
parser.add_argument("--build-dir", type=Path, default=Path("build"))
parser.add_argument("--results-dir", type=Path, default=Path("results"))
args = parser.parse_args()
if args.nx < 5 or args.ny < 5 or args.steps < 1:
parser.error("Need nx/ny >=5 and steps >=1")
if args.warmups < 0 or args.repeats < 1:
parser.error("Need warmups >=0 and repeats >=1")
marker = args.results_dir / ("check_gpu" if args.target == "gpu" else "check_cpu") / "PASS.marker"
if not marker.is_file():
parser.error(f"Numerical validation is required first: {marker}")
files = select_executables(args.build_dir, args.target)
gpu = gpu_name() if args.target == "gpu" else "CPU host fallback"
host = platform.node()
jobs = [(prec, method) for prec in PRECISIONS for method in VARIANTS]
def perform(precision: str, method: str) -> tuple[dict[str, float | str], float]:
params = [str(args.nx), str(args.ny), str(args.steps)]
if method != "cpu":
params.append(method)
values, wall_s = run_program(files[(precision, method)], params)
if values["PRECISION"].lower() != precision:
raise ValueError("Precision mismatch in solver output")
expected_method = "CPU" if method == "cpu" else f"ACC_{method}"
if values["METHOD"] != expected_method:
raise ValueError("Wrong binary or method selected")
return values, wall_s
for precision, method in jobs:
for _ in range(args.warmups):
perform(precision, method)
raw: list[dict[str, str | int | float]] = []
for trial in range(1, args.repeats + 1):
k = (trial - 1) % len(jobs)
for precision, method in jobs[k:] + jobs[:k]:
values, wall = perform(precision, method)
row: dict[str, str | int | float] = {
"target": args.target,
"host": host,
"gpu": gpu,
"precision": precision,
"method": method,
"nx": args.nx,
"ny": args.ny,
"steps": args.steps,
"trial": trial,
"process_s": wall,
}
for name in NUMERIC_KEYS:
row[name.lower()] = float(values[name])
raw.append(row)
print(
f"{args.target} {precision} {method} #{trial}: "
f"total={float(values['TOTAL_S']):.6f}s "
f"loop={float(values['LOOP_S']):.6f}s"
)
args.results_dir.mkdir(parents=True, exist_ok=True)
filename = f"bench_part3_{args.target}_{args.nx}x{args.ny}_{args.steps}steps"
raw_path = args.results_dir / f"{filename}.csv"
with raw_path.open("w", newline="", encoding="utf-8") as handle:
writer = csv.DictWriter(handle, fieldnames=list(raw[0]))
writer.writeheader()
writer.writerows(raw)
grouped: dict[tuple[str, str], list[dict[str, str | int | float]]] = defaultdict(list)
for record in raw:
grouped[(str(record["precision"]), str(record["method"]))].append(record)
summary: list[dict[str, str | int | float]] = []
for (precision, method), records in sorted(grouped.items()):
median_total = statistics.median(float(r["total_s"]) for r in records)
median_loop = statistics.median(float(r["loop_s"]) for r in records)
median_process = statistics.median(float(r["process_s"]) for r in records)
summary.append({
"target": args.target, "gpu": gpu, "precision": precision,
"method": method, "nx": args.nx, "ny": args.ny,
"steps": args.steps, "repeats": len(records),
"median_total_s": median_total, "median_loop_s": median_loop,
"median_process_s": median_process,
})
summary_path = args.results_dir / f"{filename}_median.csv"
with summary_path.open("w", newline="", encoding="utf-8") as handle:
writer = csv.DictWriter(handle, fieldnames=list(summary[0]))
writer.writeheader()
writer.writerows(summary)
print(f"RAW_CSV={raw_path}")
print(f"MEDIAN_CSV={summary_path}")
if __name__ == "__main__":
main()
compose.yaml ― Docker Compose
# Part 1/2と同じNVIDIA HPC SDK 26.9に固定して比較の再現性を保ちます。
# Docker Desktop / WSL2 / Linux Docker Engineのいずれも、
# ホストのNVIDIAドライバ・NVIDIA Container Toolkit設定が必要です。
services:
nvhpc:
build:
context: .
dockerfile: Dockerfile
working_dir: /work
volumes:
- .:/work
deploy:
resources:
reservations:
devices:
- driver: nvidia
count: 1
capabilities: [gpu]
stdin_open: true
tty: true
command: bash
Dockerfile ― 開発用Dockerfile
# 2026-11-01時点で公開済みのNVIDIA HPC SDK公式イメージ。
FROM nvcr.io/nvidia/nvhpc:26.9-devel-cuda_multi-ubuntu24.04
# CPU側の基準値をgfortranで作り、Python標準ライブラリで比較します。
# GPUドライバはホスト側で管理し、コンテナ内には入れません。
RUN apt-get update \
&& DEBIAN_FRONTEND=noninteractive apt-get install -y --no-install-recommends \
python3 gfortran \
&& rm -rf /var/lib/apt/lists/*
WORKDIR /work
付録B：最低限の起動コマンド
# ローカルCPU検証
bash scripts/check_cpu.sh
# コンテナでコンパイラとGPUを確認
docker compose build nvhpc
docker compose run --rm nvhpc nvfortran --version
docker compose run --rm nvhpc nvidia-smi
# GPU版をコンパイル・実行して全格子照合
docker compose run --rm nvhpc bash scripts/check_gpu.sh
# 数値比較をパスしてからGPUの速度を記録
docker compose run --rm nvhpc \
python3 scripts/benchmark.py --target gpu --nx 1024 --ny 1024 \
--steps 400 --warmups 1 --repeats 5