type bound procedure で OO 風
少しバグがあって ビット列長が 255byte の倍数になっている時、余分な 00 が出力されていました。
それが丁度、intel compiler の associate の関わる bug と思われるものと重なってちょっと苦しみました。associate のポインタがらみのバグも中々取れてません。
図画用の変数はスカラーの allocatable として宣言して、コンストラクタで割り付けと初期化を行います。fig = t_gif('mandel.gif', nx, ny, n_global_color, global_color_table)
FINALIZE を使って、ファイルのクローズやメモリ解放処理を行おうとしたのですが、一方で関数でコンストラクタを作ったために、代入・自動割り付け後のコンストラクタの作った右辺のオブジェクトが解放される時にファイナライズ処理が一度なされてしまい、ファイルがクローズされるので、FINALIZE はやめて、明示的に終了処理をすることにしました。気長に解決法を考えたいと思います。
(H30.5.9)
二次元配列を一次元にするのに、Fortran2008 で拡張された pointer assignment を用いました。(MFE §15.6) これによってテンポラリ配列への代入の必要がなくなりました。
(H30.5.11)
2色の時の扱いが間違っていたので修正。
出力結果
アニメでサイケ調のマンデルブロ集合の収束を眺める
2色
Modern Fortran Explained (Numerical Mathematics and Scientific Computation)
- 作者: Michael Metcalf,John Reid,Malcolm Cohen
- 出版社/メーカー: Oxford University Press, U.S.A.
- 発売日: 2011/05/08
- メディア: ペーパーバック
- この商品を含むブログを見る
Introduction to Programming with Fortran: With Coverage of Fortran 90, 95, 2003, 2008 and 77
- 作者: Ian Chivers,Jane Sleightholme
- 出版社/メーカー: Springer
- 発売日: 2015/08/17
- メディア: ペーパーバック
- この商品を含むブログを見る
Scientific Software Design: The Object-Oriented Way
- 作者: Damian Rouson,Jim Xia,Xiaofeng Xu
- 出版社/メーカー: Cambridge University Press
- 発売日: 2011/04/29
- メディア: Kindle版
- この商品を含むブログを見る
ソース・プログラム
gfortran には findloc が無いようなので、代替品をば。あと forall(integer:: ...) にも対応していないようです。
pure integer function my_findloc(dict, k) result(ires) integer, intent(in) :: dict(:), k do ires = 1, size(dict) if (dict(ires) == k) exit end do if (ires > size(dict)) ires = 0 ! not found end function my_findloc
本体
module m_gif_types use, intrinsic :: iso_fortran_env implicit none private public :: t_gif_header, t_image_block, t_graphic_control_extension, t_application_extension ! constructors interface t_gif_header procedure :: init_gif_header end interface interface t_image_block procedure :: init_image_block end interface interface t_graphic_control_extension procedure :: init_graphic_control_extension end interface interface t_application_extension procedure :: init_application_extension end interface ! type :: t_gif_header sequence character(3) :: signature = 'GIF' character(3) :: version = '89a' ! '87a' integer(int16) :: width integer(int16) :: height integer(int8) :: pck = int(B'10001010', int8) !1:Global Col. 3:Color Res. 1:Sort Flg. 3:Size of Global Col. integer(int8) :: background_color_index = 0 integer(int8) :: pixel_aspect_ratio = 0 end type type :: t_image_block sequence integer(int8) :: image_separator = int(Z'2C', int8) integer(int16) :: image_left_position = 00 integer(int16) :: image_top_position = 00 integer(int16) :: image_width integer(int16) :: image_height integer(int8) :: pck = int(B'00000000', int8) ! local color table !integer(int8) :: LZW_minimum_code_size !integer(int8) :: block_size !integer(int8) :: image_data (:) !integer(int8) :: block_terminator = int(Z'00', int8) end type t_image_block type :: t_graphic_control_extension sequence integer(int8) :: extention_introducer = int(Z'21', int8) integer(int8) :: graphic_control_label = int(Z'F9', int8) integer(int8) :: lock_size = int(Z'04', int8) integer(int8) :: pck = int(Z'00', int8) integer(int16) :: delay_time = 50 ! msec integer(int8) :: transparent_color_index = 0 integer(int8) :: block_terminator = int(Z'00', int8) end type t_graphic_control_extension type :: t_application_extension ! for animation gif sequence integer(int8) :: extension_intrducer = int(Z'21', int8) integer(int8) :: extension_label = int(Z'FF', int8) integer(int8) :: block_size_01 = int(Z'0B', int8) character(len = 8) :: application_identifier = 'NETSCAPE' character(len = 3) :: application_authentication_code = '2.0' integer(int8) :: block_size_02 = int(Z'03', int8) ! 0:block terminator integer(int8) :: n = int(Z'01', int8) !application_data integer(int16) :: nloop = 0 ! 0:unlimited loop integer(int8) :: block_terminator = int(Z'00', int8) end type t_application_extension contains ! constructors pure type(t_gif_header) function init_gif_header(nx, ny, ngcol) result(res) integer, intent(in) :: nx, ny, ngcol res%width = nx res%height = ny res%pck = iand(res%pck, int(B'11111000', int8)) + int(mod(ngcol, 8), int8) !least 3bits:size of global color:2**(ngcol + 1) 2..256 end function init_gif_header pure type(t_image_block) function init_image_block(nx, ny) result(res) integer, intent(in) :: nx, ny res%image_width = nx res%image_height = ny end function init_image_block pure type(t_graphic_control_extension) function init_graphic_control_extension(it) result(res) integer, intent(in) :: it res%delay_time = int(it, int16) ! msec end function init_graphic_control_extension pure type(t_application_extension) function init_application_extension(nloop) result(res) integer, intent(in) :: nloop res%nloop = int(nloop, int16) end function init_application_extension end module m_gif_types module m_gif_lzw use, intrinsic :: iso_fortran_env use :: m_gif_types implicit none private public :: encoder type :: t_enc integer :: kbits = 0, id = 0 end type contains subroutine enc(irgb, nbits, code) ! gif LZW integer, intent(in) :: irgb(:), nbits type (t_enc), intent(out) :: code(:) integer :: dict(0:2**12 - 1) ! 0:4095 integer :: kbits, kd, m0 ! <== m0 state variable integer :: ip, ic, ienc call clear_dict() ! clear dictionary ic = 1 code(ic) = t_enc(kbits, 2**nbits) ! clear_code m0 = irgb(1) do ip = 2, size(irgb) call subenc(ienc, irgb(ip)) if (ienc /= -1) then ! new word ic = ic + 1 code(ic) = t_enc(kbits, ienc) end if if (kd == 2**12 - 1) then ! dictionary full ic = ic + 1 code(ic) = t_enc(kbits, 2**nbits) ! clear_code call clear_dict() ! clear dictionary end if end do ic = ic + 1 code(ic) = t_enc(kbits, m0) ic = ic + 1 code(ic) = t_enc(kbits, dict(2**nbits + 1)) ! end_code contains subroutine subenc(ienc, m1) integer, intent(out) :: ienc integer, intent(in ) :: m1 integer :: k, id k = ishft(m0 + 1, 16) + m1 ! m0 + 1 to avoid 00...0 degeneracy id = findloc(dict, k, dim = 1) ! dictionary is 0-based if (id == 0) then ! not found in the dictionary kd = kd + 1 dict(kd) = k if (kd == 2**kbits + 1) kbits = kbits + 1 ienc = m0 m0 = m1 else ! found in the dictionary ienc = -1 m0 = id - 1 ! because dictionary is 0-based, m0 must be shifted by 1 end if end subroutine subenc subroutine clear_dict() kbits = nbits + 1 kd = 2**nbits + 1 ! dictinary defined dict = 0 forall(integer :: i = 0:kd) dict(i) = i ! dict(nbit):clear code; dict(nbit+1):end code end subroutine clear_dict end subroutine enc subroutine encoder(irgb, nbits, icode) integer, intent(in ) :: irgb(:), nbits integer(int8), allocatable, intent(out) :: icode(:) type(t_enc), allocatable :: code(:) integer :: i, j, k, ib, ic, nb allocate(code(size(irgb) + 100)) ! ? large enough ? call enc(irgb, nbits, code) ! compress color code to LZW code nb = ceiling(sum(code(:)%kbits) / 8.0) ! required bytes allocate(icode(nb)) icode = 0 ! pack LZW code to bit string k = 0 do i = 1, size(code) do j = 1, code(i)%kbits if (btest(code(i)%id, j - 1)) then ic = k / 8 + 1 ib = mod(k, 8) icode(ic) = ibset(icode(ic), ib) end if k = k + 1 end do end do end subroutine encoder end module m_gif_lzw module m_anime_gif use, intrinsic :: iso_fortran_env use :: m_gif_types use :: m_gif_lzw implicit none private public :: t_gif interface t_gif procedure :: init_gif end interface type :: t_gif integer :: iw, nbits integer, allocatable :: irgb2(:, :) contains procedure :: wr_gif_img, gif_close end type t_gif contains impure type(t_gif) function init_gif(fn, nx, ny, n_global_color, global_color_table, loop) result(this) character(*), intent(in) :: fn integer, intent(in) :: nx, ny, n_global_color, global_color_table(:) integer, intent(in), optional :: loop type(t_gif_header) :: head type(t_application_extension) :: anime integer :: nloop = 0 if (present(loop)) nloop = loop associate(iw => this%iw, nbits => this%nbits) nbits = max(2, n_global_color + 1) ! 1..8 -> 2,2,3,4,5,6,7,8 allocate(this%irgb2(nx, ny)) head = t_gif_header(nx, ny, n_global_color) anime = t_application_extension(nloop) open(newunit = iw, file = fn, access = 'stream') write(iw) head write(iw) int(global_color_table, int8) write(iw) anime end associate end function init_gif subroutine wr_gif_img(this, itime) class(t_gif), target, intent(in) :: this integer, intent(in) :: itime type(t_graphic_control_extension), allocatable :: ext type(t_image_block), allocatable :: img integer(int8), allocatable :: icode(:) integer :: i, nbyte integer, pointer :: irgb(:) associate(iw => this%iw, nbits => this%nbits, irgb2 => this%irgb2) ext = t_graphic_control_extension(itime) img = t_image_block(size(irgb2, 1), size(irgb2, 2)) irgb(1:size(irgb2)) => irgb2 ! 1D => 2D call encoder(irgb, nbits, icode) ! color code to GIF LZW code write(iw) ext ! write bit strings to file write(iw) img write(iw) int(nbits, int8) ! LZW_minimum_code_size = 2,2,3,4,5,6,7,8bits do i = 1, size(icode), 255 nbyte = min(i + 255, size(icode)) - i write(iw) achar(nbyte) write(iw) icode(i:i+nbyte-1) end do write(iw) achar(00) ! block_terminator end associate end subroutine wr_gif_img subroutine gif_close(this) class(t_gif), intent(in out) :: this associate(iw => this%iw) write(iw) achar(int(Z'3B')) close(iw) end associate end subroutine gif_close end module m_anime_gif program mandel use :: m_anime_gif implicit none integer, parameter :: nx = 500, ny = 500, maxiter = 254 real , parameter :: x0 = -2.0, x1 = 1.0, y0 = -1.5, y1 = 1.5 complex, allocatable :: c(:, :), z(:, :) type(t_gif), allocatable :: fig allocate(c(nx, ny), z(nx, ny)) ! make 2D-mesh: c(:, :) forall(integer :: ix = 1:nx, iy = 1:ny) & c(ix, iy) = cmplx((x1 - x0) / (nx - 1) * (ix - 1), (y1 - y0) / (ny - 1) * (iy - 1)) + cmplx(x0, y0) block ! set global color table integer :: n_global_color ! 0..7 ! 2..256 colors (2^(n_global_color + 1)) integer, allocatable :: global_color_table(:) n_global_color = 7 allocate(global_color_table(3 * 2**(n_global_color + 1))) global_color_table = 0 forall(integer :: i = 0:254) global_color_table(3*i+mod(i, 3)+1) = 255 ! global_color_table(3*i+1:3*i+3) = i fig = t_gif('mandel.gif', nx, ny, n_global_color, global_color_table) end block block ! Mandelbrot Set integer :: iter fig%irgb2 = 0 do iter = 0, maxiter where (abs(z) <= 2.0) z = z**2 + c fig%irgb2 = fig%irgb2 + 1 end where if (iter < 10 .or. mod(iter, 10) == 1) call fig%wr_gif_img(max(80 - iter, 10)) end do call fig%wr_gif_img(500) end block call fig%gif_close() deallocate(fig) end program mandel
(H30.5.9/10 微改善, 11 コメント修正)