CULA Sparse (S5, CUDA 5.5) - 超音波流体屋のプログラム備忘録

CULA Sparse はCUDA上で動作する疎行列用の線形システム反復解法ツール。
おまじないが20行くらいあって非常に煩雑であるため、
疎行列でも1000x1000くらいかそれ以下ならDense使ったほうがラク。

アカデミックならフリーで使える。ダウンロードは以下からID登録してから。
http://www.culatools.com/downloads/sparse/
linuxは落としてきた.runを管理者権限でbashで実行。

使用例: 一次元ポアソン方程式を解く

CULA Sparse 使い方の一例として、疎行列の代表例であるラプラシアンを含む微分方程式を解く。
行列の導出はこのページの最下部を見てもらうとして、こんな形状の行列を解く。
結果のyは上に凸の放物線になる。
$\left[ \begin{array}{cccccc} 2 & -1 & & & & O \\ -1 & 2 & -1 & & & \\ & -1 & 2 & & & \\ & & & \ddots && \\ & & & & 2 & -1 \\ O & & & & -1& 2 \end{array} \right] \{y\} = b_0 \left\{ \begin{array}{c} 1 \\ 1 \\ 1 \\ \vdots \\ 1 \\ 1\end{array} \right\} $
端的に言えば対角項だけ2、対角項の上下左右に-1、残りは0の行列である。

行列の行数をNとして、Nに依らない同じ結果のカーブにするには
$b_0 = \frac{8}{(N+1)^2}$とすると解は$y = 1 - x^2$になる。

スパース行列のデータ形式

こんな感じで非ゼロの行番号と列番号、値を並べてく形式をCoordinate Format (COO)と呼びます。
その他の形式はculatoolsにて。
普通の配列で問題ありませんが、非ゼロの数が予めわからない問題などは下記ソースのようにvectorクラス使うと楽です。

$int\ row[N_{nz}],\ col[N_{nz}];\ double\ A[N_{nz}];$ ($N_{nz}$は非ゼロ要素の数)として、
$row = \left\{ \begin{array}{c} 0\\ 0\\ 1\\ 1 \\ 1\\ \vdots \\ N-1 \\ N-1 \end{array} \right\}, col = \left\{ \begin{array}{c} 0\\ 1\\ 0\\ 1\\ 2 \\ \vdots\\ N-2 \\ N-1 \end{array} \right\},A = \left\{ \begin{array}{c} 2\\ -1\\ -1\\ 2\\-1 \\ \vdots \\ -1 \\ 2 \end{array} \right\}$

CSRデータ形式 (2015-07-02 追記)

CULA SparseはCOOを与えられた場合、解く前にCSRに変換して処理する。
仕様メモリギリギリの大規模行列ではこの変換作業でメモリを消費し、結果計算が失敗する場合がある。
計算時間やメモリの節約なため、予めCSR形式で作成したほうが良い。

CSR(compressed sparse row)は上記の$col, A$はCOOと同一だが、
$row$は前行までの非ゼロ要素の数を記録する。
領域は$N+1$行分用意し、最終要素は非ゼロ要素の数を加えるのに注意する。

今回の場合では$int\ row[N+1];$として、
$row =\left\{ \begin{array}{c} 0\\ 2\\ 5\\ 8 \\ 11 \\ \vdots \\ 3N-2\ (=N_{nz}) \end{array} \right\}$
のようになる。

ソース

ソルバや前処理法の選び方はCULA Toolsのフローチャートを参照されたし。
対称行列(symmetric)、正定値(positive definite)、対角優位(diagonally dominant)、帯行列(banded)

#include <cuda_runtime.h>
#include <helper_functions.h>
#include <helper_cuda.h>
#include <cula_sparse.h>
 
#define USE_CUDA 1 //実行プロセッサの切り替え (CUDA/Host)
#define USE_JACOBI 1  //前処理法の切り替え (ヤコビ/不完全LU分解)
 
int main(int argc, char** argv)
{
int devID = findCudaDevice(argc, (const char **)argv);
FILE *fout;
if ((fout = fopen("result.dat", "w")) == NULL){
  printf("file open error!!\n"); exit(EXIT_FAILURE);
}
 
//Ay = b を解くスパース行列を準備する, COO型の格納方式とする
//非ゼロ要素数の扱いが煩雑なので vector を使うのが良い
const int N = 99999;
int nNonZero;
std::vector<double> A;
std::vector<int> rowIdx;
std::vector<int> colIdx;
 
//一次元ラプラシアンを作る d2y/dx2 = -C, 定義域は[-x0, x0]
double C = 2.0;
double x0 = 1.0;
double dx = 2.0*x0/(N+1.0);
double b0 = C*dx*dx;
 
//自身が2, 隣接が-1の行列
for(int ri=0; ri<N; ri++){
  A.push_back(2.0); rowIdx.push_back(ri); colIdx.push_back(ri);
  if(ri > 0) { A.push_back(-1.0); rowIdx.push_back(ri); colIdx.push_back(ri-1);}
  if(ri < N -1){ A.push_back(-1.0); rowIdx.push_back(ri); colIdx.push_back(ri+1);}
}
nNonZero = A.size();
 
 
//解格納用ベクトル
std::vector<double> y(N);
std::vector<double> b(N, b0);
 
/////////////////////////////////////////////////////////////
//CULA SPARSE 設定
/////////////////////////////////////////////////////////////
culaSparseHandle handle; //sparse本体
culaSparseStatus status; //各関数の成否格納用
culaSparsePlan plan;     //sparse解法などの設定
culaSparseConfig config; //反復解法設定用の構造体
culaSparseResult result; //解以外の実行結果情報
 
// cula_sparseが出力するメッセージの格納用
const int msgsize = 512;
char csMsg[msgsize];
 
//初期化
culaSparseCreate(&handle);
culaSparseCreatePlan(handle, &plan);
 
//行列Aが座標指定型(COO式)スパース行列であることを指定する
culaSparseSetDcooData(handle, plan, 0, N, nNonZero, &A[0], &rowIdx[0], &colIdx[0], &y[0], &b[0]);
 
//反復解法の設定
culaSparseConfigInit(handle, &config);
config.relativeTolerance = 1e-6;
config.maxIterations = 1e5;
 
culaSparseGetConfigString(handle, &config, csMsg, msgsize);
printf("%s\n", csMsg); //設定内容の表示
 
if(USE_CUDA){
  //Hostメモリにある行列AをCUDAで解く。
  culaSparseCudaOptions platformOpts; //解法を解くデバイス指定
  culaSparseCudaOptionsInit(handle, &platformOpts);
  culaSparseSetCudaPlatform(handle, plan, &platformOpts);
}else{
  //Hostメモリにある行列AをCPUで解く。
  culaSparseHostOptions platformOpts; //解法を解くデバイス指定
  culaSparseHostOptionsInit(handle, &platformOpts);
  culaSparseSetHostPlatform(handle, plan, &platformOpts);
}
 
//反復解法ソルバーの指定 ラプラシアンは正定値なので共役勾配法が使える
culaSparseSetCgSolver(handle, plan, 0);
 
//前処理方法の指定 フローチャートではヤコビ前処理を使えとのこと　
if(USE_JACOBI) culaSparseSetJacobiPreconditioner(handle, plan, 0);
else culaSparseSetIlu0Preconditioner(handle, plan, 0); //不完全LU分解前処理
 
/////////////////////////////////////////////////////////////
//連立方程式を解く
/////////////////////////////////////////////////////////////
//実行・タイマー測定
StopWatchInterface *timer = 0;
sdkCreateTimer(&timer);
sdkStartTimer(&timer);
 
//連立方程式を解く
status =  culaSparseExecutePlan(handle, plan, &config, &result);
printf("Status: %d, NoError: %d\n", status, culaSparseNoError);
 
//タイマーストップ
sdkStopTimer(&timer);
printf("計算時間 =%f(ms)\n", sdkGetTimerValue(&timer));
 
//結果表示
culaSparseGetResultString(handle, &result, csMsg, msgsize);
printf("%s\n", csMsg);
for(int i=0;i<N;i++){
  fprintf(fout, "%f %f\n", (i+1)*dx-x0, y[i]);
  if(i == N/2+1) printf("%f %f\n", i*dx-x0, y[i]); //中心のy,1.0なら成功
}
/////////////////////////////////////////////////////////////
//終了処理
/////////////////////////////////////////////////////////////
fclose(fout);
sdkDeleteTimer(&timer);
culaSparseDestroyPlan(plan);
culaSparseDestroy(handle);
}
 

環境変数

.bashrcあたりに以下のSparseの記述を突っ込んどく

export CUDA_ROOT=/usr/local/cuda-5.5
export CULA_ROOT=/usr/local/cula
export CULA_INC_PATH=$CULA_ROOT/include
export CULA_LIB_PATH_64=$CULA_ROOT/lib64
export CULASPARSE_ROOT=/usr/local/culasparse
export CULASPARSE_INC_PATH=$CULASPARSE_ROOT/include
export CULASPARSE_LIB_PATH_64=$CULASPARSE_ROOT/lib64
export CUDA_INCPATH=$CUDA_ROOT/include:$CUDA_ROOT/samples/common/inc:$CULA_INC_PATH:$CULASPARSE_INC_PATH
 
export PATH=$CUDA_ROOT/bin:$PATH
export LD_LIBRARY_PATH=$CUDA_ROOT/lib64:$CULA_LIB_PATH_64:$CULASPARSE_LIB_PATH_64:$LD_LIBRARY_PATH
export LIBRARY_PATH=$LD_LIBRARY_PATH
export C_INCLUDE_PATH=$CUDA_INCPATH:$C_INCLUDE_PATH
export CPLUS_INCLUDE_PATH=$CUDA_INCPATH:$CPLUS_INCLUDE_PATH
 

コンパイル

Nsightが便利なので、そっちを使いましょう。
CULAのライブラリの登録だけ必要。プロジェクト右クリックで、
Properties →　Build → Settings
Tool Settingsタブ → NVCC Linker → Libraries
Addのアイコン → でてきたダイアログに cula_sparse

Makefileならこんな感じ。

CC=nvcc
ARCH=-arch=sm_13
LIB=-lcula_sparse
SRC=sample.cu
OUT=sample.out
all:
	${CC} ${ARCH} ${SRC} -o ${OUT} ${LIB} 
 

結果

N=99999, relativeTolerance=1e-6で、67行のconfig.maxIterationsをあれこれ変えた結果。
許容誤差の設定次第ではあるが、結局のとこ行数程度のイタレーションがいる模様。

計算時間

行数Nをあれこれ変えた計算時間結果。
今回は不完全LU分解前処理が特異的にうまく行く問題だったみたい。
Dense関係の線が途切れてるのは、メモリのオーバーフローで計算できなかったもの。
CPU使用率からみて、Host下のCULA Sparseは4コア、Octaveは1コアでの実行である模様。

CPU: Core i7-3770, 3.40 GHz (8 Core), Memory 16 GB
GPU: Nvidia Tesla K20, 706 MHz (2496 Core), GDDR 6 GB

結論・補足

$N < 10^3$ Denseを使う。
$N > 10^3$ Sparseを使う。その際小さいNでどのソルバ・前処理がいいか試してみること。

初期値が使えるのが反復法の良さ。コンフィグで以下を入れて、
yに初期値を入れとくと初期値次第で計算が数倍速く終わります。

config.useInitialResultVector = 1

過渡解析的な使い方をする時にどうぞ。

スパース行列をDeviceメモリ上で確保して72-82行のプラットホームをculaSparseCudaDeviceOptionsにすれば、
CULA同様デバイス上メモリに対しても計算できます。

付録1: 検証用 Matlab/Octave用プログラム

C = 2.0; x0 = 1.0;
N=99999;  dx = 2*x0/(N+1); b0 = C*dx*dx;
 
I=eye(N);
%I=speye(N); %スパース単位行列を使うと上のCG-ILU0-Hostくらい速い
A=2*I - circshift(I, [1 0]) - circshift(I, [-1 0]);
A(1,N)=0; A(N,1) = 0;
b=b0*ones(N,1);
 
tic
y=A\b;
toc
 

付録2: ポアソン方程式から線形システムまで

$\nabla^2 y = -C$
今回は単純のため一次元、$x=\pm x_0$で0となる境界条件とする。
$\frac{d^2y}{dx^2} = -C$　解は $ y = \frac{C}{2} ( x_0^2 - x^2) $

k番目のy周りで、二階微分を離散化すると
$\frac{y(k+1)+y(k-1)-2y(k)}{\Delta x^2} = -C$
整理して、区間分割数をN+1、つまり$\Delta x = \frac{2x_0}{N+1} $とすると、
$2y(k) - y(k+1) -y(k-1) = C \left( \frac{2x_0}{N+1} \right)^2 $
これを行列の形にして$C=2, x_0=1$とすると上の行列になります。

参考

CULA Sparse Reference Manual
(pdf) CULA Sparse Manual

正直上のリファレンスは手抜きな部分が多いので
cula_sparse.h などヘッダを読むか、
/usr/local/culasparse/example のサンプルコードを読む、実行してみるのが早い