フジタ技術研究報告第 48 号 2012 年

CUDA を用いた並列数値解析手法の一考察

仲 沢 武 志概要

従来、大規模な問題や複雑な形状が計算対象の場合では、スーパーコンピュータのような高価な計算リソースを使用する必要があった。これに対し、コンピュータのハード的な進化によって、現在では市販のパーソナルコンピュータ(PC)でも計算機能は目覚しく向上し、並列演算すらも比較的容易に可能となってきている。
マルチコアな PC においては、OpenMP や MPI などの CPU による並列化が可能である。また、最近では、Graphic で使用される計算機能を数値解析に適用する GP-GPU（General Purpose on Graphic Processor Unit）という手法を適用することも実施され適用が拡大されつつある。
本文では、この GP-GPU におけるフレームワークの中で Nvidia が提唱している CUDA を採用し、いくつかの例題を試行してその適用性を検討する。

Evaluation of parallel numerical calculation method using CUDA system
AbstractPreviously large size memory problems and some complex shape problems required use of asupercomputer for analysis. However recently, personal computers have achieved sufficiently highperformance to undertake these activities, and some of these personal computers can perform parallel
calculations in numerical analysis.
Personal computers with multi core CPU(s) can perform as a system for open MPI or MPI to parallel
calculation. At the same time, other components within the system such as high-performance graphics
cards can be used for a parallel calculation procedure called General Purpose on Graphic Processor Unit
(GP-GPU).
We used Nvidia’s CUDA framework for GP-GPU systems, and applied it to some problems to checkperformance of these systems.

キーワード：有限要素法 高速化
並列計算 GP-GPU CUDA

－31－フジタ技術研究報告 第 48 号

§１．はじめに

本文では、比較的簡単に並列化が実現できる GP-GPUのうち Nvidia が提唱している CUDA を使った試行に主眼

有限要素法などの数値解析は、結局のところ微分方程を置き、いくつかの例題でその適用性を検討する なお、

式を離散化し連立１次方程式に帰着させて、所要の境界

ここでの議論では CPU による並列化は対象外としてい条件の下に解を求める計算方法であり、使われ始めてかる。

らかなりの年月が経過している。この間、計算機の性能な§２．PC のハード構成

どの計算リソースの開発と関連しながら数々の計算アルゴリズムが発展して計算実効の効率化が追及されてきた。
しかしながら、最近におけるこの分野の技術は、CPU 単

本文では、計算効率の向上にグラフィック機能で並列演体の性能を飛躍的に進歩させる技術が頭打ちとなり、ま算を実施する GP-GPU を考えることとしている。

た計算アルゴリズムも一定の手法に落ち着くなど、停滞しPC の標準的なハード構成と GP-GPU における拡張例をつつあった。

図１に示す。通常の PC の構成を右側の破線内に示していところが、ここ最近、計算リソースの方向性は、いわゆるる。PC には、通常描画用のグラフィックボードが搭載されてPC と言われる個人で使うコンピュータでもマルチコアは当いる。このグラフィックボード内の複数のプロセッサーを並たり前になってきている。このような状況から計算手法も複列に使って計算を効率化するシステムが GP-GPU である。

数ある CPU を合理的に使用した並列アルゴリズムを試行

１枚のグラフィックカードを描画と同時に演算に使用することして、計算効率の向上を望める状況が復活してきた。さらも可能であるが、システムに対応したグラフィックカードを複に 、 GP-GPU(General Purpose on Graphic Processor

数搭載して、両者の役割を分担した方が効率は良いものとUnit)という概念で、グラフィックエンジンを数値解析などの予想される。

一般的用途にも応用しようという動きも活発になっている。

そこで本文では、図 1 の左側に示したように、グラフィック

GP-GPU によって並列化を実装するフレームワークにはボードを通常の描画用のものとは別にさらに１枚実装したハOpenCL および Nvidia の CUDA や AMD の ATIStreamード環境としている。などがある。

図１ PC の標準構成と CUDA での拡張例－32－

CUDA を用いた並列数値解析手法の一考察§４．適用例

以下にグラフィック以外の全体的なハード構成も示しておく。
CPU：Intel Core i-7-2600K

時間発展する物理的化学的現象の予測は、一般に多大

GPU-1：Nvidia Ge-Force GTX580な計算時間を必要とすることが多い。そこで、CUDA による（ディスプレイ表示用）

並列演算が可能となれば、計算時間が短縮され、作業が効

GPU-2：Nvidia Tesla C 2070

率化する。

（演算用）

有限要素法のような数値解析手法は、結局微分方程式をMB：ASUS SABERTOOTH P67

連立１次方程式に帰着させるわけであるが、連立１次方程チップセット：Intel P67 Express

式の演算時間は、その次数（自由度）の増加に連れて長く実装メモリー総計：16GByteなることは周知である。並列演算が最大限に活用できた場

OS：Windows7 Professional

合の効用は、グラフィックカード上の MSP の数にもよるが、図３の模式図が示すように、連立１次方程式の次数によら§３．GP-GPUずほぼ一定の時間で計算を実施できることにある。
実際には、演算チップ数の制約から、このような効用を得3.1 グラフィックカードの構成ることができる最大の自由度が存在するが、本文ではそれグラフィックボードの内部を模式的に示したものが図２でに関しては議論の対象とせず、2000 自由度程度での並列ある。演算最小単位はストリーミングプロセッサ（SM）であり、演算の適用性を調べる。

8 個の SM で構成されるセグメントがマルチストリーミングプ

時間発展問題の例として、まず地下水汚染などの状況をロセッサ（MSP）と言われている。MSP はハードのグレードに予測する移行拡散方程式に CUDA による並列演算を適用

より異なるがグラフィックボードに 256 個あるいは 512 個搭載する。拡散方程式への適用は比較的事例情報も多く、手始されている。これらが、CPU やメインメモリと通信しながら並めの検討に適切な例と考えたためである。

列演算を実施することになる。

その他に、建設分野で当然遭遇する変形解析および、Laplace 方程式への適用を試みる。この２つの例への適用

3.2 プログラミング技法

は、そもそも間接手法による連立１次方程式における係数プログラミングにおいては、基本的に並列演算させたい行列の性質の違いを調べるためである。

部分をサブルーチンで記述する。Ｃ言語での記述を例示す
ると、所要のメモリーを確保した上で、4.1 移行拡散解析への適用(sub_name)<<gd,BSZ>>(argument）

移行拡散方程式は次式のように与えられている。などとして呼び出す。
またサブルーチン側は__global__ void sub_name(argument)などと宣言して受ける。
文法の詳細は参考文献 2）～4)を参照。

図３ 並列化による計算時間短縮のイメージ
図２ グラフィックカード内の演算機能 1)－33－

フジタ技術研究報告 第 48 号

図４ 計算対象（移行拡散解析）∂C∂ 2C∂C+ vk= Dij∂t∂xk∂xi ∂x j

図６ 計算対象（変形解析）（１）

これを 3.2 に示した技法を織り交ぜて有限要素法に基づきここでの計算は時間ステップ１回分を示しているが、物質プログラミングする。離散化した行列演算に対して、CUDA

移行拡散解析は解析メッシュの最小幅や時間刻みに対での並列演算と CPU での単一演算との計算時間を比較すして、クーラン数あるいはペクレ数で表される指標の制約る。計算対象は図４のような単純な濃度拡散問題を設定すを受ける。
例えば実際に物質移行拡散解析を地下水の汚染拡散る。
得られた計算時間の比較を図５に示す。の予測に適用する場合などでは、空間メッシュは工場がここでは、離散化について詳しく記述はしないが、式（１）立地されるような規模の土地面積であり、時間的には戦を空間的には有限要素法、時間に関しては前進差分を採後から 60 年オーダを対象とするようなこともある。このよう用している。また、時間変化に対応する左辺第１項に影響な場合には、単一ステップの計算時間は短くても、時系する質量マトリックスに関しては対角化(lumping)を施してい列な計算が不可能ではないにしても事実上相当な作業る。これにより、各行が独立するので完全な並列演算を実負荷を発生する可能性もある。これを、ここで示した計算施できることになる。

時間のスピードアップによって、作業負荷が飛躍的に軽結果は図 5 に示すように、CPU での演算時間(20ms)に比減することとなる。これにより、対策工法も含めた十分な検べて 1/500 の結果(0.04ms)となった。また、計算時間の節点討資料を用意できるので、並列演算の効用は極めて意

数依存性は予想した成果を示す図４そのものとなっていて、義が深いものと考えられる。

節点数に依存しない一定値であり、CUDA による並列化の4.2 変形解析への適用効用は十分に発揮されていると思われる。

変形解析の支配方程式は下式で与えられている。

25

計算時間(ms)20∂u k ∂ 
=0Cijkl∂x j ∂xl lump(CPU)lump(GPU)15(2)有限要素法で離散化した連立１次方程式を種々の解法で計算した 5)６）。変形解析における計算対象は図６に示す単10純引張問題とした。
得られた計算時間の比較を図７に示す。この問題は対角優位性 (kii ≥ ∑ kij ) がないため、収束性が悪い。これは5i≠ j

CUDA を適用したことによる結果ではなく、そもそも解法が0持つ性質であることに注意する。つまり、一般に言われてい0100020003000

るように、対角優位性をもたない係数行列の場合、並列化

節点数の恩恵はあまり受けないと認識するのが正しい理解と思わ

図５ 計算時間の比較（移行拡散解析）れる。なお、ここでの掃き出し法では、一旦逆行列を計算し、その後の行列積に並列化を施している。

－34－

CUDA を用いた並列数値解析手法の一考察

8.0

6.0
5.0

掃き出し法（CPU)Jacobi法（CPU)
CG法（CPU)
掃き出し法（GPU)Jacobi法（GPU)
CG法（GPU)

1.6
計算時間[×E05]（ms)

計算時間[×E06]（ms)1.8

掃き出し法（CPU)Jacobi法（CPU)
CG法（CPU)
掃き出し法（GPU)Jacobi法（GPU)
CG法（GPU)

7.0

4.0
3.0
2.0

1.4
1.2
1CPU

0.8
0.6
0.4

GPU

0.2

1.0

0

0.000500

1000

1500 2000
節点数

2500

500

1000

3000

1500 2000
節点数

2500

3000

図９ 計算時間の比較（ラプラス方程式）

図７ 計算時間の比較（変形解析）§５．まとめ
4.3 ラプラス方程式への適用
一般に言われるように、対角優位性を持った係数行列の本文での知見を以下にまとめる。

場合、間接解法の収束性は良好であるとすれば CUDA に

１） 自明ではあるが、同じアルゴリズムであれば CUDA による並列計算の恩恵を得ることができる。このことを確認す

よって並列化を実装した方が計算時間は短い。

るために Laplace 方程式への適用を試みた。
Laplace 方程式は次式である。

２） 並列演算の効用が大きく期待できる問題は、例えば∂  ∂φ 
=0∂xi  ∂xi (3)係数行列を対角化して各行を独立に扱えるような問題である。

これを有限要素法で離散化し、得られた係数行列を種々のこの条件に適合する計算対象は、時間発展の問題を解法で計算した。計算対象を図８に示す。

前進差分によって時間積分した離散化である。これ

なお、蛇足ながら、図８における Direcret BC とはディリクは、移行拡散解析に CUDA を適用したことで得られた

レ境界条件であって、目的変数 φ 自体の値を境界条件とし知見である。て与えるものである。
これらから得られた計算時間は図９のようになった。
成果の予想イメージ図３のような傾向を示している。

３） 対角項以外にも非ゼロ数が要素として入る一般の係数行列を間接手法で計算する場合には、対角優位性

これより、移行拡散方程式ほどではないにしろ対角優位性をもった係数行列の場合には CUDA の恩恵を得ることがのように、扱う係数行列の性質によって収束性の良否が場合分けされる。

可能であることが確認できた。なお、ここでも掃き出し法は、一旦逆行列を計算し、その後の行列積に並列化を施してい対角優位性とは、係数行列の対角項以外の要素の総計るのは変形解析と同様である。

より対角項の要素が大きいという性質である。本文に示した

Direcret B.C.

ように、ラプラス方程式などの対角優位性を持つ問題では、

Direcret B.C.

計算手続きも簡単で収束性は良好であるため、CUDA の適用が効果的である。しかしながら、変形解析のように、対角優位性を持たない係数行列の場合には、収束性が悪い。このような場合には、バンドマトリックス法やスカイライン法などによる直接法の方が効率的である可能性がある。

図８ 計算対象（ラプラス方程式）－35－

フジタ技術研究報告 第 48 号

参考文献
1)

青木尊之、額田彰：はじめての CUDA プログラミング、工学社

2)

小山田耕二、岡田賢治：CUDA 高速 GPU プログラミング、秀和システム

3)

株式会社クイープ訳：CUDA BY EXAMPLE、インプレスジャパン

4)David B. Kirk, Wen-mei W. Hwu:CUDA：CUDA プログ
ラミング実践講座、株式会社ボーンデジタル

5)戸川隼人：マトリクスの数値計算、オーム社

6)

戸川隼人：共役勾配法、教育出版

ひ と こ とグラフィック機能を使って数値演算を並列化する GP-GPU は、比較的安価にハードウェアを装備することが可能であるし、その効果も問題によっては大きい。今後は、プリ処理も含めた効率化を進めることで、作業の合理化をシステムとして進展させることができ仲沢武志れば良いと考えている。

－36－