coarray でエラトステネスのふるい

#エラトステネスのふるいで素数を求める
先頃のトランプ大統領の有名なツイート covfefe は、「CAF をもっとやれと」ということだと天啓を受けました。ココアさんに倣って素数を数えたいと思います。
![DBQJbvYW0AUmwSh.jpg](https://qiita-image-store.s3.amazonaws.com/0/44522/d36d848a-3185-3f59-c840-aeef89eed7bb.jpeg)

coarray fortran (CAF) を利用して素数を求めます。

##シリアル実行
まずは、並列化する前のシリアル実行のプログラムを用意します。

サイズの小さい場合は以下の通りになります。数の表にチェックを入れるところは順序に依らないので、itab(i**2::i) = 0 のような三連符を用います。チェックの入らなかった数（すなわち素数）を取り出すのはフィルター関数に当たる pack 関数を用い、Fortran2003 の配列の代入による自動割り付けによって ip = pack(itab, itab /= 0) のように記述できます。 

```fortran
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 以下の素数の数
ここでは、素数の数を数えることに専念して、エラトステネスの表を１バイト論理型にして利用メモリーを減らします。

```fortran
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 関数を用意して、正負全域で同じ方向に切り下げ（切り上げ）する定義にしています。

```fortran
   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 ：素数の保存
次に素数を単にコンソールに表示するのではなく、ファイルに書き出して保存することを考えます。

各イメージが順番に、同じファイルの末尾にストリーム形式でダラダラと付け足してゆくことにします。

```fortran
   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 単精度で単精度整数が出力され死亡していましたｗ）

```fortran
    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
```