IMSL_CSC_FORMAT, int *col_ptr, int *row_ind, float *values, IMSL_MEMORY_BLOCK_SIZE, int block_size,
0)
プ ン引数
IMSL_RETURN_SPARSE_LU_FACTOR, Imsl_f_sparse_lu_factor *lu_factor 。出力)
Imsl_f_sparse_lu_factor型 構造体 ア こ 構造体 内部 インタ
imsl_f_lin_sol_gen_coordinate によ LU 分解 指すた に初期化さ ます IMSL_SUPPLY_SPARSE_LU_FACTOR, Imsl_f_sparse_lu_factor *lu_factor 。入力)
Imsl_f_sparse_lu_factor型 構造体 ア こ 構造体 プ ン
IMSL_RETURN_SPARSE_LU_FACTOR imsl_f_lin_sol_gen_coordinate によ 計算さ た入力行列 LU 分解 含 ます
IMSL_FREE_SPARSE_LU_FACTOR,
戻 前に A LU 分解 含 ンクさ た ータ構造体 解放します 因子が 必要 く た場合 こ プ ン 使用します
IMSL_RETURN_SPARSE_LU_IN_COORD, Imsl_f_sparse_elem **lu_coordinate, int **row_pivots, int **col_pivots 。出力)
LU 分解が 標対応形式 lu_coordinate 長さ nz 配列に返さ ます こ
Imsl_f_sparse_lu プ 化さ た内部表現よ 簡潔 す こ 利 点
SOLVE_ONLY 呼び出し 間 因子 内部表現が再構成さ け いこ す
し し こ 因子がプログ 終了 後 格納さ 再度 後 実行時にロー さ あ IMSL_RETURN_LU_IN_COORD IMSL_SUPPLY_LU_IN_COORD 組 合わ せ 格納 込 が単純 形式 最良 選択に ます
IMSL_SUPPLY_SPARSE_LU_IN_COORD, int nzlu, Imsl_f_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, float 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_NUMBER_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_STABILITY_FACTOR, float s_factor 。入力)
ッ 要素 絶対値 そ 行 s_factor 除算さ た絶対値 最大 要素よ 大
く く いけませ
フ : s_factor = 10.0
IMSL_GROWTH_FACTOR_LIMIT, float gf_limit 。入力)
成長因子が gf_limit 越え 場合 計算 停 します
フ : gf_limit = 1.0e16
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, float *values 。入力)
係数行列 縮さ た疎 列形式(CSC) す
こ 格納方式 説明 イン ロダク ン 参照し く さい
IMSL_MEMORY_BLOCKSIZE, int blocksize 。入力)
領域が fill-in た に割 当 け い場合 blocksize 個 新しい非
ロ 要素に対し 十分 領域 割 当 ます
フ : blocksize = nz
説明
関数 imsl_f_lin_sol_gen_coordinate A が疎 あ 線形方程式 Ax = b 解 ます フ 最初に改良さ た一般化対称 Markowitz ッ 選択法 使用し A LU 分解 行うこ に
よ いわゆ one-off (一回限 ) 問題 解 ます 因子 L 陽的に 保存さ ませ こ 行
交換に加え 消去演算 に行 わ saxpy 演算が右辺にも拡張さ た す 従 式 Ly
= b
陰的に解 ます 因子 U Ux = y 解 x 計算す 角解法に渡さ ます同 A 複数 方程式 Ax = b 解く場合 分解 一回 け計算し 種々 右辺に対し 何回も 前進
後退解法 す こ が効率的 す こ 場合 因子 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_f_lin_sol_gen_coordinate 分解 間 密 分解法に え プ ン 提供し
います こ IMSL_HYBRID_FACTORIZATION す こ プ ンが必要 2 パ ータ
１ density フ ー ッ え 前に アク 小行列 最小密度 指定します 1.0
密度 完全 fill-in 示します もう一 パ ータ order_bound 密 フ ー ッ に変換さ
アク 小行列 次数 限 設定します こ えが早す 密 分解 O (n3) 演算 た に性能が劣化す こ 防 た に使用さ ます こ プ ン ヒープ領域 増大 させ こ に注意し 下さい
例題 1
例 し 次 行列 考えます
6 0 0 0 2 1
3 1 5 0 0 1
0 1 10 0 0 2
0 0 1 3 10 0
0 0 0 0 0 10 A
xT = (1, 2, 3, 4, 5, 6) し Ax = (10, 7, 45, 33, 34, 31)T ます A 非 ロ 数 nz = 15 す
#include <imsl.h>
int main() {
Imsl_f_sparse_elem a[] = {0, 0, 10.0,
1, 1, 10.0, 1, 2, -3.0, 1, 3, -1.0, 2, 2, 15.0, 3, 0, -2.0, 3, 3, 10.0, 3, 4, -1.0, 4, 0, -1.0, 4, 3, -5.0, 4, 4, 1.0, 4, 5, -3.0, 5, 0, -1.0, 5, 1, -2.0, 5, 5, 6.0};
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
int nz = 15;
float *x;
x = imsl_f_lin_sol_gen_coordinate (n, nz, a, b, 0);
imsl_f_write_matrix ("solution", 1, n, x, 0);
imsl_free (x);
}
出力
solution
1 2 3 4 5 6
1 2 3 4 5 6
例題 2
こ 例題 A = E(1000, 10) します １ 方程式 解 LU 分解が返さ ます そ 分解さ
た同 係数行列 A 使用し 2番目 方程式 解 ます 最大 絶対誤差 実行時間 比率 プ ン し 前進 後退解法が 分解 解 計算時間 約 10パー ン あ こ 示します こ 比率 係数行列 元数 初 非 ロ 数 特に消去 際に生成さ fill-in 量によ 非常に変 化します 時間測定 結果 ンに強く依存す こ に注意し く さい
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
int main() {
Imsl_f_sparse_elem *a;
Imsl_f_sparse_lu_factor lu_factor;
float *b;
float *x;
float *mod_five;
float *mod_ten;
float error_factor_solve;
float error_solve;
double time_factor_solve;
double time_solve;
int n = 1000;
int c = 10;
int i;
int nz;
int index;
/* Get the coefficient matrix */
a = imsl_f_generate_test_coordinate (n, c, &nz, 0);
/* Set two different predetermined solutions */
mod_five = (float*) malloc (n*sizeof(*mod_five));
mod_ten = (float*) malloc (n*sizeof(*mod_ten));
for (i=0; i<n; i++) {
mod_five[i] = (float) (i % 5);
mod_ten[i] = (float) (i % 10);
}
/* Choose b so that x will approximate mod_five */
b = (float *) imsl_f_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_f_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 absolute error */
error_factor_solve = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
free (mod_five);
imsl_free (b);
imsl_free (x);
/* Get new right hand side -- b = A * mod_ten */
b = (float *) imsl_f_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_f_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_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_ten,
IMSL_INF_NORM, &index, 0);
free (mod_ten);
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) = 9.179115e-05 absolute error (solve) = 2.160072e-04 time_solve/time_factor_solve = 0.093750
例題 3
こ 例題 A = E (500, 50) 線形方程式 Ax = b 解 ます そ 後 同 方程式 大 ロ
ップ許容値 使用し 解 ます 最後に 計算した 分解 使用し 同 線形方程式 反 復改良法 解 ます 測定時間結果 ンに強く依存す こ に注意し く さい
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
int main() {
Imsl_f_sparse_elem *a;
Imsl_f_sparse_lu_factor lu_factor;
float *b;
float *x;
float *mod_five;
float error_zero_drop_tol;
float error_nonzero_drop_tol;
float error_nonzero_drop_tol_IR;
double time_zero_drop_tol;
double time_nonzero_drop_tol;
double time_nonzero_drop_tol_IR;
int nz_nonzero_drop_tol;
int nz_zero_drop_tol;
int n = 500;
int c = 50;
int i;
int nz;
int index;
/* Get the coefficient matrix */
a = imsl_f_generate_test_coordinate (n, c, &nz, 0);
for (i=0; i<nz; i++) a[i].val *= 0.05;
/* Set a predetermined solution */
mod_five = (float*) malloc (n*sizeof(*mod_five));
for (i=0; i<n; i++)
mod_five[i] = (float) (i % 5);
/* Choose b so that x will approximate mod_five */
b = imsl_f_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_zero_drop_tol = imsl_ctime();
x = imsl_f_lin_sol_gen_coordinate (n, nz, a, b, IMSL_NUM_NONZEROS_IN_FACTOR, &nz_zero_drop_tol, 0);
time_zero_drop_tol = imsl_ctime() - time_zero_drop_tol;
/* Compute max abolute error */
error_zero_drop_tol = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
imsl_free (x);
/* Solve the same problem, with drop tolerance = 0.005 */
time_nonzero_drop_tol = imsl_ctime();
x = imsl_f_lin_sol_gen_coordinate (n, nz, a, b, IMSL_RETURN_SPARSE_LU_FACTOR, &lu_factor, IMSL_DROP_TOLERANCE, 0.005,
IMSL_NUM_NONZEROS_IN_FACTOR, &nz_nonzero_drop_tol, 0);
time_nonzero_drop_tol = imsl_ctime() - time_nonzero_drop_tol;
/* Compute max abolute error */
error_nonzero_drop_tol = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
imsl_free (x);
/* Solve the same problem with IR, use last factorization */
time_nonzero_drop_tol_IR = imsl_ctime();
x = imsl_f_lin_sol_gen_coordinate (n, nz, a, b, IMSL_SUPPLY_SPARSE_LU_FACTOR, &lu_factor, IMSL_SOLVE_ONLY,
IMSL_ITERATIVE_REFINEMENT, 0);
time_nonzero_drop_tol_IR = imsl_ctime() - time_nonzero_drop_tol_IR;
/* Compute max abolute error */
error_nonzero_drop_tol_IR = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five,
IMSL_INF_NORM, &index, 0);
imsl_free (x);
imsl_free (b);
/* Print errors and ratio of execution times */
printf ("drop tolerance = 0.0\n");
printf ("\tabsolute error = %e\n", error_zero_drop_tol);
printf ("\tfillin = %d\n\n", nz_zero_drop_tol);
printf ("drop tolerance = 0.005\n");
printf ("\tabsolute error = %e\n", error_nonzero_drop_tol);
printf ("\tfillin = %d\n\n", nz_nonzero_drop_tol);
printf ("drop tolerance = 0.005 (with IR)\n");
printf ("\tabsolute error = %e\n", error_nonzero_drop_tol_IR);
printf ("\tfillin = %d\n\n", nz_nonzero_drop_tol);
printf ("time_nonzero_drop_tol/time_zero_drop_tol = %f\n", time_nonzero_drop_tol/time_zero_drop_tol);
printf ("time_nonzero_drop_tol_IR/time_zero_drop_tol = %f\n", time_nonzero_drop_tol_IR/time_zero_drop_tol);
}
出力
drop tolerance = 0.0
absolute error = 3.814697e-06 fillin = 9530
drop tolerance = 0.005
absolute error = 2.699481e+00 fillin = 8656
drop tolerance = 0.005 (with IR) absolute error = 1.907349e-06 fillin = 8656
time_nonzero_drop_tol/time_zero_drop_tol = 1.086957 time_nonzero_drop_tol_IR/time_zero_drop_tol = 0.840580
反復改良法(IR)が使用さ い 絶対誤差に注目し く さい また 反復改良法 極 時 間が こ に注意し く さい こ 場合に IR 解 大体 分解 同 い時間が ます
こ 問題 大 い ロップ許容値 反復改良法 使用 フ 2倍にあた 時間 fill-in 10パー ン 減少す こ が ました 容量が厳しい状況下 こ よう ー
フが適 場合があ ます ロップ許容値 大 く選ぶ LU に大 い誤差 持 込 反 復改良 束 妨 こ に注意し く さい