FortranでNaNを見つける方法

はじめに

Fortranで数値計算していると、NaNがどこで発生したか見つけたいと思うことが多々あると思います。
Intel fortran (ifort)にはNaNを見つけるisnanという関数があるのですが、他のコンパイラには無いようです。

方法としては、以下のような感じになります。

1. 基本的な方法

NaNは何と比較してもfalseなので

if(a/=a)then
  write(*,*) "a is NaN"
endif

もしくは何と掛けてもNaNなので、

if(a*0.0/=0.0)then
  write(*,*) "a is NaN"
endif

とすれば良いです。

2. コンパイラ最適化の回避

じつは、上の２つのコードはコンパイラ最適化によってうまく動作しなくなる可能性が高いです。
a/=aや、a*0=0といった演算は分かりきっていて意味のない演算だとコンパイラが判断してしまうので、このif文自体が取り除かれてしまいます。
これを回避するためには、

zero = 0.0
call checknan(a,zero,flag)

subroutine checknan(a,zero,flag)
  implicit none
  double precision,intent(in) :: a,zero
  logical,intent(out) :: flag
  flag = .false.
  if(a*zero/=zero)then
    flag = .true.
  endif
end subroutine

のように外部サブルーチンに分離してやると、ifが取り除かれなくなります。
コンパイラを騙しているようなものですね…
インライン展開を行うと、これでも取り除いてしまうコンパイラがあるかもしれないので、checknanのサブルーチンは別ファイルに分離するほうが安全です。

3. IEEE関数を使う方法

@implicit_none さんに教えて頂きました(コメントに有ります)。

use, intrinsic :: ieee_arithmetic

if(ieee_isnan(a))then
  write(*,*) "a is NaN"
endif

でいけるようです。便利ですね。
ただし、gcc 5以前ではコンパイル出来ないようです。
CentOS 7やRHEL 7 (デフォルトがgcc 4.8.5)のマシンはまだまだ多いので、可搬性を考えると要注意です。