CoArray Fortran を試す。

暇つぶしに Fortran2008 で導入された並列化機能 CoArray Fortran を試してみました。

例として江戸時代の暴れん坊将軍の頃にオイラーが証明したというの値を、並列に和をとって求めることにします。

n個あるイメージのi番でiから始めてn個飛びに10**8個のデータを求めて足し合わせ、最後にマスタースレッドに当たるものでその総和を取ります。並列度が8程度なので、リダクションには素朴な方法を使いました。

この総和も収束がやや遅いことで有名で、ここでもまだ10^-9程度の誤差が残っています。実は加速法を使えばわずかな項数で収束します。

実行結果

ソース・プログラム

 program coarray
   implicit none
   real(8), parameter :: pi = 4.0d0 * atan(1.0d0)
   integer, parameter :: n = 10**8
!
   real(8), allocatable :: a(:)[:]
   real(8) :: s[*], sall
   integer :: i, npe, ime
   npe = num_images()
   ime = this_image()
   allocate( a(n)[*] )
   forall(i = 1:n) a(i) = 1.0d0 / (npe * dble(i - 1) + ime)**2
   s = sum(a)
   deallocate( a )
   print *, s, this_image()
   sync all
   
   sall = 0.0d0
   if (this_image() == 1) then
     sall = s
     do i = 2, num_images()
       sall = sall + s[i]
     end do
     print *, 'sum 1/i^2 [i = 1..', n * npe, '] pi^2 / 6 ', sall, pi**2 / 6.0d0, sall - pi**2 / 6.0d0
   end if
   stop
 end program coarray

別 Version

参考：http://www.csi.cuny.edu/cunyhpc/pdf/UPC_CAF.pdf

 program coarray
   implicit none
   real(8), parameter :: pi = 4.0d0 * atan(1.0d0)
   integer, parameter :: n = 10**8
!
   real(8), allocatable :: a(:)[:]
   real(8) :: s[*], t[*]
   integer :: i, npe, ime
   npe = num_images()
   ime = this_image()
   allocate( a(n)[*] )
   forall(i = 1:n) a(i) = 1.0d0 / (npe * dble(i - 1) + ime)**2
   s = sum(a)
   deallocate( a )
   print *, s, this_image()
   
   critical  
     t[1] = t[1] + s
   end critical 
   sync all

   if (this_image() == 1) print *,'sum 1/i^2 [i = 1..', n * npe, '] pi^2 / 6 ', t[1], pi**2 / 6.0d0, t[1] - pi**2 / 6.0d0
   stop
 end program coarray