!---------------------------------------------------------------
!****s* check_mpi_in_place
!  NAME
!    check_mpi_in_place
!  SYNOPSIS
      program check_mpi_in_place
!  PURPOSE
!    Checks the correctness of the implementation of the
!    'MPI_IN_PLACE' flag in the MPI implementation used by the
!    code. 
!  USES
        ! We only use an absolute minimum of modules here - that is, mpi_tasks .
        ! All needed variables (incl. STDOUT and STDERR) may be taken from there.
        implicit none
!  ARGUMENTS
      logical :: check_success
!  INPUTS
!    o none
!  OUTPUT
!    o check_success - indicates whether the 'MPI_IN_PLACE' flag works
!  AUTHOR
!    Volker Blum 
!  HISTORY
!    Volker Blum (July 2017)
!  SOURCE

    ! size of test data array
    integer :: n_data
      
    ! array that contains test data for MPI_IN_PLACE
    real*8, allocatable :: test_data(:)

    real*8, allocatable :: aux_test_data(:)
    
    logical :: aux_check_success, check_conventional_mpi
    
    integer :: i_data
    integer :: mpierr

    integer :: n_tasks
    integer :: myid
    integer :: mpi_comm_global
    
    integer :: STDOUT = 6

    include 'mpif.h'
    
    call MPI_INIT(mpierr)

    if (mpierr.ne.MPI_SUCCESS) then
         write(STDOUT,*) "Error in MPI initialization."
         write(STDOUT,*) "Exiting..."
         stop
    end if
      
    mpi_comm_global = MPI_COMM_WORLD

    call MPI_COMM_SIZE(mpi_comm_global, n_tasks, mpierr)
    call MPI_COMM_RANK(mpi_comm_global, myid, mpierr)

    ! initialize
    check_success = .true.
    
    n_data = 10
    
    allocate(test_data(n_data),stat=mpierr)

    ! seed test data array for allreduce call below
    if (myid.eq.0) then
       test_data(:) = 1.d0
    else
       test_data(:) = 0.d0
    end if

    ! Sum the test_data array over all MPI tasks
    call MPI_ALLREDUCE(MPI_IN_PLACE, &
         test_data(:), &
         n_data, &
         MPI_DOUBLE_PRECISION, &
         MPI_SUM, &
         mpi_comm_global, &
         mpierr )

    ! The value of all entries of test_data should now be 1.d0 on all MPI tasks.
    ! If that is not the case, then the MPI_IN_PLACE flag may be broken.
    ! (NB. Understood that this test can be done more elegantly. I spelled it out.)
    do i_data = 1, n_data, 1
       if (test_data(i_data).ne.1.d0) then
          check_success = .false.
       end if
    end do

    ! Synchronize check_success variable across all MPI tasks - if one is false,
    ! aux_check_success will be false. Obviously, do not use MPI_IN_PLACE .
    call MPI_ALLREDUCE(check_success, aux_check_success, 1, MPI_LOGICAL, &
         MPI_LAND, mpi_comm_global, mpierr )

    check_success = aux_check_success

    if (.not.check_success) then
       ! Not good news - MPI_IN_PLACE did not work.
       ! For good measure, check whether the intended synchronization works at all
       ! (even without MPI_IN_PLACE)

       if (myid.eq.0) then
          write(STDOUT,'(1X,A)') "* MPI_IN_PLACE does not appear to work as intended."
          write(STDOUT,'(1X,A)') "* Checking whether MPI_ALLREDUCE works at all."
       end if
       
       ! seed test data array for allreduce call below
       if (myid.eq.0) then
          test_data(:) = 1.d0
       else
          test_data(:) = 0.d0
       end if

       allocate(aux_test_data(n_data), stat=mpierr)

       aux_test_data(:) = 0.d0
       call MPI_ALLREDUCE(test_data(:), &
         aux_test_data(:), &
         n_data, &
         MPI_DOUBLE_PRECISION, &
         MPI_SUM, &
         mpi_comm_global, &
         mpierr )
       test_data(:) = aux_test_data(:)

       ! The value of all entries of test_data should now be 1.d0 on all MPI tasks.
       ! If that is not the case, then the entire MPI implementation may be broken.
       ! (NB. Understood that this test can be done more elegantly. I spelled it out.)
       check_conventional_mpi = .true.
       do i_data = 1, n_data, 1
          if (test_data(i_data).ne.1.d0) then
             check_conventional_mpi = .false.
          end if
       end do

       ! Synchronize check_conventional_mpi variable across all MPI tasks - if one is false,
       ! aux_check_success will be false. 
       call MPI_ALLREDUCE(check_conventional_mpi, aux_check_success, 1, MPI_LOGICAL, &
            MPI_LAND, mpi_comm_global, mpierr )

       check_conventional_mpi = aux_check_success

       if (.not.check_conventional_mpi) then
       ! This is unpleasant, since it would indicate a failure of a standard
       ! MPI_ALLREDUCE call. In this (unlikely) case, the code should be stopped.
          if (myid.eq.0) then
             write(STDOUT,'(1X,A)') "* It appears that a conventional MPI_ALLREDUCE call also does not work."
          end if
       else
          if (myid.eq.0) then
             write(STDOUT,'(1X,A)') "* Without MPI_IN_PLACE, MPI_ALLREDUCE appears to work."         
          end if
       end if
          
    end if

    if (allocated(test_data)) then
      deallocate(test_data)
    end if

    if (allocated(aux_test_data)) then
      deallocate(aux_test_data)
    end if

    call MPI_FINALIZE(mpierr)
    
  end program check_mpi_in_place
