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のマニュアルページを参照)1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
このデータをダウンロード |
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
- 1行目はタイトル行で読み飛ばされます。
- 2行目は係数行列の大きさ(n=5, 5x5なので5)を与えます。 数字より後ろは行末まで読み飛ばされます。
- 3行目では係数行列(スパース)に含まれる非ゼロ要素の数(nnz=16)を与えます。数字より後ろは行末まで読み飛ばされます。
- 4行目は用いる計算手法 (method=Nag_SparseNsym_RGMRES) とプリコンディショナー (precon=Nag_SparseNsym_SSORPrec) を指定しています。 プリコンディショナーの指定より後ろは行末まで読み飛ばされます。
※ここで他の手法(Nag_SparseNsym_CGS もしくは Nag_SparseNsym_BiCGSTAB)やプリコンディショナー (Nag_SparseNsym_JacPrec もしくは Nag_SparseNsym_NoPrec) の指定も可能です。 - 5行目ではSSORプリコンディショナーのリラックスパラメータであるω (omega=1.05) を与えています。SSOR以外のプリコンディショナーを用いる場合にはこの数値は無視されます。 数字より後ろは行末まで読み飛ばされます。
- 6行目では RGMRES法の再開サブスペース空間の次元(m=1) 、計算精度(tol=1.e-10) 、最大反復数(maxitn=1000) をそれぞれ与えています。 3つの数値より後ろは行末まで読み飛ばされます。
- 7行目~22行目にはスパース行列の非ゼロ係数の値を指定しています。それぞれの行は値(a)、行番号(irow)、列番号(icol)で全部で16個(3行目で与えた値)の非ゼロ係数が指定されています。行番号及び列番号は1から最大5(2行目で指定した値)までの値を取ります。 例えば行列の一番左の一番上の値が3.5であった場合には
3.5 1 1
と与えます。 - 23行目と24行目は右辺ベクトルの値(b)を指定しています。全部で5個(2行目で指定した値)の数値が指定されています。
- 25行目と26行目は解の初期値(x)を与えています。 全部で5個(2行目で指定した値)の数値が指定されています。 ここでは初期値としてすべて 0 を与えています。
出力結果
(本関数の詳細はnag_sparse_nsym_solのマニュアルページを参照)1 2 3 4 5 6 7 8 9
この出力例をダウンロード |
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
- 1行目はタイトルです
- 2行目に解を求める際の反復数が示されます。(この例では13回)
- 3行目は求まった解の残差ノルムが出力されます。
- 4行目は解のタイトルである x が出力されます。
- 5行目~9行目に求まった解が表示されます。
ソースコード
(本関数の詳細はnag_sparse_nsym_solのマニュアルページを参照)
※本サンプルソースコードはnAG数値計算ライブラリ(Windows, Linux, MAC等に対応)の関数を呼び出します。
サンプルのコンパイル及び実行方法
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
このソースコードをダウンロード |
/* 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; }