3次元FDM用CG法NO.1(FORTRAN)プログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(p,r)/(p,Ap),β=-(r,Ap)/(p,Ap)を使用した、直交性に比較的強いバージョン
2. プログラム
C=================================================================C
SUBROUTINE CG3D1(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
C=================================================================C
C Solve Ax=b by CG NO.1 with 3 dimensional FDM C
C Alpha=(P,R)/(P,AP), Beta=-(R,AP)/(P,AP) 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/17 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 P=R=B-A*X
IERR = 0
BN = 0.0
do k=1,NZ-1
do j=1,NY-1
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)
P(i,j,k) = R(i,j,k)
end do
end do
end do
C Main Loop
do kk=1,ITER
C Q=A*P, C=(P,Q), Alpha=(P,R)/C
C = 0.0
Alpha = 0.0
do k=1,NZ-1
do j=1,NY-1
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)
C = C + P(i,j,k)*Q(i,j,k)
Alpha = Alpha + P(i,j,k)*R(i,j,k)
end do
end do
end do
Alpha = Alpha/C
C X=X+Alpha*P, R=R-Alpha*Q
RN = 0.0
Beta = 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)
RN = RN + R(i,j,k)**2
Beta = Beta + R(i,j,k)*Q(i,j,k)
end do
end do
end do
C if(ERR <= EPS) Convergent, Beta=-(R,Q)/(P,Q)
ERR = SQRT(RN/BN)
if(ERR.le.EPS) go to 100
Beta = -Beta/C
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