非対称実スパースソルバー (Nonsymmetric Sparse Solver)

C言語によるサンプルソースコード : 使用関数名:nag_sparse_nsym_sol (f11dec)

ホーム > 製品 > nAG数値計算ライブラリ > サンプルソースコード集 > 非対称実スパースソルバー (C言語/C++)

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;
}


関連情報
© 譌・譛ャ繝九Η繝シ繝。繝ェ繧ォ繝ォ繧「繝ォ繧エ繝ェ繧コ繝�繧コ繧ー繝ォ繝シ繝玲�ェ蠑丈シ夂、セ 2024
Privacy Policy  /  Trademarks