Mandelbrot 図形

##Mandelbrot 図形
Mandelbrot 図形は、各座標点での計算が独立に実行できるので、並列化が効きやすくよく並列プログラミングの例題に用いられます。Fortran の場合大抵コンパイラのスイッチさえ入れれば自動並列化・ベクトル化が効いてくれます。

##Fortran
以下のプログラムの場合、Intel Fortran では、自動並列化・AVX利用のコンパイラ・オプションを指定することで、CPU使用時間は増えますが、実計算時間は約半分強に減ります。

宣言文やI/Oを除いた、実際の計算をする部分は以下の部分で、ほぼ数学的な定義式のまま記述できます。反復の部分はループに順序があるので DO LOOP を用います。一方、各座標点に関する計算には順序は関係ないので、並列性を意味する全配列演算を用いています。漸化式の部分 z = z**2 + c は、スカラー変数の計算のように書いていますが、実際は配列演算でピクセル上の全点に渡って一度に計算しています。

ここで、数学的には発散が決まった点 (|z| > 2.0) については、それ以上の計算は必要ないので、z=z*z+c は WHERE 構文の中に入れていいのですが、並列・ベクトル計算をする場合には、むしろメモリー上の連続番地で一気に計算して、あとから不要なデータを捨てた方が効率がいいので、外に出してあります。無駄な計算はしますが　ターンアラウンド的な意味で、早く終わります。

```fortran
 ! main iteration        
 !     
      z = (0.0, 0.0)
      do i = 1, maxiter
        z = z * z + c
        where (abs(z) < 2.0) niter = i
      end do
```

本来の where 構文は、６０、７０年代の並列計算機の実装に基づいて、そういう全計算＋マスクによる取捨選択を前提として導入されたと思うのですが、Intel Fortran の場合は、自動並列化をかけてみると WHERE 構文の mask
 処理はそのようにはなっていないようでした。


```fortran
 ! main iteration        
 !     
      z = (0.0, 0.0)
      do i = 1, maxiter
        where (abs(z) < 2.0) 
          z = z * z + c
          niter = i
        end where  
      end do
```
本来の形


###実行結果
![Mandel.png](https://qiita-image-store.s3.amazonaws.com/0/44522/742d0f67-8a2a-ee2f-0e11-e4b4447a60d1.png)
図は BMP を PNG に変換しています。

###プログラム
BMP 書き出しプログラム

```fortran
    module m_bmp
      implicit none
      type :: t_bmp_file_header
        sequence  
        integer(2) :: bfType = transfer('BM', 0_2, 1) ! BitMap
        integer(4) :: bfSize          ! file size in bytes
        integer(2) :: bfReserved1 = 0 ! always 0
        integer(2) :: bfReserved2 = 0 ! always 0
        integer(4) :: bfOffBits
      end type t_bmp_file_header
      ! 
      type :: t_bmp_info_header
        sequence
        integer(4) :: biSize     = Z'28' ! size of bmp_info_header ; 40bytes 
        integer(4) :: biWidth
        integer(4) :: biHeight
        integer(2) :: biPlanes   = 1 ! always 1
        integer(2) :: biBitCount
        integer(4) :: biCompression = 0 !0:nocompression,1:8bitRLE,2:4bitRLE,3:bitfield
        integer(4) :: biSizeImage
        integer(4) :: biXPelsPerMeter = 3780 ! 96dpi
        integer(4) :: biYPelsPerMeter = 3780 ! 96dpi 
        integer(4) :: biClrUsed      = 0
        integer(4) :: biClrImportant = 0 
      end type t_bmp_info_header
      !
      type :: t_rgb
        sequence
        character :: b, g, r  ! order is b g r 
      end type t_rgb 
      !
      type :: t_bmp
        type(t_rgb), allocatable :: rgb(:, :)
      contains 
        procedure :: wr => wr_bmp
      end type
    contains   
      subroutine wr_bmp(bmp, fn)
        class(t_bmp), intent(in) :: bmp 
        character(len = *), intent(in) :: fn
        type(t_bmp_file_header) :: bmp_file_header
        type(t_bmp_info_header) :: bmp_info_header
        associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2))
          bmp_file_header%bfSize      = 14 + 40 + 0 + nx * ny * 3
          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 = nx * ny * 3
        end associate
        open(9, file = fn//'.bmp', form = 'binary', status = 'unknown')
        write(9) bmp_file_header
        write(9) bmp_info_header
        write(9) bmp%rgb
        close(9)
        return
      end subroutine wr_bmp
 ! 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
```

メイン・プログラム

```fortran
    program Mandel
      use m_bmp
      implicit none
      integer, parameter :: nx = 1920, ny = 1920, maxiter = 255
      real   , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5
      complex :: c(nx, ny), z(nx, ny)
      integer :: i, niter(nx, ny)
!
! make 2D-mesh :  c(:, :)
!      
      block
        real  :: x(nx), y(nx)
        integer :: i, j  
        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:ny) c(i, j) = cmplx(x(i), y(j)) 
      end block
 !   
 ! main iteration : niter(:, :)        
 !     
      z = (0.0, 0.0)
      do i = 1, maxiter
        z = z * z + c
        where (abs(z) < 2.0) niter = i
      end do
 !
 ! make bmp file  : Mandel.bmp   
 !
      block 
        type (t_bmp) :: bmp
        allocate(bmp%rgb(nx, ny))
        bmp%rgb = to_rgb(255 - niter,255 - niter, 255 - niter)
        call bmp%wr('Mandel')
      end block  
      stop
    end program Mandel
```