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

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

1. 概要

 3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求めるプログラムのテスト用プログラム
CG法は比較的直交性に強いバージョンと演算量の少ないバージョンの両方をテストする。

2. プログラム

C=================================================================C
C  Test Program of CG for Dimensional FDM  with symmetric         C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2003/08/17                     C
C        ( Hitachi Ltd. and Waseda University )                   C
C=================================================================C
      PARAMETER(NDX=50, NDY=50, NDZ=50, NDXY=(NDX+1)*(NDY+1)*(NDZ+1))
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(4*NDXY), B(NDXY)
      DIMENSION X(NDXY), R(NDXY), P(NDXY), Q(NDXY)
C  Initial Data
      OPEN(UNIT=1,FILE='F-CG3D.txt') 
      CALL XCLOCK(T1,3)
      EPS = 1.0e-8
      WRITE(6,*) 'Type In  NX,NY,NZ'
      READ(5,*) NX,NY,NZ
      if(NX.gt.NDX) NX = NDX
      if(NY.gt.NDY) NY = NDY
      if(NZ.gt.NDZ) NZ = NDZ
      WRITE(6,100) NX,NY,NZ
      WRITE(1,100) NX,NY,NZ
      do k=1,2
C  Set A,B,X
        CALL SET3DS(A,B,X,NX,NY,NZ)
C  Solve Ax=b by CG
        ITER = NX*NY*NZ
        CALL XCLOCK(T1,5)
        if(k.eq.1) then
          CALL CG3D1(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
          F1 = 14
          F2 = 27
        else
          CALL CG3D2(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
          F1 = 16
          F2 = 23
        end if
        CALL XCLOCK(T2,5)
C  Write
        CPU   = T2 - T1
        FLOP  = (NX-1)*(NY-1)*(NZ-1)*(F1+F2*ITER)*1.0D-6
        FLOPS = 0.0
        if(CPU.ne.0.0) FLOPS = FLOP/CPU
        WRITE(6,200) k,ITER,ERR,CPU,FLOPS
        WRITE(1,200) k,ITER,ERR,CPU,FLOPS
      end do
C  Output center-line
      CALL OUT3D(X,NX,NY,NZ)
C
      STOP
  100 FORMAT(1H ,'CG Method for 3-Dimensional FDM,  NX,NY,NZ=',3I4,
     1 /,' Type  Loop   Error  Time(s)  MFLOPS')
  200 FORMAT(1H ,I3,I5,E10.2,F8.2,F8.1)
      END
C=================================================================C
      SUBROUTINE SET3DS(A,B,X,NX,NY,NZ)
C=================================================================C
C  Set A and B by 3 Dimensional FDM  with symmetric               C
C    -div(K.grad(X)) = f,   K=1.0, f=10,0 and 1x1 square          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)
C  Clear A,B,X
      do k=0,NZ
        do j=0,NY
          do i=0,NX
            A(i,j,k,1) = 0.0
            A(i,j,k,2) = 0.0
            A(i,j,k,3) = 0.0
            A(i,j,k,4) = 0.0
            B(i,j,k)   = 0.0
            X(i,j,k)   = 0.0
          end do
        end do
      end do
C  Initial Data
      DX = NX
      DY = NY
      DZ = NZ
      HX = 1.0/DX
      HY = 1.0/DY
      HZ = 1.0/DZ
      F  = 100.0
C  Set A,B
      do k=1,NZ-1
        do j=1,NY-1
          do i=1,NX-1
            A(i,j,k,1) = -HX*HY*DZ
            if(k.eq.1) A(i,j,k,1) = 0.0
            A(i,j,k,2) = -HZ*HX*DY
            if(j.eq.1) A(i,j,k,2) = 0.0
            A(i,j,k,3) = -HY*HZ*DX 
            if(i.eq.1) A(i,j,k,3) = 0.0
            A(i,j,k,4) = 2.0*(HY*HZ*DX+HZ*HX*DY+HX*HY*DZ)
            B(i,j,k)   = HX*HY*HZ*F
          end do
        end do
      end do
C
      RETURN
      END
C=================================================================C
      SUBROUTINE OUT3D(X,NX,NY,NZ)
C-----------------------------------------------------------------C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(0:NX,0:NY,0:NZ)
C  Output Y-Z center Line
      WRITE(1,100)
      DX = NX
      do i=0,NX
        XV = i/DX
        WRITE(1,110) XV,X(i,NY/2,NZ/2)
      end do
C
  100 FORMAT(1H ,/'Output on Y,Z-center line'/'  X-axis  Solution(X)')
  110 FORMAT(1H ,2F10.6)
      RETURN
      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