Homec4science
Diffusion SPClibs (feature/fft-split) Browse

spclibs/pputils2/pputils2.f90
03540012fedefeature/fft-split

Last Change
Raw File

pputils2.f90
View Options

!>
!> @file pputils2.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> spclibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
MODULE pputils2
!
! PPUTILS2: Some MPI utilities.
!
! T.M. Tran, CRPP-EPFL
! September 2010
! September 2013: add exchange, norm2
! November 2013: add timera, hostlist
!
USE iso_fortran_env, ONLY : rkind => real64
USE mpi
IMPLICIT NONE
PRIVATE
PUBLIC :: pptransp, dist1d, exchange, ppnorm2, timera, hostlist
!
INTERFACE pptransp
MODULE PROCEDURE pptransp2_r, pptransp3_r, pptransp4_r
MODULE PROCEDURE pptransp2_c, pptransp3_c, pptransp4_c
END INTERFACE
INTERFACE exchange
MODULE PROCEDURE exchange_2d, exchange_2d_new
END INTERFACE exchange
INTERFACE ppnorm2
MODULE PROCEDURE norm2_para_2d
END INTERFACE ppnorm2
!
CONTAINS
!=======================================================================
SUBROUTINE pptransp2_r(comm, a, b)
!
! Handles double precision-type matrices.
!
! Transpose of rank 2 matrix A:
! A(n1,n2/P) -> B(n2,n1/P)
!
INTEGER, INTENT(in) :: comm
REAL(rkind), DIMENSION(:,:), INTENT(in) :: a !*** dim dependant ***!
REAL(rkind), DIMENSION(:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, PARAMETER :: lastdim = 2, mpitype=MPI_DOUBLE_PRECISION !*** dim dependant ***!
REAL(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
!
INCLUDE 'pptransp2.tpl'
!
END SUBROUTINE pptransp2_r
!=======================================================================
SUBROUTINE pptransp2_c(comm, a, b)
!
! Same as pptransp2_r, but for double complex-type matrices.
!
INTEGER, INTENT(in) :: comm
COMPLEX(rkind), DIMENSION(:,:), INTENT(in) :: a !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, PARAMETER :: lastdim = 2, mpitype=MPI_DOUBLE_COMPLEX !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
!
INCLUDE 'pptransp2.tpl'
!
END SUBROUTINE pptransp2_c
!=======================================================================
SUBROUTINE pptransp3_r(comm, a, b, dim1, dim2)
!
! Handles double precision-type matrices.
!
! Transpose dimensions dim1 and dim2 of rank 3 matrix A.
! dim1 and dim2 are such that 1 <= dim1 < dim2 <= 3.
! At input, matrix A is partitioned along dimension dim2 of matrix A.
! At exit, B = transpose(A), and B is partitioned along dimension dim1 of matrix A.
!
! For example:
! dim1 = 1, dim2 = 2 : A(n1,n2/P,n3) -> B(n2,n1/P,n3)
! dim1 = 1, dim2 = 3 : A(n1,n2,n3/P) -> B(n3,n2,n1/P)
! dim1 = 2, dim2 = 3 : A(n1,n2,n3/P) -> B(n1,n3,n2/P)
!
INTEGER, INTENT(in) :: comm
REAL(rkind), DIMENSION(:,:,:), INTENT(in) :: a !*** dim dependant ***!
REAL(rkind), DIMENSION(:,:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, INTENT(in) :: dim1, dim2
INTEGER, PARAMETER :: lastdim = 3, mpitype=MPI_DOUBLE_PRECISION !*** dim dependant ***!
REAL(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
INTEGER :: recv_order(lastdim)
!
INCLUDE 'pptransp3.tpl'
!
END SUBROUTINE pptransp3_r
!=======================================================================
SUBROUTINE pptransp3_c(comm, a, b, dim1, dim2)
!
! Same as pptransp3_r, but for double complex-type matrices.
!
INTEGER, INTENT(in) :: comm
COMPLEX(rkind), DIMENSION(:,:,:), INTENT(in) :: a !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:,:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, INTENT(in) :: dim1, dim2
INTEGER, PARAMETER :: lastdim = 3, mpitype=MPI_DOUBLE_COMPLEX !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
INTEGER :: recv_order(lastdim)
!
INCLUDE 'pptransp3.tpl'
!
END SUBROUTINE pptransp3_c
!=======================================================================
SUBROUTINE pptransp4_r(comm, a, b, dim1, dim2)
!
! Handles double precision-type matrices.
!
! Transpose dimensions dim1 and dim2 of rank 4 matrix A.
! dim1 and dim2 are such that 1 <= dim1 < dim2 <= 4.
! At input, matrix A is partitioned along dimension dim2 of matrix A.
! At exit, B = transpose(A), and B is partitioned along dimension dim1 of matrix A.
!
! For example:
! dim1 = 1, dim2 = 2 : A(n1,n2/P,n3 ,n4 ) -> B(n2,n1/P,n3 ,n4 )
! dim1 = 1, dim2 = 3 : A(n1,n2 ,n3/P,n4 ) -> B(n3,n2 ,n1/P,n4 )
! dim1 = 1, dim2 = 4 : A(n1,n2 ,n3 ,n4/P) -> B(n4,n2 ,n3 ,n1/P)
! dim1 = 2, dim2 = 3 : A(n1,n2 ,n3/P,n4 ) -> B(n1,n3 ,n2/P,n4 )
! dim1 = 2, dim2 = 4 : A(n1,n2 ,n3 ,n4/P) -> B(n1,n4 ,n3 ,n2/P)
! dim1 = 3, dim2 = 4 : A(n1,n2 ,n3 ,n4/P) -> B(n1,n2 ,n4 ,n3/P)
!
INTEGER, INTENT(in) :: comm
REAL(rkind), DIMENSION(:,:,:,:), INTENT(in ) :: a !*** dim dependant ***!
REAL(rkind), DIMENSION(:,:,:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, INTENT(in) :: dim1, dim2
INTEGER, PARAMETER :: lastdim = 4, mpitype=MPI_DOUBLE_PRECISION !*** dim dependant ***!
REAL(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
INTEGER :: recv_order(lastdim)
!
INCLUDE 'pptransp4.tpl'
!
END SUBROUTINE pptransp4_r
!=======================================================================
SUBROUTINE pptransp4_c(comm, a, b, dim1, dim2)
!
! Same as pptransp4_r, but for double complex-type matrices
!
INTEGER, INTENT(in) :: comm
COMPLEX(rkind), DIMENSION(:,:,:,:), INTENT(in) :: a !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:,:,:,:), INTENT(out) :: b !*** dim dependant ***!
INTEGER, INTENT(in) :: dim1, dim2
INTEGER, PARAMETER :: lastdim = 4, mpitype=MPI_DOUBLE_COMPLEX !*** dim dependant ***!
COMPLEX(rkind), DIMENSION(:), ALLOCATABLE :: s_buf, r_buf
INTEGER :: recv_order(lastdim)
!
INCLUDE 'pptransp4.tpl'
!
END SUBROUTINE pptransp4_c
!=======================================================================
SUBROUTINE dist1d(comm, s0, ntot, s, nloc)
!
! 1d distribute ntot elements, returns offset s and local number of
! elements nloc.
!
INTEGER, INTENT(in) :: s0, ntot
INTEGER, INTENT(out) :: s, nloc
INTEGER :: comm, me, npes, ierr, naver, rem
!
CALL MPI_COMM_SIZE(comm, npes, ierr)
CALL MPI_COMM_RANK(comm, me, ierr)
naver = ntot/npes
rem = MODULO(ntot,npes)
s = s0 + MIN(rem,me) + me*naver
nloc = naver
IF( me.LT.rem ) nloc = nloc+1
!
END SUBROUTINE dist1d
!=======================================================================
SUBROUTINE exchange_2d_new(comm, u, garea)
!
! Exhange ghost cells with (west,east,south,north) neighbors.
! Assume same ghost cells on each dimension:
! garea(1) : number of ghost cells on west and east boundaries
! garea(2) : number of ghost cells on south and north boundaries
! Both are equal to 1 by default.
!
INTEGER, INTENT(in) :: comm
REAL(rkind), ALLOCATABLE, INTENT(inout) :: u(:,:)
INTEGER, OPTIONAL, INTENT(in) :: garea(2)
INTEGER :: neighs(4), ierr
!
CALL mpi_cart_shift(comm, 0, 1, neighs(1), neighs(2), ierr)
CALL mpi_cart_shift(comm, 1, 1, neighs(3), neighs(4), ierr)
CALL exchange_2d(comm, neighs, u, garea)
END SUBROUTINE exchange_2d_new
!=======================================================================
SUBROUTINE exchange_2d(comm, neighs, u, garea)
!
! Exhange ghost cells with (west,east,south,north) neighbors.
! Assume same ghost cells on each dimension:
! garea(1) : number of ghost cells on west and east boundaries
! garea(2) : number of ghost cells on south and north boundaries
! Both are equal to 1 by default.
!
INTEGER, INTENT(in) :: comm
INTEGER, INTENT(in) :: neighs(4)
REAL(rkind), ALLOCATABLE, INTENT(inout) :: u(:,:)
INTEGER, OPTIONAL, INTENT(in) :: garea(2)
!
INTEGER :: cols, rows
INTEGER :: ierr
INTEGER, PARAMETER :: ndim=2
INTEGER, DIMENSION(ndim) :: g, lb, ub, s, e, n
!
g = [1,1]
IF(PRESENT(garea)) g = garea
lb = LBOUND(u)
ub = UBOUND(u)
s = lb + g
e = ub - g
n = ub - lb + 1 ! include ghost cells
!
! g(2) matrix full rows with stride n(1)
CALL mpi_type_vector(n(2), g(2), n(1), MPI_DOUBLE_PRECISION, rows, ierr)
CALL mpi_type_commit(rows, ierr)
!
! g(1) contiguous matrix full columns
CALL mpi_type_contiguous(n(1)*g(1), MPI_DOUBLE_PRECISION, cols, ierr)
CALL mpi_type_commit(cols, ierr)
!
! Exchange along first dimension
CALL mpi_sendrecv(u(s(1), lb(2)), 1, rows, neighs(1), 0, &
& u(e(1)+1,lb(2)), 1, rows, neighs(2), 0, &
& comm, MPI_STATUS_IGNORE, ierr)
CALL mpi_sendrecv(u(e(1)-g(1)+1,lb(2)), 1, rows, neighs(2), 0, &
& u(lb(1), lb(2)), 1, rows, neighs(1), 0, &
& comm, MPI_STATUS_IGNORE, ierr)
!
! Exchange along second dimension
CALL mpi_sendrecv(u(lb(1),s(2)), 1, cols, neighs(3), 0, &
& u(lb(1),e(2)+1), 1, cols, neighs(4), 0, &
& comm, MPI_STATUS_IGNORE, ierr)
CALL mpi_sendrecv(u(lb(1),e(2)-g(2)+1), 1, cols, neighs(4), 0, &
& u(lb(1),lb(2)), 1, cols, neighs(3), 0, &
& comm, MPI_STATUS_IGNORE, ierr)
END SUBROUTINE exchange_2d
!=======================================================================
FUNCTION norm2_para_2d(x, comm, root, garea) RESULT(res)
!
! Vector norm of 2d distributed array with ghost cells
!
USE mpi
REAL(rkind), ALLOCATABLE, INTENT(in) :: x(:,:)
INTEGER, INTENT(in) :: comm
INTEGER, INTENT(in), OPTIONAL :: root
INTEGER, INTENT(in), OPTIONAL :: garea(:)
REAL(rkind) :: res
INTEGER, PARAMETER :: ndim=2
INTEGER, DIMENSION(ndim) :: g, s, e
REAL(rkind) :: res_loc
INTEGER :: r, me, ierr
!
CALL mpi_comm_rank(comm, me, ierr)
g = [1,1]
IF(PRESENT(garea)) g = garea
r = 0
IF(PRESENT(root)) r = root
s = LBOUND(x) + g
e = UBOUND(x) - g
res_loc = SUM(x(s(1):e(1),s(2):e(2))**2)
CALL mpi_reduce(res_loc, res, 1, MPI_DOUBLE_PRECISION, MPI_SUM, r, comm, ierr)
if(me.eq.r) res = SQRT(res)
END FUNCTION norm2_para_2d
!=======================================================================
SUBROUTINE timera(cntrl, str, eltime, comm)
!
! Timers (cntrl=0/1 to Init/Update)
!
USE mpi
INTEGER, INTENT(in) :: cntrl
CHARACTER(len=*), INTENT(in) :: str
DOUBLE PRECISION, OPTIONAL, INTENT(out) :: eltime
INTEGER, OPTIONAL, INTENT(in) :: comm
!
INTEGER, PARAMETER :: ncmax=128, maxlen=32
!
INTEGER, SAVE :: icall=0, nc=0
DOUBLE PRECISION, SAVE :: startt0=0.0
DOUBLE PRECISION, DIMENSION(ncmax), SAVE :: startt = 0.0, endt = 0.0
CHARACTER(len=maxlen), SAVE :: which(ncmax)
!
DOUBLE PRECISION, DIMENSION(ncmax) :: endtmin, endtmax
INTEGER :: comm0, me, lstr, found, i, ierr
!________________________________________________________________________________
IF(PRESENT(comm)) THEN
comm0 = comm
ELSE
comm0 = MPI_COMM_WORLD
END IF
CALL mpi_comm_rank(comm0, me, ierr)
CALL mpi_barrier(comm0, ierr)
!________________________________________________________________________________
!
IF( icall .EQ. 0 ) THEN
icall = icall+1
startt0 = mpi_wtime()
END IF
lstr = MIN(LEN_TRIM(str),maxlen)
IF( lstr .GT. 0 ) found = loc(str)
!________________________________________________________________________________
!
SELECT CASE (cntrl)
!
CASE(-1) ! Current wall time
IF( PRESENT(eltime) ) THEN
eltime = mpi_wtime() - startt0
ELSE IF (me .EQ. 0 ) THEN
WRITE(*,'(/a,a,1pe10.3/)') "++ ", ' Wall time used so far = ', &
& mpi_wtime() - startt0
END IF
!
CASE(0) ! Init Timer
IF( found .EQ. 0 ) THEN ! Called for the 1st time for 'str'
nc = nc+1
which(nc) = str(1:lstr)
found = nc
END IF
startt(found) = mpi_wtime()
!
CASE(1) ! Update timer
endt(found) = mpi_wtime() - startt(found)
IF( PRESENT(eltime) ) THEN
eltime = endt(found)
ELSE IF (me .EQ. 0 ) THEN
WRITE(*,'(/a,a,1pe10.3/)') "++ "//str, ' wall clock time = ', &
& endt(found)
END IF
!
CASE(2) ! Update and reset timer
endt(found) = endt(found) + mpi_wtime() - startt(found)
startt(found) = mpi_wtime()
IF( PRESENT(eltime) ) THEN
eltime = endt(found)
END IF
!
CASE(9) ! Display all timers
IF( nc .GT. 0 ) THEN
CALL mpi_reduce(endt, endtmin, nc, MPI_DOUBLE_PRECISION, MPI_MIN, 0, comm0, ierr)
CALL mpi_reduce(endt, endtmax, nc, MPI_DOUBLE_PRECISION, MPI_MAX, 0, comm0, ierr)
IF( me .EQ. 0 ) THEN
WRITE(*,'(a)') "Minmax Timer Summary"
WRITE(*,'(a)') "===================="
DO i=1,nc
WRITE(*,'(a20,2x,2(1pe12.3))') TRIM(which(i))//":", endtmin(i), endtmax(i)
END DO
END IF
END IF
!
END SELECT
!
CONTAINS
INTEGER FUNCTION loc(str)
CHARACTER(len=*), INTENT(in) :: str
INTEGER :: i, ind
loc = 0
DO i=1,nc
ind = INDEX(which(i), str(1:lstr))
IF( ind .GT. 0 .AND. LEN_TRIM(which(i)) .EQ. lstr ) THEN
loc = i
EXIT
END IF
END DO
END FUNCTION loc
END SUBROUTINE timera
!=======================================================================
SUBROUTINE hostlist(comm)
!
! Print list of hostnames in comm
!
USE mpi
INTEGER, OPTIONAL, INTENT(in) :: comm
!
CHARACTER(len=MPI_MAX_PROCESSOR_NAME) :: procname
CHARACTER(len=MPI_MAX_PROCESSOR_NAME), ALLOCATABLE :: procnames(:)
INTEGER :: comm0, me, nprocs, ierr, i, l
!
IF(PRESENT(comm)) THEN
comm0 = comm
ELSE
comm0 = mpi_comm_world
END IF
CALL MPI_COMM_RANK(comm0, me, ierr)
CALL MPI_COMM_SIZE(comm0, nprocs, ierr)
CALL MPI_GET_PROCESSOR_NAME(procname, l, ierr)
ALLOCATE(procnames(0:nprocs-1))
CALL mpi_gather(procname,MPI_MAX_PROCESSOR_NAME,mpi_character, &
& procnames,MPI_MAX_PROCESSOR_NAME,mpi_character,0, &
& comm0,ierr)
IF(me.EQ.0) THEN
WRITE(*,'(a/(10(1x,a)))') 'Host list:', &
& (TRIM(procnames(i)),i=0,nprocs-1)
END IF
END SUBROUTINE hostlist
!=======================================================================
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Private routines/functions !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE partners(comm, ids, idr)
!
! Compute ranks of send and receive procs.
!
IMPLICIT NONE
INTEGER, INTENT(in) :: comm
INTEGER, INTENT(out) :: ids(:), idr(:)
INTEGER :: me, npes, ierr, i
!
CALL mpi_comm_rank(comm, me, ierr)
CALL mpi_comm_size(comm, npes, ierr)
IF( ispower2(npes) ) THEN
DO i=0,npes-1
ids(i+1) = IEOR(me, i)
idr(i+1) = ids(i+1)
END DO
ELSE
DO i=0,npes-1
ids(i+1) = MODULO(me+i, npes)
idr(i+1) = MODULO(me-i, npes)
END DO
END IF
END SUBROUTINE partners
!=======================================================================
LOGICAL FUNCTION ispower2(n)
INTEGER, INTENT(in) :: n
INTEGER :: l
l=2
DO WHILE ( l .LT. n )
l = 2*l
END DO
ispower2 = l .EQ. n
END FUNCTION ispower2
!=======================================================================
END MODULE pputils2
c4science · Help