Fortran で zip 的連結 その2の巻

Fortran で zip 的連結

二つの等長な配列 x と y があるとき、これを x, y 成分の対とする操作を zip と呼びます。リスト処理系の言語にはよくあります。

zip を行う関数があると、たとえば入出力などで、

do i = 1, n
 print *, x(i), y(i)
end do

が、

print *, zip(x, y)

のように１行で書けるようになり便利です。

以前、同じ型の配列間での zip 的作用が transpose で出来ることを書きました。
https://qiita.com/cure_honey/items/6d17a1a6696435abd241
今回は、配列の型が異なる場合を考えます。

考え方

今、zip 関数に対して任意の型を取り得るようにするには、入力型として無制限クラス class(*) を用いるしかありません。また zip 関数での出力型は、この入力を対にした class(*) の派生型の配列を出力します。

さらに class を使うと入出力用にユーザー定義派生型 I/O ルーチンを用意してやる必要が出てきます。書式 (format) の扱いなどが煩雑で面倒です。

Fortran は実行時のコストが高い命令は、わざと使いにくくしてあるようなので、class(*) のように、実行時にならないと型が分からず、型と値の両方の情報を運ばねばならないものは扱いがペナルティー的に面倒になります。つまり実行時性能を優先した発想になっています。バランスのとり方として、ソースプログラムの書きやすさを優先する近年の流れとは異なった価値観になっていると思われます。

ソースプログラム　intel fortran v.18

はじめに、zip が吐き出す派生型配列の元の型を定義します。配列としての処理は elemental な要素毎処理に任せることにして、基本はスカラー要素で定義します。

スカラー要素出力のところは、select type によって、各型毎に場合分けして処理する必要があり煩雑です。もう少し簡潔に出来るかもしれませんが、詳しく考えていません。

ユーザー定義派生型 I/O は、書式なしの * 出力 (List-directed I/O) の場合と、書式付きの場合に分けて処理します。書式付きのばあい DT(i1,i2,i3,i4) の４つの整数パラメータを取ることにして、それが zip する x, y 配列それぞれに対して、g フォーマットの gi1.i2、gi3.i4 と展開されることとします。

ここで面倒になるのは、複素数の取り扱いで、複素数の場合 fortran では、実変数二つが並んでいるとみなした書式を与えるので困ってしまいます。ここでは解決策として、大は小を兼ねる式に、あらゆる場合に複素数用の変数２個並びの書式を構成して、複素数以外の場合は : により、書式を中断することで対応することとします。

    module m_mod
      implicit none
      type :: t_zip
        class(*), allocatable :: a, b  
      contains
        procedure :: wr
        generic   :: write(formatted) => wr
      end type t_zip    
    contains
      function zip(a, b) 
        class(*), intent(in) :: a(:), b(:)
        type(t_zip), allocatable :: zip(:)
        integer :: i
        allocate(zip(size(a)))
        do i = 1, size(zip)
          allocate(zip(i)%a, source = a(i))
          allocate(zip(i)%b, source = b(i))
        end do
      end function zip

      subroutine wr(dtv, unit, iotype, vlist, io, iomsg)
        class(t_zip), intent(in) :: dtv
        integer, intent(in) :: unit
        character(*), intent(in) :: iotype
        integer, intent(in ) :: vlist(:)
        integer, intent(out) :: io
        character(*), intent(in out) :: iomsg
        character(20) :: fmt
        if (iotype == 'LISTDIRECTED') then
          io = 0
          fmt = '(g0:(1x, g0))'
          call sub_wr(unit, fmt, dtv%a)
          write(unit, '(1x)', advance = 'no')  
          call sub_wr(unit, fmt, dtv%b)
          write(unit, '(1x)', advance = 'yes')  
        else if (iotype == 'DT') then  
          write(fmt, '(a, g0, a, g0, a, g0, a, g0, a)') '(g', vlist(1), '.', vlist(2), ',:1x, g', vlist(1), '.', vlist(2), ')'
          call sub_wr(unit, fmt, dtv%a)
          write(unit, '(1x)', advance = 'no') 
          write(fmt, '(a, g0, a, g0, a, g0, a, g0, a)') '(g', vlist(3), '.', vlist(4), ',:1x, g', vlist(3), '.', vlist(4), ')'
          call sub_wr(unit, fmt, dtv%b)
          write(unit, '(1x)', advance = 'yes')  
        else   
          io = 999
          write(unit, *) '********************'
        end if    
      end subroutine wr

      subroutine sub_wr(unit, fmt, a)
        use, intrinsic :: iso_fortran_env
        class(*), intent(in) :: a
        integer , intent(in) :: unit
        character(*), intent(in) :: fmt
        select type (a)
        type is (integer(int32))
          write(unit, fmt, advance = 'no') a  
        type is (integer(int64))
          write(unit, fmt, advance = 'no') a  
        type is (real)
          write(unit, fmt, advance = 'no') a  
        type is (real(real64))
          write(unit, fmt, advance = 'no') a  
        type is (real(real128))
          write(unit, fmt, advance = 'no') a
        type is (complex(real32))
          write(unit, fmt, advance = 'no') a
        type is (complex(real64))
          write(unit, fmt, advance = 'no') a
        type is (complex(real128))
          write(unit, fmt, advance = 'no') a
        type is (character(*))
          write(unit, fmt, advance = 'no') a
        type is (logical(int32))
          write(unit, fmt, advance = 'no') a
        type is (logical(int64))
          write(unit, fmt, advance = 'no') a
        class default  
          write(unit, *) 'type not defined: subrotuine sub_wr'      
        end select    
      end subroutine sub_wr
    end module m_mod

    program Zipper
      use m_mod
      implicit none
      type (t_zip) :: tmp
      type (t_zip), allocatable :: t(:)     
      integer :: i
      t = zip([1, 2, 3, 4, 5, 6], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
      !
      write(6, '( *((DT(0, 0, 0, 0)):/) )') t 
      !
      print *, zip(['a', 'b', 'c'], sin([1.0d0, 2.0d0, 3.0d0]))
      !
      print *, zip([.true., .false.], [cmplx(0.0, 4 * atan(1.0)), (1.0, 3.0)])
      !
      t =  zip([.true., .false.], [cmplx(0.0, 4 * atan(1.0)), (1.0, 3.0)])
      print '( DT(0, 0, 20, 15) / DT(0, 0, 10, 5))',  t
    end program Zipper

註

      write(6, '( *((DT(0, 0, 0, 0)):/) )') t 

ここで書式の * は無限反復、: はデータ無いときここで中断、/ は改行、DT はユーザー定義派生型 I/O に渡す書式のパラメータのための命令記号です。(DT(0, 0, 0, 0)) と括弧でくくってあるのは、こうしないと余計な改行が最初に入ってしまうからで、その理由はよく分かりません。コンパイラの問題かもしれませんが、規格仕様が複雑なのでもう少し勉強したいと思います。

また、DT 書式を与えた場合、変数を経由せずに print 文に直接 zip 関数の出力を渡すと、文法エラーで叱られました。

実行結果

1 1.000000
2 2.000000
3 3.000000
4 4.000000
5 5.000000
6 6.000000
 a .8414709848078965 b .9092974268256817 c .1411200080598672
 T .000000 3.141593 F 1.000000 3.000000
T 0.00000000000000     3.14159274101257
F 1.0000     3.0000
続行するには何かキーを押してください . . .

様々な型の組み合わせに対して、zip 処理がなされています。また書式も指定通りに動作できています。

なお、おしまいの方の円周率は単精度で求めているので、7~8桁目以降は単に二進表現の展開による無意味な数字です。