エラトステネスのふるいで素数を求める
先頃のトランプ大統領の有名なツイート covfefe は、「CAF をもっとやれと」ということだと天啓を受けました。ココアさんに倣って素数を数えたいと思います。
coarray fortran (CAF) を利用して素数を求めます。
シリアル実行
まずは、並列化する前のシリアル実行のプログラムを用意します。
サイズの小さい場合は以下の通りになります。数の表にチェックを入れるところは順序に依らないので、itab(i**2::i) = 0 のような三連符を用います。チェックの入らなかった数（すなわち素数）を取り出すのはフィルター関数に当たる pack 関数を用い、Fortran2003 の配列の代入による自動割り付けによって ip = pack(itab, itab /= 0) のように記述できます。
program eratos
integer, parameter :: n = 10**3
integer :: i, itab(n)
integer, allocatable :: ip(:)
itab = 0
forall (i = 2:n) itab(i) = i
do i = 2, int(sqrt(real(n)))
if (itab(i) /= 0) itab(i**2::i) = 0
end do
ip = pack(itab, itab /= 0)
print *, ip
print *, 'number of primes', size(ip)
end program eratos
実行結果
千以下の素数リストと個数
2 3 5 7 11 13
17 19 23 29 31 37
41 43 47 53 59 61
67 71 73 79 83 89
97 101 103 107 109 113
127 131 137 139 149 151
157 163 167 173 179 181
191 193 197 199 211 223
227 229 233 239 241 251
257 263 269 271 277 281
283 293 307 311 313 317
331 337 347 349 353 359
367 373 379 383 389 397
401 409 419 421 431 433
439 443 449 457 461 463
467 479 487 491 499 503
509 521 523 541 547 557
563 569 571 577 587 593
599 601 607 613 617 619
631 641 643 647 653 659
661 673 677 683 691 701
709 719 727 733 739 743
751 757 761 769 773 787
797 809 811 821 823 827
829 839 853 857 859 863
877 881 883 887 907 911
919 929 937 941 947 953
967 971 977 983 991 997
number of primes 168
スタック利用を回避する
しかしながら、pack 関数などは暗黙にスタック領域にテンポラリな配列を確保しているため、サイズを大きくしてゆくとスタックオーバーフローが起きます。スタックサイズを増やしてもいいのですが、それはみっともないのでやらないでおきます。また、itab(i**2::i) のような三連符による添え字は intel fortran の場合 32bit で計算しているらしく 32bit 整数の範囲を超えるような場合問題が生じます。
そこで、これらを明示的に書くことで、問題を回避します。
10**n 以下の素数の数
ここでは、素数の数を数えることに専念して、エラトステネスの表を１バイト論理型にして利用メモリーを減らします。
program eratos_main
implicit none
integer(8), parameter :: np = 10
integer(8) :: i
do i = 1, np
call eratos(i)
end do
contains
subroutine eratos(np)
integer(8), intent(in) :: np
integer(8) :: i, k, n
logical(1), allocatable :: itab(:)
n = 10**np
allocate(itab(n))
itab = .true.
itab(1) = .false.
do i = 2, int(sqrt(real(n)))
if (itab(i) /= 0) forall(integer(8)::k = i*i:n:i) itab(k) = .false.
end do
k = 0
do i = 1, n
if (itab(i)) k = k + 1
end do
print '(a, i2, a, i12)', '10^', np, ': number of primes', k
end subroutine eratos
end program eratos_main
実行結果
ここで求めた素数の数は既存の表の値と一致します。
10^ 1: number of primes 4
10^ 2: number of primes 25
10^ 3: number of primes 168
10^ 4: number of primes 1229
10^ 5: number of primes 9592
10^ 6: number of primes 78498
10^ 7: number of primes 664579
10^ 8: number of primes 5761455
10^ 9: number of primes 50847534
10^10: number of primes 455052511
coarray fortran 1 ：素朴な CAF 化
まず、素朴に coarray を用いて、エラトステネスの表を分割して並列にチェックを入れてゆくことにします。ここで、種となる素数は image 1 の範囲で求まると仮定しました。この条件は image の総数が、各イメージごとに分配される表のサイズを越えなければ満たされるので、まず問題ない条件だと思います。
各イメージごとに、表に等間隔にチェックを入れてゆくときの始まりの値ですが、それは koffset = -modulo(ioffset + 1, -k) + 1 で求めました。modulo 関数は mod 関数と少し挙動が異なります。mod 関数は 原点 0 に対して折り返すように対称になっていますが、これは負の数がからむ場合に数学的に扱いにくい定義です。そこで Fortran では、別に modulo 関数を用意して、正負全域で同じ方向に切り下げ（切り上げ）する定義にしています。
program CAFefe
implicit none
integer, parameter :: n = 10**2
integer :: k[*] ! coarray
integer :: i, nm, im, imax, ioffset, koffset, itab(n)
integer, allocatable :: ip(:)
nm = num_images()
im = this_image()
imax = int( sqrt(real(n)) * sqrt(real(nm)) )
if (imax > n) error stop 'too many images'
!
ioffset = n * (im - 1)
forall (i = 1:n) itab(i) = i + ioffset
if (im == 1) itab(1) = 0 ! 1 is not a prime number
do i = 1, imax
if (im == 1) then
if (itab(i) == 0) cycle
k = i
sync images(*) ! coarray
koffset = k**2
else
sync images(1) ! coarray
k = k[1] ! coarray
if (k == 0) exit
koffset = -modulo(ioffset + 1, -k) + 1
end if
itab(koffset::k) = 0
sync all ! coarray
end do
!
if (im == 1) then
k = 0
sync images(*) ! coarray
end if
! output results
ip = pack(itab, itab /=0)
if (im > 1) sync images(im - 1) ! coarray do from 1 to num_images()
print *, ip
if (im < nm) sync images(im + 1) ! coarray end do
end program CAFefe
実行結果
ここでは、実行時の指定パラメータで１０個のイメージを動かすことにして、各イメージが 100 個分の表を担当して、1000までの素数を求めてコンソールに出力しています。
2 3 5 7 11 13
17 19 23 29 31 37
41 43 47 53 59 61
67 71 73 79 83 89
97
101 103 107 109 113 127
131 137 139 149 151 157
163 167 173 179 181 191
193 197 199
211 223 227 229 233 239
241 251 257 263 269 271
277 281 283 293
307 311 313 317 331 337
347 349 353 359 367 373
379 383 389 397
401 409 419 421 431 433
439 443 449 457 461 463
467 479 487 491 499
503 509 521 523 541 547
557 563 569 571 577 587
593 599
601 607 613 617 619 631
641 643 647 653 659 661
673 677 683 691
701 709 719 727 733 739
743 751 757 761 769 773
787 797
809 811 821 823 827 829
839 853 857 859 863 877
881 883 887
907 911 919 929 937 941
947 953 967 971 977 983
991 997
coarray fortran 2 ：素数の保存
次に素数を単にコンソールに表示するのではなく、ファイルに書き出して保存することを考えます。
各イメージが順番に、同じファイルの末尾にストリーム形式でダラダラと付け足してゆくことにします。
program CAFefe
implicit none
integer, parameter :: ntab = 10**8
integer :: k[*] ! coarray
integer :: i, nm, im, imax, ioffset, koffset, ki, np
integer, allocatable :: ip(:), itab(:)
nm = num_images()
im = this_image()
imax = int( sqrt(real(ntab)) * sqrt(real(nm)) )
if (imax > ntab) error stop 'too many images'
!
allocate(itab(ntab))
ioffset = ntab * (im - 1)
forall (i = 1:ntab) itab(i) = i + ioffset
if (im == 1) itab(1) = 0 ! 1 is not a prime number
do i = 1, imax
if (im == 1) then
if (itab(i) == 0) cycle
k = i
sync images(*) ! coarray
koffset = k**2
else
sync images(1) ! coarray
k = k[1] ! coarray
if (k == 0) exit
koffset = -modulo(ioffset + 1, -k) + 1
end if
itab(koffset::k) = 0
sync all ! coarray
end do
!
if (im == 1) then
k = 0
sync images(*) ! coarray
open(9, file = 'primes.bin', access='stream') ! file preparation
endfile(9)
close(9)
end if
!
! ip = pack(itab, itab /=0)
block
integer :: i, j, m
m = count(itab /= 0)
allocate(ip(m))
j = 0
do i = 1, ntab
if (itab(i) == 0) cycle
j = j + 1
ip(j) = itab(i)
end do
end block
sync all
if (im > 1) sync images(im - 1) ! coarray do from 1 to num_images()
open(9, file = 'primes.bin', access='stream', position='append')
write(9) ip
close(9)
if (im < nm) sync images(im + 1) ! coarray end do
!
sync all
if (im == 1) then
open(9, file = 'primes.bin', access='stream')
inquire(9, size = np)
ki = storage_size(ip(1)) / 8 ! integer storage size in bits !f2008
np = np / ki
deallocate(ip)
allocate(ip(np)) ! coarray
read(9) ip
close(9)
! print '(10i10)', ip
print *, 'number of primes', np
end if
!
end program CAFefe
実行結果
ここでは１０個のイメージで、各１０＾８個の表をチェックしてゆくことにして、合計１０＾９までの素数の個数を求めています。
number of primes 50847534
**64bit での CAF
32bit 整数では、10^9 を少し超えたところまでしか計算できないので、10^10 を超えた計算には 64bit 整数を用いて計算する必要があります。
Fortran では、デフォルトの整数と単精度実数が同じビット幅であることが求められているので、今のところ 64bit 整数と32 bit 整数を混在させると微妙な不整合にであったりするので、注意が必要です。
ここでは、シリアル版の例に習って、１バイトの論理型配列を用いて素数の数だけを求めることにします。
（10*i i 単精度で単精度整数が出力され死亡していましたｗ）
module m_CAFefe
contains
subroutine eratos(ip)
implicit none
integer(8), intent(in) :: ip
integer(8) :: ntab
integer(8),save :: k[*], np[*] ! coarray
integer(8) :: i, j, nm, im, imax, ioffset, koffset, ki
logical(1), allocatable :: itab(:)
ntab = 10**ip
nm = num_images()
im = this_image()
imax = int( sqrt(real(ntab)) * sqrt(real(nm)) )
if (imax > ntab) error stop 'too many images'
!
allocate(itab(ntab))
ioffset = ntab * (im - 1)
itab = .true.
if (im == 1) itab(1) = .false. ! 1 is not a prime number
do i = 1, imax
if (im == 1) then
if (.not. itab(i)) cycle
k = i
sync images(*) ! coarray
koffset = k**2
else
sync images(1) ! coarray
k = k[1] ! coarray
if (k == 0) exit
koffset = -modulo(ioffset + 1, -k) + 1
end if
forall (integer(8)::j = koffset:ntab:k) itab(j) = .false.
sync all ! coarray
end do
!
if (im == 1) then
k = 0
sync images(*) ! coarray
np = 0
end if
sync all
critical
np[1] = np[1] + count(itab)
end critical
sync all
if (im == 1) then
print '(a,i2,a,i11)', '10^', (ip + 1), ' number of primes', np
end if
end subroutine eratos
end module m_CAFefe
program CAFefe
use m_CAFefe!
implicit none
integer(8) :: i
do i = 1, 9
call eratos(i)
end do
end program CAFefe
***実行結果
10^ 2 number of primes 25
10^ 3 number of primes 168
10^ 4 number of primes 1229
10^ 5 number of primes 9592
10^ 6 number of primes 78498
10^ 7 number of primes 664579
10^ 8 number of primes 5761455
10^ 9 number of primes 50847534
10^10 number of primes 455052511