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