BiCGSTAB法
双共役勾配(BiCG)法の残差を減少させ安定化させたのがBiCGSTAB法。BiCGのCGとの違いは非対称行列の反復計算に使用できること。(詳しく解説しようとすると、とても難しいので割愛)
とりあえず、サンプルプログラムを作ってみた。(用いた方程式は以前CG法のプログラムで使ったものと同様)
/*==================================================*/ // Bi-CGSTAB法テストプログラム /*==================================================*/ #include <cstdio> #include <cstdlib> #include <cmath> using namespace std; const int N = 10; const int ITRMAX = 100; const double EPS = 1.0e-6; typedef double Vector[N]; // ベクトル typedef double Matrix[N][N]; // 行列 // ベクトルに行列を作用 y = Ax void matrix_x_vector(Vector y, const Matrix a, const Vector x) { for(int i = 0; i < N; i++) { double mxv = 0.0; for(int j = 0; j < N; j++) { mxv += a[i][j]*x[j]; } y[i] = mxv; } } // 内積を計算 double inner_product(const Vector x, const Vector y) { double inner_p = 0.0; for(int i = 0; i < N; i++) { inner_p += x[i]*y[i]; } return inner_p; } // ベクトルノルムを計算 // ベクトルノルム := sgm(0〜N-1)|x[i]| double vector_norm(const Vector x) { double norm = 0; for(int i = 0; i < N; i++){ norm += fabs(x[i]); } return norm; } // BiCGSTAB法 /* [Algorithm] Ax = b を BiCGSTAB法で解く 1) Compute r0 = b - Ax (Choose r0* arbitary, r0* = r0, (r0, r0*) != 0) p = r = r0 (c1 = (r0*, r0) = (r0, r0)とする) 2) Loop... c2 = (r0*, Ap) alpha = c1 / c2 e = r - alpha*Ap c3 = (e, Ae) / (Ae, Ae) x_new = x + alpha*p + c3*e r_new = e - c3*Ae if(EPS > ||r||/||b||) break c1 = (r_new, r0*) = (r_new, r0) beta = c1 / (c2*c3) p_new = r_new + beta*(p - c3*Ap) */ void bicgstab(const Matrix a, Vector x, const Vector b) { static Vector r0; // 初期残差ベクトル static Vector r; // 残差ベクトル static Vector p; // 探索方向ベクトル static Vector ax; matrix_x_vector(ax, a, x); // r0,r,pを計算 r0 := b - Ax // p = r = r0 for(int i = 0; i < N; i++) { r0[i] = b[i] - ax[i]; p[i] = r[i] = r0[i]; } // (r0, r0*) != 0, r0* = r0 double c1 = inner_product(r0, r0); if(c1 == 0.0) { fprintf(stderr, "(r0, r0*) == 0!!\n"); exit(1); } for(int iter = 1; iter < ITRMAX; iter++) { static Vector ap; matrix_x_vector(ap, a, p); double c2 = inner_product(r0, ap); double alpha = c1 / c2; static Vector e; for(int i = 0; i < N; i++) { e[i] = r[i] - alpha*ap[i]; } static Vector ae; matrix_x_vector(ae, a, e); double e_dot_ae = inner_product(e, ae); double ae_dot_ae = inner_product(ae, ae); double c3 = e_dot_ae / ae_dot_ae; for(int i = 0; i < N; i++) { x[i] += alpha*p[i] + c3*e[i]; r[i] = e[i] - c3*ae[i]; } double err = vector_norm(r)/vector_norm(b); printf("LOOP : %d\t Error : %g\n", iter, err); if(EPS > err) break; c1 = inner_product(r, r0); double beta = c1 / (c2*c3); for(int i = 0; i < N; i++) { p[i] = r[i] + beta*(p[i] - c3*ap[i]); } } } int main() { // 連立方程式 Ax = b Matrix a = { {5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0, 2.0}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 5.0} }; Vector b = {3.0, 1.0, 4.0, 0.0, 5.0, -1.0, 6.0, -2.0, 7.0, -15.0}; // 初期値は適当 Vector x = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // CG法でAx=bを解く bicgstab(a, x, b); printf("###Calc End.###\n"); for(int i = 0; i < N; i++) { printf("x[%d] = %2g\n", i, x[i]); } return 0; }
実行結果は次のようになる。
LOOP : 1 Error : 0.416594
LOOP : 2 Error : 0.126537
LOOP : 3 Error : 0.0364201
LOOP : 4 Error : 0.0123132
LOOP : 5 Error : 0.00372894
LOOP : 6 Error : 0.0011802
LOOP : 7 Error : 0.00023157
LOOP : 8 Error : 4.34888e-05
LOOP : 9 Error : 4.5223e-06
LOOP : 10 Error : 7.16945e-20
###Calc End.###
x[0] = 1
x[1] = -1
x[2] = 2
x[3] = -2
x[4] = 3
x[5] = -3
x[6] = 4
x[7] = -4
x[8] = 5
x[9] = -5
解は得られているのでBiCGSTABのコードとしては多分あってるはず。
| 固定リンク
| コメント (0)
| トラックバック (0)
最近のコメント