Pacheco section 4.2

以前の続きで、C の例題を MPI08 と CAF に、ゆるく翻案してみます。４章第２節の台形積分の例題をやってみます。

fortran66.hatenablog.com

元々が MPI の初歩用の例題なので計算の細かい所に気になる点がありますが、できるだけそのままゆるく翻案します。被積分関数は f(x) = x³ としました。積分範囲 [0,1] で解析的には 1/4=0.25 となるはずです。

MPI-08

program mpiex42
    use mpi_f08
    implicit none
    integer :: my_rank
    integer :: np
    real :: a = 0.0
    real :: b = 1.0
    integer :: n = 1024
    real :: h, x
    real :: local_a, local_b
    integer :: local_n
    real :: s, s_total
    integer :: isource
    integer :: idest = 0, itag = 0
    integer :: ierr
    type(MPI_Status) :: status

    ! Initialization
    call MPI_Init(ierr)

    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr)

    call MPI_Comm_size(MPI_COMM_WORLD, np, ierr)

    h = (b - a) / n
    local_n = n / np

    local_a =       a + my_rank * local_n * h
    local_b = local_a +           local_n * h
    s = Trap(local_a, local_b, local_n, h)

    if (my_rank == 0) then
        s_total = s
        do isource = 1, np - 1
            call MPI_Recv(s, 1, MPI_FLOAT, isource, itag, MPI_COMM_WORLD, status, ierr)
            s_total = s_total + s
        end do
    else
        call MPI_Send(s, 1, MPI_FLOAT, idest, itag, MPI_COMM_WORLD, ierr)
    end if
    !
    if (my_rank == 0) then
        print *, 'With n = ', np, ' trapezoids, our estimae'
        print *, 'of the integral from', a, ' to', b, ' =', s_total

    end if

    ! Finalization

    call MPI_Finalize(ierr)
contains

    pure real function f(x)
        import, none
        real, intent(in) :: x
        f = x**3
    end function f

    pure real function Trap(a, b, n, h)
        import, only : f
        real   , intent(in) :: a, b, h
        integer, intent(in) :: n
        real :: s, x
        integer :: i
        s = (f(a) + f(b)) / 2.0
        x = a
        do i = 1, n - 1
            x = x + h
            s = s + f(x)
        end do
        Trap = s * h
    end function Trap
end program mpiex42

dyna@dyna:~/CAF$ mpiifx mpiex42.f90
dyna@dyna:~/CAF$ mpirun -n 8 ./a.out
 With n =            8  trapezoids, our estimae
 of the integral from  0.0000000E+00  to   1.000000      =  0.2500002

CAF

program caf_ex42
     implicit none
     integer :: me, ne
     real :: a = 0.0
     real :: b = 1.0
     integer :: n = 1024
     real :: h, x
     real :: local_a, local_b
     integer :: local_n
     real :: s[*], s_total
     integer :: isource
     integer :: idest = 1

     me = this_image()
     ne = num_images()

     h = (b - a) / n
     local_n = n / ne

     local_a =       a + (me - 1) * local_n * h
     local_b = local_a +            local_n * h
     s = Trap(local_a, local_b, local_n, h)

     !  call co_sum(s, 1)! s_total <-- s
     if (me == 1) then
         s_total = s
         do isource = 2, ne
             s_total = s_total + s[isource]
         end do
     end if
     !
     if (me == 1) then
         print *, 'With n = ', ne, ' trapezoids, our estimate'
         print *, 'of the integral from', a, ' to', b, ' =', s_total
     end if
contains

    pure real function f(x)
        import, none
        real, intent(in) :: x
        f = x**3
    end function f

    pure real function Trap(a, b, n, h)
        import, only : f
        real   , intent(in) :: a, b, h
        integer, intent(in) :: n
        real :: s, x
        integer :: i
        s = (f(a) + f(b)) / 2.0
        x = a
        do i = 1, n - 1
            x = x + h
            s = s + f(x)
        end do
        Trap = s * h
    end function Trap
end program caf_ex42

dyna@dyna:~/CAF$ ifx cafex42.f90 -coarray
dyna@dyna:~/CAF$ ./a.out
 With n =           16  trapezoids, our estimate
 of the integral from  0.0000000E+00  to   1.000000      =  0.2500002