3次元FDM用SOR法(FORTRAN)プログラム

    2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してSOR法で
反復解xを求める。
加速係数(ω)を与える必要がある。収束速度はωにより大きく左右される。

2. プログラム

C=================================================================C
      SUBROUTINE SOR3D(A,B,NX,NY,NZ,X,Omega,EPS,ITER,ERR,IERR)
C=================================================================C
C  Solve Ax=b by SOR with 3 dimensional FDM                       C
C    Given Omega ( Acceleration factor )                          C
C-----------------------------------------------------------------C
C    A(0:NX,0:NY,0:NZ,7)  R*8, In, A Coefficient Matrix           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    Omega             R*8, In,  SOR Acceleration factor          C
C    EPS              R*8, In,  if ||r||/||x|| <= EPS --> return  C
C    ITER              I*4, I/O, Number of Iteration              C
C    ERR               R*8, Out, ERR=||r||/||x||                  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,7), B(0:NX,0:NY,0:NZ)
      DIMENSION X(0:NX,0:NY,0:NZ)
C  Initial
      IER = 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
          end do
        end do
      end do
C  Main Loop
      do kk=1,ITER
        RN = 0.0
        do k=1,NZ-1
         do j=1,NY-1
          do i=1,NX-1
            R = (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,5)*X(i+1,j,k) - A(i,j,k,6)*X(i,j+1,k)
     3        - A(i,j,k,7)*X(i,j,k+1) )/A(i,j,k,4) - X(i,j,k)
            X(i,j,k) = X(i,j,k) + Omega*R
            RN       = RN + R**2
          end do
         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 = kk
C
      RETURN
      END