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のコードとしては多分あってるはず。

C++入門サイト

C++を勉強するのに役立ちそうなサイトがあったのでメモ。

Link：Programing Place(旧サイト)

新サイトのコンテンツにはC++のページが移植されてないので旧サイト。解説ページの他にも参考書籍のページにはC++本のレビューが書かれている。

インクルードガード

ヘッダーを複数インクルードするとき、例えばheader1.hppとheader2.hppがあって、header2.hppが間接的にheader1.hppをインクルードしていると、header1.hpp内で定義されていることが重複して定義されてしまい、コンパイルエラーを引き起こしてしまう。

これを防ぐためにはインクルードガードしておく必要がある。。(コンパイルエラーを引き起こすのはクラス定義の重複だけだと思うけど、それ以外の定義は重複してもいいからガードしないのか？それとも気持ち悪いからちゃんとガードしておくのか？)

具体的には次のようにするようだ。

#ifndef ___XYZ
#define ___XYZ

//定義を記述

#endif

XYZはもちろんヘッダーごとに違う名前にしておく。header1.hppなら___HEADER1_HPPみたいに。先頭にアンダースコア(_)を3つつけているのは、アンダースコア2つで始まる名前がC++に予約されているためだとか。

(___XYZの部分だけど、これを__XYZ__としているものも見かける。後ろにもアンダースコアを書く書かないというのに何か意味はあるのだろうか？)

C++のヘッダーファイルの拡張子

「.h」？それとも「.hpp」？

どちらでもコンパイルできるから自分の趣味でOKなのだろうか？

C++を勉強しよう

C++の勉強を始めようと思う。

とりあえずHello,Worldプログラムを書くと

#include <iostream>

using namespace std;

int main()
{
  cout << "Hello,World!\n";

  return 0;
}

となる。全体的にはC言語っぽいけど、よく見ると違うなぁ。

名前空間ってなんだろう？

printf()じゃなくてcout(シーアウト)なのか。この「<<」って

$ echo "Hello,World!" >> hello.txt

と似てる。