Keyword: スパース, 不完全コレスキー分解
概要
本サンプルは実スパース対称行列の不完全コレスキー分解を行うC言語によるサンプルプログラムです。 本サンプルは実スパース対称行列Aを読み込み、不完全コレスキー分解を計算し、行列Aと以下に示される下三角行列Cの非ゼロ要素を出力します。
※本サンプルはnAG Cライブラリに含まれる関数 nag_sparse_sym_chol_fac() のExampleコードです。本サンプル及び関数の詳細情報は nag_sparse_sym_chol_fac のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照）| このデータをダウンロード |
nag_sparse_sym_chol_fac (f11jac) Example Program Data 7 n 16 nnz 0 0.0 lfill, dtol Nag_SparseSym_UnModFact 0.0 mic, dscale Nag_SparseSym_MarkPiv pstrat 4. 1 1 1. 2 1 5. 2 2 2. 3 3 2. 4 2 3. 4 4 -1. 5 1 1. 5 4 4. 5 5 1. 6 2 -2. 6 5 3. 6 6 2. 7 1 -1. 7 2 -2. 7 3 5. 7 7 a[i-1], irow[i-1], icol[i-1], i=1,...,nnz
- １行目はタイトル行で読み飛ばされます。
- ２行目に行列Aの次数(n)を指定しています。
- ３行目に行列Aの下三角部分の非ゼロ要素の数(nnz)を指定しています。
- ４行目の１番目のパラメータ(lfill)は値が0以上の場合、コレスキー分解で許されたフィル（fill）の最大レベルを指定しています。２番目のパラメータ(dtol)は、１番目のパラメータの値が負の値の場合、フィルイン（fill-in）を制御するためのドロップ許容値を指定しています。１番目のパラメータの値が0以上の場合は２番目のパラメータの値は参照されません。
- ５行目には行の合計を保持するためにコレスキー分解が修正されるかどうかを示すフラグ(mic)と対角スケーリング引数(dscale)を指定しています。 "Nag_SparseSym_UnModFact"はコレスキー分解が修正されないことを意味します。
- ６行目にピボット法(pstrat)を指定しています。"Nag_SparseSym_MarkPiv"はMarkowitz法を用いて対角ピボットが実行されることを意味します。
- ７〜２２行目に行列Aの非ゼロ要素(a)、行列Aの非ゼロ要素の行インデックス(irow)、列インデックス(icol)を指定しています。
出力結果
（本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照）| この出力例をダウンロード |
nag_sparse_sym_chol_fac (f11jac) Example Program Results
Original Matrix
n = 7
nnz = 16
1 4.0000e+00 1 1
2 1.0000e+00 2 1
3 5.0000e+00 2 2
4 2.0000e+00 3 3
5 2.0000e+00 4 2
6 3.0000e+00 4 4
7 -1.0000e+00 5 1
8 1.0000e+00 5 4
9 4.0000e+00 5 5
10 1.0000e+00 6 2
11 -2.0000e+00 6 5
12 3.0000e+00 6 6
13 2.0000e+00 7 1
14 -1.0000e+00 7 2
15 -2.0000e+00 7 3
16 5.0000e+00 7 7
Factorization
n = 7
nnz = 16
npivm = 0
17 5.0000e-01 1 1
18 3.3333e-01 2 2
19 3.3333e-01 3 2
20 2.7273e-01 3 3
21 -5.4545e-01 4 3
22 5.2381e-01 4 4
23 -2.7273e-01 5 3
24 2.6829e-01 5 5
25 6.6667e-01 6 2
26 5.2381e-01 6 4
27 2.6829e-01 6 5
28 3.4788e-01 6 6
29 -1.0000e+00 7 1
30 5.3659e-01 7 5
31 -5.3455e-01 7 6
32 9.0461e-01 7 7
i ipiv(i)
1 3
2 4
3 5
4 6
5 1
6 2
7 7
- ３〜２２行目には元の行列の情報が出力されています。
- ４行目に行列Aの次数が出力されています。
- ５行目に行列Aの非ゼロ要素の数が出力されています。
- ７〜２２行目に行列Aの要素の値、行インデックスと列インデックスが出力されています。
- ２４〜４４行目にはコレスキー分解して得られた行列の情報が出力されています。
- ２５行目にコレスキー分解して得られた行列の次数が出力されています。
- ２６行目にコレスキー分解して得られた行列の非ゼロ要素の数が出力されています。
- ２７行目にコレスキー分解で修正されたピボットの数が出力されています。
- ２９〜４４行目にコレスキー分解して得られた下三角行列Cの要素の値、行インデックスと列インデックスが出力されています。
- ４７〜５３行目にピボットインデックスが出力されています。
ソースコード
（本関数の詳細はnag_sparse_sym_chol_fac のマニュアルページを参照）
※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
/* nag_sparse_sym_chol_fac (f11jac) 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 dtol;
double *a;
double dscale;
Integer *irow, *icol;
Integer *ipiv, nnzc, *istr;
Integer exit_status = 0, i, n, lfill, npivm;
Integer nnz;
Integer num;
char nag_enum_arg[40];
Nag_SparseSym_Piv pstrat;
Nag_SparseSym_Fact mic;
Nag_Sparse_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_sparse_sym_chol_fac (f11jac) Example Program Results\n\n");
/* Skip heading in data file */
scanf("%*[^\n]");
/* Read algorithmic parameters */
scanf("%ld", &n);
scanf("%*[^\n]");
scanf("%ld%*[^\n]", &nnz);
scanf("%ld%lf%*[^\n]", &lfill, &dtol);
scanf("%39s%lf%*[^\n]", nag_enum_arg, &dscale);
/* nag_enum_name_to_value (x04nac).
* Converts nAG enum member name to value
*/
mic = (Nag_SparseSym_Fact) nag_enum_name_to_value(nag_enum_arg);
scanf("%39s%*[^\n]", nag_enum_arg);
pstrat = (Nag_SparseSym_Piv) nag_enum_name_to_value(nag_enum_arg);
/* Allocate memory */
num = 2 * nnz;
ipiv = nAG_ALLOC(n, Integer);
istr = nAG_ALLOC(n + 1, Integer);
irow = nAG_ALLOC(num, Integer);
icol = nAG_ALLOC(num, Integer);
a = nAG_ALLOC(num, double);
if (!ipiv || !istr || !irow || !icol || !a) {
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]);
/* Calculate incomplete Cholesky factorization */
/* nag_sparse_sym_chol_fac (f11jac).
* Incomplete Cholesky factorization (symmetric)
*/
nag_sparse_sym_chol_fac(n, nnz, &a, &num, &irow, &icol, lfill, dtol, mic,
dscale, pstrat, ipiv, istr, &nnzc, &npivm, &comm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_sparse_sym_chol_fac (f11jac).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Output original matrix */
printf(" Original Matrix \n");
printf(" n = %6ld\n", n);
printf(" nnz = %6ld\n\n", nnz);
for (i = 1; i <= nnz; ++i)
printf(" %8ld%16.4e%8ld%8ld\n", i, a[i - 1],
irow[i - 1], icol[i - 1]);
printf("\n");
/* Output details of the factorization */
printf(" Factorization\n n = %6ld \n nnz = %6ld\n", n,
nnzc);
printf(" npivm = %6ld\n\n", npivm);
for (i = nnz + 1; i <= nnz + nnzc; ++i)
printf(" %8ld%16.4e%8ld%8ld\n", i, a[i - 1],
irow[i - 1], icol[i - 1]);
printf("\n i ipiv(i) \n");
for (i = 1; i <= n; ++i)
printf(" %8ld%8ld\n", i, ipiv[i - 1]);
END:
nAG_FREE(irow);
nAG_FREE(icol);
nAG_FREE(a);
nAG_FREE(istr);
nAG_FREE(ipiv);
return exit_status;
}