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