2次元FDM用SOR法(FORTRAN)テストプログラム

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

1. 概要

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

2. プログラム

C=================================================================C
C  Test Program of SOR for Dimensional FDM                        C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2003/08/17                     C
C        ( Hitachi Ltd. and Waseda University )                   C
C=================================================================C
      PARAMETER(NDX=300, NDY=300, NDXY=(NDX+1)*(NDY+1))
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(5*NDXY), B(NDXY), X(NDXY)
C  Initial Data
      OPEN(UNIT=1,FILE='F-SOR2D.txt')
      CALL XCLOCK(T1,3)
      EPS = 1.0e-6
      WRITE(6,*) 'Type In  NX,NY,Omega'
      READ(5,*) NX,NY,Omega
      if(NX.gt.NDX) NX = NDX
      if(NY.gt.NDY) NY = NDY
      WRITE(6,100) NX,NY,Omega
      WRITE(1,100) NX,NY,Omega
C  Set A,B,X
      CALL SET2D(A,B,X,NX,NY)
C  Solve Ax=b by SOR
      ITER = 5*NX*NY
      CALL XCLOCK(T1,5)
      CALL SOR2D(A,B,NX,NY,X,Omega,EPS,ITER,ERR,IERR)
      F1 = 2
      F2 = 14
      CALL XCLOCK(T2,5)
C  Write
      CPU   = T2 - T1
      FLOP  = (NX-1)*(NY-1)*(F1+F2*ITER)*1.0D-6
      FLOPS = 0.0
      if(CPU.ne.0.0) FLOPS = FLOP/CPU
      WRITE(6,200) ITER,ERR,CPU,FLOPS
      WRITE(1,200) ITER,ERR,CPU,FLOPS
C  Output center-line
      CALL OUT2D(X,NX,NY,0)
C
      STOP
  100 FORMAT(1H ,'SOR Method for 2-Dimensional FDM,  NX,NY,Omega=',
     12I4,F6.3,/,'   Loop    Error   Time(s)  MFLOPS')
  200 FORMAT(1H ,I6,E10.2,F8.2,F8.1)
      END
C=================================================================C
      SUBROUTINE SET2D(A,B,X,NX,NY)
C=================================================================C
C  Set A and B by 2 Dimensional FDM                               C
C    -div(K.grad(X)) = f,   K=1.0, f=10,0 and 1x1 square          C
C-----------------------------------------------------------------C
C    A(0:NX,0:NY,5)  R*8, Out, A Coefficient Matrix               C
C    B(0:NX,0:NY) R*8, Out, A Right-hand Vector(b)                C
C    X(0:NX,0:NY) R*8, Out, Initial X, X=0                        C
C    NX           I*4, In,  Grid Numbers on X-axis                C
C    NY           I*4, In,  Grid Numbers on Y-axis                C
C    ID           I*4, In,  ID=0 ; NO output Un-Symmetric         C
C                             =1 ; Output Un-symmetric            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), X(0:NX,0:NY)
C  Clear A,B,X
      do j=0,NY
        do i=0,NX
          A(i,j,1) = 0.0
          A(i,j,2) = 0.0
          A(i,j,3) = 0.0
          A(i,j,4) = 0.0
          A(i,j,5) = 0.0
          B(i,j)   = 0.0
          X(i,j)   = 0.0
        end do
      end do
C  Initial Data
      DX = NX
      DY = NY 
      HX = 1.0/DX
      HY = 1.0/DY
      F  = 100.0
C  Set A,B
      do j=1,NY-1
        do i=1,NX-1
          A(i,j,1) = -HX*DY
          if(j.eq.1) A(i,j,1) = 0.0
          A(i,j,2) = -HY*DX
          if(i.eq.1) A(i,j,2) = 0.0
          A(i,j,3) = 2.0*(HY*DX+HX*DY)
          A(i,j,4) = -HY*DX
          if(i.eq.NX-1) A(i,j,4) = 0.0
          A(i,j,5) = -HX*DY
          if(j.eq.NY-1) A(i,j,5) = 0.0
          B(i,j)   = HX*HY*F
        end do
      end do
C
      RETURN
      END
C=================================================================C
      SUBROUTINE OUT2D(X,NX,NY,ID)
C-----------------------------------------------------------------C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(0:NX,0:NY)
C  Output Y
      WRITE(1,100)
      DX = NX
      do i=0,NX
        XV = i/DX
        WRITE(1,110) XV,X(i,NY/2)
      end do
C  Output Un-symmetric Solution
      if(ID.eq.1) then
        WRITE(1,*) ' ** Un-symmetrix Solution **'
        WRITE(1,*) '   i   j        X(i,j)         X(i,NY-j)'
        NO = 0
        do j=1,NY/2-1
          do i=1,NX/2-1
            if(X(i,j).ne.X(i,NY-j)) then
              WRITE(1,200) i,j,X(i,j),X(i,NY-j)
              NO = NO + 1
            end if   
          end do
        end do
        WRITE(6,*) 'Un-symmetric solution numbers=',NO
        WRITE(1,*) 'Un-symmetric solution numbers=',NO
      end if
C
      RETURN
  100 FORMAT(1H ,/'Output on Y-center line'/'  X-axis  Solution(X)')
  110 FORMAT(1H ,2F10.6)
  200 FORMAT(1H ,2I4,4(2X,Z16))
      END
C=================================================================C
      SUBROUTINE XCLOCK(CPU,ID)
C=================================================================C
C   CPU time Subroutine                                           C
C     CPU     R*8  Out, CPU Time                                  C
C     ID      I*4  In,  Dummy ( Same Hitachi FORTRAN XCLOCK )     C
C-----------------------------------------------------------------C
C  Written by Y. Ushiro ( Waseda and Hitachi )  2002.10.29        C
C=================================================================C
      REAL*8 CPU
      INTEGER*2 I1, I2, I3, I4
C
      IF(ID.GE.1) THEN
        CALL GETTIM(I1,I2,I3,I4)
        CPU = ( I1*60.0 + I2 )*60.0 + I3 + I4*0.01
      END IF
C
      RETURN
      END