【Fortran】総称名を使わずに手続きをまとめる

# 問題意識

複数の次元や変数の型を統一的に扱う1つの方法は，それぞれの場合にサブルーチンや関数（手続き）を定義した後に`interface`ブロックで総称名を付けるというものです．
例：`read_file_int_1d`, `read_file_int_2d`, `read_file_real_2d`をそれぞれ定義し，`read_file`という名前を付ける

僕自身，1つの場合に定義した手続きをコピペ→適宜編集（`(:)`→`(:,:)`，`real(4)`→`integer`等）でしのいできましたが，更新作業が煩雑になるため，総称名を使わずに済む方法をまとめました．
主に使う機能はassumed-rank引数，class(*)になります．

# Assumed-rank引数

引数となる配列の大きさを指定する方法は，これまで`a(5)`（形状明示配列explicit-shape array）や`b(:)`（形状引継ぎ配列assumed-shape array）等がありましたが（例えば[参考](https://www.softek.co.jp/SPG/Pgi/OpenACC/fortran_array.html)），
`c(..)`と表記することで配列の次元を明示せずに良くなります（Assumed-rank引数）．
細かい使い勝手は良くないですが，`size`や`shape`といった配列関数は使うことができます．

### 1次元，2次元，…の配列をまとめる

```Fortran:dim.f90
module mymod
    implicit none
contains

subroutine write_rank_shape( x )
    integer, intent(in) :: x(..)
    write(*, *) rank( x ), ':', shape( x )
end subroutine write_rank_shape

end module mymod

program main
    use mymod, only: write_rank_shape
    implicit none
    integer, allocatable :: a(:), b(:,:)

    allocate( a(5) ); a(:) = 0
    allocate( b(2,3) ); b(:,:) = 0
    call write_rank_shape( 1 )
    call write_rank_shape( a(:) )
    call write_rank_shape( b(:,:) )
end program main
```

備考：`dimension(:)`は1次元のみ，`dimension(:,:)`は2次元のみしか扱えません．
　　　この書法であれば`integer, dimension(..)`になります．

参考：
[Fortran 2018 : スカラーと配列両受け](http://fortran66.hatenablog.com/entry/2019/05/30/022957)
[Intel Fortran Compiler 16.0　での新機能](http://fortran66.hatenablog.com/entry/2015/04/09/212902)

### 1度assumed-rankで受けると，それを次元を明示した手続きに与えることはNG

一旦`arr(..)`として引数を受け取ってしまうと，それを`call mysub( arr(..) )`や`call mysub( arr )`として使用するには，呼び出す手続きでもassumed-rankとして定義している必要があります．
従って以下のコード（2次元配列`a`をassumed-rankを用いた手続き`wrapper`に与え，その手続きが2次元用の手続き`sum_dim1`を呼び出す）はコンパイルエラーになります．

```Fortran
subroutine wrapper( arr )
    integer, intent(in) :: arr(..)
    call sum_dim1( arr )
end subroutine wrapper

subroutine sum_dim1( arr )
    integer, intent(in) :: arr(:,:)
    write(*, *) sum( arr(:,1) )
end subroutine sum_dim1

program main
    use mymod
    integer :: a(3,2)
    a(:,:) = 1
    call wrapper( a(:,:) )
end program main
```

# 万能クラスunlimited polymorphic class

Fortran2003から変数の型を明示しなくても良い万能クラス`class(*)`が導入されました．
関数やサブルーチンの引数に何でも取れるようになり，加えて主プログラムでも使うことが出来ます（例えば[Qiita: 万能クラス Class(*) にいろんな型を突っ込む](https://qiita.com/cure_honey/items/e6ca856788ee6c3e42bb)）．
しかし，主プログラム内で使うには`allocatable`属性を付与する必要がある等文法が特殊である他，単純に`class(*)`に書き換えても動作しない場合が多々あります．


### 複数の変数の型をまとめる

`class(*)`で引数を受けた後に，`select type`文で変数の型について場合分けします．
通常のスイッチ構文`select case`とは異なり，

- `is`が必要である
- デフォルトは`class default`になる
- `type is (a, b)`と複数を併記することは出来ない

ことに注意が必要です．

ここでは変数のバイト数を返す関数を紹介します．

```Fortran:class.f90
module mymod
    implicit none

contains

function get_byte_size( x ) result( byte )
    class(*), intent(in) :: x
    integer :: byte
    select type ( x )
        type is ( real(8) )
            byte = 8
        type is ( integer(1) )
            byte = 1
        class default
            byte = 4
    end select
end function get_byte_size

end module mymod

program main
    use mymod, only: get_byte_size
    implicit none
    write(*, *) get_byte_size( 1.0  ) ! 4
    write(*, *) get_byte_size( 1.d0 ) ! 8
    write(*, *) get_byte_size( 1_1  ) ! 1
    write(*, *) get_byte_size( 1    ) ! 4
end program main
```

備考：`1_1`で1バイト整数の1という意味です．

参考：[Fortranにおける整数型・実数型・複素数型変数の宣言方法](https://qiita.com/implicit_none/items/ab0632357bab0f549cf2)

#### 未解決注意：上2つの合わせ技（複数次元×複数の型）

- `class(*) :: x(..)`として定義すると，`error #8793: This assumed-rank variable must not appear in a designator or expression in this context.`とコンパイルエラー．
- 別の関数`get_byte_size_arr`で`x(..)`のランクで場合分けしても，`get_byte_size( x )`で類似のエラー．
- 引数がスカラーの場合を諦めて大きさ引継ぎ行列（assumed-size array）`x(*)`を使ってもエラー．

そもそもassumed-rankで受けた引数に対して，（shapeやsumといった配列関数ではなく）各要素へのアクセスの仕方が分からず，頓挫しました．

#### 未解決注意：入力された引数の変数の型に応じた型を返す

例えば「変数の型にかかわらず，与えられた1次元配列の先頭の要素を返す」という関数を単純に書くと，

```Fortran:class_out.f90
module mymod
    implicit none

contains

function get_first_value( x ) result( x1 )
    class(*), intent(in)  :: x(:)
    class(*), allocatable :: x1
    x1 = x(1)
end function get_first_value

end module mymod

program main
    use mymod, only: get_first_value
    implicit none
    class(*), allocatable :: r ! どちらでも機能しない
!    real(4) :: r              ! どちらでも機能しない
    r = get_first_value( [1.0, 2.0, 3.0] )
    write(*, *) r
end program main
```

のようになると思います．ここでは大きさ3の単精度実数の配列から，先頭である1.0を取り出そうとしています．
しかし上記コードは機能しません．なぜなら，

- 返り値を`class(*)`として定義しているとき，それを呼び出す側の変数でも`class(*)`でなくてはならない
    - このとき`write`のところで型情報がないためエラー
- 単精度実数を入力したのだから単精度実数が返ってくるはず（上でコメントアウトしてある方）は型の不一致でコンパイルエラー
- 関数の返り値を`real(4)`としてもコンパイルエラー

つまり，プログラムのどこかで`select type`によって型の明示を行う必要があるのです．
これを関数の中に含めると以下の通りになり，動きます．

```Fortran:class_out2.f90
module mymod
    implicit none

contains

function get_first_value( x ) result( x1 )
    class(*), intent(in)  :: x(:)
    real(4) :: x1
    select type ( x )
        type is ( real(4) )
            x1 = x(1)
    end select
end function get_first_value

end module mymod

program main
    use mymod, only: get_first_value
    implicit none
    real(4) :: r
    r = get_first_value( [1.0, 2.0, 3.0] )
    write(*, *) r
end program main
```

しかし[こちらのページ](https://qiita.com/cure_honey/items/e6ca856788ee6c3e42bb)でも指摘されている通り，上記の関数は返り値が`real(4)`に限定されているため，結局はそれ以外の型を返す他の関数を定義する（＋総称名でまとめる）必要があります．これだと，引数を`class(*)`で受ける必要性を感じない…，というのが現在の印象です．

またサブルーチンにしても（例`subroutine get_first_value( x, x1 )`）同様の問題は起こります．期待通りの動作をしてくれる数少ない例は，バイナリファイルの読み書きです．


### `class(*)`を使ってもバイナリファイルの読み書きには問題無し

```Fortran:io.f90
module mymod
    implicit none
    integer, parameter :: funit = 11

contains

function get_byte_size( x ) result( byte )
    ! 上と同じなので省略
end function get_byte_size

subroutine open_bin( unit, path, arr, status )
    integer,          intent(in) :: unit
    character(len=*), intent(in) :: path, status
    class(*),         intent(in) :: arr(:)
    integer :: recl
    recl = get_byte_size( arr(1) ) * size( arr )
    open( unit, file=trim(path), form='unformatted', access='direct', recl=recl, status=status )
    write(*, '(2a,i0,a,i0)') trim( path ), ' ', get_byte_size( arr(1) ), ' ', size( arr )
end subroutine open_bin

subroutine write_file( path, arr )
    character(len=*), intent(in) :: path
    class(*),         intent(in) :: arr(:)
    call open_bin( funit, path, arr(:), 'replace' )
    write( funit, rec=1 ) arr(:)
    close( funit )
end subroutine write_file

subroutine read_file( path, arr )
    character(len=*), intent(in)  :: path
    class(*),         intent(out) :: arr(:)
    call open_bin( funit, path, arr(:), 'old' )
    read( funit, rec=1 ) arr(:)
    close( funit )
end subroutine read_file

end module mymod

program main
    use mymod, only: write_file, read_file
    implicit none
    real(4) :: a(3) = [ 1.0, 2.0, 3.0 ]
    real(8) :: b(2) = [ 5.0, 4.0 ]
    integer :: c(4) = [ 6, 7, 8, 9 ]
    integer(1) :: d(5) = [ 0, 1, 2, 3, 4 ]

    call write_file( './test1.bin', a(:) ) ! ./test1.bin 4 3
    call write_file( './test2.bin', b(:) ) ! ./test2.bin 8 2
    call write_file( './test3.bin', c(:) ) ! ./test3.bin 4 4
    call write_file( './test4.bin', d(:) ) ! ./test4.bin 1 5
    a(:) = 0.0
    b(:) = 0.0
    c(:) = 0
    d(:) = 0
    call read_file( './test1.bin', a(:) )
    call read_file( './test2.bin', b(:) )
    call read_file( './test3.bin', c(:) )
    call read_file( './test4.bin', d(:) )
    write(*, *) a(:) ! 1.000000       2.000000       3.000000
    write(*, *) b(:) ! 5.00000000000000        4.00000000000000
    write(*, *) c(:) ! 6           7           8           9
    write(*, *) d(:) ! 0    1    2    3    4
end program main
```

### `class(*)`変数のコピー

`class(*)`で受けた引数をもう1つの`class(*)`変数に代入した場合も，やはり`select type`が必要になります．

```Fortran
subroutine copy_class( x )
    class(*), intent(in)  :: x
    class(*), allocatable :: y
    y = x ! allocate( y, source=x )でも動作しました
!    write(*, *) y ! コンパイルエラー
    select type ( y )
        type is ( real(4) )
            write(*, *) y
        type is ( real(8) )
            write(*, *) y
    end select
end subroutine copy_class

program main
    use mymod, only: copy_class
    implicit none
    call copy_class( 1.0  ) ! 1.000000
    call copy_class( 2.d0 ) ! 2.00000000000000
end program main
```

しかし，`class(*), intent(out)`である引数に代入しようとすると，面倒なことになります．例えば

```Fortran
subroutine copy_class2( x )
    class(*), intent(inout) :: x
    class(*), allocatable   :: y
    y = x
    x = y ! これがコンパイルエラー
end subroutine copy_class2
```

では"error #8304: In an intrinsic assignment statement, variable shall not be a non-allocatable polymorphic."としてコンパイルエラーになります．
これを解決するためには，代入元（`y`）と代入先（`x`）について，それぞれ`select type`をする必要があります．

```Fortran
subroutine copy_class( x )
    class(*), intent(inout) :: x
    class(*), allocatable   :: y
    y = x
    select type ( x )
        type is ( real(4) ) ! xがreal(4)であると明示
            select type ( y )
                type is ( real(4) ) ! yがreal(4)であると明示
                    x = y
            end select
    end select
end subroutine copy_class
```

### `class(*)`配列の割り付け

`class(*)`型の配列を割り付ける際は，具体的な型を指定する必要があります．そのためには，`allocate`の引数である`source=`か`mold=`を使います．前者では与えた配列の値が代入されますが，後者では代入されません．

```Fortran
program main
    implicit none
    class(*), allocatable :: a(:)
    allocate( a(3), mold=0.0 ) ! mold=の引数は変数でもOK
    a(:) = 1.0
    select type ( a )
        type is ( real(4) )
            write(*, *) sum( a ) ! 3.000000
    end select
end program main
```

参考：

[stackoverflow: Use of unlimited polymorphic type for array operation in Fortran 03/08 (gfortran compiler)](https://stackoverflow.com/questions/27165854/use-of-unlimited-polymorphic-type-for-array-operation-in-fortran-03-08-gfortran)

#### `source=`と`mold=`の違い

`source=`は与えた配列の値が代入され，`mold=`は代入されないという違いの他に，`mold=`では与えた配列と異なる大きさにも割り付けることが可能だという違いがあります．

```Fortran
    class(*), allocatable :: a(:)
    real(4) :: b(5) = [1.0, 2.0, 3.0, 4.0, 5.0]
    allocate( a,    mold=b       ) ! aは大きさ5に割り付けられる
    allocate( a(3), mold=b       ) ! 3，つまりmold=よりa(3)が優先される
    allocate( a,    source=b     ) ! 5
    allocate( a,    source=b(:3) ) ! 3
    allocate( a(3), source=b     ) ! コンパイルエラー
```

従って，「与えられた配列と同じ型で，それよりも大きい配列が欲しい」という場合には，`mold=`は型を指定するためだけに使い，配列の大きさは明示的に指定するのが便利です．
※但し，`allocate`の元になる配列（ここでは`b`）を`class(*)`で受けていた場合，コンパイルは出来るもののエラーになります．このときは`select type (b)`→`allocate( a(3), mold=1.0 )`等で解決しました．

# まとめ

この記事の到達度ですと，痒い所に手が届いていない現状です．
ただ個人的には2次元or3次元のファイルの読み書きが喫緊の課題なので，

- 変数の型は`class(*)`でまとめる
- 次元は2次元と3次元で別々に定義する
- get_byte_sizeは_2d, _3dと別々に定義して総称名にしても良い？
- 領域毎に分かれたファイルを読み込んで1つの配列に統合するために，`copy_area`というサブルーチンを作っておき，2重の`select type`と代入を任せる

という方法で当分はしのごうと思います．