KOBE HPC サマースクール（初級）2019 講義資料
2019/08/26
兵庫県立大学シミュレーション学研究科
安田 修悟

連立一次方程式の解法1• LU分解
行列Aを下三角行列L（対角成分は１）と上三角行列Uに分解する。

A = LU
Aが3X3行列の場合：! a# 11# a21##" a31a12a22
a32

a13 $ ! 10& #a23 & = # b21 1&& ##a33 % " b31 b32$!
0 &# c11 c120 &# 0 c22&&##
1 %" 00c13 $&c23 &&&
c33 %連立一次方程式の解法1! !
連立一次方程式 Ax = y! !
• LU分解を行う。 LU x = y※LU分解の方法については補足資料参照．

!
!
! !
• z = U x と置くと、 Lz = y と書ける。
! z $ ! cc
# 1 & # 11 12# z2 & = # 0 c22##&& ##z
0" 3 % " 0c13 $! x1 $&#&&#c23
x2 &&#&c33 &%#" x3 &%!
0# 1# b21 1##
" b31 b32

$ !
$$!
0 &# z1 & # y1 &0 &# z2 & = # y2 && #&&#
1 &%#" z3 &% #" y3 &%連立一次方程式の解法1! !
連立一次方程式 Ax = y
解法の手順! !
!
1. Lz = y を解いて、 z を求める。
!
! !
2. U x = z を解いて、元の方程式の解 x を求める。

連立一次方程式の解法1! !
!
1. Lz = y を解いて、 z を求める。
!
0# 1# b21 1##
" b31 b32$! z $ ! y $0 &# 1 & # 1 &0 &# z2 & = # y2 && #&&#
1 &%#" z3 &% #" y3 &%(1) z1 = y1
(2) z2 = y2 − b21z1(3) z3 = y3 − b31z1 − b21z2

z1から順次，z2, z3と求められる。

連立一次方程式の解法12.

! !
U x = z を解いて、元の方程式の解! cc# 11 12# 0 c22##0" 0!
x を求める。

c13 $! x1 $ ! z1 $&#& #&c23 &# x2 & = # z2 &&#& #&c33 &%#" x3 &% #" z3 &%(1) x3 = z3 / c33(2) x2 = (z2 − c23 x3 ) / c22(3) x1 = (z1 − c12 x2 − c13 x3 ) / c11x3，x2, x1と順次，求められる。

連立一次方程式の解法2• ヤコビの反復法（Jacobi method)• ヤコビの反復法の形式

xi = c0 x0 + c1 x1 +.... + ci−1 xi−1 + ci+1 xi+1 +.... + cN x N
1. x のi成分について，成分i以外のベクトル要素の線形結合で表す．
2. i=0,…,NのN+1個の連立一次方程式を反復法で解く．

• ヤコビ反復法 （３×３行列）
非対角成分を除いて右辺に移行• ヤコビの反復法（３×３行列）

適当な初期値を使って反復計算を実行．

解が収束するまで反復する．
(x,y,z)が収束する場合には、その値は連立一次方程式の解を与える．

• ヤコビの反復法 （３×３行列）– 収束判定の例

元の方程式との残差を求めて十分小さくなるまで反復する．
サンプルコード yabobi_3by3.c

連立一次方程式の解法• サンプルコード
次の連立一次方程式をヤコビ法とLU分解を使ってそれぞれ解くプログラムがyacobi_solver.c，lu_solver.cです．
" 5$$ 0$ 1$ 1$
$ 1
#

1 2 0
2 −1 1
1 3 00 0 10 −1 0

" x % "%−1 $ 1 '5 %'$$''x2 '2$ 0 '$'
1 '$ x3 ' = $ −2 '0 '$ x4 ' $ 0 ''$'5 '&$$ x '' $# 4 '&# 5 &

熱伝達問題の差分解法• 温度場を計算する。（ラプラス方程式）−Δ𝑇 =%&'%( &+%&'%* &= 10

T=0T=0T=0yxT=0• 温度場を計算する。（ラプラス方程式）−Δ𝑇 =%&'%( &+%&'%* &= 10• 差分化する。
𝑇 𝑥. , 𝑦1 = 𝑇.,1 ,𝑖𝑥. = 𝑦. =𝑖 = 0, … , 𝑁 ,𝑁
Δℎ = 1/𝑁• 温度場を計算する。（ラプラス方程式）−Δ𝑇 =%&'%( &+%&'%* &= 10• 差分化する。

1𝑇.,1 = 10Δℎ8 + 𝑇.9:,1 + 𝑇.;:,1 + 𝑇.,19: + 𝑇.,1;:
4(𝑖, 𝑗 = 1, … , 𝑁 − 1)• 𝑇.,1 の線形連立方程式

𝑇?,1 = 𝑇@,1 = 0, 𝑗 = 0, … , 𝑁
𝑇.,? = 𝑇.,@ = 0, (𝑖 = 0, … , 𝑁)

境界条件
として固定

ヤコビの反復法
1(A)
(A)
(A)
(A)
8= 10Δℎ + 𝑇.9:,1 + 𝑇.;:,1 + 𝑇.,19: + 𝑇.,1;:
4(𝑖, 𝑗 = 1, … , 𝑁 − 1)
(A;:)
𝑇.,1

サンプルコード laplace.c

演習４• 熱伝達問題のコードで以下のことについて調べよ．
– gnuplotで温度分布を可視化せよ．
gnuplot>splot “laplace.dat”>set pm3d map>replot– 水平参照と垂直参照（2重ループの順番）を入れ替えて、計算時間を比較せよ．
– 次頁で紹介するSIMD命令の効果について確認せよ．

ベクトル化• SIMD（Single Instruction Multiple Data）– 1命令で複数要素の演算を行う．

for(i=0;i<n;i++){C[i]=A[i]+B[i];}• [SSE]: float4要素を一度に計算• [AVX]: float8要素• [AVX-512]: float 16要素
(兵県大の計算機は未対応）

A[0]

A[1]

A[2]

A[3]

+ B[0]

B[1]

B[2]

B[3]

C[0]

C[1]

C[2]

C[3]

コンパイラーによるベクトル化• ループの自動ベクトル化（-xまたは-axオプション）– 最適化レポート（-qopt-reportオプション）でベクトル化状
況を確認
icc -xhost laplace.c -qopt-report– -xhost: コンパイルを行うホスト・プロセッサーで利用可能な最上位の命令セット向けのコードを生成.
– laplace.optrptにベクトル化状況を報告.
– -no-vec: コンパイラーによるベクトル化を無効にする．

（補足）LU分解の方法

A = LU
Aが3X3行列の場合：! a# 11# a21##" a31a12a22a32

a13 $ ! 10& #a23 & = # b21 1& #a33 &% #" b31 b32$!
0 &# c11 c120 &# 0 c22&#
1 &%#" 00c13 $&
c23 &&
c33 &%

具体的に書くと、! a# 11# a21##" a31a12
a22a32

a13 $ ! c11c12& #a23 & = # b21c11
b21c12 + c22&& ##a33 % " b31c11 b31c12 + b32 c22$&&
b21c13 + c23&&
b31c13 + b32 c23 + c33 %c13

（補足）LU分解の方法
1. 行列Aを次のようにA’に変換する。

# aa1211%
A! = % a21 / a11 a22 − a!21a12%%! a12$ a31 / a11 a32 − a31

L, U成分での表現は？

& # cc12( % 11a23 − a!21a13 ( ⇒ % b21c22( %! a13 (' %$ b31 b32 c22
a33 − a31a13&((c23(b32 c23 + c33 ('c13

2. 行列A’をさらに次のようにA’’変換するとL，Uの全ての要素が決定する。

# a!
% 11
A!! = % a!21%%!
$ a31!
a12a!22! / a!22a32& # cc( % 11 12( ⇒ % b21 c22a!23(( %%! − a32!! a!23 ' $ b31 b32a33!
a13c13 &(c23 (((c33 '（補足）LU分解の方法! a! a1N# 11# " # "## a! aNN
" N1$&&&&%" cc12$ 11$ b21 c22$
$ " #$ bN1 !
#!
#bNN−1

c1N %'
c2 N ''" 'cNN '&行列AをLU分解の要素行列に変換するプログラム。
for(k=0;k<N-1;k++){w=1/a[k][k];for(i=k+1;i<N;i++){a[i][k] *= w;
for(j=k+1;j<N;j++){a[i][j] -= a[i][k]*a[k][j];}}}サンプルコード lu_decomp.c※ U行列の対角成分cii にゼロや小さな値が出てくる場合には、更に、行や列の入れ替えを伴う複雑な手順を考える必要がある。

アクセラレータとは？
• CPUの処理を代替して処理の効率を向上させる装置（ディバイス）。
– GPGPU (NVIDIA), Xeon Phi (Intel), APU (AMD)など
ホスト

コ
ア
０

コ
ア
１

コ
ア
2

ホスト側メモリ

アクセラレータコ
ア３

制御

Many coresMulti Threads

デバイスメモリ
PCIバス

アクセラレータの実行イメージ
ハイブリット構成 （CPU＋アクセラレータ）

ホスト

重たい計算部分をオフロード

CPU

GPU・Phi等
PCIバス

Mainメモリ

アクセラレータ

デバイスメモリ※データのオフロード通信がボトルネックとなる• CPUのメモリ空間とディバイスのメモリ空間が異なる• 一般的なメモリ帯域幅に比べてPCIバスでの低速な転送OpenACC• Accelerators用のAPIの一つ。
• CPUからディバイスに処理をオフロードする際に使用する。
• Fortran/C/C++言語でディレクティブを指定（OpenMPと同様）• CUDA等の専用プログラムに比べて、導入が容易。

OpenACCプログラミング• 基本フォーマット
#pragma␣acc␣ディレクティブ␣節１␣節2␣・・・{
構造化ブロック（途中で抜け出したり、終了しない）}
例#paragma ␣acc ␣kernels␣copy(a)for(i=0;i<imax;i++){a[i]=a[i]*a[i];・・・・}

OpenACCプログラミング• 基本フォーマット (Fortranの場合)!$acc␣ディレクティブ␣節１␣節2␣・・・
構造化ブロック!$acc end ディレクティブ!$acc ␣kernels␣copy(a)do i=0, imax
a[i]=a[i]*a[i];・・・・
end do!$acc end kernels

基本的なディレクティブ構文① Accelerator compute構文– オフロードするループ対象部分を指定．
– kernels構文とparallel構文．

② DATA構文– ホストとアクセラレータのデータ転送を明示的に指示．
– 高速化の鍵を握る．

③ Loop構文–①のオフロード領域のループのベクトル長や並列分割の詳細を明示的に指示することができる．

PGIコンパイラー• Module環境の設定
module␣load␣pgi– module␣listでintelが設定されている場合
module␣unload␣intelでIntel環境をunloadする。

• OpenACCのコンパイルpgcc␣-acc␣<prog.c> ␣-o␣acc.out
pgfortran –acc <prog.f90>– “-acc”でOpenACCディレクティブを有効にする．

PGIコンパイラー• OpenMPのコンパイルpgcc␣-mp␣program.c ␣-o␣mp.out“-mp”はOpenMPディレクティブを有効にするためのオプション。
• 最適化オプション– “-fastsse” : 一般的な最適化オプション– “-O3”, “-O4”: より高度な最適化オプション（”-fastsse”より後に入れる
例 pgcc␣-fastsse␣-O3␣program.c）

PGIコンパイラー• 環境変数（実行前に設定）– PGI_ACC_TIME• 実行後に簡易プロファイル情報を標準出力に出力。
• 0（ゼロ）を指定すると機能を抑制（デフォルト）、ゼロ以外の整数値で機能する。

– PGI_ACC_NOTIFY• ディバイスの実行イベントを表示する• 整数値１：Kernel launch のイベントを出力、整数値２：データ転送のイベントの出力、など• 0（ゼロ）を指定すると機能を抑制（デフォルト）GPUキューの実行• インタラクティブ実行
qsub␣-I␣-q␣G

GPUを搭載したノードにてプログラムを実行。

• ホストプログラムの並列化に応じてジョブスクリプ
トを作成． (以下の例は，ホストが並列化無しの場合）#!/bin/bash#PBS␣-q␣G#PBS␣-I␣select=1:ncpus=1:mem=1gb#PBS␣-N␣job_gpu#PBS␣-j␣oesource␣/etc/profile.d/modules.shmodule␣purge
・現在のモジュール環境を全て破棄
module␣load␣pgi
・PGIコンパイラの環境設定
cd␣${PBS_O_WORKDIR}
dplace␣./a.out

OpenACCの演習1• pi1.cとpi2.cのループにOpenACC指示文#pragma acc kernels
を入れて書き換え，CPU計算とGPU計算の時間を比較せよ．
• pi1.cとpi2.cのオフロード処理の違いに着目して、GPUでの計算時間の違いについて考察しなさい．

DATA構文• #pragma␣acc␣data␣節• 節には以下のようなものがある．（Kernel構文の節にも使える．）– copyin(a,b): 配列a, bをホストからアクセラレータにコピーする．
– copyout(c): 配列cをアクセラレータからホストにコピーする．
– copy(A)：copyinとcpyoutの双方を行う．
– create(B): アクセラレータでローカルに使用する配列Bの領域を作成．
– present(C): 対象となる処理に入ったときにすでにデバイスメモリにある配列Cのデータを使用する．

DATA構文• pi2.cのOpenACCプログラムに対して，1. GPU内だけで必要な変数を予め作成する．
data create(f,x,y){.......
}
2. Kernel構文でデータ属性を指示する．
– kernel present(f) copyin(dh,nmax)– kernel present(f) copy(sum) copyin(nmax)

OpenACCの演習2• ラプラス方程式のOpenACCプログラム
1. laplace.cの主要なループをkernel構文を使って書き換える．
2. データ構文を使ってオフロード処理を出来るだけ抑制するコードを作成する．
3. シリアル計算とGPU計算の時間を比較する．

OpenACCの演習３• セルオートマトンのOpenACCプログラム(セルオートマトンについては補足資料を参照）
1. r184.cの主要なループをkernel構文を使って書き換える．
2. データ構文を使ってオフロード処理を出来るだけ抑制するコードを作成する．
3. シリアル計算とGPU計算の時間を比較する．

（補足）セルオートマトン• 格子状のセルと単純な規則による，離散的計算モデルである．（ウィキペディアより）• 交通流モデルの一つであるルール184を扱う．
– １つのセルに入ることのが出来る車両は最大1台．
– 車両は１タイムステップで右に１セル分だけ進める．
– １つ先のセルに車両が存在すれば，次のタイムステップでは先のセルに進めず，そのセルに留まる．
– １つ先のセルに車両が存在しなければ，次のタイムステップに先のセルに進める．

（補足）セルオートマトン• 境界条件– 周期境界条件・・・右端のセルから出て行く車両は左端のセルに移動する．円周上に並んだセルを考えることができる．
– 開放境界条件・・・左端のセルに車両が入ってくる確率と，右端のセルから車両が出て行く確率を予め設定する．

全体の流れ（補足）フローチャート

初期密度r0 入力
初期配置作成状態出力
状態更新
Cell_old から Cell
状態コピー
Cell_old := Cell

Yesn<NMAX
NoFin

IMAX: セルサイズ
NMAX: タイムステップ数
Cell[IMAX]; 現在の状態
Cell_old[IMAX]; 過去の状態

状態更新の方法は？