CoArray Fortran で Mandelbrot (再々)

CoArray Fortran で Mandelbrot

以前もやりましたが、しつこくやります。
https://qiita.com/cure_honey/items/3af353390fe913ae86a9

本質的には変わりはありませんが、Fortran 2008 の submodule を使ったりしています。またcoarray の使い方を少し変えています。Mandelbrot 図形計算も、すこし変えてあります。

parameterized type は使い方が難しく、OO の type-bound procedure にし難く、coarray とのつながりも分かりにくいです。parameterized type は、型種 KIND を実行時パラメータに出来て動的に変えられる点で他に無い機能を持ちますが、配列の大きさを変えるのは、素直に allocate でやった方がいいかも知れません。まだ判断つきかねます。

実行結果

intel fortran ver. 18.0.1 で実行しています。submodule は、まだバグ？が残っている感じです。optional 変数を認識しないとか、type-bound procedure を procedure :: wr => wr_bmp のようにポインタで結ぼうとすると、Linker がうまく見つけられなくて Linker エラーが出ます。

image 1 writing..
image 2 writing..
image 3 writing..
image 4 writing..
image 5 writing..
image 6 writing..
image 7 writing..
image 8 writing..
続行するには何かキーを押してください . . .

ソース・プログラム

fortran2008 の新知識をありったけぶち込んでみました。明示的に coarray 命令を使っているところにはコメントをつけておきました。抜けてたらごめんなさいｗ　

追記:文法上の警告を減らすようにニ三か所修正しました。内容的には変化ありません。

    module m_bmp
      implicit none
      private
      public :: t_bmp, rd_bmp, to_rgb
      type :: t_rgb
        sequence
        character :: b = char(0), g = char(0), r = char(0)
      end type t_rgb 

      type :: t_bmp(nx, ny)
        integer, len :: nx = 0, ny = 0  
        type (t_rgb) :: rgb(nx, ny)
      contains 
        generic :: wr => wr_bmp
        procedure :: wr_bmp
        procedure :: wr_header
      end type  
!
      interface assignment(=)
        module procedure :: assign_bmp
      end interface assignment(=) 

      interface t_bmp
        module procedure :: construct_bmp
      end interface t_bmp
!
      interface  
        module function construct_bmp(nx, ny) result(res)
          integer, intent(in) :: nx, ny
          type (t_bmp(:, :)), allocatable :: res
        end function construct_bmp    

        module subroutine assign_bmp(dest, src) 
          type (t_bmp(:, :)), intent(out), allocatable :: dest
          type (t_bmp(*, *)), intent(in ) :: src
        end subroutine assign_bmp

        module subroutine wr_bmp(bmp, fn)
          class (t_bmp(*, *)), intent(in) :: bmp
          character (len = *), intent(in) :: fn
        end subroutine wr_bmp

        module subroutine wr_header(bmp, fn)
          class (t_bmp(*, *)), intent(in) :: bmp
          character (len = *), intent(in) :: fn
        end subroutine wr_header

        module function rd_bmp(fn) result(bmp)
          character (len = *), intent(in) :: fn
          type (t_bmp(:, :)), allocatable :: bmp
        end function rd_bmp
      end interface  
    contains
    ! convert to t_RGB    
      pure elemental type(t_rgb) function to_rgb(ir, ig, ib)
        integer, intent(in) :: ir, ig, ib
        to_rgb = t_rgb(achar(ib), achar(ig), achar(ir))
      end function to_rgb  
    end module m_bmp

    submodule (m_bmp) m_bmp_body
      use, intrinsic :: iso_fortran_env
      implicit none
      type :: t_bmp_file_header
        sequence  
        integer(int16) :: bfType = transfer('BM', 0_int16) ! BitMap
        integer(int32) :: bfSize          ! file size in bytes
        integer(int16) :: bfReserved1 = 0 ! always 0
        integer(int16) :: bfReserved2 = 0 ! always 0
        integer(int32) :: bfOffBits
      end type t_bmp_file_header
      !
      type :: t_bmp_info_header
        sequence
        integer(int32) :: biSize          = Z'28' ! size of bmp_info_header ; 40bytes 
        integer(int32) :: biWidth
        integer(int32) :: biHeight
        integer(int16) :: biPlanes        = 1 ! always 1
        integer(int16) :: biBitCount
        integer(int32) :: biCompression   = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
        integer(int32) :: biSizeImage
        integer(int32) :: biXPelsPerMeter = 3780 ! 96dpi
        integer(int32) :: biYPelsPerMeter = 3780 ! 96dpi 
        integer(int32) :: biClrUsed       = 0
        integer(int32) :: biClrImportant  = 0 
      end type t_bmp_info_header
    contains   
      module procedure construct_bmp
        allocate(t_bmp(nx, ny)::res)
      end procedure construct_bmp    

      module procedure assign_bmp
        allocate(dest, source = src)
      end procedure assign_bmp    

      module procedure wr_header
        type (t_bmp_file_header) :: bmp_file_header
        type (t_bmp_info_header) :: bmp_info_header
        integer :: i, j, iw
        associate(nx => bmp%nx, ny => bmp%ny * num_images())                          ! coarray
          bmp_file_header%bfSize      = 14 + 40 + 0 + (3 * nx + mod(nx, 4)) * ny 
          bmp_file_header%bfOffBits   = 14 + 40
          bmp_info_header%biWidth     = nx
          bmp_info_header%biHeight    = ny 
          bmp_info_header%biBitCount  = 24 
          bmp_info_header%biSizeImage = (3 * nx + mod(nx, 4)) * ny 
          open(newunit = iw, file = fn//'.bmp', access = 'stream', status = 'replace')
          write(iw) bmp_file_header
          write(iw) bmp_info_header
          close(iw)
        end associate
      end procedure wr_header

      module procedure wr_bmp
        integer :: i, j, iw
        associate(nx => bmp%nx, ny => bmp%ny)
          open(newunit = iw, file = fn//'.bmp', access = 'stream', status = 'unknown', position = 'append')
          write(iw) (bmp%rgb(:, i), (achar(0), j = 1, mod(nx, 4)), i = 1, ny)
          close(iw)
        end associate  
      end procedure wr_bmp

      module procedure rd_bmp
        type (t_bmp_file_header) :: bmp_file_header
        type (t_bmp_info_header) :: bmp_info_header
        integer :: i, j, ir
        character :: dummy
        associate(nx => bmp_info_header%biWidth, ny => bmp_info_header%biHeight)
          open(newunit = ir, file = fn//'.bmp', access = 'stream', status = 'old')
          read(ir) bmp_file_header
          read(ir) bmp_info_header
          bmp = t_bmp(nx, ny)
          read(ir) (bmp%rgb(:, i), (dummy, j = 1, mod(nx, 4)), i = 1, ny)
          close(ir)
        end associate
      end procedure rd_bmp
    end submodule m_bmp_body

    program mandel
      implicit none
      integer, parameter :: nx = 1279, ny = 1280, maxiter = 255
      real   , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5
      complex, allocatable :: c(:, :)[:], z(:, :)[:]                ! coarray
      integer :: i, j, me, ni, my, iter
      integer, allocatable :: niter(:, :)[:]                        ! coarray
      ni = num_images()
      me = this_image()
      if (mod(ny, ni) /= 0) stop 'ny must be a multiple of num_image()'
      my = ny / ni
      allocate(c(nx, my)[*], z(nx, my)[*]) ! implicitly sync all    ! coarray
!
! make 2D-mesh :  c(:, :)
!      
      block
        real  :: x(nx), y(ny)
        x = [( (x1 - x0) / (nx - 1) * (i - 1) + x0, i = 1, nx )]
        y = [( (y1 - y0) / (ny - 1) * (i - 1) + y0, i = 1, ny )]
        forall (i = 1:nx, j = 1:my) c(i, j) = cmplx(x(i), y((me - 1) * my + j)) 
      end block
 !   
 ! main iteration : niter(:, :)        
 !           
      allocate(niter(nx, my)[*]) ! implicitly sync all              ! coarray    
      niter = 0
      do iter = 0, maxiter
        where (abs(z) <= 2.0) 
          z = z * z + c
          niter = niter + 1
        end where  
      end do
 !
 ! make bmp file: Mandel.bmp   
 !    
      block 
        use m_bmp
        type (t_bmp(:, :)), allocatable :: bmp[:]                   ! coarray
        allocate(bmp[*], source = t_bmp(nx, my))
        bmp%rgb = to_rgb(256 - niter, 256  - niter, 256 - niter)
        if (me == 1) call bmp%wr_header('mandel') ! write header
        if (me > 1 ) sync images(me - 1)                            ! coarray
          print '(3g0)', 'image ', me, ' writing..' 
          call bmp%wr('mandel')                   ! write body
        if (me < ni) sync images(me + 1)                            ! coarray  
      end block
    end program mandel