IMSL_SMALLEST_PIVOT, float *small_pivot
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros,
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, f_complex *values, IMSL_MEMORY_BLOCK_SIZE, int block_size,
0)
プ ン引数
IMSL_RETURN_SPARSE_LU_FACTOR,Imsl_c_sparse_lu_factor *lu_factor (出力)
Imsl_c_sparse_lu_factor型 構造体 ア 構造体内 インタ
imsl_c_lin_sol_gen_coordinate によ LU 分解 指すように初期化さ ます IMSL_SUPPLY_SPARSE_LU_FACTOR,Imsl_c_sparse_lu_factor *lu_factor (入力)
Imsl_c_sparse_lu_factor型 構造体 ア
こ 構造体 IMSL_RETURN_SPARSE_LU_FACTOR プ ン
imsl_c_lin_sol_gen_coordinate によ 計算さ た入力行列 LU分解 含 ま す
IMSL_FREE_SPARSE_LU_FACTOR,
戻 前に A LU 分解 含 ータ構造体 ンク 解放します こ プ ン 使用す こ 因子が必要 い場合 け す
IMSL_RETURN_SPARSE_LU_IN_COORD, Imsl_c_sparse_elem **lu_coordinate, int **row_pivots, int **col_pivots (出力)
LU 分解が lu_coordinate に長さ nz 標対応形式 返さ ます こ
Imsl_c_sparse_lu_factor に プ 化さ た内部表現よ コンパク す 点
SOLVE_ONLY 呼び出し 間 こ 因子 内部表現が再構成さ け いこ
す し し 因子がプログ 終了後に格納さ 再度 後 実行時にロー さ も あ IMSL_RETURN_LU_IN_COORD IMSL_SUPPLY_LU_IN_COORD 組 合わせ 格納や 込 が単純 フ ー ッ 最良 選択に ます IMSL_SUPPLY_SPARSE_LU_IN_COORD, int nzlu, Imsl_c_sparse_elem *lu_coordinate,
int *row_pivots, int *col_pivots (入力) 標対応形式 LU 分解 提供します
こ 説明 IMSL_RETURN_SPARSE_LU_IN_COORD 参照し く さい IMSL_FACTOR_ONLY,
入力行列 LU 分解 計算します 引数 b 無視さ ます IMSL_SOLVE_ONLY,
A LU 分解 用い Ax = b 解 ます こ プ ン プ ン引数
IMSL_SUPPLY_SPARSE_LU_FACTOR 又 IMSL_SUPPLY_SPARSE_LU_IN_COORD 使用 す 必要があ ます
IMSL_RETURN_USER, f_complex x[] (出力) 解 x 含 長さ n ー 割 当 配列 IMSL_TRANSPOSE
,
問題 ATx = b 解 ます こ プ ン 分解 行うい プ ン 組 合わせ
使用す こ が ます IMSL_CONDITION, float *condition,
A L1 条件数 推定し 変数 condition に返します IMSL_PIVOTING_STRATEGY, Imsl_pivot method (入力)
以下 う 一 設定方法によ ッ 選択 選択します IMSL_ROW_MARKOWITZ IMSL_COLUMN_MARKOWITZ 又 IMSL_SYMMETRIC_MARKOWITZ
フ ： IMSL_SYMMETRIC_MARKOWITZ
IMSL_NUM_OF_SEARCH_ROWS
,
int num_search_row(入力)
ッ 要素 し 検索さ 非 ロ 要素 最小数 持 行数 フ ： num_search_row = 3
IMSL_ITERATIVE_REFINEMENT,
反復改良法が必要 場合 こ プ ン 選択します
IMSL_DROP_TOLERANCE, float tolerance (入力)
可能 fill-in がこ 許容値に対し チ ックさ ます 新しい要素 絶対値が
tolerance よ 小さけ そ 捨 ます
フ ： tolerance = 0.0
IMSL_HYBRID_FACTORIZATION, float density, int order_bound, (入力)
アク 小行列 密度が 0.0 ≤ density ≤ 1.0 に アク 小行列 次数
が order_bound 以下 場合 密 分解法に え ます
IMSL_GROWTH_FACTOR_LIMIT, float gf_limit (入力)
成長因子が gf_limit 越え 場合 計算 停 します フ ： gf_limit = 1.e16
IMSL_GROWTH_FACTOR, float *gf (出力)
gf ウ 消去法 任意 過程 絶対値が最大 要素 A 絶対値が最大 要素 除したも す
IMSL_SMALLEST_PIVOT, float *small_pivot (出力) 大 さが最小 ッ 要素 値へ インタ
IMSL_NUM_NONZEROS_IN_FACTOR, int *num_nonzeros (出力) 因子にあ 非 ロ数 合計 含 へ インタ
IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, f_complex *values (入力)
縮さ た疎 列形式 (CSC) 係数行列 扱います こ 格納方式 説明 イン ロダ
ク ン 参照し く さい
IMSL_MEMORY_BLOCK_SIZE, int blocksize (入力)
fill-in た に割 当 領域が必要 場合 blocksize 個 新しい非 ロ 要素
に対し 充分 領域 割 当 ます フ ： blocksize = nz
説明
関数 imsl_c_lin_sol_gen_coordinate A が疎 あ 線形方程式 Ax = b 解 ます フ 改良さ た一般化対称Markowitz ッ 選択法 使用し A LU 分解 最初に実行す こ
によ いわゆ one-off (一回限 ) 問題 解 ます 消去演算 に実行さ saxpy 演算が
行交換に加え 右辺に拡張さ 因数 L 陽的に 保存さ ませ 従 式 Ly = b 陰
的に解 ます そ 因子 U Ux =y 解 x 計算す 角解法に渡さ ます
複数 方程式 Ax = b A 変更せ に解く場合 通常 分解 一回 け計算し 種々 右辺に対し
複数 前進 後退解法 実行す こ が効率的 す こ 場合 因子 L 陽的に保存さ 全 列 行 交換が記録さ ます そ 角方程式 Ly = b Ux = y 解 ます ー
IMSL_RETURN_SPARSE_LU_FACTOR IMSL_RETURN_LU_IN_COORD い プ ン 指定 し 分解 し IMSL_SOLVE_ONLY 組 合わせ IMSL_SUPPLY_SPARSE_LU_FACTOR 又 IMSL_SUPPLY_SPARSE_LU_IN_COORD い 使用し こ 分解 戻し 異 右辺に対し 順に こ 関数 呼びます IMSL_RETURN_SPARSE_LU_FACTOR 使用す 場合
imsl_lin_sol_gen_coordinate 最終 呼び出しに L U 保存す た に使用さ たヒープ 領域 解放す IMSL_FREE_SPARSE_LU_FACTOR 含 け ませ
ATx = b 解が必要 あ IMSL_TRANSPOSE プ ン 指定します こ キーワー 前進消去
後退代入 変え UTy = b LTx = y 解 得 ようにします 従 分解 行 う一回 呼び出し
Ax = b ATx = b 両方 解 求 こ が ます
プ ン IMSL_CONDITION A L1 条件数 推定値 計算す に使用さ ます こ ア
Higham に基 ます IMSL_CONDITION 指定す one-off (一回限 ) 問題が解 にし
も 完全 L が計算さ 保存さ こ に ます こ Higham法が Az = r ATz = r 問題 解 必要 す す
フ ッ 交換法 対称Markowitz法 す 行優先 あ い 列優先 問題 場合
IMSL_ROW_MARKOWITZ あ い IMSL_COLUMN_MARKOWITZ い 選ぶこ によ fill-in
低減に ます Markowitz法 ッ 候補 し あ 選 数 行 又 列 検索します
フ 数 3 すが IMSL_NUM_OF_SEARCH_ROWS 使用し 変更す こ が可能 す
プ ン IMSL_DROP_TOLERANCE fill-in 減少させ こ が 許容値 設定す た に使用 ます こ 指定した ロップ許容値よ も大 さが小さい新しい fill-in 要素が分解に追加さ こ
防 す す こ 誤差 分解に持 込 最終的 解 よ 正確にす た に
IMSL_ITERATIVE_REFINEMENT 使用す こ 勧 します こ 場合 ー フ ロップ
許容値 領域 大 さ 節約 こ 精度改良 た に必要 反復解 ップに余分 時 間が こ す
関数 imsl_c_lin_sol_gen_coordinate 分解 間 あ 点 密 分解法に え プ ン
提供し います IMSL_HYBRID_FACTORIZATION 選択す こ によ こ プ ンが有効
に ます こ プ ンが必要 2 パ ータ １ density フ ー ッ え
前に アク 小行列 た 最小密度 指定します 1.0 密度 完全 fill-in 示します も
う一 パ ータ order_bound 密 フ ー ッ に変換さ アク 小行列 階数 限
設定します こ 密 分解 O (n3) 演算が性能劣化 原因 に えが早す い ようにす た に使用さ ます こ プ ン ヒープ領域 増大す こ に注意し 下さい
例題 1
次 行列 考えます：
10 7 0 0 0 0 0
0 3 2 3 1 2 0 0
0 0 4 2 0 0 0
2 4 0 0 1 6 1 3 0
5 4 0 0 5 12 2 7 7
1 12 2 8 0 0 0 3 7
i
i i
A i
i i i
i i i
i i i
xT = (1 + i, 2 + 2i, 3 + 3i, 4 + 4i, 5 + 5i, 6 + 6i) す 次 ように ます
Ax = (3 + 17i, 19 + 5i, 6 + 18i, 38 + 32i, 63 + 49i, 57 + 83i)T
#include <imsl.h>
int main() {
static Imsl_c_sparse_elem a[] = {0, 0, {10.0, 7.0},
1, 1, {3.0, 2.0}, 1, 2, {-3.0, 0.0}, 1, 3, {-1.0, 2.0}, 2, 2, {4.0, 2.0}, 3, 0, {-2.0, -4.0}, 3, 3, {1.0, 6.0}, 3, 4, {-1.0, 3.0}, 4, 0, {-5.0, 4.0}, 4, 3, {-5.0, 0.0}, 4, 4, {12.0, 2.0}, 4, 5, {-7.0, 7.0}, 5, 0, {-1.0, 12.0}, 5, 1, {-2.0, 8.0}, 5, 5, {3.0, 7.0}};
static f_complex b[] =
{{3.0, 17.0}, {-19.0, 5.0}, {6.0, 18.0}, {-38.0, 32.0}, {-63.0, 49.0}, {-57.0, 83.0}};
int n = 6;
int nz = 15;
f_complex *x;
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b, 0);
imsl_c_write_matrix ("solution", n, 1, x, 0);
imsl_free (x);
}
出力
solution
1 ( 1, 1) 2 ( 2, 2) 3 ( 3, 3) 4 ( 4, 4) 5 ( 5, 5) 6 ( 6, 6)
例題 2
こ 例題 A = E(1000, 10) します １ 線形方程式 解 LU 分解が返さ ます そ
分解さ た同 係数行列 A 使用し 2番目 方程式 解 ます 最大 絶対誤差 実行時間 比 率 プ ン し 前進 後退解法が 分解 解 計算時間 約 10パー ン あ こ 示します こ 比率 係数行列 元数 初 非 ロ 数 特に消去 際に生成さ fill-in 量によ 非常 に変化します 時間測定 結果 ンに強く依存す こ に注意し く さい
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
int main() {
Imsl_c_sparse_elem *a;
Imsl_c_sparse_lu_factor lu_factor;
f_complex *b, *x, *mod_five, *mod_ten;
float error_factor_solve, error_solve;
double time_factor_solve, time_solve;
int n = 1000, c = 10, i, nz, index;
/* Get the coefficient matrix */
a = imsl_c_generate_test_coordinate (n, c, &nz, 0);
/* Set two different predetermined solutions */
mod_five = (f_complex*) malloc (n*sizeof(*mod_five));
mod_ten = (f_complex*) malloc (n*sizeof(*mod_ten));
for (i=0; i<n; i++) {
mod_five[i] = imsl_cf_convert ((float)(i % 5), 0.0);
mod_ten[i] = imsl_cf_convert ((float)(i % 10), 0.0);
}
/* Choose b so that x will approximate mod_five */
b = imsl_c_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_five,
0);
/* Time the factor/solve */
time_factor_solve = imsl_ctime();
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b, IMSL_RETURN_SPARSE_LU_FACTOR, &lu_factor, 0);
time_factor_solve = imsl_ctime() - time_factor_solve;
/* Compute max abolute error */
error_factor_solve = imsl_c_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
imsl_free (b);
imsl_free (x);
/* Get new right hand side -- b = A * mod_ten */
b = imsl_c_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, mod_ten, 0);
/* Use the previously computed factorization to solve Ax = b */
time_solve = imsl_ctime();
x = imsl_c_lin_sol_gen_coordinate (n, nz, a, b, IMSL_SUPPLY_SPARSE_LU_FACTOR, &lu_factor, IMSL_SOLVE_ONLY,
0);
time_solve = imsl_ctime() - time_solve;
error_solve = imsl_c_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_ten,
IMSL_INF_NORM, &index, 0);
imsl_free (b);
imsl_free (x);
/* Print errors and ratio of execution times */
printf ("absolute error (factor/solve) = %e\n", error_factor_solve);
printf ("absolute error (solve) = %e\n", error_solve);
printf ("time_solve/time_factor_solve = %f\n", time_solve/time_factor_solve);
}
出力
absolute error (factor/solve) = 2.389053e-06 absolute error (solve) = 7.656095e-06 time_solve/time_factor_solve = 0.070313