program maximum
    IMPLICIT NONE
    include 'mpif.h'
    integer :: my_rank, num_procs, root
    integer :: rc
    !
    real(8) :: x, xmn=0.0, xmx=10.0, my_valmx, my_locmx, my_xmn, my_xmx
    integer, parameter :: MxFuncVals=10000
    real(8), dimension(:), allocatable :: fvals
    real(8), dimension(2) :: inmx, outmx
    integer :: i, j, my_imn, my_imx, per_proc, my_flen

    call MPI_Init(rc)
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, rc)
    call MPI_Comm_size(MPI_COMM_WORLD, num_procs, rc)
    !
    ! Define function and interval per each process
    !
    ! local dimensions:
    per_proc=(MxFuncVals+num_procs-1)/num_procs
    my_imn=1+my_rank*per_proc; my_imx=min(MxFuncVals,(my_rank+1)*per_proc)
    my_flen = (my_imx-my_imn+1)
    ! local array:
    allocate(fvals(my_flen),stat=rc)
    if ( rc .NE. 0 ) then
       print*,"Process ",my_rank," can't allocate FVALS(",my_flen,") !"
    end if
    ! function values:
    DO i = my_imn, my_imx
       x=xmn+(xmx-xmn)*(i-1)/(MxFuncVals-1)
       j = i-my_imn+1
       fvals(j) = (5.0-x)*sin(x*7.0)*cos(x*15.0)
    END DO
    !
    ! Compute global maximum value and location
    !
    ! maximum value and location per each process:
    my_valmx=maxval(fvals,my_flen)
    my_locmx=my_imn-1+maxloc(fvals,my_flen)
    inmx(1) = my_valmx; inmx(2) = my_locmx   ! copy into buffer
    root=0
    call MPI_Reduce( inmx, outmx, 1, MPI_2DOUBLE_PRECISION, MPI_MAXLOC, root, MPI_COMM_WORLD, rc );
    if ( rc .NE. 0 ) then
       print*,"For process ",my_rank," MPI_Reduce rc = ",rc
    end if
    if ( my_rank .EQ. root ) then
       i=outmx(2); x=xmn+(xmx-xmn)*(i-1)/(MxFuncVals-1)
       print*,"Maximum value is ",outmx(1)," in x = ",x
    end if
    print*," Bye bye from ",my_rank
    deallocate(fvals)
    call MPI_Finalize(rc)
    stop
end program maximum