« BiCGSTAB法 | トップページ | GKrellM »

2008/12/09

CG法

まずはCG法(Conjugate Gradient Method)で連立方程式を解く練習。

CG方のアルゴリズムはWikipediaに任せることにして、CG法のプログラムは以下のようになった。

/*==================================================*/
// CG法テストプログラム
/*==================================================*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N    10
#define TMAX 100
#define EPS  (1.0e-6)

typedef double Vector[N];       // ベクトル
typedef double Matrix[N][N];    // 行列

// ベクトル初期化
void init_vector(Vector x)
{
  int i;
  for(i = 0; i < N; i++){
    x[i] = 0;
  }
}

// ベクトルに行列を作用 y = Ax
void vector_x_matrix(Vector y, Matrix a, Vector x)
{
  int    i, j;
  double vxm;
  for(i = 0; i < N; i++){
    vxm = 0;
    for(j = 0; j < N; j++){
      vxm += a[i][j]*x[j];
    }
    y[i] = vxm;
  }
}

// 内積を計算
double dot_product(Vector x, Vector y)
{
  int    i;
  double dot_p = 0;
  for(i = 0; i < N; i++){
    dot_p += x[i]*y[i];
  }
  return dot_p;
}

// ベクトルノルムを計算
// ベクトルノルム := sgm(0〜N-1)|x[i]|
double vector_norm(Vector x)
{
  int    i;
  double norm = 0;
  for(i = 0; i < N; i++){
    norm += fabs(x[i]);
  }
  return norm;
}

// CG法
void cg_method(Matrix a, Vector x, Vector b)
{
  static Vector p;   // 探索方向ベクトル
  static Vector r;   // 残差ベクトル
  static Vector ax;
  static Vector ap;
  int           i, iter;

  // Axを計算
  vector_x_matrix(ax, a, x);

  // pとrを計算 p = r := b - Ax
  for(i = 0; i < N; i++){
    p[i] = b[i] - ax[i];
    r[i] = p[i];
  }

  // 反復計算
  for(iter = 1; iter < TMAX; iter++){
    double alpha, beta, err = 0;

    // alphaを計算
    vector_x_matrix(ap, a, p);
    alpha = +dot_product(p, r)/dot_product(p, ap);
   
    for(i = 0; i < N; i++){
      x[i] += +alpha*p[i];
      r[i] += -alpha*ap[i];
    }
   
    err = vector_norm(r);   // 誤差を計算
    printf("LOOP : %d\t Error : %g\n", iter, err);
    if(EPS > err) break;

    // EPS < err ならbetaとpを計算してループを継続
    beta = -dot_product(r, ap)/dot_product(p, ap);
    for(i = 0; i < N; i++){
      p[i] = r[i] + beta*p[i];
    }
  }
}

int main(void)
{
  // 連立方程式 Ax = b
  // 行列Aは正定値対象行列
  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 = {1.0, 1.0, 1.0, 1.0, 1.0,  1.0, 1.0,  1.0, 1.0,   1.0};
  int    i;

  // CG法でAx=bを解く
  cg_method(a, x, b);
 
  printf("###Calc End.###\n");
  for(i = 0; i < N; i++){
    printf("x[%d] = %2g\n", i, x[i]);
  }
 
  return 0;
}

これを実効すると以下のようになる。

$ ./a.out
LOOP : 1         Error : 33.625
LOOP : 2         Error : 19.5236
LOOP : 3         Error : 8.64938
LOOP : 4         Error : 3.7894
LOOP : 5         Error : 1.49629
LOOP : 6         Error : 0.757903
LOOP : 7         Error : 0.348423
LOOP : 8         Error : 0.141565
LOOP : 9         Error : 0.0454147
LOOP : 10        Error : 1.17553e-16
###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法 | トップページ | GKrellM »

C」カテゴリの記事

Numerical Computing」カテゴリの記事

コメント

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: CG法:

« BiCGSTAB法 | トップページ | GKrellM »