Commit ffadaf89 authored by Marco Govoni's avatar Marco Govoni
Browse files

Added parallel interpolation

parent 95fe6f17
...@@ -17,7 +17,6 @@ MODULE fft_at_gamma ...@@ -17,7 +17,6 @@ MODULE fft_at_gamma
! Everything is done following dffts ! Everything is done following dffts
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE gvect, ONLY : nl,nlm
USE gvecs, ONLY : nls,nlsm USE gvecs, ONLY : nls,nlsm
USE fft_interfaces, ONLY : fwfft,invfft USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor USE fft_types, ONLY : fft_type_descriptor
......
...@@ -17,8 +17,6 @@ MODULE fft_at_k ...@@ -17,8 +17,6 @@ MODULE fft_at_k
! Everything is done following dffts ! Everything is done following dffts
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE gvect, ONLY : nl,nlm
USE gvecs, ONLY : nls,nlsm
USE fft_interfaces, ONLY : fwfft,invfft USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor USE fft_types, ONLY : fft_type_descriptor
! !
...@@ -32,6 +30,8 @@ MODULE fft_at_k ...@@ -32,6 +30,8 @@ MODULE fft_at_k
! !
! !
SUBROUTINE single_invfft_k(dfft,n,nx,a1,b,cdriver,igk) SUBROUTINE single_invfft_k(dfft,n,nx,a1,b,cdriver,igk)
!
USE gvecs, ONLY : nls
! !
! INVFFT : G ---> R ! INVFFT : G ---> R
! !
...@@ -73,6 +73,8 @@ MODULE fft_at_k ...@@ -73,6 +73,8 @@ MODULE fft_at_k
! !
! !
SUBROUTINE single_fwfft_k(dfft,n,nx,a,b1,cdriver,igk) SUBROUTINE single_fwfft_k(dfft,n,nx,a,b1,cdriver,igk)
!
USE gvecs, ONLY : nls
! !
! FWFFT : R ---> G ! FWFFT : R ---> G
! !
......
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
! This file is part of WEST. ! This file is part of WEST.
! !
! Contributors to this file: ! Contributors to this file:
! He Ma, Marco Govoni ! He Ma, Marco Govoni, Han Yang
! !
MODULE fourier_interpolation MODULE fourier_interpolation
! ----------------------------------------------------------------- ! -----------------------------------------------------------------
...@@ -17,23 +17,107 @@ MODULE fourier_interpolation ...@@ -17,23 +17,107 @@ MODULE fourier_interpolation
! !
CONTAINS CONTAINS
! !
SUBROUTINE get_G2R_mapping (n1, n2, n3, ng, ngx, ndim, nl) !SUBROUTINE get_G2R_mapping (n1, n2, n3, ng, ngx, ndim, nl)
! !
! ! INPUT : n1, n2, n3 = dimension of arbitrary R grid
! ! ng, ngx = actual number of PW
! ! ndim = dimensions of nl, 1: (n), 2:(n,-n)
! ! OUTPUT : nl = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
! !
! USE kinds, ONLY : DP
! USE gvect, ONLY : g
! USE cell_base, ONLY : at, bg, tpiba
! !
! IMPLICIT NONE
! !
! ! I/O
! !
! INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim
! INTEGER, INTENT(OUT) :: nl(ngx,ndim) ! mapping from [1, npw] to [1, n1*n2*n3]
! !
! ! Workspace
! !
! INTEGER :: ig, i_dim
! INTEGER :: h, k, l, hmax, kmax, lmax, hidx, kidx, lidx
! !
! nl = 0
! !
! ! Maximum Miller indexes associated to the R grid
! !
! hmax = (n1 - 1) / 2
! kmax = (n2 - 1) / 2
! lmax = (n3 - 1) / 2
! !
! DO i_dim = 1, ndim
! DO ig = 1, ng
! !
! ! Get the current Miller index
! !
! SELECT CASE (i_dim)
! CASE(1)
! h = NINT(SUM(g(:,ig) * at(:,1)))
! k = NINT(SUM(g(:,ig) * at(:,2)))
! l = NINT(SUM(g(:,ig) * at(:,3)))
! CASE(2)
! h = NINT(SUM(-g(:,ig) * at(:,1)))
! k = NINT(SUM(-g(:,ig) * at(:,2)))
! l = NINT(SUM(-g(:,ig) * at(:,3)))
! CASE DEFAULT
! CALL errore("","Bad mapping given: only 1 .OR. 2",1)
! END SELECT
! !
! IF( ABS(h) <= hmax .AND. ABS(k) <= kmax .AND. ABS(l) <= lmax ) THEN
! !
! ! derive position of G vector in (n1, n2, n3) grid
! !
! IF ( h >= 0 ) THEN
! hidx = h + 1
! ELSE
! hidx = h + n1 + 1
! ENDIF
! !
! IF ( k >= 0 ) THEN
! kidx = k + 1
! ELSE
! kidx = k + n2 + 1
! ENDIF
! !
! IF ( l >= 0 ) THEN
! lidx = l + 1
! ELSE
! lidx = l + n3 + 1
! ENDIF
! !
! nl(ig,i_dim) = hidx + (kidx - 1) * n1 + (lidx - 1) * n1 * n2
! !
! ENDIF
! !
! ENDDO
! ENDDO
! !
!END SUBROUTINE
!
!
SUBROUTINE set_nl (dfft, ng, ngx, ndim, nl, igk)
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : dfft = FFT descriptor
! ng, ngx = actual number of PW ! ng, ngx = actual number of PW, max
! ndim = dimensions of nl, 1: (n), 2:(n,-n) ! ndim = dimensions of nl, 1: (n), 2:(n,-n)
! OUTPUT : nl = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] ) ! OUTPUT : nl = computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE gvect, ONLY : g USE gvect, ONLY : g
USE cell_base, ONLY : at, bg, tpiba USE cell_base, ONLY : at, bg, tpiba
USE fft_types, ONLY : fft_type_descriptor
! !
IMPLICIT NONE IMPLICIT NONE
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim TYPE ( fft_type_descriptor ), INTENT(IN) :: dfft
INTEGER, INTENT(OUT) :: nl(ngx,ndim) ! mapping from [1, npw] to [1, n1*n2*n3] INTEGER, INTENT(IN) :: ng, ngx, ndim
INTEGER, INTENT(OUT) :: nl(ndim,ngx) ! mapping from [1, npw] to [1, n1*n2*n3]
INTEGER,INTENT(IN),OPTIONAL :: igk(ng)
! !
! Workspace ! Workspace
! !
...@@ -42,29 +126,44 @@ MODULE fourier_interpolation ...@@ -42,29 +126,44 @@ MODULE fourier_interpolation
! !
nl = 0 nl = 0
! !
! Maximum Miller indexes associated to the R grid ! Maximum Miller indexes associated to the R grid
! !
hmax = (n1 - 1) / 2 hmax = (dfft%nr1 - 1) / 2
kmax = (n2 - 1) / 2 kmax = (dfft%nr2 - 1) / 2
lmax = (n3 - 1) / 2 lmax = (dfft%nr3 - 1) / 2
! !
DO i_dim = 1, ndim DO i_dim = 1, ndim
DO ig = 1, ng DO ig = 1, ng
! !
! Get the current Miller index ! Get the current Miller index
! !
SELECT CASE (i_dim) IF(PRESENT(igk)) THEN
CASE(1) SELECT CASE (i_dim)
h = NINT(SUM(g(:,ig) * at(:,1))) CASE(1)
k = NINT(SUM(g(:,ig) * at(:,2))) h = NINT(SUM(g(:,igk(ig)) * at(:,1)))
l = NINT(SUM(g(:,ig) * at(:,3))) k = NINT(SUM(g(:,igk(ig)) * at(:,2)))
CASE(2) l = NINT(SUM(g(:,igk(ig)) * at(:,3)))
h = NINT(SUM(-g(:,ig) * at(:,1))) CASE(2)
k = NINT(SUM(-g(:,ig) * at(:,2))) h = NINT(SUM(-g(:,igk(ig)) * at(:,1)))
l = NINT(SUM(-g(:,ig) * at(:,3))) k = NINT(SUM(-g(:,igk(ig)) * at(:,2)))
CASE DEFAULT l = NINT(SUM(-g(:,igk(ig)) * at(:,3)))
CALL errore("","Bad mapping given: only 1 .OR. 2",1) CASE DEFAULT
END SELECT CALL errore("","Bad mapping given: only 1 .OR. 2",1)
END SELECT
ELSE
SELECT CASE (i_dim)
CASE(1)
h = NINT(SUM(g(:,ig) * at(:,1)))
k = NINT(SUM(g(:,ig) * at(:,2)))
l = NINT(SUM(g(:,ig) * at(:,3)))
CASE(2)
h = NINT(SUM(-g(:,ig) * at(:,1)))
k = NINT(SUM(-g(:,ig) * at(:,2)))
l = NINT(SUM(-g(:,ig) * at(:,3)))
CASE DEFAULT
CALL errore("","Bad mapping given: only 1 .OR. 2",1)
END SELECT
ENDIF
! !
IF( ABS(h) <= hmax .AND. ABS(k) <= kmax .AND. ABS(l) <= lmax ) THEN IF( ABS(h) <= hmax .AND. ABS(k) <= kmax .AND. ABS(l) <= lmax ) THEN
! !
...@@ -73,23 +172,28 @@ MODULE fourier_interpolation ...@@ -73,23 +172,28 @@ MODULE fourier_interpolation
IF ( h >= 0 ) THEN IF ( h >= 0 ) THEN
hidx = h + 1 hidx = h + 1
ELSE ELSE
hidx = h + n1 + 1 hidx = h + dfft%nr1 + 1
ENDIF ENDIF
! !
IF ( k >= 0 ) THEN IF ( k >= 0 ) THEN
kidx = k + 1 kidx = k + 1
ELSE ELSE
kidx = k + n2 + 1 kidx = k + dfft%nr2 + 1
ENDIF ENDIF
! !
IF ( l >= 0 ) THEN IF ( l >= 0 ) THEN
lidx = l + 1 lidx = l + 1
ELSE ELSE
lidx = l + n3 + 1 lidx = l + dfft%nr3 + 1
ENDIF ENDIF
! !
nl(ig,i_dim) = hidx + (kidx - 1) * n1 + (lidx - 1) * n1 * n2 IF (hidx>dfft%nr1 .OR. kidx>dfft%nr2 .OR. lidx>dfft%nr3) CALL errore('setnl','Mesh too small?',ig)
! !
#if defined (__MPI) && !defined (__USE_3D_FFT)
nl(i_dim,ig) = lidx + ( dfft%isind (hidx + (kidx - 1) * dfft%nr1x) - 1) * dfft%nr3x
#else
nl(i_dim,ig) = hidx + (kidx - 1) * dfft%nr1x + (lidx - 1) * dfft%nr1x * dfft%nr2x
#endif
ENDIF ENDIF
! !
ENDDO ENDDO
...@@ -97,150 +201,327 @@ MODULE fourier_interpolation ...@@ -97,150 +201,327 @@ MODULE fourier_interpolation
! !
END SUBROUTINE END SUBROUTINE
! !
SUBROUTINE single_interp_invfft_k(dfft,n,nx,a1,b,cdriver,nl)
!
USE kinds, ONLY : DP
USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor
!
! INVFFT : G ---> R
!
! INPUT : n = actual number of PW
! a1 = ONE COMPLEX arrays containing ONE COMPLEX functions in G space
! lda = leading dimension of a1
! ldb = leading dimension of b
! OUTPUT : b = ONE COMPLEX array containing ONE REAL functions in R space + 0
!
IMPLICIT NONE
!
! I/O
!
TYPE(fft_type_descriptor), INTENT(IN) :: dfft
INTEGER,INTENT(IN) :: n, nx
COMPLEX(DP),INTENT(IN) :: a1(nx)
COMPLEX(DP),INTENT(OUT) :: b(dfft%nnr)
CHARACTER(LEN=*),INTENT(IN) :: cdriver
INTEGER,INTENT(IN) :: nl(1,nx)
!
! Workspace
!
INTEGER :: ig
!
b=0.0_DP
!
DO ig=1,n
b(nl(1,ig))=a1(ig)
ENDDO
!
CALL invfft(cdriver,b,dfft)
!
END SUBROUTINE
!
! !
SUBROUTINE single_fwfft_fromArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk) SUBROUTINE single_interp_fwfft_k(dfft,n,nx,a,b1,cdriver,nl)
! !
! Note that the interpolation needs to be called by all processors within one band group USE kinds, ONLY : DP
USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor
! !
! FWFFT : R ---> G ! FWFFT : R ---> G
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : n = actual number of PW
! ng = actual number of PW ! a = ONE COMPLEX array containing ONE REAL functions in R space + 0
! ngx = leading dimendion for fg ! lda = leading dimension of a
! ndim = 1, 2 ! ldb = leading dimension of b1
! fr = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 ) DESTROYED ! OUTPUT : b1 = ONE COMPLEX array containing ONE COMPLEX functions in G space
! nl = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] ) !
! OUTPUT : fg = ONE COMPLEX array containing ONE functions in G space (note that the array is distributed ) IMPLICIT NONE
!
! I/O
!
TYPE(fft_type_descriptor), INTENT(IN) :: dfft
INTEGER,INTENT(IN) :: n,nx
COMPLEX(DP),INTENT(INOUT) :: a(dfft%nnr)
COMPLEX(DP),INTENT(OUT) :: b1(nx)
CHARACTER(LEN=*),INTENT(IN) :: cdriver
INTEGER,INTENT(IN) :: nl(1,nx)
!
! Workspace
!
INTEGER :: ig
!
CALL fwfft(cdriver, a, dfft)
! !
USE kinds, ONLY : DP b1=0.0_DP
USE fft_scalar, ONLY : cfft3d !
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp DO ig=1,n
USE mp, ONLY : mp_bcast b1(ig) = a(nl(1,ig))
ENDDO
!
END SUBROUTINE
!
!
SUBROUTINE single_interp_invfft_gamma(dfft,n,nx,a1,b,cdriver,nl)
!
USE kinds, ONLY : DP
USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor
!
! INVFFT : G ---> R
!
! INPUT : n = actual number of PW
! a1 = ONE COMPLEX arrays containing ONE COMPLEX functions in G space
! lda = leading dimension of a1
! ldb = leading dimension of b
! OUTPUT : b = ONE COMPLEX array containing ONE REAL functions in R space + 0
!
IMPLICIT NONE
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim TYPE(fft_type_descriptor), INTENT(IN) :: dfft
INTEGER, INTENT(IN) :: nl(ngx,ndim) INTEGER,INTENT(IN) :: n, nx
COMPLEX(DP), INTENT(INOUT) :: fr(n1*n2*n3) COMPLEX(DP),INTENT(IN) :: a1(nx)
COMPLEX(DP), INTENT(OUT) :: fg(ngx) COMPLEX(DP),INTENT(OUT) :: b(dfft%nnr)
INTEGER, INTENT(IN), OPTIONAL :: igk(ng) CHARACTER(LEN=*),INTENT(IN) :: cdriver
INTEGER,INTENT(IN) :: nl(2,nx)
! !
! Workspace ! Workspace
! !
INTEGER :: ig, idx INTEGER :: ig
! !
! FFT is serial... !
! !$OMP PARALLEL private(ig)
IF ( me_bgrp == 0 ) THEN !$OMP DO
! DO ig=1,dfft%nnr
CALL cfft3d( fr, n1, n2, n3, n1, n2, n3, 1, -1) b(ig)= (0.0_DP,0.0_DP)
! ENDDO
ENDIF !$OMP ENDDO
! !$OMP DO
! ... result is broadcast DO ig=1,n
! b(nl(1,ig))= a1(ig)
CALL mp_bcast( fr, 0, intra_bgrp_comm ) b(nl(2,ig))= DCONJG(a1(ig))
!
IF( PRESENT(igk) ) THEN
DO ig=1, ng
!
idx = nl(igk(ig),1)
IF (idx > 0) fg(ig) = fr(nl(igk(ig),1))
!
ENDDO
ELSE
DO ig=1, ng
!
idx = nl(ig,1)
IF (idx > 0) fg(ig) = fr(nl(ig,1))
!
ENDDO
ENDIF
!
DO ig = (ng+1), ngx
!
fg(ig) = (0.0_DP,0.0_DP)
!
ENDDO ENDDO
!$OMP ENDDO
!$OMP END PARALLEL
! !
ENDSUBROUTINE CALL invfft(cdriver, b, dfft)
!
END SUBROUTINE
! !
! !
SUBROUTINE single_invfft_toArbitraryRGrid (fr, n1, n2, n3, ng, ngx, ndim, nl, fg, igk) SUBROUTINE single_fwfft_gamma(dfft,n,nx,a,b1,cdriver,nl)
! !
! Note that the interpolation needs to be called by all processors within one band group USE kinds, ONLY : DP
USE fft_interfaces, ONLY : fwfft,invfft
USE fft_types, ONLY : fft_type_descriptor
! !
! INVFFT : R <--- G ! FWFFT : R ---> G
! !
! INPUT : n1, n2, n3 = dimension of arbitrary R grid ! INPUT : n = actual number of PW
! n = actual number of PW ! a = ONE COMPLEX array containing ONE REAL functions in R space + 0
! nx = leading dimendion for fg ! lda = leading dimension of a
! ndim = 1,2 ! ldb = leading dimension of b1
! fg = ONE COMPLEX array containing ONE function in R space (note that the array is not distributed, i.e. dimension = n1*n2*n3 ) ! OUTPUT : b1 = ONE COMPLEX array containing ONE COMPLEX functions in G space
! nl = pre-computed mapping from G to R space (i,e. from [1,n] to [1, n1*n2*n3] )
! OUTPUT : fr = ONE COMPLEX array containing ONE functions in G space (note that the array is distributed )
! !
USE kinds, ONLY : DP IMPLICIT NONE
USE fft_scalar, ONLY : cfft3d
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp
USE mp, ONLY : mp_bcast, mp_sum
! !
! I/O ! I/O
! !
INTEGER, INTENT(IN) :: n1, n2, n3, ng, ngx, ndim TYPE(fft_type_descriptor), INTENT(IN) :: dfft
INTEGER, INTENT(IN) :: nl(ng,ndim) INTEGER,INTENT(IN) :: n, nx
COMPLEX(DP), INTENT(OUT) :: fr(n1*n2*n3) COMPLEX(DP),INTENT(INOUT) :: a(dfft%nnr)
COMPLEX(DP), INTENT(IN):: fg(ngx) COMPLEX(DP),INTENT(OUT) :: b1(nx)
INTEGER, INTENT(IN), OPTIONAL :: igk(ng) CHARACTER(LEN=*),INTENT(IN) :: cdriver
INTEGER,INTENT(IN) :: nl(2,nx)
! !
! Workspace ! Workspace
! !
INTEGER :: ig, idx INTEGER :: ig
! !
fr = 0._DP !
! CALL fwfft(cdriver, a, dfft)
IF( PRESENT(igk) ) THEN ! Keep only G>=0
DO ig=1, ng !
! !$OMP PARALLEL private(ig)
idx = nl(igk(ig),1) !$OMP DO
IF (idx > 0) THEN DO ig=1,n
fr(nl(igk(ig),1)) = fg(ig) b1(ig) = a(nl(1,ig))
ENDIF ENDDO
! !$OMP ENDDO
ENDDO !$OMP END PARALLEL
ELSE !
DO ig=1, ng DO ig = (n+1), nx
! b1(ig) = (0.0_DP,0.0_DP)
idx = nl(ig,1) ENDDO
IF (idx > 0) THEN !
fr(nl(ig,1)) = fg(ig) END </