IMSL_FCN_W_DATA, void amultp (float *p, float *z, void *data), void *data (入力)
z = Ap 計算す ー 提供 関数 こ ー によ 提供さ た ータへ イン
タも受け入 ます data ー 提供 関数に ータ 受け渡す インタ す 詳 細 本書 イン ロダク ン ー 提供 関数へ ータ 受け渡し 参照し く さい
IMSL_PRECOND_W_DATA, void precond (float *r, float *z, void *data), void *data, 入力
z = M -1 r 設定す ー 提供 関数 M 前処理行列 す こ ー によ
提供さ ータへ インタも受け入 ます data ー 定義関数に ータ 受け渡す インタ す 詳細 本書 イン ロダク ン ー 提供 関数へ ータ 受け渡し 参照し く さい
説明
関数imsl_f_lin_sol_def_cg プ ン 前処理 持 共役勾配法 使 対称定値線形方程
式 Ax = b 解 ます こ 手法 Golub Van Loan (1983年 第10章) Hageman Young (1981年 第7章) に詳しく説明さ います
前処理行列 M A 近似す 行列 そ 線形方程式 Mz = r が簡単に解け も す こ 2 性質 相反します こ バ ン こ が 現在 研究 ックに います
imsl_f_lin_sol_def_cg フ M = I す プ ン IMSL_JACOBI が選択さ
M A 対角項に ッ さ ます
必要 反復回数 行列や誤差 許容度によ ます 大ま に 次 ように ます：
maxit =
n
n >>1 に対し 詳細 記 参考文献 参照し く さいM 前処理行列 b p r x z ベク τ 望ま 相対誤差 します 使用さ ア 以下 通 す
ここに 反復行列 G = I M -1
A
最大固有値 max(G)
推定 す 停 基準 次に基 ます：(Hageman Young 1981年 ージ 148-151)
ここに
又，次式が成 立 ます：
ここ Tn 対称 重対角行列
1 2
2 2 3
3 3
Tn
また
1
= 1
1/
1k = 1 k
/
k-1 1/
k,
/ 1k Bk k
通常固有値 少し 反復 得 ます
例題 1
こ 例題 ，線形方程式 解 ます 係数行列 フ 行列 し 格納さ ます
#include <imsl.h>
static void amultp (float*, float*);
int main() {
int n = 3;
float b[] = {27.0, -78.0, 64.0};
float *x;
x = imsl_f_lin_sol_def_cg (n, amultp, b, 0);
imsl_f_write_matrix ("x", 1, n, x, 0);
}
static void amultp (float *p, float *z) {
static float a[] = {1.0, -3.0, 2.0,
-3.0, 10.0, -5.0, 2.0, -5.0, 6.0};
int n = 3;
imsl_f_mat_mul_rect ("A*x", IMSL_A_MATRIX, n, n, a, IMSL_X_VECTOR, n, p, IMSL_RETURN_USER, z, 0);
}
出力
x
1 2 3 1 -4 7
例題 2
こ 例題 異 前処理行列 用い 規則的 c c グ ッ c = 100 プ 方程式 有限差分解 求 線形方程式 解 ます 行列 A = E (c2
, c)
す 最初 コビ前処理 選択し M = diag (A) します 実行さ た反復回数 最大絶対誤差がプ ン さ ます
次に A 対称 重対角部分によ よ 複雑 前処理行列 M 使用します
対称正定値バン 解法プログ 一回 M 分解す に使用さ 続い 前進後退解法 実行す こ に注意し 下さい 実行さ た反復回数 最大絶対誤差がプ ン さ ます 反復 大幅 減 少に注意し 下さい
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
static void amultp (float*, float*);
static void precond (float*, float*);
static Imsl_f_sparse_elem *a;
static int n = 2500;
static int c = 50;
static int nz;
int main() {
int maxit = 1000;
int i;
int index;
float *b;
float *x;
float *mod_five;
float *diagonal;
float norm;
n = c*c;
mod_five = (float*) malloc (n*sizeof(*mod_five));
diagonal = (float*) malloc (n*sizeof(*diagonal));
b = (float*) malloc (n*sizeof(*b));
/* Generate coefficient matrix */
a = imsl_f_generate_test_coordinate (n, c, &nz,
0);
/* Set a predetermined answer and diagonal */
for (i=0; i<n; i++) {
mod_five[i] = (float) (i % 5);
diagonal[i] = 4.0;
}
/* Get right hand side */
amultp (mod_five, b);
/* Solve with jacobi preconditioning */
x = imsl_f_lin_sol_def_cg (n, amultp, b, IMSL_MAX_ITER, &maxit,
IMSL_JACOBI, diagonal, 0);
/* Find max absolute error, print results */
norm = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five, IMSL_INF_NORM, &index,
0);
printf ("iterations = %d, norm = %e\n", maxit, norm);
imsl_free (x);
/* Solve same system, with different preconditioner */
x = imsl_f_lin_sol_def_cg (n, amultp, b, IMSL_MAX_ITER, &maxit,
IMSL_PRECOND, precond, 0);
norm = imsl_f_vector_norm (n, x, IMSL_SECOND_VECTOR, mod_five, IMSL_INF_NORM, &index,
0);
printf ("iterations = %d, norm = %e\n", maxit, norm);
}
/* Set z = Ap */
static void amultp (float *p, float *z) {
imsl_f_mat_mul_rect_coordinate ("A*x", IMSL_A_MATRIX, n, n, nz, a,
IMSL_X_VECTOR, n, p,
IMSL_RETURN_USER_VECTOR, z, 0);
}
/* Solve Mz = r */
static void precond (float *r, float *z) {
static float *m;
static float *factor;
static int first = 1;
float *null = (float*) 0;
if (first) {
/* Factor the first time through */
m = imsl_f_generate_test_band (n, 1, IMSL_SYMMETRIC_STORAGE,
0);
imsl_f_lin_sol_posdef_band (n, m, 1, null, IMSL_FACTOR, &factor,
IMSL_FACTOR_ONLY, 0);
first = 1;
}
/* Perform the forward and back solves */
imsl_f_lin_sol_posdef_band (n, m, 1, r, IMSL_FACTOR_USER, factor,
IMSL_SOLVE_ONLY, IMSL_RETURN_USER, z, 0);
}
出力
iterations = 115, norm = 1.382828e-05 iterations = 75, norm = 7.319450e-05
重大 ー
IMSL_STOP_USER_FCN ー 定義関数 計算 停 が要求さ た ー フ グ= "#"