From e282f210f6032dd7a1707ddeb53037a67bf40562 Mon Sep 17 00:00:00 2001 From: ignis Date: Sat, 6 Jul 2019 02:43:46 +0900 Subject: [PATCH] split ludcmp alloc and calc, explicit real type size --- code/Compact.f90 | 178 +++++++++++++++++++++++++++++++++++------------ 1 file changed, 135 insertions(+), 43 deletions(-) diff --git a/code/Compact.f90 b/code/Compact.f90 index e75fcf1..6c4ff64 100644 --- a/code/Compact.f90 +++ b/code/Compact.f90 @@ -1,15 +1,12 @@ MODULE Compact IMPLICIT NONE - PRIVATE - REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & + REAL*8, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & lyf,lys,wyf,wys, & lzf,lzs,wzf,wzs ! lyzf,lyzs,wyzf,wyzs INTEGER :: nxc,nyc,nzc - REAL, PARAMETER :: ezero = 1.0e-14 + REAL*8, PARAMETER :: ezero = 1.0e-14 - PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp - CONTAINS SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp) @@ -20,6 +17,22 @@ nxc=nx nyc=ny nzc=nz + + CALL ludcmp_allocate(nx,ny,nz,xp,yp,zp) + + CALL ludcmp_calculate(nx,ny,nz,xp,yp,zp) + + END SUBROUTINE ludcmp + + SUBROUTINE ludcmp_allocate(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz + ! IF(nyc /= nzc) PRINT*,'ny should be equal nz' ! xp, yp, zp = 0 : periodic @@ -32,9 +45,6 @@ IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' ALLOCATE(wxs(nxc),STAT=ierr) IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' - CALL p_lud(1,nxc) - ELSE - CALL nonp_lud(1,nxc) ENDIF ALLOCATE(lyf(nyc),STAT=ierr) @@ -46,9 +56,6 @@ IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' ALLOCATE(wys(nyc),STAT=ierr) IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' - CALL p_lud(2,nyc) - ELSE - call nonp_lud(2,nyc) ENDIF ALLOCATE(lzf(nzc),STAT=ierr) @@ -60,19 +67,95 @@ IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' ALLOCATE(wzs(nzc),STAT=ierr) IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ENDIF + + END SUBROUTINE ludcmp_allocate + + SUBROUTINE ludcmp_deallocate(xp,yp,zp) + INTEGER, INTENT(IN) :: xp,yp,zp + +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + DEALLOCATE(lxf) + DEALLOCATE(lxs) + IF(xp.eq.0) THEN + DEALLOCATE(wxf) + DEALLOCATE(wxs) + ENDIF + + DEALLOCATE(lyf) + DEALLOCATE(lys) + IF(yp.eq.0) THEN + DEALLOCATE(wyf) + DEALLOCATE(wys) + ENDIF + + DEALLOCATE(lzf) + DEALLOCATE(lzs) + IF(zp.eq.0) THEN + DEALLOCATE(wzf) + DEALLOCATE(wzs) + ENDIF + + END SUBROUTINE ludcmp_deallocate + + SUBROUTINE ludcmp_testalloc + + IF (.not. ALLOCATED(lxf)) print *, "lxf not allocated" + IF (.not. ALLOCATED(lxs)) print *, "lxs not allocated" + IF (.not. ALLOCATED(wxf)) print *, "wxf not allocated" + IF (.not. ALLOCATED(wxs)) print *, "wxs not allocated" + + IF (.not. ALLOCATED(lyf)) print *, "lyf not allocated" + IF (.not. ALLOCATED(lys)) print *, "lys not allocated" + IF (.not. ALLOCATED(wyf)) print *, "wyf not allocated" + IF (.not. ALLOCATED(wys)) print *, "wys not allocated" + + IF (.not. ALLOCATED(lzf)) print *, "lzf not allocated" + IF (.not. ALLOCATED(lzs)) print *, "lzs not allocated" + IF (.not. ALLOCATED(wzf)) print *, "wzf not allocated" + IF (.not. ALLOCATED(wzs)) print *, "wzs not allocated" + + END SUBROUTINE ludcmp_testalloc + + SUBROUTINE ludcmp_calculate(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz + + CALL ludcmp_testalloc + +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + IF(xp.eq.0) THEN + CALL p_lud(1,nxc) + ELSE + CALL nonp_lud(1,nxc) + ENDIF + + IF(yp.eq.0) THEN + CALL p_lud(2,nyc) + ELSE + call nonp_lud(2,nyc) + ENDIF + + IF(zp.eq.0) THEN CALL p_lud(3,nzc) ELSE call nonp_lud(3,nzc) ENDIF -! CALL x_lud -! CALL yz_lud - - END SUBROUTINE ludcmp + END SUBROUTINE ludcmp_calculate SUBROUTINE nonp_lud(xyz,xx) INTEGER :: i,xyz,xx - REAL, DIMENSION(xx) :: aa + REAL*8, DIMENSION(xx) :: aa aa=3. aa(1)=0.5 ; aa(2)=4. aa(xx-1)=4. ; aa(xx)=0.5 @@ -92,7 +175,7 @@ SUBROUTINE p_lud(xyz,xx) INTEGER :: i,xyz,xx - REAL :: a + REAL*8 :: a a=3. ! first derivative IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction @@ -105,8 +188,8 @@ SUBROUTINE stdlu(a,n,l) INTEGER :: n - REAL :: a(n),l(n) - REAL :: d + REAL*8 :: a(n),l(n) + REAL*8 :: d INTEGER :: i l(1)=1.0/a(1) DO i=2,n @@ -117,9 +200,9 @@ SUBROUTINE ptdlu(a,n,l,w) INTEGER :: n - REAL :: a,l(n),w(n) + REAL*8 :: a,l(n),w(n) INTEGER :: i - REAL :: aa(n),d + REAL*8 :: aa(n),d DO i=1,n-1 aa(i)=a @@ -143,11 +226,14 @@ SUBROUTINE dfnonp(n,h,x,dx,nd,dir) INTEGER,INTENT(IN) :: n,nd,dir - REAL,INTENT(IN) :: h - REAL,INTENT(IN),DIMENSION(nd,n) :: x - REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx INTEGER :: i,j - REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + REAL*8 :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + + + ! print *, "dfnonp received (nd,n)", nd, n h1=1./h @@ -185,11 +271,14 @@ SUBROUTINE dfp(n,h,x,dx,nd,dir) INTEGER,INTENT(IN) :: n,nd,dir - REAL,INTENT(IN) :: h - REAL,INTENT(IN),DIMENSION(nd,n) :: x - REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx INTEGER :: i,j - REAL :: r1,r2,h1 + REAL*8 :: r1,r2,h1 + + + ! print *, "dfnonp received (nd,n)", nd, n h1=1./h r1=7./3. @@ -221,10 +310,10 @@ SUBROUTINE ptdslv(r,n,l,w,nd) INTEGER,INTENT(IN) :: n,nd - REAL,INTENT(INOUT),DIMENSION(nd,n) :: r - REAL,INTENT(IN),DIMENSION(:) :: l,w + REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL*8,INTENT(IN),DIMENSION(:) :: l,w INTEGER i,j - REAL, DIMENSION(nd) :: sum + REAL*8, DIMENSION(nd) :: sum DO j=1,nd sum(j)=w(1)*r(j,1) r(j,1)=r(j,1)*l(1) @@ -249,11 +338,13 @@ SUBROUTINE d2fp(n,h,x,dx,nd,dir) INTEGER,INTENT(IN) :: n,nd,dir - REAL,INTENT(IN) :: h - REAL,INTENT(IN),DIMENSION(nd,n) :: x - REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx INTEGER :: i,j - REAL :: h2,r1,r2,t1,t2 + REAL*8 :: h2,r1,r2,t1,t2 + + h2=1./(h*h) r1=6. r2=3./8. @@ -312,10 +403,10 @@ SUBROUTINE tdslv(r,n,l,nd) INTEGER,INTENT(IN) :: n,nd - REAL,INTENT(INOUT),DIMENSION(nd,n) :: r - REAL,INTENT(IN),DIMENSION(:) :: l + REAL*8,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL*8,INTENT(IN),DIMENSION(:) :: l INTEGER i,j - REAL t1 + REAL*8 t1 DO j=1,nd r(j,1)=r(j,1)*l(1) ENDDO @@ -334,11 +425,12 @@ SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) INTEGER,INTENT(IN) :: n,nd,dir - REAL,INTENT(IN) :: h - REAL,INTENT(IN),DIMENSION(nd,n) :: x - REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + REAL*8,INTENT(IN) :: h + REAL*8,INTENT(IN),DIMENSION(nd,n) :: x + REAL*8,INTENT(OUT),DIMENSION(nd,n) :: dx INTEGER :: i,j - REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2 + REAL*8 :: h2,r1,r2,r3,a,b,c,e,t1,t2 + h2=1./(h*h) r1=6.