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