2次元FDM用CG法(FORTRAN)テストプログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式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=300, NDY=300, NDXY=(NDX+1)*(NDY+1))
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(3*NDXY), B(NDXY)
DIMENSION X(NDXY), R(NDXY), P(NDXY), Q(NDXY)
C Initial Data
OPEN(UNIT=1,FILE='F-CG2D.txt')
CALL XCLOCK(T1,3)
EPS = 1.0e-8
WRITE(6,*) 'Type In NX,NY'
READ(5,*) NX,NY
if(NX.gt.NDX) NX = NDX
if(NY.gt.NDY) NY = NDY
WRITE(6,100) NX,NY
WRITE(1,100) NX,NY
do k=1,2
C Set A,B,X
CALL SET2DS(A,B,X,NX,NY)
C Solve Ax=b by CG
ITER = NX*NY
CALL XCLOCK(T1,5)
if(k.eq.1) then
CALL CG2D1(A,B,NX,NY,X,R,P,Q,EPS,ITER,ERR,IERR)
F1 = 10
F2 = 23
else
CALL CG2D2(A,B,NX,NY,X,R,P,Q,EPS,ITER,ERR,IERR)
F1 = 12
F2 = 19
end if
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) k,ITER,ERR,CPU,FLOPS
WRITE(1,200) k,ITER,ERR,CPU,FLOPS
C Output center-line
CALL OUT2D(X,NX,NY,1)
end do
C
STOP
100 FORMAT(1H ,'CG Method for 2-Dimensional FDM, NX,NY=',2I4,
1 /,' Type Loop Error Time(s) MFLOPS')
200 FORMAT(1H ,I3,I5,E10.2,F8.2,F8.1)
END
C=================================================================C
SUBROUTINE SET2DS(A,B,X,NX,NY)
C=================================================================C
C Set A and B by 2 Dimensional FDM with symmetric 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,3) R*8, Out, A Coefficient Matrix C
C with symmetric 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,3), 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
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)
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