プ ン引数
IMSL_CSC_FORMAT, int col_ptr[], int row_ind[], float values[] (入力)
縮さ た疎 列格納方式 (CSC) 係数行列 入力 す こ 方式に い 本書 イン ロダク ン 参照
IMSL_PREORDERING, int preorder (入力)
行列 A 事前並べ え 用い 最小次数並べ え(MDO: Minimum Degree Ordering)
ア バ ー ン：
preorder 方法
0 George Liu 商最小次数 (Quotient Minimum Degree) ア
1 George Liu 商最小次数 変形 前処理 外部次数 用い
フ : preorder = 0
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (出力)
タイプ Imsl_snodal_symbolic_factor 構造体へ インタ 戻 入力行列 ーパ
ー ーダ ン ック因子分解 含
Imsl_snodal_symbolic_factor 構造体 詳細 次 通 ：
表 1.7 構造体 Imsl_snodal_symbolic_factor
パ ータ ータ型 説明
nzsub int ** Cholesky 因子 非 ロ非対角要素 縮さ た行指標 含 配列へ インタ
xnzsub int ** *nzsub に対す イン ック 含 長さ n+1 配列
へ インタ Cholesky 因子 列 j にあ 非 ロに対
す 行指標 (*nzsub)[(*xnzsub)[j]] 順に
保存さ
maxsub int 指標 し 用い 配列 *nzsub 要素数
*nzsub 大 さ maxsub よ 大 いこ もあ こ に注意
xlnz int ** 長さ n+1 配列 配列 *alnz 非 ロ非対角要
素 出す に用い 始 終わ イン ック
含 (alnz に い プ ン引数
IMSL_RETURN_NUMERIC_FACTOR 説明 参照). 因子
行列 列 j に対し *alnz 始 終わ イン
ック そ (*xlnz) [j] (*xlnz) [j+1] に
格納さ
maxlnz int Cholesky 因子 非 ロ非対角要素 数
perm int ** 置換ベク 含 長さ n 配列へ インタ
invp int ** 逆置換ベク 含 長さ n 配列へ インタ
multifrontal
_space int こ 変数 用い い
nsuper int Cholesky 因子 ーパー ー 数
snode int ** 長さ n 配列へ インタ 要素 (*snode)[j]
列 j が所属す 基本 ーパー ー 数 含 snode_ptr int ** 長さ nsuper+1 配列へ インタ 各 ーパー
ー 開始列 含
nleaves int 対称的に置換さ た入力行列 A 後順 削除木 葉
数
etree_leaves int ** 削除木 葉 数 含 長さ nleaves+1 配列へ インタ
こ 構造体 領域 解放す に imsl_free_snodal_symbolic_factor 使う IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (入力)
Imsl_snodal_symbolic_factor 型 構造体へ インタ こ 構造体
IMSL_RETURN_SYMBOLIC_FACTOR プ ン imsl_f_sparse_cholesky_smp によ 計算さ た入力行列 ン ック因子分解 含 こ 構造体
IMSL_RETURN_SYMBOLIC_FACTOR プ ン引数 説明さ こ 構造体に割 当 た領域 解放す に 関数 imsl_free_snodal_symbolic_factor 使う IMSL_SYMBOLIC_FACTOR_ONLY (入力)
入力行列 ン ック因子分解 計算し戻 引数 b 無視さ IMSL_RETURN_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor (出力)
Imsl_f_numeric_factor 型 構造体へ インタ 戻 入力行列 数値的因子分解
含 Imsl_f_numeric_factor 構造体 詳細 関数
imsl_f_lin_sol_posdef_coordinate プ ン
IMSL_RETURN_NUMERIC_FACTOR 説明にあ ます こ 構造体に割 当 た領域 解放す に 関数 imsl_f_free_numeric_factor 使う
IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor (入力)
タイプ Imsl_f_numeric_factor 構造体へ インタ こ 構造体
IMSL_RETURN_NUMERIC_FACTOR プ ン imsl_f_sparse_cholesky_smp 計算さ た入力行列 数値的因子分解 含 こ 構造体 説明 関数
imsl_f_lin_sol_posdef_coordinate IMSL_RETURN_NUMERIC_FACTOR プ ン 引数 説明 参照 こ 構造体に割 当 た領域 解放す に 関数 imsl_f_free_numeric_factor 使う
IMSL_NUMERIC_FACTOR_ONLY (入力)
入力行列 数値的因子分解 計算す 引数 b 無視さ IMSL_SOLVE_ONLY (入力)
Ax = b A 数値的あ い ン ック因子分解 用い 解く こ プ ン
IMSL_SUPPLY_NUMERIC_FACTOR あ い IMSL_SUPPLY_SYMBOLIC_FACTOR 使 うこ 要求します
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *smallest_element (出力)
数値的因子分解 最小 対角要素 含 へ インタ こ プ ン
imsl_f_sparse_cholesky_smp へ こ 呼び出し 数値的因子分解が計算さ 場合 け有効 す
IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element (出力)
数値的因子分解 最大 対角要素 含 へ インタ こ プ ン
imsl_f_sparse_cholesky_smp へ こ 呼び出し 数値的因子分解が計算さ 場合 け有効 す
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros (出力) 因子にあ 非 ロ 総数 含 へ インタ
IMSL_RETURN_USER, float x[] (出力)
解 x 含 ー 割 当 長さ n 配列
説明
関数 imsl_f_sparse_cholesky_smp 疎 対称正定値係数行列 A 持 線型方程式 解 ます こ 関数 フ 使い方 最初に係数行列 ン ック因子分解 計算します 次に OpenMP 並 列 数値的因子分解 行います そ 数値的因子 用い 線型方程式 解 求 ます
ン ック因子分解 最小次数並べ え 決定し Cholesky 因子 L 疎 ーパー ーダ ータ 構造体 設定します こ ップ 疎 係数行列 “パターン” け 必要 します す わ 非
ロ要素 位置 け 要素そ も 使いませ 従 Imsl_f_sparse_elem 構造体 val 項
無視さ ます アプ ー ンが同 疎 パターン 持 疎 対称正定値係数行列 いく 生成す 場 合 IMSL_RETURN_SYMBOLIC_FACTOR IMSL_SUPPLY_SYMBOLIC_FACTOR 用い ン ッ ク因子分解 一度 けや よい す
Cholesky 因子 L に対す 疎 ータ構造が ン ック因子 し 与え 数値的因子分解によ
PAPT = LLT
L が求ま ます ここ P 最小次数並べ え 決定さ 置換行列 す
数値的因子分解 並列 ーパー ーダ 解法 使います こ 左方 よび右方列計算 キー 組 合わ
せます こ ア Rauber et al. (1999) に説明があ ます
アプ ー ンが 同 係数行列 右辺が異 いく 線型方程式 解く場合 プ ン IMSL_RETURN_NUMERIC_FACTOR IMSL_SUPPLY_NUMERIC_FACTOR 用い Cholesky 因子 事前 に計算します そ IMSL_SOLVE_ONLY プ ン すべ 線型方程式 効率よく解くこ が ます
数値的因子分解が与え 解 x 次 計算 求ま ます：
Ly1 = Pb LT y2 = y1
x = PT y2
置換 情報 P 数値的因子構造体 Imsl_f_numeric_factor にあ ます
例題 1
5 × 5 係数行列 A
します A 下 角 非 ロ 数 nz = 10 す b := Ax = (55, 83, 103, 97, 82)T し 線型方 程式 Ax = b 解 xT = (5, 4, 3, 2, 1) 得ます
#include <imsl.h>
int main() {
Imsl_f_sparse_elem a[] = {0, 0, 10.0,
1, 1, 20.0, 2, 0, 1.0,
2, 2, 30.0, 3, 2, 4.0, 3, 3, 40.0, 4, 0, 2.0, 4, 1, 3.0, 4, 3, 5.0, 4, 4, 50.0};
float b[] = {55.0, 83.0, 103.0, 97.0, 82.0};
int n = 5, nz = 10;
float *x = NULL;
x = imsl_f_sparse_cholesky_smp (n, nz, a, b, 0);
imsl_f_write_matrix ("solution", 1, n, x, 0);
imsl_free (x);
}
出力
solution
1 2 3 4 5 5 4 3 2 1
例題 2
ここ ン ック因子 再利用 示します 最初に A = E(2500, 50) 線型方程式 Ax = b 解 A ン ック因子 得ます 次に C y = d, C = A+2*I 線型方程式 そ ン ック因子 用い 解 ます ここ I 恒等行列 す こ A C が同 非 ロ構造 持 同 ン ッ ク因子 使うこ が こ 利用し います
#include <imsl.h>
#include <stdlib.h>
#include <stdio.h>
int main() {
Imsl_f_sparse_elem *a = NULL, *c = NULL;
Imsl_snodal_symbolic_factor symbolic_factor;
float *b = NULL, *d = NULL, *x = NULL, *y = NULL;
float *mod_vector = NULL;
int n, ic, nz, i, index;
float error_1, error_2;
ic = 50;
n = ic * ic;
mod_vector = (float*) malloc (n * sizeof(float));
/* Build coefficient matrix A */
a = (Imsl_f_sparse_elem *) imsl_f_generate_test_coordinate (n, ic, &nz,
IMSL_SYMMETRIC_STORAGE, 0);
/* Build coefficient matrix C */
c = (Imsl_f_sparse_elem*) malloc (nz * sizeof(Imsl_f_sparse_elem));
for (i = 0; i < nz; i++) c[i] = a[i];
for (i = 0; i < n; i++) c[i].val = 6.0;
/* Form right hand side b */
for (i = 0; i < n; i++)
mod_vector[i] = (float) (i % 5);
b = (float *) imsl_f_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_vector, IMSL_SYMMETRIC_STORAGE, 0);
/* Form right hand side d */
d = (float *) imsl_f_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, c,
IMSL_X_VECTOR, n, mod_vector, IMSL_SYMMETRIC_STORAGE, 0);
/* Solve Ax = b and return the symbolic factorization */
x = imsl_f_sparse_cholesky_smp (n, nz, a, b, IMSL_RETURN_SYMBOLIC_FACTOR, &symbolic_factor, 0);
/* Compute solution error |x - mod_vector| */
error_1 = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_vector, IMSL_INF_NORM, &index,
0);
/* Solve Cy = d given the symbolic factorization */
y = imsl_f_sparse_cholesky_smp (n, nz, c, d, IMSL_SUPPLY_SYMBOLIC_FACTOR, &symbolic_factor, 0);
/* Compute solution error |y - mod_vector| */
error_2 = imsl_f_vector_norm (n, y, IMSL_SECOND_VECTOR, mod_vector, IMSL_INF_NORM, &index,
0);
printf ("Solution error |x - mod_vector| = %e\n", error_1);
printf ("Solution error |y - mod_vector| = %e\n", error_2);
/* Free allocated memory */
if (b) imsl_free(b);
if (d) imsl_free(d);
if (x) imsl_free(x);
if (y) imsl_free(y);
if (mod_vector) free(mod_vector);
if (a) imsl_free(a);
if (c) free(c);
imsl_free_snodal_symbolic_factor (&symbolic_factor);
}
出力
Solution error |x - mod_vector| = 4.529953e-005 Solution error |y - mod_vector| = 2.861023e-006
例題 3
ここ A = E(2500, 50) し 最初に Ax = b1 解 数値的因子 得ます 次にそ 数値的因子 用い
Ax = b2 解 ます 実行時間 プ ン します 時間 計算機に大 く依存す こ に注意し く
さい
#include <imsl.h>
#include <stdio.h>
#include <omp.h>
int main() {
int n, ic, nz;
float *b_1 = NULL, *b_2 = NULL, *x_1 = NULL, *x_2 = NULL;
double time_1, time_2;
Imsl_f_sparse_elem *a = NULL;
Imsl_f_numeric_factor numeric_factor;
ic = 50;
n = ic * ic;
/* Generate two right hand sides */
imsl_random_seed_set (1234567);
b_1 = imsl_f_random_uniform (n, 0);
b_2 = imsl_f_random_uniform (n, 0);
/* Build coefficient matrix a */
a = imsl_f_generate_test_coordinate (n, ic, &nz, IMSL_SYMMETRIC_STORAGE,
0);
/* Now solve Ax_1 = b_1 and return the numeric factorization */
time_1 = omp_get_wtime();
x_1 = imsl_f_sparse_cholesky_smp (n, nz, a, b_1, IMSL_RETURN_NUMERIC_FACTOR, &numeric_factor, 0);
time_1 = omp_get_wtime() - time_1;
/* Now solve Ax_2 = b_2 given the numeric factorization */
time_2 = omp_get_wtime();
x_2 = imsl_f_sparse_cholesky_smp (n, nz, a, b_2, IMSL_SUPPLY_NUMERIC_FACTOR, &numeric_factor, IMSL_SOLVE_ONLY,
0);
time_2 = omp_get_wtime() - time_2;
printf("time_2/time_1 = %lf\n", time_2/time_1);
/* Free allocated memory */
if (x_1) imsl_free(x_1);
if (x_2) imsl_free(x_2);
if (b_1) imsl_free(b_1);
if (b_2) imsl_free(b_2);
if (a) imsl_free(a);
imsl_f_free_numeric_factor (&numeric_factor);
}
出力
time_2/time_1 = 0.029411
重大 ー
IMSL_BAD_SQUARE_ROOT ロあ い 負 平方根が因子分解 発生した
係数行列 正定値 い
sparse_cholesky_smp 複素数
疎 ー 正定値行列 A Cholesky 分解 OpenMP 並列 ーパー ーダ 法 計算し 線型方程 式 Ax = b 解く
梗概
#include <imsl.h>
f_complex *imsl_f_sparse_cholesky_smp (int n, int nz, Imsl_c_sparse_elem a[],f_complex b[],
…, 0)
void imsl_free_snodal_symbolic_factor (Imsl_snodal_symbolic_factor *sym_factor) void imsl_c_free_numeric_factor (Imsl_c_numeric_factor *num_factor)
double 型 関数 imsl_z_sparse_cholesky_smp imsl_z_free_numeric_factor
必須 引数
int n (入力) 入力行列 次数 int nz (入力)
行列 下 角にあ 非 ロ 数 Imsl_c_sparse_elem a[] (入力)
長さ nz 配列 行列 下 角 各非 ロ ン 位置 値 含 f_complex b[] (入力)
長さ n 配列 右辺 含
戻 値
疎 ー 正定値線型方程式 Ax = b 解 x インタ こ 領域 解放す に imsl_free
使う 解が得 い場合 NULL が返さ
プ ン引数 梗概
#include <imsl.h>
f_complex *imsl_c_sparse_cholesky_smp (int n, int nz, Imsl_c_sparse_elem a[],f_complex b[], IMSL_CSC_FORMAT, int col_ptr[], int row_ind[],f_complex values[],
IMSL_PREORDERING, int preorder,
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor, IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor, IMSL_SYMBOLIC_FACTOR_ONLY,
IMSL_RETURN_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor, IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor, IMSL_NUMERIC_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *smallest_element, IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element, IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,
IMSL_RETURN_USER, f_complex x[],
0)
プ ン引数
IMSL_CSC_FORMAT, int col_ptr[], int row_ind[],f_complex values[] (入力)
縮さ た疎 列格納方式 (CSC) 係数行列 入力 す こ 方式に い 本書 イン ロダク ン 参照
IMSL_PREORDERING, int preorder (入力)
行列 A 事前並べ え 用い 最小次数並べ え (MDO : Minimum Degree Ordering)
ア バ ー ン：
preorder 方法
0 George Liu 商最小次数 (Quotient Minimum Degree) ア
1 George Liu 商最小次数 変形 前処理 外部次数 用い
フ : preorder = 0
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (出力)
タイプ Imsl_snodal_symbolic_factor 構造体へ インタ 戻 入力行列 ーパ
ー ーダ ン ック因子分解 含
Imsl_snodal_symbolic_factor 構造体 詳細 次 通 ：
表 1.8 構造体 Imsl_snodal_symbolic_factor
パ ータ ータ型 説明
nzsub int ** Cholesky 因子 非 ロ非対角要素 縮さ た行指標 含 配列へ インタ
xnzsub int ** *nzsub に対す イン ック 含 長さ n+1 配列
へ インタ Cholesky 因子 列 j にあ 非 ロに対
す 行指標 (*nzsub)[(*xnzsub)[j]] 順に
保存さ
maxsub int 指標 し 用い 配列 *nzsub 要素数
*nzsub 大 さ maxsub よ 大 いこ もあ こ に注意
xlnz int ** 長さ n+1 配列 配列 *alnz 非 ロ非対角要
素 出す に用い 始 終わ イン ック
含 (alnz に い プ ン引数
IMSL_RETURN_NUMERIC_FACTOR 説明 参照). 因子
行列 列 j に対し *alnz 始 終わ イン
ック そ (*xlnz) [j] (*xlnz) [j+1] に
格納さ
maxlnz int Cholesky 因子 非 ロ非対角要素 数
perm int ** 置換ベク 含 長さ n 配列へ インタ
invp int ** 逆置換ベク 含 長さ n 配列へ インタ
multifrontal
_space int こ 変数 用い い
nsuper int Cholesky 因子 ーパー ー 数
snode int ** 長さ n 配列へ インタ 要素 (*snode)[j]
列 j が所属す 基本 ーパー ー 数 含 snode_ptr int ** 長さ nsuper+1 配列へ インタ 各 ーパー
ー 開始列 含
nleaves int 対称的に置換さ た入力行列 A 後順 削除木 葉
数
etree_leaves int ** 削除木 葉 数 含 長さ nleaves+1 配列へ インタ
こ 構造体 領域 解放す に imsl_free_snodal_symbolic_factor 使う IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_snodal_symbolic_factor *sym_factor (入力)
Imsl_snodal_symbolic_factor 型 構造体へ インタ こ 構造体
IMSL_RETURN_SYMBOLIC_FACTOR プ ン imsl_c_sparse_cholesky_smp によ 計算さ た入力行列 ン ック因子分解 含 こ 構造体
IMSL_RETURN_SYMBOLIC_FACTOR プ ン引数 説明さ こ 構造体に割 当 た領域 解放す に 関数 imsl_free_snodal_symbolic_factor 使う IMSL_SYMBOLIC_FACTOR_ONLY (入力)
入力行列 ン ック因子分解 計算し戻 引数 b 無視さ IMSL_RETURN_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor (出力)
Imsl_c_numeric_factor 型 構造体へ インタ 戻 入力行列 数値的因子分解
含 Imsl_c_numeric_factor 構造体 詳細 関数
imsl_c_lin_sol_posdef_coordinate プ ン
IMSL_RETURN_NUMERIC_FACTOR 説明にあ ます こ 構造体に割 当 た領域 解放す に 関数 imsl_c_free_numeric_factor 使う
IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_c_numeric_factor *num_factor (入力)
タイプ Imsl_c_numeric_factor 構造体へ インタ こ 構造体
IMSL_RETURN_NUMERIC_FACTOR プ ン imsl_c_sparse_cholesky_smp 計算さ た入力行列 数値的因子分解 含 こ 構造体 説明 関数
imsl_c_lin_sol_posdef_coordinate IMSL_RETURN_NUMERIC_FACTOR プ ン 引数 説明 参照 こ 構造体に割 当 た領域 解放す に 関数 imsl_c_free_numeric_factor 使う
IMSL_NUMERIC_FACTOR_ONLY (入力)
入力行列 数値的因子分解 計算す 引数 b 無視さ IMSL_SOLVE_ONLY (入力)
Ax = b A 数値的あ い ン ック因子分解 用い 解く こ プ ン
IMSL_SUPPLY_NUMERIC_FACTOR あ い IMSL_SUPPLY_SYMBOLIC_FACTOR 使 うこ 要求します
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *smallest_element (出力)
数値的因子分解 最小 対角要素 含 へ インタ こ プ ン
imsl_c_sparse_cholesky_smp へ こ 呼び出し 数値的因子分解が計算さ 場合 け有効 す
IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element (出力)
数値的因子分解 最大 対角要素 含 へ インタ こ プ ン
imsl_c_sparse_cholesky_smp へ こ 呼び出し 数値的因子分解が計算さ 場合 け有効 す
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros (出力) 因子にあ 非 ロ 総数 含 へ インタ