!  Copyright (C) 2017 The authors of Py-ChemShell
!
!  This file is part of Py-ChemShell.
!
!  Py-ChemShell 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.
!
!  Py-ChemShell 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 Lesser General Public License for more details.
!
!  You should have received a copy of the GNU Lesser General Public
!  License along with Py-ChemShell.  If not, see
!  <http://www.gnu.org/licenses/>.

! WG = Workgroup
module ChemShellParallelModule

    use iso_c_binding
    use mpi
!    use IFPORT, only : MAX_HOSTNAM_LENGTH
    use DL_PY2F
    use ChemShellTaskfarmModule

    implicit none

!    include 'mpif.h'

    private

    ! public to libparallel.so procedures but invisible outside libparallel.so
    public initialise_parallel, finalise_parallel
    public create_python_object, synchronise_python_object, clean_python_object
    public getWG, getNode
    public synchronise_parallel
    public c_hostnames, c_world_hostnames, c_world_ranks
    public chemsh_parallel_intsize

    integer, parameter :: intType = sizeof(MPI_INTEGER)
    integer, parameter :: chemsh_parallel_intsize = sizeof(MPI_INTEGER)

    type parallelType

        private

        integer(kind=intType), public :: nprocs = 1
        integer(kind=intType), public :: nnodes_subset
        integer(kind=intType), public :: nnodes_wg
        integer(kind=intType), public :: rank
        integer(kind=intType), public :: nworkgroups = 1
        integer(kind=intType), public :: workgroup
        integer(kind=intType), public :: node_subset
        integer(kind=intType), public :: node
        integer(kind=intType), public :: comm_workgroup
        integer(kind=intType), public :: comm_node
        integer(kind=intType), public :: comm_subset
        integer(kind=intType), public :: ga_memory = 300
        integer(kind=intType), public :: debug  = 0
        integer(kind=intType), public :: null = MPI_COMM_NULL
        type(taskfarmType), pointer, private :: taskfarm => null()
        type(dictType)    , pointer, private :: PyObject => null()

        contains

        generic, public :: broadcast => private_broadcastScalar,    &
                                        private_broadcastOneDimInt, &
                                        private_broadcastOneDimLng, &
                                        private_broadcastOneDimDbl, &
                                        private_broadcastTwoDimInt, &
                                        private_broadcastTwoDimLng, &
                                        private_broadcastTwoDimDbl, &
                                        private_broadcastCPtr

        generic, public :: send      => private_sendScalar,    &
                                        private_sendOneDimInt, &
                                        private_sendOneDimLng, &
                                        private_sendOneDimDbl, &
                                        private_sendTwoDimInt, &
                                        private_sendTwoDimLng, &
                                        private_sendTwoDimDbl, &
                                        private_sendCPtr

        generic, public :: receive   => private_receiveScalar,    &
                                        private_receiveOneDimInt, &
                                        private_receiveOneDimLng, &
                                        private_receiveOneDimDbl, &
                                        private_receiveTwoDimInt, &
                                        private_receiveTwoDimLng, &
                                        private_receiveTwoDimDbl, &
                                        private_receiveCPtr

        generic, public :: gather    => private_gatherTwoDimDbl

        procedure,         public  :: init        => private_initialise
        procedure,         public  :: initialise  => private_initialise
        procedure,         public  :: finalise    => private_finalise
        procedure,         public  :: summarise   => private_summarise
        procedure,         public  :: sync        => private_synchronise
        procedure,         public  :: synchronise => private_synchronise
        procedure,         public  :: print       => private_print
        procedure,         public  :: is_root
        procedure,         public  :: is_subset
        procedure,         public  :: is_taskfarmed
        procedure, nopass, public  :: is_parallel
        procedure, nopass, public  :: exec_function
        procedure, nopass, public  :: get_function_handle
        procedure,         private :: private_broadcastScalar
        procedure,         private :: private_broadcastOneDimInt
        procedure,         private :: private_broadcastOneDimLng
        procedure,         private :: private_broadcastOneDimDbl
        procedure,         private :: private_broadcastTwoDimInt
        procedure,         private :: private_broadcastTwoDimLng
        procedure,         private :: private_broadcastTwoDimDbl
        procedure,         private :: private_broadcastCPtr
        procedure,         private :: private_sendScalar
        procedure,         private :: private_sendOneDimInt
        procedure,         private :: private_sendOneDimLng
        procedure,         private :: private_sendOneDimDbl
        procedure,         private :: private_sendTwoDimInt
        procedure,         private :: private_sendTwoDimLng
        procedure,         private :: private_sendTwoDimDbl
        procedure,         private :: private_sendCPtr
        procedure,         private :: private_receiveScalar
        procedure,         private :: private_receiveOneDimInt
        procedure,         private :: private_receiveOneDimLng
        procedure,         private :: private_receiveOneDimDbl
        procedure,         private :: private_receiveTwoDimInt
        procedure,         private :: private_receiveTwoDimLng
        procedure,         private :: private_receiveTwoDimDbl
        procedure,         private :: private_receiveCPtr
        procedure,         private :: private_gatherTwoDimDbl
        procedure,         private :: private_initialise
        procedure,         private :: private_finalise
        procedure,         private :: private_summarise
        procedure,         private :: private_synchronise
        procedure,         private :: private_print
        procedure,         private :: private_get_comm_from_level

    endtype parallelType

    ! the only thing Fortran application code can see from this API
    ! TODO: make this instance allocated in ChemShell rather than here to make
    !       the parallel module fully general-purpose
    type(parallelType), pointer    , public  :: chemsh_parallel

    real(kind=8)      , allocatable, private :: onedimdbl(:)    , twodimdbl(:,:)
    real(kind=8)      , allocatable, private :: onedimdblbuff(:), twodimdblbuff(:,:)

    character(len=:, kind=c_char), pointer :: c_hostnames, c_world_hostnames
    integer(kind=c_int)          , pointer :: c_world_ranks(:)

    contains

!!! API METHODS !!!
! public

    ! general-purpose API for Python
    integer(c_int) function create_python_object(objectPtr, onedimsize, twodimsizem, twodimsizen) result(ierror) bind(c)
!    integer(c_int) function create_python_object() result(ierror) bind(c)

        type(PyType)  , intent(in) :: objectPtr
        integer(c_int), intent(in) :: onedimsize
        integer(c_int), intent(in) :: twodimsizem
        integer(c_int), intent(in) :: twodimsizen

        ! result dictionary object
        if(.not.associated(chemsh_parallel%PyObject)) then
            allocate(chemsh_parallel%PyObject, source=ptr2dict(objectPtr))
        endif

        ! 1D array
        if(.not.allocated(onedimdbl).and.onedimsize.gt.0) then
            allocate(onedimdbl(onedimsize))
        endif
        if(.not.allocated(onedimdblbuff).and.onedimsize.gt.0) then
            allocate(onedimdblbuff(onedimsize*chemsh_parallel%nworkgroups))
        endif

        ! 2D array
        if(.not.allocated(twodimdbl).and.twodimsizem.gt.0.and.twodimsizen.gt.0) then
            allocate(twodimdbl(twodimsizem,twodimsizen))
        endif
        if(.not.allocated(twodimdblbuff).and.twodimsizem.gt.0.and.twodimsizen.gt.0) then
            allocate(twodimdblbuff(twodimsizem,twodimsizen*chemsh_parallel%nworkgroups))
        endif

        ierror = 0

    endfunction create_python_object
    ! FIXME: this is purpose-specific rather than API-like
    integer(c_int) function clean_python_object() result(ierror) bind(c)

        if(associated(chemsh_parallel%PyObject)) then
            deallocate(chemsh_parallel%PyObject)
        endif
        if(allocated(onedimdbl)) then
            deallocate(onedimdbl)
        endif
        if(allocated(onedimdblbuff)) then
            deallocate(onedimdblbuff)
        endif
        if(allocated(twodimdbl)) then
            deallocate(twodimdbl)
        endif
        if(allocated(twodimdblbuff)) then
            deallocate(twodimdblbuff)
        endif

        ierror = 0

    endfunction clean_python_object
    integer(c_int) function synchronise_python_object(objPtr, entry, attach) result(ierror) bind(c)

        ! A character dummy argument with length other than 1 is not interoperable
        type(PyType)                            , intent(in)           :: objPtr
        character(len=1, kind=c_char)           , intent(in)           :: entry(*)
        integer(c_int)                          , intent(in), optional :: attach

        integer(kind=intType)                                :: i, j, size, ierror_cast
        integer(kind=intType)        , parameter             :: root = 0
        type(dictType)               , pointer               :: PyObj
        character(len=16)                                    :: strbuff = ""
        integer                                              :: shp(2)
        integer                      , parameter             :: LENCHAR = 16
        integer(kind=8)              , allocatable           :: twodimint(:,:), twodimintbuff(:,:), toinserts(:,:)
        integer                      , allocatable           :: lpos(:), rpos(:)

        allocate(PyObj, source=ptr2dict(objPtr))

        ! convert c_char to Fortran character
        do i = 1, LENCHAR
            if(entry(i).ne.C_NULL_CHAR) then
                strbuff(i:i) = entry(i)
            else
                exit
            endif
        enddo
        shp = PyObj%enquireShape(strbuff)

        ! have to cast the type
        size = shp(1)*shp(2)

        allocate(twodimint(shp(1),shp(2)))
        allocate(toinserts(shp(1),shp(2)))
        allocate(twodimintbuff(shp(1),shp(2)*chemsh_parallel%nworkgroups))
        allocate(lpos(shp(2)))
        allocate(rpos(shp(2)))

        ! TODO: only integers so far
        call PyObj%get(trim(strbuff), twodimint)

        if(chemsh_parallel%is_taskfarmed()) then

            print '(3A)', " >>> merging ", trim(strbuff), " across workgroups..."

            ! workgroup 0: gather results from all other workgroups
            call MPI_GATHER(twodimint                , &
                            size                     , &
                            MPI_INTEGER8             , &
                            twodimintbuff            , &
                            size                     , &
                            MPI_INTEGER8             , &
                            root                     , &
                            chemsh_parallel%comm_node, &
                            ierror_cast)

            ! workgroup 0
            if(chemsh_parallel%rank.eq.0) then
                ! do maths
                if(present(attach)) then
                    ! attach to existing values, e.g.:
                    !     ( 1, 2,-1,-1,-1,...) merged with ( 0, 3,-1,-1,-1,...)
                    !     results in ( 1, 2, 0, 3,-1,...)
                    if(attach.eq.1) then
                        do i = 1, chemsh_parallel%nworkgroups-1
                            toinserts = twodimintbuff(1:shp(1),i*shp(2)+1:i*shp(2)+shp(2))
                            lpos = minloc(twodimint, 1)
                            rpos = minloc(toinserts, 1)
                            do j = 1, shp(2)
                                twodimint(lpos(j):lpos(j)+rpos(j)-1,j) = toinserts(:rpos(j)-1,j)
                            enddo
                        enddo
                    ! sum up
                    else
                        twodimint = twodimint + twodimintbuff(1:shp(1),i*shp(2)+1:i*shp(2)+shp(2))
                    endif
                ! sum up
                else
                    twodimint = twodimint + twodimintbuff(1:shp(1),i*shp(2)+1:i*shp(2)+shp(2))
                endif
            endif

            ! workgroup 0: broadcast to all workgroups
            call MPI_BCAST(twodimint                , &
                           size                     , &
                           MPI_INTEGER8             , &
                           root                     , &
                           chemsh_parallel%comm_node, &
                           ierror_cast)

            ! push to Python side
            if(chemsh_parallel%is_root()) then
                call PyObj%set(trim(strbuff), twodimint)
            endif

        endif

        deallocate(twodimint, toinserts, twodimintbuff)
        deallocate(lpos, rpos)

        ierror = ierror_cast

    endfunction synchronise_python_object
    integer(c_int) function initialise_parallel(nwgs, nnodes_subset, ga_memory, ChemShell_debug) &
                            result(ierror)                                                       &
                            bind(c, name="chemsh_parallel_f_initialise")

        integer(c_int), intent(in)           :: nwgs
        integer(c_int), intent(in)           :: nnodes_subset
        integer(c_int), intent(in)           :: ga_memory
        integer(c_int), intent(in), optional :: ChemShell_debug

        allocate(chemsh_parallel)

        if(present(ChemShell_debug)) then
            chemsh_parallel%debug = ChemShell_debug
        endif

        chemsh_parallel%nworkgroups = nwgs
        ! GA global memory (heap = stack = global/2): taskfarm%ga is private
        chemsh_parallel%ga_memory   = ga_memory

        ! number of processes in sub-group should not exceed total number of processes
        chemsh_parallel%nnodes_subset = nnodes_subset

        if(chemsh_parallel%initialise().ne.0) then
            print *, " >>> ERROR: failed initialising parallel module"
        endif

        ierror = chemsh_parallel%print()

        ! get a list of hostnames all over the world
        ierror = getWorldHostnames()

        ! get a list of hostnames within each workgroup
        ierror = getHostnames()

        ! get a list of world ranks within each workgroup
        ierror = getWorldRanks()

    endfunction initialise_parallel
    integer(c_int) function finalise_parallel() result(ierror) bind(c, name="chemsh_parallel_f_finalise")

        ierror = chemsh_parallel%finalise()

        deallocate(c_world_hostnames, c_hostnames, c_world_ranks)

        deallocate(chemsh_parallel)

    endfunction finalise_parallel
    integer(c_int) function broadcast(itype, buffer, lenbuf, root) result(ierror) bind(c, name="chemsh_parallel_f_broadcast")
    
        integer(c_int), intent(in) :: itype, lenbuf, root
        type(c_funptr), intent(in) :: buffer

        if(associated(chemsh_parallel)) then
            ierror = chemsh_parallel%broadcast(buffer, lenbuf, root, 0)
        else
            print *, " >>> ERROR: parallel moduled not initialised"
        endif

    endfunction broadcast
    integer(c_int) function synchronise_parallel() result(ierror) bind(c, name="chemsh_parallel_f_synchronise")

        ierror = chemsh_parallel%synchronise()
    
    endfunction synchronise_parallel
    ! myWorkgroup
    integer(c_int) function getWG() result(result) bind(c, name="chemsh_parallel_f_get_wg")
    
        result = chemsh_parallel%workgroup
    
    endfunction getWG
    ! myWorkgroupRank
    integer(c_int) function getNode() result(result) bind(c, name="chemsh_parallel_f_get_node")
    
        result = chemsh_parallel%node
    
    endfunction getNode
    ! number of world processes
    integer(c_int) function getNumProcesses() result(result) bind(c, name="chemsh_parallel_f_get_nprocs")
    
        result = chemsh_parallel%nprocs
    
    endfunction getNumProcesses
    ! number of workgroups
    integer(c_int) function getNumWorkgroups() result(result) bind(c, name="chemsh_parallel_f_get_nwgs")
    
        result = chemsh_parallel%nworkgroups
    
    endfunction getNumWorkgroups
    ! number of nodes in a workgroup
    integer(c_int) function getWGNumNodes() result(result) bind(c, name="chemsh_parallel_f_get_nnodes_wg")
    
        result = chemsh_parallel%nnodes_wg
    
    endfunction getWGNumNodes
    ! number of nodes in a workgroup's subset
    integer(c_int) function getSubsetNumNodes() result(result) bind(c, name="chemsh_parallel_f_get_nnodes_subset")
    
        result = chemsh_parallel%nnodes_subset
    
    endfunction getSubsetNumNodes

!!! PRIVATE METHODS !!!

    integer(kind=intType) function private_initialise(this) result(ierror)

        class(parallelType), intent(inout) :: this

        ! init
        call MPI_INIT(ierror)

        ! get nprocs
        call MPI_COMM_SIZE(MPI_COMM_WORLD, this%nprocs, ierror)

        ! get world rank
        call MPI_COMM_RANK(MPI_COMM_WORLD, this%rank  , ierror)

        allocate(this%taskfarm)

        ! GA global memory (heap = stack = global/2)
        this%taskfarm%ga_memory = this%ga_memory

        ierror = this%taskfarm%initialise(this%nworkgroups  , &
                                          this%nnodes_subset, &
                                          this%debug)

        this%workgroup      = this%taskfarm%colour
        this%node           = this%taskfarm%key
        this%node_subset    = this%taskfarm%key_subset
        ! YL 24/11/2018: corrected "this%nnodes_wg = this%taskfarm%ncolours" as "nkeys"
        this%nnodes_wg      = this%taskfarm%nkeys
        ! YL 24/11/2018: update nnodes_subset after initilisation
        this%nnodes_subset  = this%taskfarm%nsubsetkeys
        this%comm_workgroup = this%taskfarm%comm_colour
        this%comm_node      = this%taskfarm%comm_key
        this%comm_subset    = this%taskfarm%comm_subset
        this%null           = MPI_UNDEFINED

    endfunction private_initialise
    integer(kind=intType) function private_finalise(this) result(ierror)
    
        class(parallelType), intent(inout) :: this

!        call MPI_COMM_FREE(this%taskfarm%comm_subset, ierror)

        deallocate(this%taskfarm)

        call MPI_FINALIZE(ierror)
        
    endfunction private_finalise
    integer(kind=intType) function private_print(this) result(ierror)

        class(parallelType), intent(in) :: this

        ierror = this%taskfarm%print()

    endfunction private_print
    ! old way: bytewise broadcast of scalar value to workgroup
    integer(kind=intType) function private_broadcastCPtr(this, buffer, lenbuf, root, dummy, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        type(c_funptr)     , intent(in)           :: buffer
        class(*)           , intent(in)           :: lenbuf, root
        class(*)           , intent(in)           :: dummy
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastCPtr
    integer(kind=intType) function private_broadcastScalar(this, buffer, lenbuf, root, level) result(ierror)
    
        class(parallelType), intent(in)           :: this
        class(*)           , intent(in)           :: buffer, lenbuf, root
        character(len=*)   , intent(in), optional :: level
    
        integer(kind=intType)                       :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        ! Intel/GNU: there is no way to use iso_c_binding types (e.g., c_ptr) in selecttype procedure, otherwise "The type specification in TYPE IS or CLASS IS statement must not specify a sequence derived type or a type with the BIND attribute."
        selecttype(tmp=>buffer)
            class default
                call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)
        endselect
        call flush(6)

    endfunction private_broadcastScalar
    integer(kind=intType) function private_broadcastOneDimInt(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        integer(kind=4)    , intent(in)           :: buffer(:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastOneDimInt
    integer(kind=intType) function private_broadcastOneDimLng(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        integer(kind=8)    , intent(in)           :: buffer(:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastOneDimLng
    integer(kind=intType) function private_broadcastOneDimDbl(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        real(kind=8)       , intent(in)           :: buffer(:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastOneDimDbl
    integer(kind=intType) function private_broadcastTwoDimInt(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        integer(kind=4)    , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastTwoDimInt
    integer(kind=intType) function private_broadcastTwoDimLng(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        integer(kind=8)    , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastTwoDimLng
    integer(kind=intType) function private_broadcastTwoDimDbl(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        real(kind=8)       , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level

        integer(kind=intType)                       :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_BCAST(buffer, getCastInt(lenbuf), MPI_BYTE, getCastInt(root), communicator, ierror)

    endfunction private_broadcastTwoDimDbl
    integer(kind=intType) function private_sendCPtr(this, buffer, lenbuf, dest, tag, dummy, level) result(ierror)
 
        class(parallelType), intent(in) :: this
        type(c_funptr)     , intent(in) :: buffer
        class(*)           , intent(in) :: lenbuf, dest, tag
        class(*)           , intent(in) :: dummy
        character(len=*)   , intent(in), optional :: level
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendCPtr
    integer(kind=intType) function private_sendScalar(this, buffer, lenbuf, dest, tag) result(ierror)
    
        class(parallelType), intent(in) :: this
        class(*)           , intent(in) :: buffer
        class(*)           , intent(in) :: lenbuf, dest, tag
    
        ! Intel/GNU: there is no way to use iso_c_binding types (e.g., c_ptr) in selecttype procedure, otherwise "The type specification in TYPE IS or CLASS IS statement must not specify a sequence derived type or a type with the BIND attribute."
        selecttype(buffer)
            class default
        endselect

        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendScalar
    integer(kind=intType) function private_sendOneDimInt(this, buffer, lenbuf, dest, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=4)    , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, dest, tag
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendOneDimInt
    integer(kind=intType) function private_sendOneDimLng(this, buffer, lenbuf, dest, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=8)    , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, dest, tag
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendOneDimLng
    integer(kind=intType) function private_sendOneDimDbl(this, buffer, lenbuf, dest, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        real(kind=8)       , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, dest, tag
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendOneDimDbl
    integer(kind=intType) function private_sendTwoDimInt(this, buffer, lenbuf, dest, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=4)    , intent(in) :: buffer(:,:)
        class(*)           , intent(in) :: lenbuf, dest, tag
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendTwoDimInt
    integer(kind=intType) function private_sendTwoDimLng(this, buffer, lenbuf, dest, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=8)    , intent(in) :: buffer(:,:)
        class(*)           , intent(in) :: lenbuf, dest, tag
 
        call MPI_SEND(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(dest)   , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      ierror)

    endfunction private_sendTwoDimLng
    integer(kind=intType) function private_sendTwoDimDbl(this, buffer, lenbuf, dest, tag, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        real(kind=8)       , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, dest, tag
        character(len=*)   , intent(in), optional :: level

        integer(kind=intType)                     :: communicator
 
        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_SEND(buffer            , &
                      getCastInt(lenbuf), &
                      MPI_BYTE          , &
                      getCastInt(dest)  , &
                      getCastInt(tag)   , &
                      communicator      , &
                      ierror)

    endfunction private_sendTwoDimDbl
    integer(kind=intType) function private_receiveCPtr(this, buffer, lenbuf, source, tag, dummy, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        type(c_funptr)     , intent(in)           :: buffer
        class(*)           , intent(in)           :: lenbuf, source, tag
        class(*)           , intent(in)           :: dummy
        character(len=*)   , intent(in), optional :: level
 
        integer(kind=intType)                     :: communicator

        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        call MPI_RECV(buffer            , &
                      getCastInt(lenbuf), &
                      MPI_BYTE          , &
                      getCastInt(source), &
                      getCastInt(tag)   , &
                      communicator      , &
                      MPI_STATUS_IGNORE , &
                      ierror)

    endfunction private_receiveCPtr
    integer(kind=intType) function private_receiveScalar(this, buffer, lenbuf, source, tag) result(ierror)
    
        class(parallelType), intent(in) :: this
        class(*)           , intent(in) :: buffer
        class(*)           , intent(in) :: lenbuf, source, tag
    
        ! Intel/GNU: there is no way to use iso_c_binding types (e.g., c_ptr) in selecttype procedure, otherwise "The type specification in TYPE IS or CLASS IS statement must not specify a sequence derived type or a type with the BIND attribute."
        selecttype(buffer)
            class default
        endselect

        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveScalar
    integer(kind=intType) function private_receiveOneDimInt(this, buffer, lenbuf, source, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=4)    , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, source, tag
 
        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveOneDimInt
    integer(kind=intType) function private_receiveOneDimLng(this, buffer, lenbuf, source, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=8)    , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, source, tag
 
        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveOneDimLng
    integer(kind=intType) function private_receiveOneDimDbl(this, buffer, lenbuf, source, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        real(kind=8)       , intent(in) :: buffer(:)
        class(*)           , intent(in) :: lenbuf, source, tag
 
        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveOneDimDbl
    integer(kind=intType) function private_receiveTwoDimInt(this, buffer, lenbuf, source, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=4)    , intent(in) :: buffer(:,:)
        class(*)           , intent(in) :: lenbuf, source, tag
 
        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveTwoDimInt
    integer(kind=intType) function private_receiveTwoDimLng(this, buffer, lenbuf, source, tag) result(ierror)
 
        class(parallelType), intent(in) :: this
        integer(kind=8)    , intent(in) :: buffer(:,:)
        class(*)           , intent(in) :: lenbuf, source, tag
 
        call MPI_RECV(buffer             , &
                      getCastInt(lenbuf) , &
                      MPI_BYTE           , &
                      getCastInt(source) , &
                      getCastInt(tag)    , &
                      this%comm_workgroup, &
                      MPI_STATUS_IGNORE  , &
                      ierror)

    endfunction private_receiveTwoDimLng
    integer(kind=intType) function private_receiveTwoDimDbl(this, buffer, lenbuf, source, tag, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        real(kind=8)       , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, source, tag
        character(len=*)   , intent(in), optional :: level

        integer(kind=intType)                     :: communicator
 
        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif
        
        call MPI_RECV(buffer            , &
                      getCastInt(lenbuf), &
                      MPI_BYTE          , &
                      getCastInt(source), &
                      getCastInt(tag)   , &
                      communicator      , &
                      MPI_STATUS_IGNORE , &
                      ierror)

    endfunction private_receiveTwoDimDbl
    ! stub: not in use
    integer(kind=intType) function private_gatherTwoDimDbl(this, buffer, lenbuf, root, level) result(ierror)
 
        class(parallelType), intent(in)           :: this
        real(kind=8)       , intent(in)           :: buffer(:,:)
        class(*)           , intent(in)           :: lenbuf, root
        character(len=*)   , intent(in), optional :: level

        integer(kind=intType)                     :: communicator
 
        if(present(level)) then
            communicator = this%private_get_comm_from_level(level)
        else
            communicator = this%comm_workgroup
        endif

        ierror = 0       
 
!        call MPI_GATHER(buffer                   , &
!                        size                     , &
!                        MPI_DOUBLE_PRECISION     , &
!                        twodimdblbuff            , &
!                        size                     , &
!                        MPI_DOUBLE_PRECISION     , &
!                        getCastInt(root)         , &
!                        communicator             , &
!                        ierror_cast)
    endfunction private_gatherTwoDimDbl
    integer(kind=intType) function private_summarise(this, level) result(ierror)

        class(parallelType), intent(in) :: this
        integer            , intent(in) :: level

        ierror = this%taskfarm%summarise(level)

    endfunction private_summarise
    integer(kind=intType) function private_synchronise(this) result(ierror)
    
        class(parallelType), intent(in) :: this

        ierror = this%taskfarm%synchronise()
        
    endfunction private_synchronise
    integer(kind=intType) function private_get_comm_from_level(this, level) result(communicator)

        class(parallelType), intent(in) :: this
        character(len=*)   , intent(in) :: level

        if(trim(level).eq."global".or.trim(level).eq."world") then
            communicator = MPI_COMM_WORLD
        elseif(trim(level).eq."subworkgroup".or.trim(level).eq."subset") then
            communicator = this%comm_subset
        elseif(trim(level).eq."workgroup") then
            communicator = this%comm_workgroup
        elseif(trim(level).eq."node".or.trim(level).eq."counterpart") then
            communicator = this%comm_node
        else
            communicator = this%comm_workgroup
        endif

    endfunction private_get_comm_from_level
    character(len=6) function getWGDirname() result(iret)

        character(len=3) :: WG_index

        ! get directory name and create it
        write(WG_index, '(I0.3)') getWG()
        write(iret, '(A6)') "wg_" // WG_index

    endfunction getWGDirname
    ! will typically be called from a Python interpreter that belongs to a single work group
    integer(c_int) function getWorldRanks() result(iret) bind(c, name="chemsh_parallel_f_get_c_world_ranks")

        integer(kind=intType), parameter :: root = 0
        integer :: i, ierror
        integer                          :: status(MPI_STATUS_SIZE)

        allocate(c_world_ranks(chemsh_parallel%nnodes_wg))
        
        call MPI_COMM_RANK(MPI_COMM_WORLD, c_world_ranks(chemsh_parallel%node+1), ierror)
        if(chemsh_parallel%is_root()) then
            do i = 1, chemsh_parallel%nnodes_wg-1
                call MPI_RECV(c_world_ranks(i+1), 1, MPI_INT, i, 254, chemsh_parallel%taskfarm%comm_colour, status, ierror)
            enddo
        else
            call MPI_SEND(c_world_ranks(chemsh_parallel%node+1), &
                          1,                                     &
                          MPI_INT,                               &
                          root,                                  &
                          254,                                   &
                          chemsh_parallel%taskfarm%comm_colour,  &
                          ierror)
        endif

        iret = ierror

    endfunction getWorldRanks 
    ! will be called when the MPI is initialised and the saved string can be accessed from a Python interpreter that belongs to a single work group
    integer(c_int) function getWorldHostnames() result(iret) bind(c, name="chemsh_parallel_f_get_c_world_hostnames")

        character(len=MPI_MAX_PROCESSOR_NAME)              :: procname
        character(len=:)                     , allocatable :: allprocnames(:)
        integer(kind=intType)                , parameter   :: root = 0
        integer                                            :: i, ierror, lhostname, c_len
        integer                                            :: status(MPI_STATUS_SIZE)

        call MPI_GET_PROCESSOR_NAME(procname, lhostname, ierror)

        c_len = lhostname + 1

        allocate(character(len=lhostname) :: allprocnames(chemsh_parallel%nprocs))
        allocate(character(len=c_len*chemsh_parallel%nprocs+1, kind=c_char) :: c_world_hostnames)

        ! world root
        if(chemsh_parallel%rank.eq.0) then
            allprocnames(1) = procname(1:lhostname)
            do i = 1, chemsh_parallel%nprocs-1
                call MPI_RECV(allprocnames(i+1), lhostname, MPI_CHARACTER, i, 256, MPI_COMM_WORLD, status, ierror)
            enddo
        else
            call MPI_SEND(procname, lhostname, MPI_CHARACTER, root, 256, MPI_COMM_WORLD, ierror)
        endif

        call MPI_BARRIER(MPI_COMM_WORLD, IERROR)

        ! broadcast
        call MPI_BCAST(allprocnames                    , &
                       lhostname*chemsh_parallel%nprocs, &
                       MPI_CHARACTER                   , &
                       root                            , &
                       MPI_COMM_WORLD                  , &
                       ierror)

        ! formulate c_world_hostnames in all workgroups
        if(chemsh_parallel%is_root()) then
            do i = 1, chemsh_parallel%nprocs
                write(c_world_hostnames((i-1)*c_len+1:i*c_len-1), '(A)') allprocnames(i)
                write(c_world_hostnames( i   *c_len  :i*c_len  ), '(A)') trim(";")
            enddo
            c_world_hostnames(c_len*chemsh_parallel%nprocs+1:c_len*chemsh_parallel%nprocs+1) = C_NULL_CHAR
        endif

        deallocate(allprocnames)

        iret = ierror

    endfunction getWorldHostnames
    ! will be called when the MPI is initialised and the saved string can be accessed from a Python interpreter that belongs to a single work group
    integer(c_int) function getHostnames() result(iret) bind(c, name="chemsh_parallel_f_get_c_hostnames")

        character(len=MPI_MAX_PROCESSOR_NAME)              :: procname
        character(len=:)                     , allocatable :: allprocnames(:)
        integer(kind=intType)                , parameter   :: root = 0
        integer                                            :: i, ierror, lhostname, c_len
        integer                                            :: status(MPI_STATUS_SIZE)

        call MPI_GET_PROCESSOR_NAME(procname, lhostname, ierror)

        c_len = lhostname + 1

        allocate(character(len=lhostname) :: allprocnames(chemsh_parallel%nnodes_wg))
        allocate(character(len=c_len*chemsh_parallel%nnodes_wg+1, kind=c_char) :: c_hostnames)

        if(chemsh_parallel%is_root()) then
            allprocnames(1) = procname(1:lhostname)
            do i = 1, chemsh_parallel%nnodes_wg-1
                call MPI_RECV(allprocnames(i+1),                    &
                              lhostname,                            &
                              MPI_CHARACTER,                        &
                              i,                                    &
                              255,                                  &
                              chemsh_parallel%taskfarm%comm_colour, &
                              status,                               &
                              ierror)
            enddo
            do i = 1, chemsh_parallel%nnodes_wg
                write(c_hostnames((i-1)*c_len+1:i*c_len-1), '(A)') allprocnames(i)
                write(c_hostnames( i   *c_len  :i*c_len  ), '(A)') trim(";")
            enddo
            c_hostnames(c_len*chemsh_parallel%nnodes_wg+1:c_len*chemsh_parallel%nnodes_wg+1) = C_NULL_CHAR
        else
            call MPI_SEND(procname, lhostname, MPI_CHARACTER, root, 255, chemsh_parallel%taskfarm%comm_colour, ierror)
        endif

        deallocate(allprocnames)

        iret = ierror

    endfunction getHostnames
    character(len=10) function getWGOutName() result(iret)

!        iret = getWGDirname() // ".out"
        write(iret, '(A10)') trim(getWGDirname()) // ".out"

    endfunction getWGOutName
    character(len=10) function getWGErrName() result(iret)

!        iret = getWGDirname() // ".err"
        write(iret, '(A10)') trim(getWGDirname()) // ".err"

    endfunction getWGErrName
    character(len=14) function getBackupName(oldname) result(iret)

        character(*), intent(in) :: oldname

!        iret = trim(oldname) // ".OLD"
        write(iret, '(A14)') trim(oldname) // ".OLD"

    endfunction getBackupName
    integer(c_int) function get_function_handle(func_pointer) result(handle)
    
        interface
            integer(c_int) function chemsh_parallel_c_register_function(func_pointer) bind(c, name='chemsh_parallel_c_register_function')
                use iso_c_binding
                type(c_funptr), value :: func_pointer
            endfunction chemsh_parallel_c_register_function
        endinterface
      
        type(c_funptr), value :: func_pointer
      
        handle = chemsh_parallel_c_register_function(func_pointer)
    
    endfunction get_function_handle
    integer(c_int) function exec_function(handle, debug) result(ierror)
    
        interface
            integer(c_int) function chemsh_parallel_c_exec_function(handle, debug) bind(c, name='chemsh_parallel_c_exec_function')
                use iso_c_binding
                integer(c_int), value :: handle
                integer(c_int), value :: debug
            endfunction chemsh_parallel_c_exec_function
        endinterface
      
        integer(c_int), intent(in)           :: handle
        integer(c_int), intent(in), optional :: debug

        integer(c_int) :: c_debug
    
        if(.not.present(debug)) then
            c_debug = 0
        else
            c_debug = debug
        endif

        ierror = chemsh_parallel_c_exec_function(handle, c_debug)
    
    endfunction exec_function
    logical function is_taskfarmed(this) result(result)

        class(parallelType), intent(in) :: this

        result = this%nworkgroups.gt.1

    endfunction is_taskfarmed
    ! query build
    logical function is_parallel() result(result)

        result = .true.
    
    endfunction is_parallel
    logical function is_root(this) result(result)

        class(parallelType), intent(in) :: this

        if(associated(this%taskfarm)) then
            result = this%taskfarm%is_root()
        else
            print *, ">>> ChemShell ERROR: parallelism not yet initialised!"
            stop
        endif

    endfunction is_root
    logical function is_subset(this) result(result)

        class(parallelType), intent(in) :: this

        if(associated(this%taskfarm)) then
            result = this%taskfarm%is_subset()
        else
            print *, ">>> CHEMSH_PARALLEL ERROR: MPI not yet initialised!"
            stop
        endif

    endfunction is_subset
    ! convert any kind=4 or kind=8 integer to intType
    integer(kind=intType) function getCastInt(input) result(result)

        class(*), intent(in) :: input

        selecttype(tmp=>input)
            type is(integer(kind=4))
                result = int(tmp, kind=4)
            type is(integer(kind=8))
                result = int(tmp, kind=8)
!            type is(integer(c_int))
!                result = tmp
        endselect
        
    endfunction getCastInt

endmodule ChemShellParallelModule

