ター さ た一般化最小残差法 (GMRES : Generalized Minimum Residual) 使用し 線形方程式 Ax
= b
解 ます梗概
#include <imsl.h>
float *imsl_f_lin_sol_gen_min_residual (int n, void amultp (float *p, float *z), float *b, ..., 0) double型 関数 imsl_d_lin_sol_gen_min_residual
必須 引数
int n 。入力) 行列 行数
void amultp (float *p, float *z)
z = Ap 計算す ー 提供 関数
float *b 。入力)
右辺 含 長さ n ベク
戻 値
線形方程式 Ax = b 解 x へ インタ こ 領域 解放す に imsl_free 使用します 解
が得 い場合 NULL が返さ ます
プ ン引数 梗概
#include <imsl.h>
float *imsl_f_lin_sol_gen_min_residual (int n, void amultp (), float *b, IMSL_RETURN_USER, float x[],
IMSL_MAX_ITER, int *maxit, IMSL_REL_ERR, float tolerance, IMSL_PRECOND, void precond(),
IMSL_MAX_KRYLOV_SUBSPACE_DIM, int kdmax, IMSL_HOUSEHOLDER_REORTHOG,
IMSL_FCN_W_DATA, void amultp (), void *data, IMSL_PRECOND_W_DATA, void precond(), void *data, 0)
プ ン引数
IMSL_RETURN_USER, float x[] 。出力)
解 x 含 長さ n ー 割 当 配列 IMSL_MAX_ITER, int *maxit 。入力/出力)
整数へ インタ 最初に 許さ GMRES 反復 最大数に設定さ ます 終了時 に 使用さ た反復数が返さ ます
フ : maxit = 1000
IMSL_REL_ERR, float tolerance 。入力)
こ ア
||b Ax ||
2 ||b||2 満たす x 計算しますここ = tolerance す
フ : tolerance = sqrt(imsl_f_machine(4))
IMSL_PRECOND, void precond (float *r, float *z) 。入力)
z = M -1 r 設定す ー 提供 関数 す ここ M 前処理行列 す
IMSL_MAX_KRYLOV_SUBSPACE_DIM, int kdmax 。入力)
Krylov 部分空間 最大 次元 す す わ ター さ 前に許さ 最大
GMRES 反復数 す
フ : kdmax = imsl_i_min(n, 20)
IMSL_HOUSEHOLDER_REORTHOG,
Gram-Schmidt プロ く Householder変換によ 直交化 実行します
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 -1r ッ す ー 提供 関数 ここ M 前処理行列 す こ また ー
によ 提供さ ータへ インタも受け data ー 定義関数に渡さ ータへ インタ す 詳細 イン ロダク ン ー 提供 関数に ータ 受け渡す方法 参照し 下さい
説明
imsl_f_lin_sol_gen_min_residual GMRES 使用し 線形方程式 Ax = b 解 ます こ
コー H.F. Walker によ FORTRAN ーチン GMRES に基 ます こ 方法 Saad Schultz
(1986年) また Walker (1988年)に詳細が説明さ います
GMRES 法 近似解 x0 初期残差 r0
= b - Ax
0 ター ます 反復 m 修正 zm 次 Krylov部分空間 決定さ ます：
m( ) span( , A ,..., A
m1 )
ν= r
0 次 最小二乗問題 解 す：(zk
min
m(r0))b A ( x
0 z )
2反復 m xm = x0
+ z
mHouseholder変換によ 直交化 Gram-Schmidt法よ 必要 しませ が 計算量 多い す
し しWalker (1988年) 特に残差 減少が限界に達す に Householder法 方がよ 安定し い
こ 報告し います
例題 1
例 し 次 行列 考えます
6 0 0 0 2 1
3 1 5 0 0 1
0 1 10 0 0 2
0 0 0 15 0 0
0 0 1 3 10 0
0 0 0 0 0 10 A
Ax = (10, 7, 45, 33, 34, 31)T ように xT= (1, 2, 3, 4, 5, 6) にします 関数 imsl_f_mat_mul_rect_coordinate 用い 積 Ax 作 ます
#include <imsl.h>
void amultp (float*, float*);
int main() {
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
float *x;
x = imsl_f_lin_sol_gen_min_residual (n, amultp, b, 0);
imsl_f_write_matrix ("Solution, x, to Ax = b", 1, n, x, 0);
}
void amultp (float *p, float *z) {
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};
int n = 6;
int nz = 15;
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);
}
出力
Solution, x, to Ax = b
1 2 3 4 5 6
1 2 3 4 5 6
例題 2
こ 例題 最初 例題 同 方程式 解 ます 今回 前処理行列 使います 前処理行列
A 対角項 します
#include <imsl.h>
#include <stdio.h>
void amultp (float*, float*);
void precond (float*, float*);
int main() {
float b[] = {10.0, 7.0, 45.0, 33.0, -34.0, 31.0};
int n = 6;
float *x;
int maxit = 1000;
x = imsl_f_lin_sol_gen_min_residual (n, amultp, b, IMSL_MAX_ITER, &maxit,
IMSL_PRECOND, precond, 0);
imsl_f_write_matrix ("Solution, x, to Ax = b", 1, n, x, 0);
printf ("\nNumber of iterations taken = %d\n", maxit);
}
/* Set z = Ap */
void amultp (float *p, float *z) {
static 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};
int n = 6;
int nz = 15;
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 */
void precond (float *r, float *z) {
static float diagonal_inverse[] =
{0.1, 0.1, 1.0/15.0, 0.1, 1.0, 1.0/6.0};
int n = 6;
int i;
for (i=0; i<n; i++)
z[i] = diagonal_inverse[i]*r[i];
}
出力
Solution, x, to Ax = b
1 2 3 4 5 6
1 2 3 4 5 6
Number of iterations taken = 5
重大 ー
IMSL_STOP_USER_FCN ー 定義関数 計算 停 が要求さ た ー フ グ= "#"