LU 分解 並列化プログラム

MPI を用いた LU 分解の並列実装。逐次実装から、1次元列ブロック分割、2次元ブロック分割（同時多段多列消去法）、タイル分割まで段階的に最適化を行っている。

ファイル構成

ファイル	説明
default.c	逐次 LU 分解（ベースライン実装）
lu.c	1次元列ブロック分割による並列 LU 分解（pivoting あり）
lu_block.c	2次元ブロック分割 + ノンブロッキング通信 + ループアンローリング
lu_tile.c	タイル分割による 2次元ブロック並列 LU 分解
Makefile	ビルド用（全ターゲットをコンパイル）
run.sh	ビルド＋実行スクリプト（パラメータをコマンドラインから指定可能）
lu.bash	PBS ジョブスクリプト（2ノード、計196プロセス）
lu-debug.bash	PBS ジョブスクリプト（1ノード、デバッグ用）

各実装の詳細

default.c — 逐次 LU 分解

• MPI 環境上で動作するが、LU 分解自体は逐次処理
• MyLUsolve() で LU 分解 → 前進代入 → 後退代入を実行
• pivoting なし
• 定数: N=8960, NPROCS=224

lu.c — 1次元列ブロック分割（pivoting あり）

• 行列を列方向に NPROCS 個のブロックに分割し、各プロセスが 1 ブロックを担当
• LU 分解: 各ステップで対角要素を持つプロセスがピボット選択・L 列の計算を行い、MPI_Send/MPI_Recv で他プロセスに送信
• 前進代入・後退代入: ブロック単位でウェーブフロント的に処理。対角ブロックを持つプロセスが計算を主導し、結果を順次送信
• pivoting: 部分ピボット選択を実装。perm[] 配列で行の置換を管理
• 定数: N=2240, NPROCS=224

lu_block.c — 2次元ブロック分割 + ノンブロッキング通信

• XPROCS x YPROCS (14x14 = 196) プロセスで行列を 2次元ブロック分割
• 同時多段多列消去法: 左上から右下に向かってブロックを順次処理
• ユーティリティ関数:
  • LU() — ブロック内 LU 分解
  • MatMatSub() — A -= B * C（2x2 ループアンローリング適用）
  • BInvA() — B^{-1}A の計算
  • ABInv() — AB^{-1} の計算
  • LInvA() — L^{-1}A の計算
  • UInvA() — AU^{-1} の計算
• ノンブロッキング通信: MPI_Isend/MPI_Irecv + MPI_Waitsome で通信と計算のオーバーラップを実現
• 前進代入・後退代入: 2次元分割に対応したウェーブフロント処理
• 定数: N=2240, NPROCS=196, XPROCS=14, YPROCS=14

lu_tile.c — タイル分割

• 行列を double A[TILE_NUM][TILE_NUM][TILE_SIZE][TILE_SIZE] の 4次元配列で管理
• TILE_SIZE=64 ごとに LU 分解を実行。キャッシュ効率と MPI 通信の一括送信に有利
• アルゴリズム:
  1. PROCS x PROCS (14x14) で 2次元ブロック分割
  2. 各スーパーブロック内で MatMatSub による更新 → LU 分解
  3. L パネル・U パネルの更新と送信
  4. MPI_Bcast で計算済みタイルを全プロセスに共有
• ユーティリティ関数:
  • TileLU_Local() — タイル内 LU 分解
  • TileMatMatSub_Local() — A -= L * U（2x2 ループアンローリング適用）
  • TileUInvA_Local() — tile * U^{-1}
  • TileLInvA_Local() — L^{-1} * tile
• 前進代入・後退代入は lu_block.c と同様のロジック（A_AT マクロで 4次元配列にアクセス）
• 定数: N=8960, TILE_SIZE=64, NPROCS=4, PROCS=2, ITER_NUM=6

ビルド

make        # 全ターゲット (lu_tile, lu_block, lu, default) をコンパイル
make clean  # 実行ファイルの削除

環境変数でコンパイラを変更可能:

変数	デフォルト	説明
CC	mpicc	MPI コンパイラ
CFLAGS	-O3	コンパイルオプション

実行

run.sh（ローカル実行）

run.sh はビルドと実行を一括で行うスクリプト。問題サイズやプロセス数などの各パラメータをコマンドラインから指定できる。

./run.sh <target> [options]

オプション一覧

オプション	説明	対象	デフォルト
-n N	問題サイズ (行列の次元)	全て	ソースコードの定数
-np NP	MPI プロセス数	全て	ターゲット依存
-x XPROCS	x方向の分割数	lu_block	2
-y YPROCS	y方向の分割数	lu_block	2
-p PROCS	分割数 (PROCS*PROCS = NP)	lu_tile	2
-t TILE_SIZE	タイルサイズ	lu_tile	64
-i ITER_NUM	イテレーション回数	全て	6

デフォルトのプロセス数

ターゲット	デフォルト NP
default	224
lu	224
lu_block	196
lu_tile	4

使用例

# タイル分割: N=1024, 4プロセス, 2x2分割, タイルサイズ32, 10回イテレーション
bash run.sh lu_tile -n 2240 -np 4 -p 2 -t 32 -i 10

# 2次元ブロック分割: N=2240, 4プロセス, 2x2分割
bash run.sh lu_block -n 2240 -np 4 -x 2 -y 2

# 1次元列ブロック分割: N=2240, 4プロセス
bash run.sh lu -n 2240 -np 4

# 逐次: N=1000, 1プロセス, 3回イテレーション
bash run.sh default -n 1000 -np 1 -i 3

PBS ジョブスケジューラ（クラスタ実行）

qsub lu.bash        # 2ノード x 98プロセス = 196プロセス
qsub lu-debug.bash  # 1ノード x 100プロセス（デバッグ用）

検証（テスト）

各プログラムには組み込みの検証ルーチンが含まれている。

検証方法: Ax = b の解 x が全て 1.0 になることを検証する。

1. 行列生成: 対称行列 A を生成（MATRIX=1 の場合）
2. 右辺ベクトル: b = A * [1, 1, ..., 1]^T として設定（解が x = [1, 1, ..., 1]^T になるように）
3. LU 分解実行: LU 分解 → 前進代入 → 後退代入で x を求める
4. 誤差チェック: 各プロセスが担当分の ||x - 1||_2 を計算し、MPI_Reduce で集約
5. 判定: 集約した誤差が EPS * N * N^3（= マシンイプシロン x 問題サイズに依存する閾値）以下なら "OK! Test is passed." を出力

Pass value: 1.598930e+07      # 許容誤差
Calculated value: 2.345678e-05  # 実際の誤差
 OK! Test is passed.           # テスト合格

全実装で ITER_NUM 回（デフォルト6回）の反復実行を行い、1回目をウォームアップとして除外した ITER_NUM - 1 回の平均実行時間を報告する。run.sh の -i オプションで回数を変更可能。

性能結果

測定環境

ローカル PC

項目	値
CPU	AMD Ryzen 5 3580U (4コア / 8スレッド)
クロック	1.4 - 2.1 GHz
L1d / L1i	128 KiB / 256 KiB (4 instances)
L2	2 MiB (4 instances)
L3	4 MiB
OS	Linux 6.18.3-arch1-1 (x86_64)
コンパイラ	GCC 15.2.1 (mpicc / OpenMPI 5.0.9)
最適化	-O3

miyabi スパコン

項目	仕様
CPU	Intel Xeon Platinum 8480+ (Sapphire Rapids)
ソケット数	2
コア数/ソケット	56 (スレッド数/コア: 2)
合計論理 CPU 数	224
ベースクロック / 最大クロック	800 MHz / 3.8 GHz
L1d / L1i	5.3 MiB / 3.5 MiB (112 instances)
L2	224 MiB (112 instances)
L3	210 MiB (2 instances)
NUMA ノード	2
命令セット	AVX-512, AMX-BF16/INT8

ローカル PC: 4プロセス (2x2), ITER_NUM=6, ウォームアップ1回除く5回平均

N	default.c 1並列 (s)	lu.c 4並列 (s)	lu_block.c 2x2 (s)	lu_tile.c 2x2 (s)
1120	0.578884	0.719870	0.182842	0.081484
2240	4.36	4.81	2.94	0.38
3360	14.930323	16.305339	10.518179	1.459952

miyabi スパコン: 大規模並列

N	default.c (s)	lu.c 224並列 (s)	lu_block.c 14x14 (s)	lu_tile.c 14x14 (s)
2240	4.18	1.63	0.0331	0.177
4480	35.1	4.41	0.199	0.325
8960	258	12.12	2.25	1.13

参考文献

• Alfredo Buttari et al. "A class of parallel tiled linear algebra algorithms for multicore architectures". Parallel Computing 35.1 (2009), pp. 38-53.
• Jack J Dongarra et al. "The LINPACK Benchmark: past, present and future". Concurrency and Computation: practice and experience 15.9 (2003), pp. 803-820.