GTC Japan 2014
2014年7月16日
cuBLASとcuSPARSEを援用した有限要素法のCUDA Fortran実装
出川智啓
（沼津工業高等専門学校電子制御工学科）
（現
望月麟太郎
電気通信大学情報理工学部）
概要
実装方法の特徴
 カーネルの記述が不要
 CUDAプログラミング経験が浅くても比較
的容易に実装可能
 Fortran 2003の機能を利用してcuSPARSE
ライブラリの関数を直接呼び出し（部分的
にC言語で記述したラッパを併用）
 有限要素法をGPUに実装するにあたり，cuBLASと
cuSPARSEライブラリを利用した実装方法を提示
 有限要素法の計算過程を，ベクトルのスカラ倍，
ベクトル成分同士の和，ベクトル成分同士の積，
ベクトルの内積，行列－ベクトル積に分解し，
cuBLASおよびcuSPARSEライブラリの関数を適用
課題
 ライブラリを用いることによる処理の冗長
化への対策
 係数行列の（再）構成
 Fortran の Interface と cuSPARSE Function
呼び出しの不整合の原因究明と対策
有限要素法による流体計算
支配方程式
無次元化された非圧縮性流体の連続の式
とNavier‐Stokes方程式
n n
n n
2 n
2 n 
~


1
u
u
u
v
u
u 




n
 2 



u  u  Δt 
2 
Re  x
y
y 

 x

n n
n n
2 n
2 n 


v~  v n  Δt u v  v v  1   v   v 

2
2 


Re  x
y
y 
 x

u v

0
x y
 u
u
u
p
1   2u  2u 
 2  2 
v


 u
x
y
x Re  x
y 
 t

2
2







v
v
v
p
1
v
v
 u v
 2  2 


 t
x
y
y Re  x
y 

x, y
t
u, v
p
Re
u~, v~
有限要素方程式
流速修正法[1,2]
 2 p n 1  2 p n 1
1  u~ v~ 





2
2
x
y
Δt  x y 
~
1
 x n n
n
y
n n　n
-1
[D]u m
u  u  Δt  [G ]u u  [G ]u v 

Re



 v~  v n  Δt  [G x ]u n v n  [G y ]v n v n  1 [D]v n m -1

Re


[D] p
n 1

1
~  [G y ]v~
[G x ]u

Δt
Compute r0=−[D]p0. Set r0*=r0,−1=0.
For j=0,1,…, until ||r||/|||| < , Do
cj = rj+j−1(cj−1−j−1[D]cj−1)
j = (r0*, rj)/(r0*, [D]cj)

zj = rj−j[D]cj
j = ([D]zj, zj)/([D]zj, [D]zj)
 
 
~  Δt G x p n1m 1
u n1  u
 n1 ~
y
n 1
1


v
v
Δt
G
p
m

 n1 ~
p n1
u  u  Δt
x
n 1

p

 v n1  v~  Δt

y
: 空間方向
: 時間
: x, y方向速度
: 時間
: レイノルズ数
: x, y方向中間速度
Bi‐CGStab法[3]
Galerkin法による定式化[1,2]
pj+1 = pj+jcj+jzj
rj+1 = zj−j[D]zj
j = (r0*, rj)/{j(r0*, [D]pj)}
[1] 内山知実，Javaによる連続体力学の有限要素法，森北出版，2001．
[2] 日本計算工学会流れの有限要素法研究委員会編，続・有限要
素法による流れのシミュレーション，シュプリンガージャパン，2008．
[Gx],[Gy] : x, y方向勾配行列
[D]
: 拡散行列
m
: 集中化質量行列
EndDo
[3]藤野清次，張紹良，反復法の数理，朝倉書店，1996．
cuBLAS, cuSPARSEによる実装
計算手順の置換
cuSPARSE用Interface
有限要素方程式，Bi‐CGStab法の計算手順を
 ベクトルのスカラ倍
 ベクトル成分同士の和
 ベクトル成分同士の積
 ベクトルの内積
 行列－ベクトル積
の組合せに置換
cuBLAS/cuSPERSE
Function
Procedure
*1
ベクトルのスカラ倍
cublasDscal
ベクトル成分同士の和
cublasDaxpy
ベクトル成分同士の積*1
cublasDgbmv
ベクトルの内積
cublasDdot
行列－ベクトル積*2
cusparseDcrsmv
Fortran 2003から追加されたiso_c_bindingの機能を利用してinterfaceを記述．CUDA Fortran
用のmoduleファイルが提供されていなくてもcusparse.libをリンクするだけで呼び出し可能
中間速度計算の置換
1
 x n n
n
y
n n　n
~
u  u  Δt  [G ]u u  [G ]u v 
[ D]u m -1
Re


cuSPARSE Helper Function用interfaceの例
integer(c_int) function cuSparseCreate(handle) bind(C,name="cusparseCreate")
use iso_c_binding
implicit none
type(c_ptr) :: handle
end function cuSparseCreate
integer(c_int) function cuSparseGetVersion(handle,version) bind(C,name="cusparseGetVersion")
use iso_c_binding
implicit none
type(c_ptr),value :: handle
integer(c_int) :: version
end function cuSparseGetVersion
Du ← cusparseDcrsmv([D], un)
cuSPARSE level 1 Function用interfaceの例
cuSPARSE level 2 Function用interfaceの例
← cublasDgbmv(un, m)
integer(c_int) function cuSparseDaxpyi(handle, & N_Nonzero, alpha, x, idx, y, INDEX_BASE) &
bind(C,name="cusparseDaxpyi")
use iso_c_binding
implicit none
! y = y+alpha*x
type(c_ptr) ,value :: handle
integer(c_int),value :: N_Nonzero
real(c_double),value :: alpha
real(c_double),device :: x(:)
integer(c_int),device :: idx(:)
real(c_double),device :: y(:)
integer(c_int),value :: INDEX_BASE
end function cuSparseDaxpyi
integer(c_int) function cuSparseDcsrmv(handle,OPERATION,&
N_Row,N_Col,N_Nonzero,alpha,descr,Aval,Arowptr, &
Acolidx,x, beta, y) bind(C,name="cusparseDcsrmv")
use iso_c_binding
implicit none
! y = alpha*OPERATION(A)*x + beta*y
type(c_ptr),value
:: handle
integer(c_int),value :: OPERATION
integer(c_int),value :: N_Row
integer(c_int),value :: N_Col
integer(c_int),value :: N_Nonzero
real(c_double) :: alpha
type(c_ptr) ,value :: descr
real(c_double),device :: Aval(:)
integer(c_int),device :: Arowptr(:)
integer(c_int),device :: Acolidx(:)
real(c_double),device :: x(:)
real(c_double) :: beta
real(c_double),device :: y(:)
end function cuSparseDcsrmv
uu ← cublasDgbmv(un, un)
vu
← cublasDgbmv (vn, un)
Guu ← cusparseDcrsmv([Gx], uu)
Guv ← cusparseDcrsmv ([Gy], vu)
~
u
~
u
~
u
~
u
~
u
~
~
← cublasDaxpy(u, −t×Gvu)
~
← cublasDaxpy(u, t/Re×Du)
~ −1
u
← cublasDaxpy(u, −t×Guu)
← cublasDgbmv( , m )
速度修正の置換
 
~  Δt G x p n1m 1
u n1  u
実行して結果を正しく得ることができるが，
本来はポインタ渡しする変数を値渡しにし
ないと正しく動作しない．なぜ？
 ベクトル成分同士の積はcuBLAS, cuSPARSEに無い
 帯行列－ベクトル積cublasDgbmvを利用し，対角行列
Gp ← cusparseDcrsmv([Gx], pn+1)
Gp ← cublasDgbmv(Gp, m−1)
~ ← cublasDaxpy(u
~, −tGp)
u
~
un+1 ← u
とベクトルの積として計算
*2 疎行列の格納形式にはCRS形式を採用
現状では，Fortranの機能だけではcuSPARSEを
扱えず，C言語でラッパを作成する必要がある
cuSPARSE level 1 routineと同様にinterfaceを作成しても，
CUSPARSE_STATUS_INVALID_VALUEが発生
有限要素法プログラムの作成と移植，性能評価
有限要素法プログラムの作成




計算機環境
CPU Intel Xeon W3565
GPU NVIDIA Tesla C2075
CUDA 5.0, 5.5
PGI CUDA Fortran 14.6
計算結果
圧力，速度分布(t=600)
節点数 7414
要素数 14400
U
円柱後流に放出される
カルマン渦列がとらえ
られている．ストローハ
ル数や円柱が受ける
抗力が既存の計算結
果とよく一致した．
16d
プログラム作成，GPU移植担当
沼津高専5年生（当時）
プログラミング言語の経験
 C言語 1年
 Java
半年
 数値計算の経験無し
計算条件
62d
移植に要した時間
有限要素法プログラム
の作成（C言語）
3ヶ月
（6時間/週，うち2週間は
定期試験により中断）
GPU移植
（CUDA 5.0）
1週間
CUDA Fortran
へ移植
1日
y
x





レイノルズ数 Re=Ud/=100
円柱直径
d=1
流入流速
U=1
計算時間間隔 t=0.01
計算終了時間 t=1000
計算時間
CPU
1週間
GPU
GPU
（CUDA C）
（CUDA Fortran）
8時間
8時間
CUDAプログラミング初心者でも，C
言語で作成した各処理をカーネル
に置き換えるだけで有限要素法を
GPUへ移植可能である．