2次元FDM用CG法(FORTRAN)プログラム Y-対称NO.1
2003/08/26 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)を使用した、演算量の少ないバージョン
係数行列Aと右辺ベクトルbが解析領域に対してyの中心で対称なら、計算解xも対称に
なる工夫をしたプログラムである。但し、Cコンパイラーが同一式中の加減算の順序を
変更するときには、対称性は保証されない。
2. プログラム
C=================================================================C
SUBROUTINE CG2DS(A,B,NX,NY,X,R,P,Q,EPS,ITER,ERR,IERR)
C=================================================================C
C Solve Ax=b by CG NO.1 with 2 dimensional FDM C
C Alpha=(R,R)/(P,AP), Beta=new(R,R)/old(R,R) C
C for symmetric solution on Y-axis C
C-----------------------------------------------------------------C
C A(0:NX,0:NY,3) R*8, In, A Coefficient Matrix C
C with symmetric C
C B(0:NX,0:NY) R*8, In, A Right-hand Vector(b) C
C NX I*4, In, Grid Numbers on X-axis C
C NY I*4, In, Grid Numbers on Y-axis C
C X(0:NX,0:NY) R*8, I/O, Initial and Solution vector C
C R(0:NX,0:NY) R*8, Out, Residual vector C
C P(0:NX,0:NY) R*8, Wk, Work Vector (P) C
C Q(0:NX,0:NY) R*8, Wk, Work Vector (Q) C
C EPS R*8, In, if ||r||/||b|| <= EPS --> return C
C ITER I*4, I/O, Number of Iteration C
C ERR R*8, Out, ERR=||r||/||b|| C
C IERR I*4, Out, IERR=0, Normal Return C
C =1, No Convergent C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/08/25 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(0:NX,0:NY,3), B(0:NX,0:NY)
DIMENSION X(0:NX,0:NY), R(0:NX,0:NY)
DIMENSION P(0:NX,0:NY), Q(0:NX,0:NY)
C Even or Odd of NY
MNY = mod(NY,2)
if(MNY.eq.0) then
NY2 = NY/2 - 1
else
NY2 = (NY - 1)/2
end if
C P=R=B-A*X
IERR = 0
BN = 0.0
C0 = 0.0
do j=1,NY2
k = NY - j
do i=1,NX-1
BN = BN + B(i,j)**2 + B(i,k)**2
R(i,j) = B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
1 - A(i,j,3)*X(i,j) - A(i+1,j,2)*X(i+1,j)
2 - A(i,j+1,1)*X(i,j+1)
R(i,k) = B(i,k) - A(i,k+1,1)*X(i,k+1) - A(i,k,2)*X(i-1,k)
1 - A(i,k,3)*X(i,k) - A(i+1,k,2)*X(i+1,k)
2 - A(i,k,1)*X(i,k-1)
C0 = C0 + R(i,j)**2 + R(i,k)**2
P(i,j) = R(i,j)
P(i,k) = R(i,k)
end do
end do
C Y-center line
if(MNY.eq.0) then
j = NY/2
do i=1,NX-1
BN = BN + B(i,j)**2
R(i,j) = B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
1 - A(i,j,3)*X(i,j) - A(i+1,j,2)*X(i+1,j)
2 - A(i,j+1,1)*X(i,j+1)
C0 = C0 + R(i,j)**2
P(i,j) = R(i,j)
end do
end if
C Main Loop
do kk=1,ITER
C Q=A*P, Alpha=C0/(P,Q)
Alpha = 0.0
do j=1,NY2
k = NY - j
do i=1,NX-1
Q(i,j) = A(i,j,1)*P(i,j-1) + A(i,j,2)*P(i-1,j)
1 + A(i,j,3)*P(i,j) + A(i+1,j,2)*P(i+1,j)
2 + A(i,j+1,1)*P(i,j+1)
Q(i,k) = A(i,k+1,1)*P(i,k+1) + A(i,k,2)*P(i-1,k)
1 + A(i,k,3)*P(i,k) + A(i+1,k,2)*P(i+1,k)
2 + A(i,k,1)*P(i,k-1)
Alpha = Alpha + P(i,j)*Q(i,j) + P(i,k)*Q(i,k)
end do
end do
C Y-center line
if(MNY.eq.0) then
j = NY/2
do i=1,NX-1
Q(i,j) = A(i,j,1)*P(i,j-1) + A(i,j,2)*P(i-1,j)
1 + A(i,j,3)*P(i,j) + A(i+1,j,2)*P(i+1,j)
2 + A(i,j+1,1)*P(i,j+1)
Alpha = Alpha + P(i,j)*Q(i,j)
end do
end if
Alpha = C0/Alpha
C X=X+Alpha*P, R=R-Alpha*Q
C1 = 0.0
do j=1,NY-1
do i=1,NX-1
X(i,j) = X(i,j) + Alpha*P(i,j)
R(i,j) = R(i,j) - Alpha*Q(i,j)
C1 = C1 + R(i,j)**2
end do
end do
C if(ERR <= EPS) Convergent, Beta=C1/C0
ERR = SQRT(C1/BN)
if(ERR.le.EPS) go to 100
Beta = C1/C0
C0 = C1
C P=R+Beta*P
do j=1,NY-1
do i=1,NX-1
P(i,j) = R(i,j) + Beta*P(i,j)
end do
end do
end do
IERR = 1
C
100 continue
ITER = kk
C
RETURN
END