Keyword: スパースソルバー, sparse solver, RGMRES, 疎行列, SSOR
概要
本サンプルは以下に示す非対称実スパース連立方程式を解くC言語によるプログラムです。 連立方程式は係数行列と右辺ベクトルとして与えます。
以下のプログラム例ではRGMRES法 (restarted generalized minimal residual method) を用いていますが、その他にCGS法 (conjugate gradient squared method)、及びBICGSTAB法 (bi-conjugate gradient stabilised method) に対応しています。（入力データのNag_SparseNsym_RGMRESの部分をNag_SparseNsym_CGS、Nag_SparseNsym_BiCGSTABにそれぞれ置き換える事で手法を選ぶことが可能です）
またプリコンディショナーにはSSOR (symmetric successive-over-relaxation) を用いていますが、その他にJacobi の指定及びプリコンディショナーを使用しない設定が可能です。（入力データのNag_SparseNsym_SSORPrecの部分をそれぞれNag_SparseNsym_JacPrec、Nag_SparseNsym_NoPrecに置き換えることで指定可能です）
※本サンプルはnAG Cライブラリに含まれる関数nag_sparse_nsym_sol()のExampleコードです。本サンプル及び関数の詳細情報はnag_sparse_nsym_solのマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_sparse_nsym_solのマニュアルページを参照）| このデータをダウンロード |
nag_sparse_nsym_sol (f11dec) Example Program Data 5 n 16 nnz Nag_SparseNsym_RGMRES Nag_SparseNsym_SSORPrec method, precon 1.05 omega 1 1.e-10 1000 m, tol, maxitn 2. 1 1 1. 1 2 -1. 1 4 -3. 2 2 -2. 2 3 1. 2 5 1. 3 1 5. 3 3 3. 3 4 1. 3 5 -2. 4 1 -3. 4 4 -1. 4 5 4. 5 2 -2. 5 3 -6. 5 5 a[i-1], irow[i-1], icol[i-1], i=1,...,nnz 0. -7. 33. -19. -28. b[i-1], i=1,...,n 0. 0. 0. 0. 0. x[i-1], i=1,...,n
- １行目はタイトル行で読み飛ばされます。
- ２行目は係数行列の大きさ（n=5, 5x5なので5)を与えます。 数字より後ろは行末まで読み飛ばされます。
- ３行目では係数行列（スパース）に含まれる非ゼロ要素の数(nnz=16)を与えます。数字より後ろは行末まで読み飛ばされます。
- ４行目は用いる計算手法 (method=Nag_SparseNsym_RGMRES) とプリコンディショナー (precon=Nag_SparseNsym_SSORPrec) を指定しています。 プリコンディショナーの指定より後ろは行末まで読み飛ばされます。
※ここで他の手法（Nag_SparseNsym_CGS もしくは Nag_SparseNsym_BiCGSTAB）やプリコンディショナー （Nag_SparseNsym_JacPrec もしくは Nag_SparseNsym_NoPrec） の指定も可能です。 - ５行目ではSSORプリコンディショナーのリラックスパラメータであるω (omega=1.05) を与えています。SSOR以外のプリコンディショナーを用いる場合にはこの数値は無視されます。 数字より後ろは行末まで読み飛ばされます。
- ６行目では RGMRES法の再開サブスペース空間の次元(m=1) 、計算精度(tol=1.e-10) 、最大反復数(maxitn=1000) をそれぞれ与えています。 ３つの数値より後ろは行末まで読み飛ばされます。
- ７行目〜２２行目にはスパース行列の非ゼロ係数の値を指定しています。それぞれの行は値(a)、行番号(irow)、列番号(icol)で全部で16個（３行目で与えた値）の非ゼロ係数が指定されています。行番号及び列番号は1から最大5（２行目で指定した値）までの値を取ります。 例えば行列の一番左の一番上の値が3.5であった場合には
3.5 1 1
と与えます。 - ２３行目と２４行目は右辺ベクトルの値(b)を指定しています。全部で５個（２行目で指定した値）の数値が指定されています。
- ２５行目と２６行目は解の初期値(x)を与えています。 全部で５個（２行目で指定した値）の数値が指定されています。 ここでは初期値としてすべて 0 を与えています。
出力結果
（本関数の詳細はnag_sparse_nsym_solのマニュアルページを参照）| この出力例をダウンロード |
nag_sparse_nsym_sol (f11dec) Example Program Results
Converged in 13 iterations
Final residual norm = 5.087e-09
x
1.000000e+00
2.000000e+00
3.000000e+00
4.000000e+00
5.000000e+00
- １行目はタイトルです
- ２行目に解を求める際の反復数が示されます。（この例では１３回）
- ３行目は求まった解の残差ノルムが出力されます。
- ４行目は解のタイトルである x が出力されます。
- ５行目〜９行目に求まった解が表示されます。
ソースコード
（本関数の詳細はnag_sparse_nsym_solのマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_sparse_nsym_sol (f11dec) Example Program.
*
* CLL6I261D/CLL6I261DL Version.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.1, 2017.
*/
#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nag_string.h>
#include <nagf11.h>
int main(void)
{
double *a = 0, *b = 0, *x = 0;
double omega;
double rnorm;
double tol;
Integer exit_status = 0;
Integer *icol = 0, *irow = 0;
Integer i, m, n;
Integer maxitn, itn;
Integer nnz;
char nag_enum_arg[40];
Nag_SparseNsym_Method method;
Nag_SparseNsym_PrecType precon;
Nag_Sparse_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_sparse_nsym_sol (f11dec) Example Program Results\n");
/* Skip heading in data file */
scanf("%*[^\n]");
scanf("%ld%*[^\n]", &n);
scanf("%ld%*[^\n]", &nnz);
scanf("%39s", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
method = (Nag_SparseNsym_Method) nag_enum_name_to_value(nag_enum_arg);
scanf("%39s%*[^\n]", nag_enum_arg);
precon = (Nag_SparseNsym_PrecType) nag_enum_name_to_value(nag_enum_arg);
scanf("%lf%*[^\n]", &omega);
scanf("%ld%lf%ld%*[^\n]", &m, &tol, &maxitn);
x = nAG_ALLOC(n, double);
b = nAG_ALLOC(n, double);
a = nAG_ALLOC(nnz, double);
irow = nAG_ALLOC(nnz, Integer);
icol = nAG_ALLOC(nnz, Integer);
if (!irow || !icol || !a || !x || !b) {
printf("Allocation failure\n");
exit_status = 1;
goto END;
}
/* Read the matrix a */
for (i = 1; i <= nnz; ++i)
scanf("%lf%ld%ld%*[^\n]", &a[i - 1], &irow[i - 1],
&icol[i - 1]);
/* Read right-hand side vector b and initial approximate solution x */
for (i = 1; i <= n; ++i)
scanf("%lf", &b[i - 1]);
scanf("%*[^\n]");
for (i = 1; i <= n; ++i)
scanf("%lf", &x[i - 1]);
scanf("%*[^\n]");
/* Solve Ax = b using nag_sparse_nsym_sol (f11dec) */
/* nag_sparse_nsym_sol (f11dec).
* Solver with no Jacobi/SSOR preconditioning (nonsymmetric)
*/
nag_sparse_nsym_sol(method, precon, n, nnz, a, irow, icol, omega, b, m, tol,
maxitn, x, &rnorm, &itn, &comm, &fail);
printf("%s%10ld%s\n", "Converged in", itn, " iterations");
printf("%s%16.3e\n", "Final residual norm =", rnorm);
/* Output x */
printf(" x\n");
for (i = 1; i <= n; ++i)
printf(" %16.6e\n", x[i - 1]);
END:
nAG_FREE(irow);
nAG_FREE(icol);
nAG_FREE(a);
nAG_FREE(x);
nAG_FREE(b);
return exit_status;
}