CoArray Fortran で Mandelbrot (Google Colaboratory 編)　

#Google Colaboratory で CAF 

Google Colaboratory で Fortran が使え、あまつさえ root 権限があるので、OpenCoarray も インストールできるというので、早速使ってみました。

[【乞食速報】本が安い！Google Colaboratory で CAF、MS Azure で IJulia](http://fortran66.hatenablog.com/entry/2018/09/28/220201)

でも、google collaboratory で使える cpu 数は２個っぽいので、caf で並列化してもあまり恩恵はなく、あくまで文法学習用としての利用となると思います。以下の例では 8 並列で実行していますが、かえって遅くなっていると思われます。

##実行結果
いつものマンデルブロ集合の図を描きます。

gfortran は parameterized derived type が ver.8 以降でしか使えず、apt で入れられる caf は、ver.8 対応ではないので、これを使わないプログラムにする必要がありました。

少し古い版のプログラムを元に修正したので、出力 bmp ファイルの横幅が 4 の倍数でないといけません。（2,3 行だけ直せばいいのですが、台風のせいでだるいので。）

実行は apt install も含めて、結構待たずに進みます。

##Colaboratory への入力
- open coarray fortran の導入

セル入力

```
!apt install mpich
!apt install open-coarrays-bin
!apt install less
!caf --version
```
応答（確認用に再実行したのでインストール済みと表示されていますが、実際はダラダラメッセージが出ます）

```
Reading package lists... Done
Building dependency tree       
Reading state information... Done
mpich is already the newest version (3.3~a2-2).
0 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.
Reading package lists... Done
Building dependency tree       
Reading state information... Done
open-coarrays-bin is already the newest version (1.9.1-1).
0 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.
Reading package lists... Done
Building dependency tree       
Reading state information... Done
less is already the newest version (481-2.1ubuntu2).
0 upgraded, 0 newly installed, 0 to remove and 0 not upgraded.

OpenCoarrays Coarray Fortran Compiler Wrapper (caf version 1.9.1)
Copyright (C) 2015-2016 Sourcery, Inc.

OpenCoarrays comes with NO WARRANTY, to the extent permitted by law.
You may redistribute copies of OpenCoarrays under the terms of the
BSD 3-Clause License.  For more information about these matters, see
the file named LICENSE.

```

- fortran program の入力

セル入力

```fortran
%%writefile caf.f90
    module m_bmp
        implicit none
        type :: t_bmp_file_header
            sequence  
            integer(2) :: bfType = Z'4D42' !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_header => wr_bmp_header
            procedure :: wr => wr_bmp
        end type
    contains
! write out BMP file    
        subroutine wr_bmp_header(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
            integer :: me, ni, ir, ic
            ni = num_images()
            me = this_image()
            associate(nx => size(bmp%rgb, 1), ny => size(bmp%rgb, 2) * ni)
                bmp_file_header%bfSize      = 14 + 40 + 0 + nx * ny * 3
                bmp_file_header%bfOffBits   = 14 + 40
                bmp_info_header%biWidth     = nx    ! nx shouold be a multiple of 4
                bmp_info_header%biHeight    = ny
                bmp_info_header%biBitCount  = 24 
                bmp_info_header%biSizeImage = nx * ny * 3
            end associate
            open(9, file = fn//'.bmp', access = 'stream', status = 'replace')
            write(9) bmp_file_header
            write(9) bmp_info_header
            endfile(9)
            close(9)
        end subroutine wr_bmp_header

        
        subroutine wr_bmp(bmp, fn)
            class(t_bmp), intent(in) :: bmp 
            character(len = *), intent(in) :: fn
            open(9, file = fn//'.bmp', access = 'stream', status = 'old', position = 'append')
            write(9) bmp%rgb
            close(9)
        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

    
    program Mandel
        use m_bmp
        implicit none       ! nx must be a multiple of 4
        integer, parameter :: nx = 1920, ny = 1920, maxiter = 255
        real   , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5
        complex, allocatable :: c(:, :)[:], z(:, :)[:]
        integer, allocatable :: niter(:, :)[:]
        integer :: me, ni, my
        type (t_bmp) :: bmp
        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)[*], bmp%rgb(nx, my)[*])     
!
! make 2D-mesh :  c(:, :)
!      
        block
            integer :: i, j
            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(:, :)        
!     
        call stamp('start Mandel')
        allocate(niter(nx, my)[*]) ! implicitly sync all
        niter = imandel(c)         
        call stamp('end   Mandel')
!
! make bmp file: Mandel.bmp   
!
        bmp%rgb = to_rgb(256 - niter, 256 - niter, 256 - niter)
        call stamp('start bmp')
        if (me == 1) call bmp%wr_header('Mandel')  
        if (me > 1) sync images(me - 1) ! order 1, 2, 3...., n 
        call bmp%wr('Mandel')  
        if (me < ni) sync images(me + 1) 
        call stamp('end   bmp')
        call stamp('normal end')
        deallocate(bmp%rgb) ! implicitly sync all    
!        stop
    contains

  
        pure elemental integer function imandel(c)
            complex, intent(in) :: c
            complex :: z
            z = c
            do imandel = 1, maxiter
                if (abs(z) > 2.0) exit
                z = z * z + c
            end do    
        end function imandel
        

        subroutine stamp(text)
            character(*), intent(in) :: text
            real :: t0 = 0.0, t1      ! t0 save
            integer :: c0 = 0, c1, cr ! c0 save
            logical :: first = .true.
            if (first) then 
                first = .false.
                call cpu_time(t1)
                call system_clock(c1, cr)
                t0 = t1
                c0 = c1
            end if  
            call cpu_time(t1)
            call system_clock(c1, cr)
            print '(a, i3, a, f12.4, a, f12.4, 2a)', 'image= ', this_image(), ': time= ',&
                                     real(c1 - c0) / cr,' cpu= ', t1 - t0, ': ', text
            t0 = t1
            c0 = c1
        end subroutine stamp
    end program Mandel
```
応答

```
Overwriting caf.f90
```
- コンパイル


セル入力

```
!caf caf.f90
```
応答

```
(無応答）
```
- 実行

セル入力

```
!cafrun -np 8 --allow-run-as-root a.out
```
応答

```
image=   5: time=       0.0000 cpu=       0.0000: start Mandel
image=   1: time=       0.0000 cpu=       0.0000: start Mandel
image=   6: time=       0.0000 cpu=       0.0000: start Mandel
image=   7: time=       0.0000 cpu=       0.0000: start Mandel
image=   2: time=       0.0000 cpu=       0.0000: start Mandel
image=   8: time=       0.0000 cpu=       0.0000: start Mandel
image=   3: time=       0.0000 cpu=       0.0000: start Mandel
image=   4: time=       0.0000 cpu=       0.0000: start Mandel
image=   8: time=       0.3760 cpu=       0.0915: end   Mandel
image=   1: time=       0.3940 cpu=       0.0956: end   Mandel
image=   1: time=       0.0640 cpu=       0.0158: start bmp
image=   8: time=       0.0960 cpu=       0.0150: start bmp
image=   7: time=       0.7100 cpu=       0.1737: end   Mandel
image=   2: time=       0.7390 cpu=       0.1706: end   Mandel
image=   7: time=       0.0670 cpu=       0.0157: start bmp
image=   2: time=       0.0620 cpu=       0.0139: start bmp
image=   1: time=       0.3600 cpu=       0.0876: end   bmp
image=   1: time=       0.0000 cpu=       0.0000: normal end
image=   3: time=       3.7550 cpu=       0.8708: end   Mandel
image=   6: time=       3.7800 cpu=       0.8757: end   Mandel
image=   3: time=       0.0680 cpu=       0.0155: start bmp
image=   6: time=       0.0570 cpu=       0.0148: start bmp
image=   2: time=       3.0450 cpu=       0.7064: end   bmp
image=   2: time=       0.0000 cpu=       0.0001: normal end
image=   4: time=       7.4240 cpu=       1.8315: end   Mandel
image=   5: time=       7.5450 cpu=       1.8372: end   Mandel
image=   4: time=       0.0860 cpu=       0.0191: start bmp
image=   3: time=       3.7120 cpu=       0.9114: end   bmp
image=   3: time=       0.0010 cpu=       0.0001: normal end
image=   5: time=       0.0620 cpu=       0.0135: start bmp
image=   4: time=       0.1840 cpu=       0.0455: end   bmp
image=   4: time=       0.0040 cpu=       0.0000: normal end
image=   5: time=       0.3410 cpu=       0.0741: end   bmp
image=   5: time=       0.0080 cpu=       0.0001: normal end
image=   6: time=       4.3150 cpu=       1.0487: end   bmp
image=   6: time=       0.0070 cpu=       0.0001: normal end
image=   7: time=       7.5720 cpu=       1.8442: end   bmp
image=   7: time=       0.0050 cpu=       0.0001: normal end
image=   8: time=       8.0740 cpu=       1.9225: end   bmp
image=   8: time=       0.0000 cpu=       0.0001: normal end
```
- 作図 matplotlib

セル入力

```python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
plt.rcParams['figure.figsize'] = (35.0, 35.0)
img = Image.open("Mandel.bmp")
plt.imshow(img)
```
応答

```
<matplotlib.image.AxesImage at 0x7f62711dffd0>
 （ 図 ）
```

![mandelCAF.png](https://qiita-image-store.s3.amazonaws.com/0/44522/a38d8c25-0c61-b0eb-11ab-61903d369d90.png)

乞食は三日やったらやめられない！