2次元FDM用SOR法(FORTRAN)プログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求める。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。
2. プログラム
C=================================================================C
SUBROUTINE SOR2D(A,B,NX,NY,X,Omega,EPS,ITER,ERR,IERR)
C=================================================================C
C Solve Ax=b by SOR with 2 dimensional FDM C
C Given Omega ( Acceleration factor ) C
C-----------------------------------------------------------------C
C A(0:NX,0:NY,5) R*8, In, A Coefficient Matrix 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 Omega R*8, In, SOR Acceleration factor 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,5), B(0:NX,0:NY)
DIMENSION X(0:NX,0:NY)
C Initial
IER = 0
BN = 0.0
do j=1,NY-1
do i=1,NX-1
BN = BN +B(i,j)**2
end do
end do
C Main Loop
do k=1,ITER
RN = 0.0
XN = 0.0
do j=1,NY-1
do i=1,NX-1
R = (B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
1 - A(i,j,4)*X(i+1,j) - A(i,j,5)*X(i,j+1) )
2 / A(i,j,3) - X(i,j)
X(i,j) = X(i,j) + Omega*R
RN = RN + (R*A(i,j,3))**2
end do
end do
C if(ERR <= EPS) Convergent
ERR = SQRT(RN/BN)
if(ERR.le.EPS) go to 100
end do
IERR = 1
C
100 continue
ITER = k
C
RETURN
END