Google Colaboratory で CAF
Google Colaboratory で Fortran が使え、あまつさえ root 権限があるので、OpenCoarray も インストールできるというので、早速使ってみました。
【乞食速報】本が安い！Google Colaboratory で CAF、MS Azure で IJulia
でも、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 の入力
セル入力
%%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
セル入力
%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>
（ 図 ）
乞食は三日やったらやめられない！