C=================================================================C
C  Test Program of CG for 3-Dimensional FDM  with symmetric       C
C    for symmetric solution on Y-Z-axis                           C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
C     後 保範（東京工芸大学）                                     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-CG3DS.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
      do k=1,2
        WRITE(6,*) '----- 3D FDM with CG --------------------------'
        WRITE(1,*) '----- 3D FDM with CG---------------------------'
        if(k.eq.1) then
          WRITE(6,*) ' Symmetric solution CG for Y-axis'
          WRITE(1,*) ' Symmetric solution CG for Y-axis'
         else
          WRITE(6,*) ' Unsymmetric solution CG for Y-axis'
          WRITE(1,*) ' Unsymmetric solution CG for Y-axis'
        end if
        WRITE(6,100) NX,NY,NZ
        WRITE(1,100) NX,NY,NZ
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)
        F1 = 16
        F2 = 23
        if(k.eq.1) then
C   Symmetric solution
          CALL CG3DS(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
         else
C   Unsymmetric solution
          CALL CG3D2(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
        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
C  Output center-line
        CALL OUT3D(X,NX,NY,NZ,1)
      end do
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,   2007/06/04                     C
C        ( Tokyo Polytechnic 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,ID)
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         C
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  Output Un-symmetric Solution
      if(ID.eq.1) then
        WRITE(1,*) ' ** Un-symmetrix Solution **'
        WRITE(1,*) '  i  j  k     X(i,j,k)        X(i,NY-j,k)',
     1             '      X(i,j,NZ-k)       X(i,NY-j,NZ-k)'
        NO = 0
        do k=1,NZ/2-1
          do j=1,NY/2-1
            do i=1,NX/2-1
              XX = X(i,j,k)
              XL = X(i,NY-j,NZ-k)
              if(XX.ne.X(i,NY-j,k) .or. XX.ne.X(i,j,NZ-k) .or.
     1           XX.ne.X(i,NY-j,NZ-k) ) then 
                WRITE(1,200) i,j,k,XX,X(i,NY-j,k),X(i,j,NZ-k),XL
                NO = NO + 1
              end if
            end do
          end do
        end do
        WRITE(6,*) 'Un-symmetric solution numbers=',NO
        WRITE(1,*) 'Un-symmetric solution numbers=',NO
      end if
C
  100 FORMAT(1H ,/'Output on Y,Z-center line'/'  X-axis  Solution(X)')
  110 FORMAT(1H ,2F10.6)
  200 FORMAT(1H ,3I3,4(2X,Z16))
      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 Yasunori Ushiro,   2007/06/04                     C
C        ( Tokyo Polytechnic University )                         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