fortran66のブログ

fortran について書きます。

配列の一括操作

Fortran での配列

Fortranの配列は、等質なデータがメモリ上に連続的に等間隔に並んでいるものに対応しています。何かの次に何かがあるという順序数に対応する一次元リスト構造のイメージより、等質なものが規則的に並列に並んでいるイメージに近いと思います。

Fortran でポインタが嫌われるのは、ポインタではデータの連続性が成り立たない場合も許されるからです。連続したメモリ上にデータがないとコンパイラの最適化が阻害されてしまいます。

real, target :: a(n), c(n, n)
real, pointer :: b(:), d(:, :)
b => a(1:n:2)
d => c(2:n - 1, 2:n -1)

ここで、b,dの指し示しているメモリーは、とびとびの番地に対応していて連続していません。

allocatable 属性は pointer 属性ほど柔軟性がありませんが、連続メモリーが保障されるので、同内容のことが出来るなら allocatable を優先すべきだということになります。また allocatable はメモリ解放も自動でやってくれる利点もあります。

Fortran での MAP 演算

Fortran での配列データに対する MAP 操作は、要素型(elemental)関数や要素型サブルーチンで実現されます。これはスカラー引数に対して定義しておけば、同じ型の配列引数を取ることが出来て結果も同数のデータの配列で返すというものです。露わに MAP 等を書かなくとも同様の働きをします。

Fortran90 からはほとんどの組み込み関数・サブルーチンが要素型になっています。また Fortran95 からはユーザー定義の関数・サブルーチンを要素型に定義することが出来るようになりました。

例 組み込み関数

real :: x(n), y(n)
integer :: i
x = [(i / 10.0, i = 1, n)] ! [0.1, 0.2, ... , 1.0] 
y = sin(x) * cos(x) 

例 ユーザ定義要素型関数

      elemental pure real function hsin2x(x)
        real, intent(in) :: x
        hsin2x = sin(x) * cos(x)
        return
      end function hsin2x

Fortran での Filter 演算

Fortran での配列データに対する FILTER 演算は、pack 関数によって実現されます。
例 

real :: x(100)
reaal, allocatable :: y(:)
..
y = pack(x, mask = x > 0.0)

ここで、y には x のうち正の値のみが入ります。mask 引数は論理型配列を取ります。『mask =』は省略できます。pack の出力する配列のデータ数は不定個ですが、Fortran2003 の allocatable 配列の自動割り付け機能により、サイズを気にせず代入できます。

本来 pack 関数は、unpack 関数と組みになっていて、mask 条件を満たすデータを選び出して何らかの処理をしたのち元の配列に戻すために使われることを意図したものですが、元々の意図を離れて FILTER 演算のために使われることが多いようです。

例 元々の使い方の例

real :: x(n), y(n)
..
forall (i = 1:n, x(i) > 0.0)
  y(i) = i * f( x(i) ) 
end forall

ここで、全要素に演算が投機的に実行されて、mask 条件を満たすところだけ代入操作が行われることが暗黙に意図されています。(もともと80年代の並列計算機のモデルに適した機能として導入されているので。)しかし mask 条件が満たされる割合が小さいときは、無駄な計算が多すぎることになります。そこで pack で mask 条件を満たすものだけをあらかじめ抽出しておいて、必要な計算だけを実行し、unpack で戻すことですることで計算の無駄を省けます。

real :: x(n), y(n)
real, allocatable :: tmp1(:), tmp2(:)
tmp1 = pack(x, x > 0.0)
forall (i = 1:size(tmp1))
  tmp2(i) = i * f( tmp1(i) )
end forall
y = unpack(tmp2, x > 0.0)

Fortran でのリダクション演算

Fortran ではいくつかのリダクション演算が用意されています。リダクション演算とは、配列として複数個のデータを引数としてうけて、単一の結果に減らし(レデュースして)て返す演算です。マスクも使えます。
例 総和演算 sum

real :: x(n), s
..
s = sum(x, mask = x > 0.0)

例 総積演算 product

real :: x(n), s
..
s = product(x, mask = x /= 0.0)

例 計数 count

real :: x(n)
integer :: k
..
k = count(x > 0.0)

例 論理和・論理積 any, all

real:: x(n), y(n)
..
if ( any(x > y) ) exit
..
if ( all(x > 0.0) ) return
..

count, any, all では、引数が MASK そのものになります。