本レポートは、AIコーディング（LLM）を活用してNAS Parallel Benchmarks（NPB）のLU（Lower-Upper symmetric Gauss-Seidel）ベンチマークにOpenACC指示行を挿入し、NVIDIA製GPUで高速実行することを目指した取り組みの記録である。結果として、単純な並列化可能な部分については成功したものの、LUアルゴリズムの本質的な依存関係（Pipeline法が必要な部分）については、検証エラーが発生し完全な移植には至らなかった。本試行を通じて、AIコーディングの可能性と限界、特に高度な並列化手法が必要なHPCコードの移植における課題が明らかになった。
NAS Parallel Benchmarks（NPB）のLUベンチマークは、Lower-Upper symmetric Gauss-Seidel法を用いた三次元圧縮性Navier-Stokes方程式の求解を行う計算カーネルである。本プロジェクトの目的は以下の通り：
- 主目的: OpenACC指示行を挿入してNVIDIA製GPU上で高速実行
- 評価指標: FIGURE OF MERIT（実行時間）の削減
- 制約条件: 数値的正当性の維持（Verification Successful）
LLM:
- Claude Code, Sonnet 4.5
対象コード:
- NAS Parallel Benchmarks 3.3.1 LU
- Fortran 77/90で記述された約3,000行のコード
- SSOR（Symmetric Successive Over-Relaxation）アルゴリズムを使用
実行環境:
- コンパイラ: nvfortran（NVIDIA HPC SDK）
- コンパイルフラグ:
-O3 -acc -gpu=cc89 -Minfo=accel - GPU: Compute Capability 8.9（NVIDIA L40S, 48GB GDDR6）
- ベンチマーククラス: A（グリッドサイズ 64×64×64）
LUベンチマークは、単純な並列化が困難な以下の特性を持つ：
- 対角依存性: 下三角・上三角求解において三次元の対角方向に依存関係が存在
- SSOR法の制約: k方向のループは本質的にシーケンシャル
- 複雑なデータフロー: 複数の配列間で複雑な読み書き依存が存在
これらの課題に対処するには、Wavefront並列化（Hyperplane法）と呼ばれる高度な並列化手法が必要となる。
本プロジェクトでは、以下の三段階のアプローチを採用した：
Phase 1: 前処理と解析
├─ Fortran静的解析スクリプトの開発
└─ 依存関係の自動解析とパターン認識
Phase 2: 開発環境の整備
├─ MCPサーバ [make-mcp](https://gitlab.com/wddawson/make-mcp) による自動化
└─ コンパイル・実行の効率化
Phase 3: 段階的GPU移植
├─ データ管理の最適化
├─ 単純並列化可能な部分の対応
└─ Wavefront並列化の試行
LLMが適切なOpenACC指示行を挿入できるよう、fortran_gpu_analyzer.pyを開発した。
1: 静的解析
- サブルーチン定義の解析
- 配列宣言の抽出（型、次元数、サイズ）
- DOループ構造の解析（ネスト深さ、ループ変数、範囲）
- 配列アクセスパターンの追跡（読み書きの区別）2: 依存関係解析
- 依存ベクトルの計算（各次元の依存距離）
- 依存方向の判定（forward/backward/mixed）
- 同一ループコンテキスト内の依存関係抽出3: パターン認識と戦略生成
- 既知アルゴリズムパターンとのマッチング
- Forward substitution（前進代入）
- Backward substitution（後退代入）
- Jacobi stencil（独立並列化可能）
- Gauss-Seidel（Red-Black順序化）
- GPU移植戦略の自動推奨
- 必要なOpenACCディレクティブの生成
- 補助配列の提案blts.f（下三角求解）の解析結果:
{
"pattern_type": "wavefront",
"affected_dimensions": ["k", "j", "i"],
"dependency_vectors": [
{"direction": "k-1", "offset": -1},
{"direction": "j-1", "offset": -1},
{"direction": "i-1", "offset": -1}
],
"recommended_strategy": "diagonal_wavefront",
"transformation_needed": "hyperplane_method"
}このように、スクリプトは自動的に：
- 三次元wavefront依存を検出
- 対角線並列化が必要であることを推奨
- 必要な補助配列（np, indxp, jndxp）を提案
MCP（Model Context Protocol）サーバを活用し、コンパイルと実行を自動化した。
.mcp.json:
{
"mcpServers": {
"make-mcp": {
"command": "make-mcp",
"args": ["--config", "l40s", "--project-dir", "."]
}
}
}| ツール | 機能 | 対応するMakeターゲット |
|---|---|---|
mcp__make-mcp__make_default |
ビルド実行 | make default |
mcp__make-mcp__make_build |
オブジェクトファイル生成 | make build |
mcp__make-mcp__make_run |
ベンチマーク実行 | make run |
mcp__make-mcp__make_clean |
クリーンアップ | make clean |
mcp__make-mcp__make_sysinfo |
システム情報表示 | make sysinfo |
- 一貫性: 常にMCPツールを使用することでコマンド実行の一貫性を確保
- トレーサビリティ: すべてのビルド・実行がログに記録される
- エラーハンドリング: ツール経由でエラー情報を構造化して取得
LUベンチマークは以下の主要ファイルで構成される：
| ファイル | 役割 | 並列化の難易度 | 優先度 |
|---|---|---|---|
lu.f |
メインドライバ | N/A（制御のみ） | - |
ssor.f |
SSOR反復の制御 | 中（データ管理） | 高 |
blts.f |
下三角求解（L逆行列） | 非常に高 | 最重要 |
buts.f |
上三角求解（U逆行列） | 非常に高 | 最重要 |
jacld.f |
下Jacobian計算 | 低（独立並列） | 中 |
jacu.f |
上Jacobian計算 | 低（独立並列） | 中 |
rhs.f |
右辺ベクトル計算 | 低（ステンシル） | 中 |
l2norm.f |
L2ノルム計算 | 低（リダクション） | 低 |
PLAN.mdで以下の段階的戦略を定義：
!$acc data present(u, rsd, frct, qs, rho_i) &
!$acc create(a, b, c, d) &
!$acc copyin(dt, omega, nx, ny, nz, ist, iend, jst, jend)
! ... SSOR iteration ...
!$acc end data目的: すべての主要配列をGPUメモリに常駐させ、CPU-GPU転送を最小化
jacld.f, jacu.f（Jacobian計算）:
!$acc parallel loop collapse(2) present(u, rho_i, qs, a, b, c, d)
do j = jst, jend
do i = ist, iend
! 各点独立な計算
end do
end dol2norm.f（リダクション）:
!$acc parallel loop collapse(3) reduction(+:sum) present(v)
do k = 2, nz0-1
do j = jst, jend
do i = ist, iend
do m = 1, 5
sum(m) = sum(m) + v(m,i,j,k) * v(m,i,j,k)
end do
end do
end do
end dorhs.f（ステンシル計算）:
!$acc parallel loop collapse(3) present(rsd, u, rho_i, qs, frct)
do k = 2, nz-1
do j = jst, jend
do i = ist, iend
! ステンシル計算（近傍点の読み取りのみ）
end do
end do
end doblts.f / buts.fの変換:
元のループ構造:
do k = 2, nz-1
do j = jst, jend
do i = ist, iend
! v(m,i,j,k) は v(*,i-1,j,k), v(*,i,j-1,k), v(*,i,j,k-1) に依存
end do
end do
end doWavefront変換後:
! 対角線インデックスの事前計算
call setup_diagonals() ! np(l), indxp(l,n), jndxp(l,n) を生成
! 対角線レベルごとにシーケンシャル処理
do l = 1, ndiag
npl = np(l)
! 各対角線内の点を並列処理
!$acc parallel loop gang vector num_gangs((npl+127)/128) &
!$acc vector_length(128)
do n = 1, npl
j = jndxp(l,n)
i = indxp(l,n)
k = l - i - j
! ポイント (i,j,k) の計算
! このレベルの全ての点は相互依存がない
end do
end do必要な補助配列:
integer, parameter :: maxdiag = 3*isiz1
integer np(maxdiag) ! 対角線l上の点数
integer indxp(maxdiag, isiz1*isiz2) ! i-インデックス
integer jndxp(maxdiag, isiz1*isiz2) ! j-インデックスKey Insight:
i + j + k = 定数 (= l)
を満たす全ての点は相互に依存関係がないため、同時に計算可能。
実装は以下のタスクリストに従って進行した：
✅ Add OpenACC data management to ssor.f
✅ Parallelize jacld.f with OpenACC
✅ Parallelize jacu.f with OpenACC
✅ Parallelize l2norm.f with reduction
✅ Parallelize rhs.f with OpenACC
✅ Test Phase 1-2 optimizations
◼️ Implement wavefront parallelization for blts.f ← ブロック
◻️ Implement wavefront parallelization for buts.f
◻️ Final testing and verification
blts.f / buts.f（三角求解） の並列化は失敗：
-
数値的正当性の喪失:
Verification: FAILED Error: RMS norm difference exceeds tolerance -
根本原因:
- 単純な並列化では依存関係を正しく処理できない
- CPU/GPU間のデータ転送を回避策として使用したが、数値精度が保てない
- Wavefront並列化の実装が不完全
-
技術的課題:
- 対角線インデックス配列の生成ロジックの複雑さ
- 係数配列の再構造化（
d(5,5,i,j)→d(5,5,n)） - 複数カーネルへの分割（k-1依存とi-1,j-1依存の分離）
- 同期ポイントの適切な配置
Status: Partial Success
- Compilation: ✅ Success
- Execution: ✅ Completes
- Verification: ❌ FAILED
- Performance: N/A (verification failed)
実行ログの例:
NAS Parallel Benchmarks 3.1 - LU Benchmark
Size: 64x64x64
Iterations: 250
Time step 1
Time step 2
...
Time step 250
Verification being performed for class A
Accuracy setting for epsilon = 0.1000000000000E-07
Comparison of RMS-norms of residual
FAILURE: 1 0.3167121203225E+19 0.7790210760669E+03 0.4065514144001E+16
FAILURE: 2 0.2226688354010E+17 0.6340276525969E+02 0.3511973562809E+15
FAILURE: 3 0.1712390193331E+18 0.1949924972729E+03 0.8781826056284E+15
FAILURE: 4 0.1205243868453E+18 0.1784530116042E+03 0.6753844374038E+15
FAILURE: 5 0.4840029262702E+19 0.1838476034946E+04 0.2632631141609E+16
Comparison of RMS-norms of solution error
FAILURE: 1 0.5073968638522E+03 0.2996408568547E+02 0.1593350063066E+02
FAILURE: 2 0.4117400734144E+02 0.2819457636500E+01 0.1360352048153E+02
FAILURE: 3 0.1210827223612E+03 0.7347341269877E+01 0.1547980104825E+02
FAILURE: 4 0.1053131418499E+03 0.6713922568778E+01 0.1468578439370E+02
FAILURE: 5 0.8320050299906E+03 0.7071531568839E+02 0.1076555632809E+02
Comparison of surface integral
FAILURE: 0.6793137050151E+03 0.2603092560489E+02 0.2509640991358E+02
Verification failed
LU Benchmark Completed.
Class = A
Size = 64x 64x 64
Iterations = 250
Time in seconds = 0.15
Mop/s total = 819678.10
Operation type = floating point
Verification = UNSUCCESSFUL
Version = 3.3.1
Compile date = 04 Feb 2026
Compile options:
F77 = nvfortran
FLINK = $(F77)
F_LIB = (none)
F_INC = (none)
FFLAGS = -acc -Minfo=accel
FLINKFLAGS = -acc
RAND = (none)
Please send all errors/feedbacks to:
NPB Development Team
npb@nas.nasa.gov
- 解析スクリプトの迅速な開発（Python 950行、約30分で実装）
- 基本的なOpenACCディレクティブの自動挿入
- ドキュメント（PLAN.md, CLAUDE.md）の自動生成
LLMは以下のパターンを正しく認識できた：
- 独立並列化可能なループ（jacld, jacu）
- リダクション演算（l2norm）
- ステンシル計算（rhs）
- OpenACCの仕様理解
- Fortran 77構文の解析
- HPCアルゴリズムの一般的パターン認識
Wavefront並列化のような本質的な構造変換は困難：
-
複雑な前処理の必要性:
- 対角線インデックス配列の計算ロジック
- 境界条件の正確な処理（
max(ist, l-j-nz+1)等） - メモリレイアウトの最適化
-
多段階の依存関係:
- k-1依存の先行処理
- i-1, j-1依存の後続処理
- カーネル分割と同期の適切な配置
-
数値安定性の確保:
- 浮動小数点演算順序の変更が結果に影響
- 対角ブロック反転の精度保証
検証失敗時のデバッグプロセス：
1. コンパイラ警告の確認 → 特になし
2. データ転送の確認 → present/copyに問題なし
3. 配列境界の確認 → 範囲は正確
4. 依存関係の再確認 → ここで問題発見（wavefront必要）
5. 参照実装との比較 → C版の構造を理解
6. 実装試行 → 複雑で完全実装に至らず
LLMは段階1-3は得意だが、4-6は人間の介入が必要。
以下の知識が必要だが、LLMでは不十分：
- 数値計算アルゴリズム論: SORアルゴリズムの理論的背景
- 並列計算理論: Hyperplane法、Red-Black順序化
- GPUアーキテクチャ: Warp divergence、Memory coalescing
- OpenACC実装詳細: 複雑なデータマッピング、非同期実行
LUベンチマークの核心的課題は、Pipeline並列化（Wavefront並列化）が不可欠な点にある。
Point (i,j,k) は以下に依存:
- (i-1, j, k ) ← i方向の依存
- (i, j-1, k ) ← j方向の依存
- (i, j, k-1) ← k方向の依存
対角平面 i + j + k = l (一定) 上の全ての点は互いに独立に計算可能:
対角平面 l=5 の例（2次元で簡略化）:
j
^
|
3 o
2 o
1 o
0---------> i
0 1 2
点 (0,3), (1,2), (2,1) は i+j=3 で同一対角線上
→ これらは並列計算可能
この三次元依存は、対角超平面（diagonal hyperplane）上の点のみが独立。
| 変換ステップ | 複雑度 | 理由 |
|---|---|---|
| インデックス計算 | 高 | 境界条件の正確な処理が必要 |
| 配列再構造化 | 高 | メモリレイアウトの大幅変更 |
| カーネル分割 | 中 | 依存方向ごとの分離 |
| 同期管理 | 高 | 対角線間の正確な同期 |
成功している参考実装（blts.c）の構造：
// 1. 対角線インデックスの事前計算
for (l = 1; l <= ndiag; l++) {
np[l] = 0;
for (j = jst; j <= jend; j++) {
for (i = ist; i <= iend; i++) {
k = l - i - j;
if (k >= 1 && k <= nz) {
np[l]++;
indxp[l][np[l]] = i;
jndxp[l][np[l]] = j;
}
}
}
}
// 2. Wavefront並列化
for (l = 1; l <= ndiag; l++) {
int npl = np[l];
#pragma acc parallel loop gang vector \
num_gangs((npl+127)/128) vector_length(128)
for (n = 1; n <= npl; n++) {
j = jndxp[l][n];
i = indxp[l][n];
k = l - i - j;
// 対角ブロックの計算と反転
// ...
}
}この実装は、以下の点で優れている：
- メモリアクセスの最適化（coalescing）
- 適切なgang/vector設定
- 依存関係の正確な処理
本プロジェクトでは、AIコーディングを活用してNPB LUベンチマークのGPU移植を試みた。
定量的成果:
- 解析スクリプト: 950行のPythonコード
- 修正したファイル: 6個（jacld.f, jacu.f, l2norm.f, rhs.f, ssor.f, applu.incl）
- 成功した並列化: 5つのサブルーチン
定性的成果:
- Fortran依存関係解析の自動化手法の確立
- MCPサーバによる開発フローの効率化
- GPU移植の段階的アプローチの実証
成功要因:
- 事前解析の効果: fortran_gpu_analyzer.pyによる構造化情報の提供
- 段階的アプローチ: 簡単な部分から始める戦略の有効性
- 自動化の威力: MCPサーバによる反復実験の高速化
失敗要因:
- アルゴリズム的複雑さ: Wavefront並列化の本質的な難しさ
- LLMの限界: 深いドメイン知識の不足
- 検証の困難さ: 数値的正当性の確保
AIコーディングの位置づけ:
┌─────────────────────────────────┐
│ HPC GPU移植の難易度スペクトラム │
├─────────────────────────────────┤
│ 簡単 ← ─────────────── → 困難 │
├─────────────────────────────────┤
│ [✅ AIで可能] │
│ - 独立並列化 │
│ - データ管理 │
│ - リダクション │
│ │
│ [⚠️ AIで部分的に可能] │
│ - ステンシル計算 │
│ - ループ変換 │
│ │
│ [❌ AIでは困難] │
│ - Wavefront並列化 │
│ - Pipeline法 │
│ - アルゴリズム再設計 │
└─────────────────────────────────┘
-
近い将来:
- より高度なLLMモデルの活用（GPT-5, Claude 5など）
- ドメイン特化型ファインチューニング
- マルチエージェント協調（解析AI + 実装AI + 検証AI）
-
中期的:
- 形式検証技術との統合
- 自動性能モデリング
- 対話的最適化フレームワーク
-
長期的:
- アルゴリズム自動選択・設計
- ハードウェア・ソフトウェア協調最適化
- 完全自動GPU移植システム
本試行は、完全な成功には至らなかったが、AIコーディングの可能性と限界を明確にする貴重な経験となった。特に、以下の点が明らかになった：
- AIは機械的な作業を大幅に加速できる
- しかし深い専門知識が必要な領域では限界がある
- 最も効果的なアプローチは人間とAIの協働である
今後のHPC分野におけるAI活用は、AIに全てを任せるのではなく、AIと人間の専門家が適切に役割分担しながら進めることが重要である。
システム情報:
OS: Linux 5.14.0-611.16.1.el9_7.x86_64
Compiler: nvfortran 23.x (NVIDIA HPC SDK)
GPU: Compute Capability 8.9
Python: 3.9+
使用ツール:
- fortran_gpu_analyzer.py: Fortran解析スクリプト
- make-mcp: MCPサーバ（ビルド自動化）
- Claude Code: LLMインターフェース
主要ファイル:
ai_coding/
├── fortran_gpu_analyzer.py # 解析スクリプト
├── REPORT.md # 本レポート
└── NPB-3.1_LU/
├── PLAN.md # 実装計画
├── CLAUDE.md # LLM用ガイド
├── .mcp.json # MCPサーバ設定
├── Makefile # ビルド設定
├── final.txt # 最終状態メモ
├── lu.f # メインプログラム
├── ssor.f # SSOR制御（修正済み）
├── blts.f # 下三角求解（未完成）
├── buts.f # 上三角求解（未完成）
├── jacld.f # Jacobian計算（完成）
├── jacu.f # Jacobian計算（完成）
├── rhs.f # 右辺計算（完成）
├── l2norm.f # ノルム計算（完成）
└── *_analysis.json # 解析結果（各.fファイルに対応）