Fortranのanyとallで条件判定を簡略化する

概要

allとany関数を用いることで，複数の条件の真偽を簡潔に考慮できます．例えば，配列の形状や値の等価性を簡単に判別できるようになります．

環境

• Windows 10 64 bit
• gfortran 10.3
• Intel Fortran classic 2021.5.0
• NAG Fortran 7.1.0
• fypp 3.1

Fortranの関係演算子

all関数，any関数の使い方を説明する前に，Fortranにおける関係演算子（==, /=, <, <=, >, >=）について言及しておきます．

関係演算子は，演算子右と左の値を比較して，演算子が示す関係を満たしていれば.true.，そうでなければ.false.を返します．

print *, (1 == 2), (1 /= 2), (1 < 2), (1 > 2) ! F T T F

言語によっては，関係演算子の左右にはスカラしか置けない場合もありますが，Fortranは配列を置くことができます．その場合，関係演算子の戻値も配列になります．

integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))     ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b)) ! [[6, 5], [4, 3], [2, 1]]
print *, a < b ! [[T, T], [T, F], [F, F]]

配列要素の関係を，関係演算子を用いてまとめて評価できるのは非常に便利なのですが，その結果に基づいてif分岐をしようとすると，問題が生じます．if文の条件式は，スカラの論理値にしか対応していないからです．

program main
    use, intrinsic :: iso_fortran_env
    implicit none

    integer(int32) :: a(2, 3), b(2, 3), i
    a = reshape([(i, i=1, product(shape(a)))], shape(a))
    b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b))
    if (a < b) print *, "a is smaller than b"
end program main

このプログラムをコンパイルしようとすると，エラーが発生します．

    8 |     if (a < b) print *, "a is smaller than b"
      |        1
Error: IF clause at (1) requires a scalar LOGICAL expression

こういう場合に，any関数やall関数が役に立ちます．

all関数, any関数

Fortranにおけるall関数，any関数は，配列のreductionを行う関数の一つです．

関数	戻値	機能
all(mask [, dim])	logiacl	論理型配列maskが全要素が.true.の場合，真を返す．dimで検査する配列の次元を指定できる．
any(mask [, dim])	logiacl	論理型配列maskのいずれかの要素が.true.の場合，真を返す．dimで検査する配列の次元を指定できる．

インタフェースとしては非常に簡単で，論理型配列を渡せば，その内容に応じて.true.または.false.が返ってきます．

前節の例のように，関係演算子による比較結果に基づいてif文で判別を行う場合に，比較結果の論理型配列からさらに一つの結果（全要素が小さい，全要素が等しい等）を得る場合に利用します．

integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b))
if (all(a < b)) print *, "a is smaller than b" ! aはbより大きい要素も小さい要素も持っているのでこれは表示されない

配列の全要素の値を比較する場合は，下記のように書けるでしょう．

integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))     ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=product(shape(b)), 1, -1)], shape(b)) ! [[6, 5], [4, 3], [2, 1]]

if (all(a == b)) then
    print *, "a is equal to b"
else
    print *, "a is not equal to b."
    print *, "maximum absolute difference = ", maxval(abs(a - b)), "at ", maxloc(abs(a - b))
    print *, "minimum absolute difference = ", minval(abs(a - b)), "at ", minloc(abs(a - b))
    ! abs(a - b) == [[5, 3], [1, 1], [3, 5]]
end if

実行結果は下記のようになります．

 a is not equal to b.
 maximum absolute difference =            5 at            1           1
 minimum absolute difference =            1 at            1           2

配列要素のどれか一つでも一致していればよい場合については，下記のような例を挙げることができます．

integer(int32) :: a(2, 3), b(2, 3), i
a = reshape([(i, i=1, product(shape(a)))], shape(a))      ! [[1, 2], [3, 4], [5, 6]]
b = reshape([(i, i=1, product(shape(b))*2, 2)], shape(b)) ! [[1, 3], [5, 7], [9, 11]]

if (any(a == b)) then
    print *, "a has", count((a - b) == 0), " value(s) equal to b"
else
    print *, "a does not have the same value as b."
end if

実行結果は下記のようになります．

 a has           1  value(s) equal to b

その他のreductionを行う関数

既にいくつか現れていますが，all, any以外に，配列のreductionを行う関数を列挙しておきます．

関数	戻値	機能
count(mask [, dim])	integer	論理型配列maskにおける.true.の要素数を返す．dimで検査する配列の次元を指定できる．
sum(array, [, dim] [, mask])	arrayと同じ型	arrayの各要素の和を返す．dimで検査する配列の次元を指定できる．maskで計算対象となる要素を指定できる．
product(array, [, dim] [, mask])	arrayと同じ型	arrayの各要素の積を返す．dimで検査する配列の次元を指定できる．maskで計算対象となる要素を指定できる．
maxval(array [, dim] [, mask])	arrayと同じ型（整数型もしくは実数型）	配列array内の最大値を返す．dimで検査する配列の次元を指定できる．maskで検査する要素を指定できる．
minval(array [, dim] [, mask])	arrayと同じ型（整数型もしくは実数型）	配列array内の最小値を返す．dimで検査する配列の次元を指定できる．maskで検査する要素を指定できる．
maxloc(array [, dim] [, mask] [, kind])	整数型配列（次元はarrayと同じ．型種別はkindで指定）	配列array内で最初に最大値が現れる配列添字を返す．dimで検査する配列の次元を指定できる．maskで検査する要素を指定できる．
minloc(array [, dim] [, mask] [, kind])	整数型配列（次元はarrayと同じ．型種別はkindで指定）	配列array内で最初に最小値が現れる配列添字を返す．dimで検査する配列の次元を指定できる．maskで検査する要素を指定できる．
findloc(array, value [, dim] [, mask] [, kind] [, back])	整数型配列（次元はarrayと同じ．型種別はkindで指定）	配列array内で最初にvalueが現れる配列添字を返す．dimで検査する配列の次元を指定できる．maskで検査する要素を指定できる．back=.true.にすると，配列を後ろから検査する．

allとanyの使いどころ

配列の値の等価性比較

既に例示していますが，二つの配列の等価性を比較する際に，allやanyが活用できます．

二つの配列aとbの全要素が同じであることを判定したいのであれば，

all(a==b)

と書けばよく，どれか一つの要素でも値が等しくないことを判定したいのであれば，

any(a/=b)

と書けるでしょう．

実数の場合は，==で等価性を比較してはいけないと教えられたと思います．（私も教えられました）
実数型配列の等価性の比較は，二つの配列の差が，全要素においてあるしきい値（一般的にはマシンイプシロン）以下であれば等しいと判別することが多いと思います．
その場合でも，allを使えば

all(abs(a - b) < epsilon(a))

といった書き方ができます．

ユーザ定義演算子を用いれば，実数に対する $\approx$ 演算子が定義できます．

    interface operator(.approx.)
        procedure :: approximately_equal_real32
        procedure :: approximately_equal_real32_rank1
    end interface
contains
    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a
        real(real32), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real32

    pure function approximately_equal_real32_rank1(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: approx_eq

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank1

配列の形状チェック

二つの配列aとbを比較する場合，前提として配列の形状が同じでなければなりません．そのため，処理を安全側に振るのであれば，配列の値を比較する前に，配列の形状が一致することを確認しておく方がよいでしょう．

その場合も，allを使えば簡単です．

all(shape(a)==shape(b))

shape手続は，引数の配列形状を配列で返すので，allを使えば配列の各次元の要素数を個別に比較できます．

integer(int32) :: a(2, 3), b(3, 2)
print *, shape(a) ! [2, 3]
print *, shape(b) ! [3, 2]
print *, all(shape(a) == shape(b))

           2           3
           3           2
 F

配列の次元が異なる場合はコンパイルエラーが出て比較自体ができないので，注意してください．

複数の条件が.and.あるいは.or.で評価されている場合の論理演算の簡略化

色々な条件が関係する場合，if文で判定しようとすると，複数の条件を.and.または.or.で評価する必要がでてきて，条件式が非常に長くなり，見通しも悪くなります．

全ての条件が.and.で評価されている場合，あるいは.or.で評価されている場合は，allあるいはanyを用いて簡略化できます．

allやanyは，引数に論理型配列を渡せばよいので，配列構成子[]と組み合わせれば簡単に複数の条件を考慮できます．

例えば，二つの配列の添字上下限を比較する場合，allを使わないと，条件がかなり長くなります．

if (lbound(a, dim=1) == lbound(b, dim=1) .and. &
    lbound(a, dim=2) == lbound(b, dim=2) .and. &
    ubound(a, dim=1) == ubound(b, dim=1) .and. &
    ubound(a, dim=2) == ubound(b, dim=2)) then
    ! 添字の上下限が全て一致した場合の処理
end if

これは2次元配列の例ですが，Fortran規格では15次元配列（Coarrayを使う場合は15-corank）まで認められているので，最大で30個の条件を.and.で評価しなければなりません．

allを用いれば，非常に簡潔になります．

if (all([lbound(a) == lbound(b), &
         ubound(a) == ubound(b)])) then
    ! 添字の上下限が全て一致した場合の処理
end if

複数の条件が.or.で評価されている場合は，anyを用います．

ただし，簡単になるからと全ての条件をallやanyでまとめればよいかというと，そうは考えていません．条件の意味のまとまりごとに考えて，allでまとめる場合と.and.を用いる場合を使い分けた方がよいでしょう．

配列内にNaNが存在しているかの確認

Fortranの面白い機能にelemental属性があります．elemental属性を付与された手続は，定義上はスカラ値のみを引数に取りますが，配列を受け取れるようになります．配列の全要素に対して，手続内で定義された演算を実行し，配列を返します．

例えば，組込手続のsinはelemental functionに分類されており，リファレンス上は実数（あるいは複素数）のスカラ値を引数に取るようになっていますが，配列を渡すことも可能です．

real(real32), parameter :: PI = acos(-1.)
print *, sin(0.), sin(PI)               ! 0.0000000  -8.7422777E-08
print *, sin([0., 0.5, 1., 1.5, 2.]*PI) ! 0.0000000   1.0000000  -8.7422777E-08  -1.0000000   1.7484555E-07

実数がNaNかどうかをFortran標準に従って判別するには，ieee_arithmeticモジュールで定義されているieee_is_nan関数を用います．ieee_is_nanもelemental functionなので，配列を渡すと全要素についてNaNかどうかを調べ，結果を配列で返します．

返ってきた配列のなかに一つでも.true.があればNaNが生じているということですがら，any関数を用いてやれば簡単に検証できます．

use, intrinsic :: ieee_arithmetic
real(real32) :: a(10, 10)
a(:, :) = 0d0

a(2, 3) = ieee_value(0._real32, ieee_quiet_nan)
if (any(ieee_is_nan(a))) then
    print *, "a has ", count(ieee_is_nan(a)), " NaN at ", findloc(ieee_is_nan(a), .true.)
end if

 a has            1  NaN at            2           3

下記のような，配列を引数にとり，その中に一つでもNaNがあれば.true.を返すような手続があると，数値計算では役に立つかも知れません．

    interface has_nan
        procedure :: has_nan_real32_rank1
        procedure :: has_nan_real32_rank2
    end interface
contains
    function has_nan_real32_rank1(array) result(has_nan)
        implicit none
        real(real32), intent(in) :: array(:)
        logical :: has_nan

        has_nan = any(ieee_is_nan(array))
    end function has_nan_real32_rank1

    function has_nan_real32_rank2(array) result(has_nan)
        implicit none
        real(real32), intent(in) :: array(:, :)
        logical :: has_nan

        has_nan = any(ieee_is_nan(array))
    end function has_nan_real32_rank2

.approx.演算子およびhas_nan手続のためのメタプログラミング

allとanyの使い方をいくつか紹介してきました．その中で，演算子や関数をいくつか示しました．引数は，実数のスカラ値か配列を取ります．Fortranでは，異なる引数を取る同じような手続を，総称名を利用してまとめることができます．しかし，個別の引数の型に対する手続は書く必要があります．

実行する処理自体はFortranの機能によって，配列次元を問わず全く同一に書くことができます．一方で，Fortranにはメタプログラミングの機能が無いため，同じ処理で引数の異なる手続を何個も書く必要があります．

配列の形状に関しては，Fortran 2018から加わったassumed-rank（次元引継？）仮引数の機能である程度簡略化は可能です．assumed-rank仮引数を使うには仮引数の配列形状に..を用います．

real(real32), intent(in) :: a(..)

実引数の形状に応じた処理を記述するには，select rank文を用います．

select rank (a)
rank (0) ! 実引数がスカラの場合
    has_nan = ieee_is_nan(a)
end select

これはselect typeと類似で，なかなか融通が利きません．これまで上で見てきたように，ieee_is_nanを用いて配列内のNaNを検出するには，配列の次元がいくつであったとしても，any(ieee_is_nan(a))と書けました．assumed-rankの機能を利用する場合，rankを一つずつ書いていく必要があります．

select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1)
    has_nan = any(ieee_is_nan(a))
rank (2)
    has_nan = any(ieee_is_nan(a))
rank (3)
    has_nan = any(ieee_is_nan(a))
rank (4)
    has_nan = any(ieee_is_nan(a))
rank (5)
    has_nan = any(ieee_is_nan(a))
rank (6)
    has_nan = any(ieee_is_nan(a))
rank (7)
    has_nan = any(ieee_is_nan(a))
rank (8)
    has_nan = any(ieee_is_nan(a))
rank (9)
    has_nan = any(ieee_is_nan(a))
rank (10)
    has_nan = any(ieee_is_nan(a))
rank (11)
    has_nan = any(ieee_is_nan(a))
rank (12)
    has_nan = any(ieee_is_nan(a))
rank (13)
    has_nan = any(ieee_is_nan(a))
rank (14)
    has_nan = any(ieee_is_nan(a))
rank (15)
    has_nan = any(ieee_is_nan(a))
end select

おそらくコンパイル時に配列次元が静的に決まっていないといけないことが理由だと思いますが，もう少し簡潔に

select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1,2,3,4,...)
    has_nan = any(ieee_is_nan(a))
end select

とか

select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank (1:15)
    has_nan = any(ieee_is_nan(a))
end select

とか

select rank (a)
rank (0)
    has_nan = ieee_is_nan(a)
rank default
    has_nan = any(ieee_is_nan(a))
end select

とかが通ってくれればよいのですが，前二つは規格上対応しておらず，後者はコンパイラによってはieee_is_nanがassumed-rank仮引数に対応していないようです．Intel Fortranではビルドができましたが，gfortranやNAG Fortranでは無理でした．

今までは，関数インタフェースも含めて同じ処理を何個も書かなければいけなかったことを思うと，大分便利になったと思います．

    pure function has_nan_real64_rank1(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(:)
        logical :: has_nan
        has_nan = any(ieee_is_nan(a))
    end function has_nan_real64_rank1
    pure function has_nan_real64_rank2(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(:, :)
        logical :: has_nan
        has_nan = any(ieee_is_nan(a))
    end function has_nan_real64_rank2

とはいえ，いちいち繰り返すのも大変なので，以前紹介したfyppを用います．

fyppは，Pythonによるプリプロセッサであり，（他の言語でも使えるものの，主に）Fortranのプリプロセッサとして開発されました．Pythonを使ったことがある人には，非常に馴染みやすいと思います．

変数の定義や処理については以前紹介した記事を参考にしてもらうとして，fyppを利用すると，assumed-rankを用いたNaNの検出手続が下記のように書けます．

    public :: has_nan

    interface has_nan
    #:for real_kind in REAL_KINDS
        procedure :: has_nan_${real_kind}$
    #:endfor
    end interface
contains
    #:for real_kind in REAL_KINDS
    pure function has_nan_${real_kind}$(a) result(has_nan)
        implicit none
        real(${real_kind}$), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        #:for rank in RANKS
        rank (${rank}$)
            has_nan = any(ieee_is_nan(a))
        #:endfor
        end select
    end function has_nan_${real_kind}$

    #:endfor

この中で出ているREAL_KINDSやRANKSはfyppの機能で定義できる変数であり，ソースの先頭で下記のように定義しています．

#! nubmer of bits of real type
#:set REAL_BITS = [32, 64]

#:if defined('real128')
#:set REAL_BITS = REAL_BITS + [128]
#:endif

#! real kinds
#:set REAL_KINDS = ["real{}".format(bits) for bits in REAL_BITS]

#! target array rank
#:set MAX_RANK = 15

#! rank
#:set RANKS = [rank for rank in range(1, MAX_RANK+1)]

これをfyppで処理すると，下記のようなFortranのコードが生成されます．生成されたコードに適切なモジュール名をつけてやれば完成です．

生成されたコード

    public :: has_nan

    interface has_nan
        procedure :: has_nan_real32
        procedure :: has_nan_real64
    end interface

contains
    pure function has_nan_real32(a) result(has_nan)
        implicit none
        real(real32), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        rank (1)
            has_nan = any(ieee_is_nan(a))
        rank (2)
            has_nan = any(ieee_is_nan(a))
        rank (3)
            has_nan = any(ieee_is_nan(a))
        rank (4)
            has_nan = any(ieee_is_nan(a))
        rank (5)
            has_nan = any(ieee_is_nan(a))
        rank (6)
            has_nan = any(ieee_is_nan(a))
        rank (7)
            has_nan = any(ieee_is_nan(a))
        rank (8)
            has_nan = any(ieee_is_nan(a))
        rank (9)
            has_nan = any(ieee_is_nan(a))
        rank (10)
            has_nan = any(ieee_is_nan(a))
        rank (11)
            has_nan = any(ieee_is_nan(a))
        rank (12)
            has_nan = any(ieee_is_nan(a))
        rank (13)
            has_nan = any(ieee_is_nan(a))
        rank (14)
            has_nan = any(ieee_is_nan(a))
        rank (15)
            has_nan = any(ieee_is_nan(a))
        end select
    end function has_nan_real32

    pure function has_nan_real64(a) result(has_nan)
        implicit none
        real(real64), intent(in) :: a(..)
        logical :: has_nan

        select rank (a)
        rank (0)
            has_nan = ieee_is_nan(a)
        rank (1)
            has_nan = any(ieee_is_nan(a))
        rank (2)
            has_nan = any(ieee_is_nan(a))
        rank (3)
            has_nan = any(ieee_is_nan(a))
        rank (4)
            has_nan = any(ieee_is_nan(a))
        rank (5)
            has_nan = any(ieee_is_nan(a))
        rank (6)
            has_nan = any(ieee_is_nan(a))
        rank (7)
            has_nan = any(ieee_is_nan(a))
        rank (8)
            has_nan = any(ieee_is_nan(a))
        rank (9)
            has_nan = any(ieee_is_nan(a))
        rank (10)
            has_nan = any(ieee_is_nan(a))
        rank (11)
            has_nan = any(ieee_is_nan(a))
        rank (12)
            has_nan = any(ieee_is_nan(a))
        rank (13)
            has_nan = any(ieee_is_nan(a))
        rank (14)
            has_nan = any(ieee_is_nan(a))
        rank (15)
            has_nan = any(ieee_is_nan(a))
        end select
    end function has_nan_real64

.approx.演算子の実装としていたapproximately_equal_real32_rank1等の手続については，assume-rank仮引数を用いるのは現状では止めた方がよいでしょう．
二つの仮引数のrankを同時に調べられないので，aのrankを調べた後bのrankを調べる必要があり，あまりスマートではありません．

    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(..)
        real(real32), intent(in) :: b(..)
        logical :: approx_eq

        select rank (a)
        rank (0); select rank (b); rank (0)
                approx_eq = abs(a - b) <= epsilon(0.); end select
        rank (1); select rank (b); rank (1)
                approx_eq = all(abs(a - b) <= epsilon(a)); end select
...

また，select rank内において，総称名でまとめられた手続を呼ぶと，引数の配列次元を正しく認識できず，コンパイルエラーになる事例も確認できました．そのため，assumed-rankの使用について，下記のように自身の基準を決めました．

• assumed-rank仮引数を用いてもよい
  • 引数が一つの手続
• assumed-rank仮引数は用いない方がよい
  • 引数が二つ以上
  • 引数が一つ以上で，select rank内で総称名でまとめられた手続を呼び出している

そのため，メタプログラミングでコードを生成する場合は，仮引数の配列次元が異なる手続を生成するようにします．

    public :: operator(.approx.)
    public :: are_same_shape

    interface operator(.approx.)
    #:for real_kind in REAL_KINDS
        procedure :: approximately_equal_${real_kind}$
        #:for rank in RANKS
        procedure :: approximately_equal_${real_kind}$_rank${rank}$
        #:endfor
    #:endfor
    end interface

    interface are_same_shape
    #:for real_kind in REAL_KINDS
    #:for rank in RANKS
        procedure :: are_same_shape_${real_kind}$_rank${rank}$
    #:endfor
    #:endfor
    end interface

contains
    #:for real_kind in REAL_KINDS
    pure function approximately_equal_${real_kind}$(a, b) result(approx_eq)
        implicit none
        real(${real_kind}$), intent(in) :: a
        real(${real_kind}$), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a-b)<=epsilon(a)
    end function approximately_equal_${real_kind}$

        #:for rank in RANKS
        pure function approximately_equal_${real_kind}$_rank${rank}$(a, b) result(approx_eq)
            implicit none
            real(${real_kind}$), intent(in) :: a${ranksuffix(rank)}$
            real(${real_kind}$), intent(in) :: b${ranksuffix(rank)}$
            logical :: approx_eq

            if(.not.are_same_shape(a, b))then
                approx_eq = .false.
                return
            end if

            approx_eq = all(abs(a-b)<=epsilon(a))
        end function approximately_equal_${real_kind}$_rank${rank}$

        #:endfor

        #:for rank in RANKS
        pure function are_same_shape_${real_kind}$_rank${rank}$(a, b) result(are_same_shape)
            implicit none
            real(${real_kind}$), intent(in) :: a${ranksuffix(rank)}$
            real(${real_kind}$), intent(in) :: b${ranksuffix(rank)}$
            logical :: are_same_shape

            are_same_shape = all(shape(a) == shape(b))
        end function are_same_shape_${real_kind}$_rank${rank}$

        #:endfor

    #:endfor

同様にfyppで処理すると，Fortranのコードが生成されます．異様に長くなるので，MAX_RANKを3として生成したコードを示します．

生成されたコード

    public :: are_same_shape
    public :: has_nan

    interface operator(.approx.)
        procedure :: approximately_equal_real32
        procedure :: approximately_equal_real32_rank1
        procedure :: approximately_equal_real32_rank2
        procedure :: approximately_equal_real32_rank3
        procedure :: approximately_equal_real64
        procedure :: approximately_equal_real64_rank1
        procedure :: approximately_equal_real64_rank2
        procedure :: approximately_equal_real64_rank3
    end interface

    interface are_same_shape
        procedure :: are_same_shape_real32_rank1
        procedure :: are_same_shape_real32_rank2
        procedure :: are_same_shape_real32_rank3
        procedure :: are_same_shape_real64_rank1
        procedure :: are_same_shape_real64_rank2
        procedure :: are_same_shape_real64_rank3
    end interface

contains
    pure function approximately_equal_real32(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a
        real(real32), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real32

    pure function approximately_equal_real32_rank1(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank1

    pure function approximately_equal_real32_rank2(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:, :)
        real(real32), intent(in) :: b(:, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank2

    pure function approximately_equal_real32_rank3(a, b) result(approx_eq)
        implicit none
        real(real32), intent(in) :: a(:, :, :)
        real(real32), intent(in) :: b(:, :, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real32_rank3

    pure function are_same_shape_real32_rank1(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:)
        real(real32), intent(in) :: b(:)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank1

    pure function are_same_shape_real32_rank2(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:, :)
        real(real32), intent(in) :: b(:, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank2

    pure function are_same_shape_real32_rank3(a, b) result(are_same_shape)
        implicit none
        real(real32), intent(in) :: a(:, :, :)
        real(real32), intent(in) :: b(:, :, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real32_rank3

    pure function approximately_equal_real64(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a
        real(real64), intent(in) :: b
        logical :: approx_eq

        approx_eq = abs(a - b) <= epsilon(a)
    end function approximately_equal_real64

    pure function approximately_equal_real64_rank1(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:)
        real(real64), intent(in) :: b(:)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank1

    pure function approximately_equal_real64_rank2(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:, :)
        real(real64), intent(in) :: b(:, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank2

    pure function approximately_equal_real64_rank3(a, b) result(approx_eq)
        implicit none
        real(real64), intent(in) :: a(:, :, :)
        real(real64), intent(in) :: b(:, :, :)
        logical :: approx_eq

        if (.not. are_same_shape(a, b)) then
            approx_eq = .false.
            return
        end if

        approx_eq = all(abs(a - b) <= epsilon(a))
    end function approximately_equal_real64_rank3

    pure function are_same_shape_real64_rank1(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:)
        real(real64), intent(in) :: b(:)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank1

    pure function are_same_shape_real64_rank2(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:, :)
        real(real64), intent(in) :: b(:, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank2

    pure function are_same_shape_real64_rank3(a, b) result(are_same_shape)
        implicit none
        real(real64), intent(in) :: a(:, :, :)
        real(real64), intent(in) :: b(:, :, :)
        logical :: are_same_shape

        are_same_shape = all(shape(a) == shape(b))
    end function are_same_shape_real64_rank3

まとめ

非常に便利なのにあまり使われていないallとanyの使い方について説明しました．

特に配列を多用して数値計算を行うFortranでは，配列の等価性の比較や配列形状の比較，配列内のNaNの検出といったところでの利用価値が高いと考え，それらの例も示しました．

Fortran 2018の機能であるassumed-rankについても言及しました．Fortranの強みである配列操作とは非常に相性が良さそうですが，コンパイラの対応状況なのか，実装の都合なのか，もう一息頑張って欲しいと思います．