OOPSE/libmdtools/mpiSimulation_module.F90



!! MPI support for long range forces using force decomposition 
!! on a square grid of processors.
!! Corresponds to mpiSimunation.cpp for C++
!! mpi_module also contains a private interface for mpi f90 routines.
!!
!! @author Charles F. Vardeman II
!! @author Matthew Meineke
!! @version $Id: mpiSimulation_module.F90,v 1.9 2003-11-07 21:46:56 chuckv Exp $, $Date: 2003-11-07 21:46:56 $, $Name: not supported by cvs2svn $, $Revision: 1.9 $

module mpiSimulation  
  use definitions
#ifdef IS_MPI
  use oopseMPI
  implicit none
  PRIVATE


!! PUBLIC  Subroutines contained in this module
!! gather and scatter are a generic interface
!! to gather and scatter routines
  public :: gather, scatter
  public :: setupSimParallel
  public :: replanSimParallel
  public :: getNcol
  public :: getNrow
  public :: isMPISimSet
  public :: printComponentPlan
  public :: getMyNode

!! PUBLIC  Subroutines contained in MPI module
  public :: mpi_bcast
  public :: mpi_allreduce
!  public :: mpi_reduce
  public :: mpi_send
  public :: mpi_recv
  public :: mpi_get_processor_name
  public :: mpi_finalize

!! PUBLIC mpi variables
  public :: mpi_comm_world
  public :: mpi_character
  public :: mpi_integer
  public :: mpi_double_precision
  public :: mpi_sum
  public :: mpi_max
  public :: mpi_status_size
  public :: mpi_any_source


!! Safety logical to prevent access to ComponetPlan until
!! set by C++.
  logical, save :: ComponentPlanSet = .false.


!! generic mpi error declaration.
  integer, public :: mpi_err

#ifdef PROFILE
  public :: printCommTime

  real(kind = dp ), save :: commTime = 0.0_dp
  real(kind = dp ) :: commTimeInitial,commTimeFinal
#endif

!! Include mpiComponentPlan type. mpiComponentPlan is a 
!! dual header file for both c and fortran.
#define __FORTRAN90
#include "mpiComponentPlan.h"


!! Tags used during force loop for parallel simulation
  integer, allocatable, dimension(:) :: tagLocal
  integer, public, allocatable, dimension(:) :: tagRow
  integer, public, allocatable, dimension(:) :: tagColumn

!! Logical set true if mpiSimulation has been initialized
  logical, save :: isSimSet = .false.


  type (mpiComponentPlan), save :: mpiSim

!! gs_plan contains plans for gather and scatter routines
  type, public :: gs_plan
     private
     type (mpiComponentPlan), pointer :: gsComponentPlan => NULL()
     integer :: gsPlanSize !! size of this plan (nDim*nComponents)
     integer :: globalPlanSize !! size of all components in plan
     integer, dimension(:), pointer :: displs !! Displacements array for mpi indexed from 0.
     integer, dimension(:), pointer :: counts !! Counts array for mpi indexed from 0.
     integer :: myPlanComm  !! My communicator for this plan
     integer :: myPlanRank  !! My rank in this plan
     integer :: planNprocs  !! Number of processors in this plan
  end type gs_plan

! plans for different decompositions
  type (gs_plan), public, save :: plan_row
  type (gs_plan), public, save :: plan_row3d
  type (gs_plan), public, save :: plan_col
  type (gs_plan), public, save :: plan_col3d
  type (gs_plan),  public, save :: plan_row_Rotation
  type (gs_plan),  public, save :: plan_col_Rotation

  type (mpiComponentPlan), pointer :: simComponentPlan

! interface for gather scatter routines
!! Generic interface for gather.
!! Gathers an local array into row or column array
!! Interface provided for integer, double and double
!! rank 2 arrays.
  interface gather
     module procedure gather_integer
     module procedure gather_double
     module procedure gather_double_2d
  end interface

!! Generic interface for scatter.
!! Scatters a row or column array, adding componets
!! and reducing them to a local nComponent array.
!! Interface provided for double and double rank=2 arrays.

  interface scatter
     module procedure scatter_double
     module procedure scatter_double_2d
  end interface


contains

!! Sets up mpiComponentPlan with structure passed from C++.
  subroutine setupSimParallel(thisComponentPlan,ntags,tags,status)
!  Passed Arguments
!    integer, intent(inout) :: nDim !! Number of dimensions
    !! mpiComponentPlan struct from C
    type (mpiComponentPlan), intent(inout) :: thisComponentPlan
!! Number of tags passed, nlocal  
    integer, intent(in) :: ntags
!! Result status, 0 = normal, -1 = error
    integer, intent(out) :: status
    integer :: localStatus
!! Global reference tag for local particles
    integer, dimension(ntags),intent(inout) :: tags 

    write(*,*) 'mpiSim_mod thinks node', thisComponentPlan%myNode, ' has tags(1) = ', tags(1)


    status = 0
    if (componentPlanSet) then
       return
    endif
    componentPlanSet = .true.
    
 !! copy c component plan to fortran   
    mpiSim = thisComponentPlan 
    write(*,*) "Seting up simParallel"

    call make_Force_Grid(mpiSim,localStatus)
    if (localStatus /= 0) then
       write(default_error,*) "Error creating force grid"
       status = -1
       return
    endif

    call updateGridComponents(mpiSim,localStatus)
    if (localStatus /= 0) then
       write(default_error,*) "Error updating grid components"
       status = -1
       return
    endif
    

    !! initialize gather and scatter plans used in this simulation
    call plan_gather_scatter(1,mpiSim%myNlocal,&
         mpiSim,mpiSim%rowComm,plan_row)
    call plan_gather_scatter(nDim,mpiSim%myNlocal,&
         mpiSim,mpiSim%rowComm,plan_row3d)
    call plan_gather_scatter(9,mpiSim%myNlocal,&
         mpiSim,mpiSim%rowComm,plan_row_Rotation)
    call plan_gather_scatter(1,mpiSim%myNlocal,&
         mpiSim,mpiSim%columnComm,plan_col)
    call plan_gather_scatter(nDim,mpiSim%myNlocal,&
         mpiSim,mpiSim%columnComm,plan_col3d)
   call plan_gather_scatter(9,mpiSim%myNlocal,&
         mpiSim,mpiSim%columnComm,plan_col_Rotation)


!  Initialize tags    
    call setTags(tags,localStatus)
    if (localStatus /= 0) then
       status = -1
       return
    endif
    isSimSet = .true.

!    call printComponentPlan(mpiSim,0)
  end subroutine setupSimParallel

  subroutine replanSimParallel(thisComponentPlan,status)
!  Passed Arguments
    !! mpiComponentPlan struct from C
    type (mpiComponentPlan), intent(inout) :: thisComponentPlan   
    integer, intent(out) :: status
    integer :: localStatus
    integer :: mpierror
    status = 0

    call updateGridComponents(thisComponentPlan,localStatus)
    if (localStatus /= 0) then
       status = -1
       return
    endif
    
    !! Unplan Gather Scatter plans
    call unplan_gather_scatter(plan_row)
    call unplan_gather_scatter(plan_row3d)
    call unplan_gather_scatter(plan_row_Rotation)
    call unplan_gather_scatter(plan_col)
    call unplan_gather_scatter(plan_col3d)
    call unplan_gather_scatter(plan_col_Rotation)

    !! initialize gather and scatter plans used in this simulation
    call plan_gather_scatter(1,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%rowComm,plan_row)
    call plan_gather_scatter(nDim,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%rowComm,plan_row3d)
    call plan_gather_scatter(9,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%rowComm,plan_row_Rotation)
    call plan_gather_scatter(1,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%columnComm,plan_col)
    call plan_gather_scatter(nDim,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%rowComm,plan_col3d)
    call plan_gather_scatter(9,thisComponentPlan%myNlocal,&
         thisComponentPlan,thisComponentPlan%rowComm,plan_col_Rotation)


  end subroutine replanSimParallel

!! Updates number of row and column components for long range forces.
  subroutine updateGridComponents(thisComponentPlan,status)
    type (mpiComponentPlan) :: thisComponentPlan !! mpiComponentPlan

!! Status return
!! -  0 Success
!! - -1 Failure
    integer, intent(out) :: status 
    integer :: nComponentsLocal
    integer :: nComponentsRow = 0
    integer :: nComponentsColumn = 0
    integer :: mpiErrors

    status = 0
    if (.not. componentPlanSet) return

    if (thisComponentPlan%myNlocal == 0 ) then
       status = -1
       return
    endif
    
    nComponentsLocal = thisComponentPlan%myNlocal

    write(*,*) "UpdateGridComponents: myNlocal ", nComponentsLocal
    call mpi_allreduce(nComponentsLocal,nComponentsRow,1,mpi_integer,&
         mpi_sum,thisComponentPlan%rowComm,mpiErrors)
    if (mpiErrors /= 0) then
       status = -1
       return
    endif

    call mpi_allreduce(nComponentsLocal,nComponentsColumn,1,mpi_integer, &
         mpi_sum,thisComponentPlan%columnComm,mpiErrors)    
    if (mpiErrors /= 0) then
       status = -1
       return
    endif

    thisComponentPlan%nComponentsRow = nComponentsRow
    thisComponentPlan%nComponentsColumn = nComponentsColumn
    write(*,*) "UpdateGridComponents: myNRow ",&
         thisComponentPlan%nComponentsRow
    write(*,*) "UpdateGridComponents: myNColumn ",&
         thisComponentPlan%nComponentsColumn

  end subroutine updateGridComponents


!! Creates a square force decomposition of processors into row and column
!! communicators.
  subroutine make_Force_Grid(thisComponentPlan,status)
    type (mpiComponentPlan) :: thisComponentPlan
    integer, intent(out) :: status !! status returns -1 if error
    integer ::  nColumnsMax !! Maximum number of columns
    integer ::  nWorldProcessors !! Total number of processors in World comm.
    integer :: rowIndex !! Row for this processor.
    integer :: columnIndex !! Column for this processor.
    integer :: nRows !! Total number of rows.
    integer :: nColumns !! Total number of columns.
    integer :: mpiErrors !! MPI errors.
    integer :: rowCommunicator !! MPI row communicator.
    integer :: columnCommunicator
    integer :: myWorldRank
    integer :: i
    

    if (.not. ComponentPlanSet) return
    status = 0
  
!! We make a dangerous assumption here that if numberProcessors is
!! zero, then we need to get the information from MPI.
    if (thisComponentPlan%numberProcessors == 0 ) then 
       call mpi_comm_size( MPI_COMM_WORLD, nWorldProcessors,mpiErrors)
       if ( mpiErrors /= 0 ) then
          status = -1 
          return
       endif
       call mpi_comm_rank( MPI_COMM_WORLD,myWorldRank,mpiErrors)
       if ( mpiErrors /= 0 ) then
          status = -1 
          return
       endif

    else
       nWorldProcessors = thisComponentPlan%numberProcessors
       myWorldRank = thisComponentPlan%myNode
    endif


    nColumnsMax = nint(sqrt(real(nWorldProcessors,kind=dp)))

    do i = 1, nColumnsMax
       if (mod(nWorldProcessors,i) == 0) nColumns = i
    end do

    nRows = nWorldProcessors/nColumns

    rowIndex = myWorldRank/nColumns


    call mpi_comm_split(mpi_comm_world,rowIndex,0,rowCommunicator,mpiErrors)
    if ( mpiErrors /= 0 ) then
       write(default_error,*) "MPI comm split failed at row communicator"
       status = -1 
       return
    endif

    columnIndex = mod(myWorldRank,nColumns)
    call mpi_comm_split(mpi_comm_world,columnIndex,0,columnCommunicator,mpiErrors)
    if ( mpiErrors /= 0 ) then
       write(default_error,*) "MPI comm split faild at columnCommunicator"
       status = -1 
       return
    endif

! Set appropriate components of thisComponentPlan
    thisComponentPlan%rowComm = rowCommunicator
    thisComponentPlan%columnComm = columnCommunicator
    thisComponentPlan%rowIndex = rowIndex
    thisComponentPlan%columnIndex = columnIndex
    thisComponentPlan%numberRows = nRows
    thisComponentPlan%numberColumns = nColumns


  end subroutine make_Force_Grid


  !! initalizes a gather scatter plan
  subroutine plan_gather_scatter( nDim,nComponents,thisComponentPlan, &
       thisComm, this_plan,status)  
    integer, intent(in) :: nDim !! Number of dimensions for gather scatter plan
    integer, intent(in) :: nComponents
    type (mpiComponentPlan), intent(in), target :: thisComponentPlan
    type (gs_plan), intent(out) :: this_plan !! MPI Component Plan
    integer, intent(in) :: thisComm !! MPI communicator for this plan

    integer :: arraySize !! size to allocate plan for
    integer, intent(out), optional :: status
    integer :: ierror
    integer :: i,junk

    if (present(status)) status = 0
    
   
!! Set gsComponetPlan pointer 
!! to the componet plan we want to use for this gather scatter plan.
!! WARNING this could be dangerous since thisComponentPlan was origionally
!! allocated in C++ and there is a significant difference between c and 
!! f95 pointers....  
    this_plan%gsComponentPlan => thisComponentPlan

! Set this plan size for displs array.
    this_plan%gsPlanSize = nDim * nComponents

! Duplicate communicator for this plan
    call mpi_comm_dup(thisComm,this_plan%myPlanComm,mpi_err)
    if (mpi_err /= 0) then
       if (present(status)) status = -1
       return
    end if
    call mpi_comm_rank(this_plan%myPlanComm,this_plan%myPlanRank,mpi_err)
    if (mpi_err /= 0) then
       if (present(status)) status = -1
       return
    end if

    call mpi_comm_size(this_plan%myPlanComm,this_plan%planNprocs,mpi_err)

    if (mpi_err /= 0) then
       if (present(status)) status = -1
       return
    end if

    !! counts and displacements arrays are indexed from 0 to be compatable
    !! with MPI arrays.
    allocate (this_plan%counts(0:this_plan%planNprocs-1),STAT=ierror)
    if (ierror /= 0) then
       if (present(status)) status = -1
       return
    end if

    allocate (this_plan%displs(0:this_plan%planNprocs-1),STAT=ierror)
    if (ierror /= 0) then
       if (present(status)) status = -1
       return
    end if

   !! gather all the local sizes into a size # processors array.
    call mpi_allgather(this_plan%gsPlanSize,1,mpi_integer,this_plan%counts, &
         1,mpi_integer,thisComm,mpi_err)

    if (mpi_err /= 0) then
       if (present(status)) status = -1
       return
    end if
  

    !! figure out the total number of particles in this plan
    this_plan%globalPlanSize = sum(this_plan%counts)
  
 
    !! initialize plan displacements.
    this_plan%displs(0) = 0
    do i = 1, this_plan%planNprocs - 1,1
       this_plan%displs(i) = this_plan%displs(i-1) + this_plan%counts(i-1) 
    end do


  end subroutine plan_gather_scatter


  subroutine unplan_gather_scatter(this_plan)
    type (gs_plan), intent(inout) :: this_plan
    
    
    this_plan%gsComponentPlan => null()
    call mpi_comm_free(this_plan%myPlanComm,mpi_err)

    deallocate(this_plan%counts)
    deallocate(this_plan%displs)

  end subroutine unplan_gather_scatter

  subroutine gather_integer( sbuffer, rbuffer, this_plan, status)

    type (gs_plan), intent(inout) :: this_plan
    integer, dimension(:), intent(inout) :: sbuffer
    integer, dimension(:), intent(inout) :: rbuffer
    integer :: noffset
    integer, intent(out), optional :: status
    integer :: i


    if (present(status)) status = 0
    noffset = this_plan%displs(this_plan%myPlanRank)

!    if (getmyNode() == 1) then
!       write(*,*) "Node 0 printing allgatherv vars"
!       write(*,*) "Noffset: ", noffset
!       write(*,*) "PlanSize: ", this_plan%gsPlanSize
!       write(*,*) "PlanComm: ", this_plan%myPlanComm
!    end if

    call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_integer, &
         rbuffer,this_plan%counts,this_plan%displs,mpi_integer, &
         this_plan%myPlanComm, mpi_err)

    if (mpi_err /= 0) then
      if (present(status)) status  = -1
    endif

  end subroutine gather_integer

  subroutine gather_double( sbuffer, rbuffer, this_plan, status)

    type (gs_plan), intent(in) :: this_plan
    real( kind = DP ), dimension(:), intent(inout) :: sbuffer
    real( kind = DP ), dimension(:), intent(inout) :: rbuffer
    integer :: noffset
    integer, intent(out), optional :: status


    if (present(status)) status = 0
    noffset = this_plan%displs(this_plan%myPlanRank)
#ifdef PROFILE
    commTimeInitial = mpi_wtime()
#endif
    call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
         rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
         this_plan%myPlanComm, mpi_err)
#ifdef PROFILE
    commTimeFinal = mpi_wtime()
    commTime = commTime + commTimeFinal - commTimeInitial
#endif

    if (mpi_err /= 0) then
      if (present(status)) status  = -1
    endif

  end subroutine gather_double

  subroutine gather_double_2d( sbuffer, rbuffer, this_plan, status)

    type (gs_plan), intent(in) :: this_plan
    real( kind = DP ), dimension(:,:), intent(inout) :: sbuffer
    real( kind = DP ), dimension(:,:), intent(inout) :: rbuffer
    integer :: noffset,i,ierror
    integer, intent(out), optional :: status
    
    external  mpi_allgatherv
    
   if (present(status)) status = 0

!    noffset = this_plan%displs(this_plan%me)
#ifdef PROFILE
   commTimeInitial = mpi_wtime()
#endif

    call mpi_allgatherv(sbuffer,this_plan%gsPlanSize, mpi_double_precision, &
        rbuffer,this_plan%counts,this_plan%displs,mpi_double_precision, &
        this_plan%myPlanComm, mpi_err)

#ifdef PROFILE
    commTimeFinal = mpi_wtime()
    commTime = commTime + commTimeFinal - commTimeInitial
#endif

    if (mpi_err /= 0) then
      if (present(status)) status = -1
    endif

 end subroutine gather_double_2d

  subroutine scatter_double( sbuffer, rbuffer, this_plan, status)

    type (gs_plan), intent(in) :: this_plan
    real( kind = DP ), dimension(:), intent(inout) :: sbuffer
    real( kind = DP ), dimension(:), intent(inout) :: rbuffer
    integer, intent(out), optional :: status
    external mpi_reduce_scatter

   if (present(status)) status = 0

#ifdef PROFILE
   commTimeInitial = mpi_wtime()
#endif
    call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
         mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
#ifdef PROFILE
    commTimeFinal = mpi_wtime()
    commTime = commTime + commTimeFinal - commTimeInitial
#endif

    if (mpi_err /= 0) then
     if (present(status))  status = -1
    endif

  end subroutine scatter_double

  subroutine scatter_double_2d( sbuffer, rbuffer, this_plan, status)

    type (gs_plan), intent(in) :: this_plan
    real( kind = DP ), dimension(:,:), intent(inout) :: sbuffer
    real( kind = DP ), dimension(:,:), intent(inout) :: rbuffer
    integer, intent(out), optional :: status
    external mpi_reduce_scatter
 
   if (present(status)) status = 0
#ifdef PROFILE
   commTimeInitial = mpi_wtime()
#endif

    call mpi_reduce_scatter(sbuffer,rbuffer, this_plan%counts, &
         mpi_double_precision, MPI_SUM, this_plan%myPlanComm, mpi_err)
#ifdef PROFILE
    commTimeFinal = mpi_wtime()
    commTime = commTime + commTimeFinal - commTimeInitial
#endif

    if (mpi_err /= 0) then
      if (present(status)) status = -1
    endif

  end subroutine scatter_double_2d


  subroutine setTags(tags,status)
    integer, dimension(:) :: tags
    integer :: status

    integer :: alloc_stat
    
    integer :: ncol
    integer :: nrow

    status = 0
! allocate row arrays
    nrow = getNrow(plan_row)
    ncol = getNcol(plan_col)

    if (.not. allocated(tagRow)) then
       allocate(tagRow(nrow),STAT=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1 
          return
       endif
    else
       deallocate(tagRow)
       allocate(tagRow(nrow),STAT=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1 
          return
       endif

    endif
! allocate column arrays
    if (.not. allocated(tagColumn)) then
       allocate(tagColumn(ncol),STAT=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1 
          return
       endif
    else
       deallocate(tagColumn)
       allocate(tagColumn(ncol),STAT=alloc_stat)
       if (alloc_stat /= 0 ) then
          status = -1 
          return
       endif
    endif
    
    call gather(tags,tagRow,plan_row)
    call gather(tags,tagColumn,plan_col)

  end subroutine setTags

!  pure function getNcol(thisplan) result(ncol)
  function getNcol(thisplan) result(ncol)
    type (gs_plan), intent(in) :: thisplan
    integer :: ncol
    ncol = thisplan%gsComponentPlan%nComponentsColumn
  end function getNcol

!  pure function getNrow(thisplan) result(nrow)
  function getNrow(thisplan) result(nrow)
    type (gs_plan), intent(in) :: thisplan
    integer :: nrow
    nrow = thisplan%gsComponentPlan%nComponentsRow
  end function getNrow

  function isMPISimSet() result(isthisSimSet)
    logical :: isthisSimSet
    if (isSimSet) then 
       isthisSimSet = .true.
    else
       isthisSimSet = .false.
    endif
  end function isMPISimSet
  

  subroutine printComponentPlan(this_plan,printNode)

    type (mpiComponentPlan), intent(in) :: this_plan
    integer, optional :: printNode
    logical :: print_me = .false.

    if (present(printNode)) then
       if (printNode == mpiSim%myNode) print_me = .true.
    else
       print_me = .true.
    endif

    if (print_me) then
       write(default_error,*) "SetupSimParallel: writing component plan"
       
       write(default_error,*) "nMolGlobal: ", mpiSim%nMolGlobal
       write(default_error,*) "nAtomsGlobal: ", mpiSim%nAtomsGlobal
       write(default_error,*) "nBondGlobal: ", mpiSim%nBondsGlobal
       write(default_error,*) "nTorsionsGlobal: ", mpiSim%nTorsionsGlobal
       write(default_error,*) "nSRIGlobal: ", mpiSim%nSRIGlobal
       write(default_error,*) "myNlocal: ", mpiSim%myNlocal
       write(default_error,*) "myNode: ", mpiSim%myNode
       write(default_error,*) "numberProcessors: ", mpiSim%numberProcessors
       write(default_error,*) "rowComm: ", mpiSim%rowComm
       write(default_error,*) "columnComm: ", mpiSim%columnComm
       write(default_error,*) "numberRows: ", mpiSim%numberRows
       write(default_error,*) "numberColumns: ", mpiSim%numberColumns
       write(default_error,*) "nComponentsRow: ", mpiSim%nComponentsRow
       write(default_error,*) "nComponentsColumn: ", mpiSim%nComponentsColumn
       write(default_error,*) "rowIndex: ", mpiSim%rowIndex
       write(default_error,*) "columnIndex: ", mpiSim%columnIndex
    endif
  end subroutine printComponentPlan

  function getMyNode() result(myNode)
    integer :: myNode
    myNode = mpiSim%myNode
  end function getMyNode

#ifdef PROFILE
  subroutine printCommTime()

    write(*,*) "MPI communication time is: ", commTime

  end subroutine printCommTime
#endif

#endif // is_mpi
end module mpiSimulation


Revision:	858
Committed:	Fri Nov 7 21:46:56 2003 UTC (20 years, 8 months ago) by chuckv
File size:	22426 byte(s)
Log Message:	Added support for compiling fortran without use of mpich modules. We use mpif.h instead.: