2013-04-01

配列を intent(out) で返したつもりが一部しか返っていなかった

Fortran 90で配列の上下限を省略してコロンを使うことについて にやや関連した話題。

まず、呼び出し側の配列として

  real(8) :: x(201,3)

というものがあって、subroutine 側に対応する

  real(8), intent(out) :: c(51,3)

というのがあった。んで、subroutine を

  call sub( c = x )   !! intent(out)

という風に call していた。このとき、

  x(1:51,1:3) = c(1:51,1:3)

のように値が入ってくれるのだろう、と、勝手に思っていたが、恐るべきことに実際には、

  x(1:51,1) = c(1:51,1)

となっていた…(戦慄)。これがもし

  x(1,1) = c(1,1)

だけだったらさすがにすぐ気づいたんだろうけど、なぜか第1次元は全部埋まっていたので発見が遅れた。


解決方法を4つ考えた。

1. ダミー配列を用意する
つまり
  real(8) :: tmp(51,3)
というのを用意しておき、
  call sub( c = tmp )   !! intent(out)
  x(1:51,1:3) = tmp(1:51,1:3)
とする。わかりやすいが、計算リソースを食う。

2. call 時に上下限を指定する。
つまり
  call sub( c = x(1:51,1:3) )   !! intent(out)
ということ。ちょっとわかりにくい。

3. subroutine 側の配列のサイズを呼び出し側に合わせる
つまり
  real(8), intent(out) :: c(201,3)
とする。使わない分まで確保するため、メモリの無駄。だが、今まで一番多くやってたのはたぶんこれ…。実際には 201 という数値ではなくて、変数になっていて、一緒に intent(in) して割付、というのを最近はけっこうやっていた。

4. subroutine 側の配列を assumed shape にする。下限を指定する。
つまり
  real(8), intent(out) :: c(:,:)
としておく。これが一番一般的なのかな。パフォーマンスに若干影響するとかしないとか見かけた気がするけど気のせいかもしれない。この方法は楽だけど、下限が1でない配列を使う時は下限がずれる可能性があり、危険。下限だけ指定するのがいいかもしれない。

参考:
http://d.hatena.ne.jp/rontas/20130113/1358047540
このブログの Fortran カテゴリは自分が知らなかったことばかりでためになった。


これまでの自分は基本的に 1 か 3 を使っていたが、徐々に 4 に移行すべきなのかもしれない……。というか、実際にちょっと使い始めてみたけど、2 と 4 は併用するものかも。

No comments: