疎 実対称正定値 線形方程式 Ax = b 解 ます プ ン引数 使用し 関連す 計算 行うこ が ます こ に A ン ック分解 A 数値的分解 ン ック 又 数値的分解
い が与え たAx = b 解 計算があ ます
梗概
#include <imsl.h>
float *imsl_f_lin_sol_posdef_coordinate (int n, int nz, Imsl_f_sparse_elem *a, float *b, ..., 0) void imsl_free_symbolic_factor (Imsl_symbolic_factor *sym_factor)
void imsl_f_free_numeric_factor (Imsl_f_numeric_factor *num_factor) double型 関数 imsl_d_lin_sol_posdef_coordinate
imsl_d_free_numeric_factor
必須 引数
int n 。入力) 行列 行数
int nz 。入力)
行列 下 角にあ 非 ロ 数
Imsl_f_sparse_elem *a 。入力)
行列 下 角にあ 非 ロ入力値 位置 値 含 長さ nz ベク float *b 。入力)
右辺 含 長さn ベク
戻 値
疎対称正定値線形方程式 Ax = b 解 x へ インタ こ 領域 解放す に imsl_free 使
用します 解が得 い場合 NULL が返さ ます
プ ン引数 梗概
#include <imsl.h>
float *imsl_f_lin_sol_posdef_coordinate (int n, int nz, Imsl_f_sparse_elem *a, float *b,
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor, IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor, IMSL_SYMBOLIC_FACTOR_ONLY,
IMSL_RETURN_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor, IMSL_SUPPLY_NUMERIC_FACTOR, Imsl_f_numeric_factor *num_factor, IMSL_NUMERIC_FACTOR_ONLY,
IMSL_SOLVE_ONLY,
IMSL_MULTIFRONTAL_FACTORIZATION, IMSL_RETURN_USER, float x[],
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *small_element, IMSL_LARGEST_DIAGONAL_ELEMENT, float *largest_element, IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, float *values, 0)
プ ン引数
IMSL_RETURN_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor 。出力)
戻 入力行列 ン ック分解 含 Imsl_symbolic_factor型 構造体へ イン
タ
Imsl_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 フロンタ 行列 タック た に必要 作
業用領域 大 さ チフロンタ 因子分 解が用い い場合 こ 値 ロに ッ さ
こ 構造体内 アロ ー さ た 解放す に 関数 imsl_free_symbolic_factor 使うこ
IMSL_SUPPLY_SYMBOLIC_FACTOR, Imsl_symbolic_factor *sym_factor 。入力) Imsl_symbolic_factor型 構造体へ インタ
こ 構造体 IMSL_RETURN_SYMBOLIC_FACTOR プ ン
imsl_f_lin_sol_posdef_coordinate によ 計算さ 入力行列 ン ック 分解 含 ます こ 構造体 IMSL_RETURN_SYMBOLIC_FACTOR プ ン引数に説 明があ ます こ 構造体内 アロ ー さ た 解放す に 関数
imsl_free_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構造体 詳細 次 通
パ ータ ータ型 説明
nzsub int ** Cholesky因子 非 ロ 非対角要素に対す 行指標
含 配列へ インタ こ 配列 長さ nz 割 当 が 配列 全 要素 使うこ い
xnzsub int ** nzsub に対す 指標 含 長さ n + 1 配列へ
インタ Cholesky因子 列 j にあ 非 ロ 要
素に対す 行指標 nzsub[xnzsub[j]]
連続領域に格納さ
xlnz int ** 長さ n + 1 配列へ インタ こ 配列
alnz 非 ロ 非対角要素 抽出す に使
わ 始 終わ 指標 含 因子行列 列 j
に対し alnz 始 終わ 指標 そ
xlnz[j] xlnz[j + 1] に格納さ
alnz float ** Cholesky因子 非 ロ 非対角要素 含
perm int ** 置換ベク 含 長さ n 配列へ インタ
diag float ** Cholesky因子 対角要素 含 長さ n 配列へ インタ
L A コ キー因子 し num_nonzeros L 非 ロ 数 します に述べた構
造体 L 対角要素 diag に格納さ ます L 非 ロ 非対角要素 alnzに格納さ
ます 列 j に対し alnz L 非 ロ要素 抽出す に使わ 始 終わ
指標 そ xlnz[j] xlnz[j+1] に格納さ ます L 非 ロ 要素 列イン
ック nzsub に含ま ます xnzsub[i] 列 i に対し L 行イン ック 抽出
す 始 nzsub イン ック 含 ます こ 次 コー 説明さ ます ここ
構造体 成分 因子行列L 下 角 再構築す や 方 示します：
Imsl_f_numeric_factor numfctr;
. . .
for (i = 0; i < n; i++){
L[i][i] = (*numfctr.diag)[i];
if ((*numfctr.xlnz)[i] > (num_nonzeros-n)) continue;
start = (*numfctr.xlnz)[i]-1;
stop = (*numfctr.xlnz)[i+1]-1;
k = (*numfctr.xnzsub)[i]-1;
for (j = start; j < stop; j++){
L[(*numfctr.nzsub)[k]-1][i] = (*numfctr.alnz)[j];
k++;
} }
こ 構造体内 アロ ー さ た 解放す に 関数 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_lin_sol_posdef_coordinate によ 計算さ た入力行列 数値的分解 含 います こ 構造体 プ ン引数 IMSL_RETURN_NUMERIC_FACTOR 説明さ
います
IMSL_NUMERIC_FACTOR_ONLY,
入力行列 数値的分解 計算し 戻 ます 引数b 無視さ ます IMSL_SOLVE_ONLY,
A 数値的 又 ン ック分解が与え た Ax = b 解 ます こ プ ン IMSL_SUPPLY_NUMERIC_FACTOR 又 IMSL_SUPPLY_SYMBOLIC_FACTOR 使 用 必要 します
IMSL_MULTIFRONTAL_FACTORIZATION,
チフロンタ 法 使用し 数値的分解 します フ 疎 縮格納方式に 基 い 標準 分解 します
IMSL_RETURN_USER, float x[] 。出力) 解 x 含 長さ n ー 割 当 配列
IMSL_SMALLEST_DIAGONAL_ELEMENT, float *small_element 。出力)
数値的分解 発生す 最小 対角要素 含 へ インタ こ プ ン 数値的分解が imsl_f_lin_sol_posdef_coordinate 呼び出し 計算さ 場合
有効 す
IMSL_LARGEST_DIAGONAL_ELEMENT, float *large_element 。出力)
数値的分解 発生す 最大 対角要素 含 へ インタ こ プ ン 数値的分解が imsl_f_lin_sol_posdef_coordinate 呼び出し 計算さ 場合
有効 す
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros 。出力) 分解因子 非 ロ 総数 含 へ インタ
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, float *values 。入力)
縮した疎 列形式。CSC) 係数行列 使います こ 格納方式 説明 イン ロダク ン 参照し く さい
説明
関数 imsl_f_lin_sol_posdef_coordinate 疎 対称正定値係数行列 A 持 線形方程式 解 ます こ 関数 フ 使用法 最初に係数行列 置換 ン ック分解が計算さ そ 数値的分解が計算さ ます 線形方程式 解 こ 数値的分解因子 使用し 求 こ が ます
ン ック分解 ップ 最小次数 並び え 決定 コ キー因子L 疎 ータ構造体 設 定 ます こ ップ 疎 係数行列 パターン す わ 非 ロ要素 位置 け 必要
し 要素そ も 扱いませ 従 Imsl_f_sparse_elem 構造体 val 欄 無視さ ます
あ アプ ー ンが 同 疎 パターン 持 要素 値が異 複数 疎 対称正定値係数行列 扱う場合 IMSL_RETURN_SYMBOLIC_FACTOR IMSL_SUPPLY_SYMBOLIC_FACTOR 使用す こ
によ ン ック分解 一度 け計算す よい す
ン ック分解因子によ コ キー因子L 疎 ータ構造体が与え 数値的分解 次式 が成 立 ように Lに ン 値 生成します：
PAP T = LL T ここ P 最小次数 並び えによ 決定さ 置換行列 す
数値的分解 2 方法 い 実行す こ が ます フ 疎 縮格納方式に基
い 標準 分解 行 います こ 方法 George Liu (1981年) に記述さ います プ ン
チフロンタ 法 使用す こ が ます チフロンタ 法 よ 多く 必要 し ますが 場合によ に よ 高速に ます チフロンタ 法 Liu (1987年) ーチンに基 い
います こ 方法 詳細 説明 Liu 1990年 よび Duff Reid 1983年 1984年 Ashcraft
1987年 Ashcraftそ 他 1987年 Liu 1986年 1989年 参照し く さい
アプ ー ンが 同 係数行列 右辺が異 いく 線形方程式 解く必要があ 場合 プ ン IMSL_RETURN_NUMERIC_FACTOR IMSL_SUPPLY_NUMERIC_FACTOR 使 コ キー分
解因子 前処理 す こ が ます そ IMSL_SOLVE_ONLY プ ン 使い 全 方程式
効率良く解くこ が ます
数値分解が与え そ 解 x 次によ 得 ます：
L y1
= Pb
LT y2= y
1x = PT y2
置換情報 P 数値因子構造体に保持さ ます
例題 1
例 し 次 5 × 5 係数行列 考えます
50 5 0 3 2
5 40 4 0 0
0 4 30 0 1
3 0 0 20 0
2 0 1 0 10 a
Ax = (55, 83, 103, 97, 82)T に ように xT = (5, 4, 3, 2, 1) します A 下 角 非 ロ 数 nz
= 10 す 下 角 疎 標対応形式 次 ように ます：
行 0 1 2 2 3 3 4 4 4 4 列 0 1 0 2 2 3 0 1 3 4 値 10 20 1 30 4 40 2 3 5 50 こ 表現 唯一 く 次も同等 す：
行 3 4 4 4 0 1 2 2 3 4 列 3 0 1 3 0 1 0 2 2 4 値 40 2 3 5 10 20 1 30 4 50
#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;
int nz = 10;
float *x;
x = imsl_f_lin_sol_posdef_coordinate (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 A = E(2500, 50) します 方程式 Ax = b1 解 そこ 数値的分解 返しま
す 計算さ たこ 数値的分解 使用し 方程式 Ax = b2 解 ます また 実行時間 比率がプ ン さ ます こ 時間 測定結果 ンに強く依存す こ に注意し 下さい
#include <imsl.h>
#include <stdio.h>
int main() {
Imsl_f_sparse_elem *a;
Imsl_f_numeric_factor numeric_factor;
float *b_1;
float *b_2;
float *x_1;
float *x_2;
int n;
int ic;
int nz;
double time_1;
double time_2;
ic = 50;
n = ic*ic;
/* Generate two right hand sides */
b_1 = imsl_f_random_uniform (n*sizeof(*b_1), 0);
b_2 = imsl_f_random_uniform (n*sizeof(*b_2), 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 = imsl_ctime ();
x_1 = imsl_f_lin_sol_posdef_coordinate (n, nz, a, b_1, IMSL_RETURN_NUMERIC_FACTOR, &numeric_factor,
0);
time_1 = imsl_ctime () - time_1;
/* Now solve Ax_2 = b_2 given the numeric factorization */
time_2 = imsl_ctime ();
x_2 = imsl_f_lin_sol_posdef_coordinate (n, nz, a, b_2, IMSL_SUPPLY_NUMERIC_FACTOR, &numeric_factor,
IMSL_SOLVE_ONLY, 0);
time_2 = imsl_ctime () - time_2;
printf("time_2/time_1 = %lf\n", time_2/time_1);
}
出力
time_2/time_1 = 0.037037