【連立一次方程式】双共役勾配法(BiCG法) C++コード付き
今回は正定値対称行列にしか適用できない共役勾配法(Conjugate Gradient法)をパワーアップした手法である双共役勾配法(Bi-Conjugate Gradient法)のC++コードを公開します。CG法についてはこちらをご覧下さい。
しくみはおいおい説明するとして、とりあえずコードをおいておきます。というより現在進行形で勉強中です!お待ちください。
今回解く連立一次方程式は
とあらわされます。ここで、 は の行列、 と は 次元のベクトルです。 と が既知で、 が未知です。具体的には
であり、これを解くと となります。コードを書く際に、
- 作者: 二宮市三
- 出版社/メーカー: 共立出版
- 発売日: 2006/02
- メディア: 単行本
- クリック: 2回
- この商品を含むブログ (2件) を見る
の「第4章 大型線型方程式の反復解法」を参考しています。ではコードを以下に示します。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <stdio.h> #include <time.h> #include <stdlib.h> using namespace std; int i,j,k,y,x; //matix vector multiplication Ax=b inline void mvm(double A[],double x[],double b[], int size) { for(i=0;i<size;i++) { b[i]=0.; for(j=0;j<size;j++) { b[i]+=A[i*size+j]*x[j]; } } } //matix vector multiplication Atx=b (transpose) inline void tmvm(double At[],double x[],double b[], int size) { for(j=0;j<size;j++) { b[j]=0.; for(i=0;i<size;i++) { b[j]+=At[j*size+i]*x[i]; } } } //inner product// inline double ip(double a[], double b[], int size) { double value=0.; for(i=0;i<size;i++) value+=a[i]*b[i]; return value; } main() { int size=3; double alpha,beta; double r0,rk,rk1; double eps=pow(2.0,-50); double *r=new double[size]; double *rs=new double[size]; double *p=new double[size]; double *b=new double[size]; double *ps=new double[size]; double *Atps=new double[size]; double *x=new double[size]; double *Ax=new double[size]; double *Ap=new double[size]; double *A=new double[size*size]; A[0*size+0]=3.;A[0*size+1]=1.;A[0*size+2]=1.; A[1*size+0]=1.;A[1*size+1]=4.;A[1*size+2]=1.; A[2*size+0]=1.;A[2*size+1]=1.;A[2*size+2]=5.; b[0]=5.;b[1]=6.;b[2]=7.; //IC// for(i=0;i<size;i++) x[i]=0.; mvm(A,x,Ax,size); for(i=0;i<size;i++) r[i]=b[i]-Ax[i]; for(i=0;i<size;i++) p[i]=r[i]; for(i=0;i<size;i++) ps[i]=r[i]; for(i=0;i<size;i++) rs[i]=r[i]; r0=sqrt(ip(r,r,size)); //main loop// for(k=0;k<10000;k++) { mvm(A,p,Ap,size); alpha=ip(rs,r,size)/ip(ps,Ap,size); rk=ip(rs,r,size); for(i=0;i<size;i++) x[i]=x[i]+alpha*p[i]; for(i=0;i<size;i++) r[i]=r[i]-alpha*Ap[i]; rk1=sqrt(ip(r,r,size)); if(rk1/r0<=eps) break; tmvm(A,ps,Atps,size); for(i=0;i<size;i++) rs[i]=rs[i]-alpha*Atps[i]; beta=ip(rs,r,size)/rk; for(i=0;i<size;i++) p[i]=r[i]+beta*p[i]; for(i=0;i<size;i++) ps[i]=rs[i]+beta*ps[i]; } for(i=0;i<size;i++) cout<<x[i]<<endl; delete [] r,rs,p,b,ps,Atps,x,Ax,Ap,A; return 0; }
Bi-CG法を安定化した、安定化双共役勾配法(Bi-conjugate gradient stabilized法)はこちらをご覧ください。