superLU using

直接法のソルバー
複素値にも対応

公式サイト
http://crd.lbl.gov/~xiaoye/SuperLU/
(英語版マニュアル有り)

とにかく使う事だけできるようになる覚書．

前提

• makefileでコンパイル．
• 行列はCRS形式で記憶．
• SuperLUディレクトリは"/home/tlab"においてある．

解くべき方程式

nn:総節点数，
A：全体行列 (nn×nn)，
b：右辺ベクトル（nn×1）
として
Ax = b
を解く．

CRS構造の中身

として，CRS構造を持つ全体行列 A は

• A.nnz //非零要素数
• A.val //非零要素値ベクトル（length:nnz）
• A.row_ptr //各行のもっとも左にある非零要素がA.valにおける何番目の要素かを収めるベクトル(length:nn+1)
• A.col_ind //各非零要素の列番号を収めたベクトル(length:nnz)

からなる．
（当たり前ですが，配列は0から始まる前提です．）

使い方

複素行列の場合

"slu_zdefs.h"をインクルード．
"SuperLU/SRC/slu_dcomplex.h"を参考にして，A.val と b を"doublecomplex"型で宣言，行列生成．
(参考：doublecomplex型の宣言)

typedef struct { 
     double r, i; 
 } doublecomplex;
 

ソースコードの書き方

実行列	⇔	複素行列
slu_ddefs.h		slu_zdefs.h

をインクルードしておく．
行列 A (A.nnz, A.val, A.col_ind, A.row_ptr)と右辺ベクトル b が準備してある状態で
以下のハイライト部分を貼り付ければOK．

/* 1. 必要な変数の宣言*/
	SuperMatrix P,L,U,B;
	int nrhs=1, info, permc_spec;
	superlu_options_t options;
	SuperLUStat_t stat;
 	int      *perm_r; /* row permutations from partial pivoting */
 	int      *perm_c; /* column permutation vector */
	if ( !(perm_r = intMalloc(nn)) ) ABORT("Malloc fails for perm_r[].");
	if ( !(perm_c = intMalloc(nn)) ) ABORT("Malloc fails for perm_c[].");
 
/* 2. オプションの設定*/
	set_default_options(&options);
	options.ColPerm = NATURAL;
 
/* 3. ここは正直なにをしてるのはわかりません．*/
	StatInit(&stat);
 
/* 4. SuperLu用に改めて全体行列 P と右辺ベクトル B を生成*/
	zCreate_CompRow_Matrix(&P, nn, nn, A.nnz, A.val, A.col_ind, A.row_ptr, SLU_NR, SLU_Z, SLU_GE);
	zCreate_Dense_Matrix(&B, nn, nrhs, b, nn, SLU_DN, SLU_Z, SLU_GE);
 
/* 5. 求解の実施*/
	zgssv(&options, &P, perm_c, perm_r, &L, &U, &B, &stat, &info);
 
/* 6. LU分解した際の行列 L,U の廃棄*/
	Destroy_SuperNode_Matrix(&L);
	Destroy_CompCol_Matrix(&U);
 
/* 7. 全体行列 P の廃棄と何かの配列解法*/
	Destroy_CompRow_Matrix(&P);
	free(perm_r);
	free(perm_c);
 
 

• 1. 必要な変数の宣言

ここはとくに何もなし．そのままでOK.

• 2．オプションの設定

普通に解くならそのままでOK．前処理など使いたいときはマニュアル．

• 3. ここは正直なにをしてるのはわかりません．

文字通り．

• 4. SuperLu用に改めて全体行列 P と右辺ベクトル B を生成
• 5. 求解の実施

先に作った A と b を"SuperMatrix"という型にはめる操作．
行列用メモリを2重に取っているわけではない．
ちなみに，実行列と複素行列の場合は以下の部分を変える．

実行列	⇔	複素行列
SLU_D		SLU_Z
dgssv		zgssv

ベクトル b に解が入ってるので出力．

• 6. LU分解した際の行列 L,U の廃棄．

タイムステップごとに連立方程式を解いてる場合などはその都度実施しないとメモリがいっぱいになる．（上書きしてくれない．）

• 7. 全体行列 P の廃棄と何かの配列解放．

この"Destroy_***"はCRSセットの配列(A.val,A.col_ind,A.row_ptr)も解放するので注意．

makefileの書き方

"main.c"で方程式を作り，サブルーチン"superlu.c"で求解を行う"sample"という実行ファイルを生成するときの例

CC = gcc
SLU = home/tlab/Superlu_4.2    (自分のバージョンで.ここでは ver.4.2)
 
sample: main.c superlu.o
	$(CC) -o sample -L /$(SLU)/lib -I /$(SLU)/SRC main.c　superlu.o /$(SLU)/lib/libsuperlu_4.2.a /$(SLU)/lib/libblas.a -lm
superlu.o: superlu.c
	$(CC) -c -I /$(SLU)/SRC superlu.c -lm
 

• "-L /$(SLU)/lib"

SuperLUライブラリファイルの場所指定(ライブラリファイルの検索場所追加)
静的ライブラリなので，先頭に書く．

• "-I /$(SLU)/SRC"

SuperLU用ヘッダーファイルの置き場指定(ヘッダーファイルの検索場所追加)
"slu_ddefs.h"や"slu_zdefs.h"をインクルードしているファイルのコンパイルには全部つける．
（"superlu.c"はインクルードしているので，"superlu.o"の生成にこのオプションが必要．）

• "/$(SLU)/lib/libsuperlu_4.2.a"
• "/$(SLU)/lib/libblas.a"

フルパスで指定してるけど，これらはフルパス書かなくてもよさそうな気もする．

• "-lm"

linuxでコンパイルするとき用のオプション．
cygwinなら不要．