3次元FDM用CG法(FORTRAN)テストプログラム
2003/08/26 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求めるプログラムのテスト用プログラム。
解析対象場がy,z軸対称なら、CG法による計算解xも対称性が保たれるかのテストプログラム。
2. プログラム
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, 2003/08/25 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-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,3
WRITE(6,*) '----- 3D FDM with CG --------------------------'
WRITE(1,*) '----- 3D FDM with CG---------------------------'
if(k.eq.1) then
WRITE(6,*) 'Modified symmetric solution CG for Y-axis'
WRITE(1,*) 'Modified symmetric solution CG for Y-axis'
else if(k.eq.2) 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 Modified symmetric solution
CALL CG3DM(A,B,NX,NY,NZ,X,R,P,Q,EPS,ITER,ERR,IERR)
else if(k.eq.2) 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, 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,ID)
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 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