Fortranの配列は一味違う？

「天体の位置計算 増補版」という書籍がある。この中で、位置の計算式を行列を使って計算する方法が紹介されている。行列計算は配列同士の計算なので配列への入出力といったことを理解しておかないとうまくいかない。Fortranの配列は列優先という考え方である。これはどういうことかというと、以下の３×３の配列を見てみよう。行列的に表現してある。

\begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}

１〜９までの数値を配列に代入したと考えてみると、１行目に１〜３、２行目に４〜６、３行目に７〜９の数値が入ると考える。普通の考え方であり問題ない。では、これをFortranで書いて実行してみよう。

テストコード matrix-1.f90

program matrix
    implicit none

    integer :: a(3,3)
    integer :: i,j
    
    data a /1, 2, 3, &
            4, 5, 6, &
            7, 8, 9/

    do i = 1, 3
        do j = 1, 3
            write(*, fmt='(i2)', advance='no') a(i, j)
        end do
        write(*,*)
    end do
    stop
end program matrix

結果予想

 1 2 3
 4 5 6
 7 8 9

実行結果

 1 4 7
 2 5 8
 3 6 9

考えたとおりにならず違った結果になる。出力された値の並び方を見てみるとよくわかるが、これが列優先ということになる。これを意識して配列への代入を考えないと、行列計算したときに正しくない結果をもたらすことになる。総和をとる場合や配列の値を使って行列計算などを行う場合は、data文を使わずに配列への代入を次のように記述する。

a(1, 1:3) = (/1, 2, 3/)
a(2, 1:3) = (/4, 5, 6/)
a(3, 1:3) = (/7, 8, 9/)

data文を使う場合は、先に述べたように列優先を意識してdata文に値を記述しないと結果が正しくならない。Fortranは配列の考え方がちょっと違うのである。

計算例

「天体の位置計算 増補版」55頁に以下の行列計算式が掲載されている。この行列計算式の数値を配列を指定して代入する場合とdata文を使って配列へ代入した場合の計算結果の違いを見てみよう。

\[ \begin{eqnarray} \begin{pmatrix} L_{2} \\ M_{2} \\ N_{2} \end{pmatrix} &=& \begin{pmatrix} \begin{array}{rrr} 0.99997 60 & -0.00635 67 & -0.00276 33 \\ 0.00635 67 & 0.99997 98 & 0.00000 88 \\ 0.00276 33 & -0.00000 88 & 0.99999 62 \end{array} \end{pmatrix}\begin{pmatrix} \begin{array}{rrr} 0.56473 73 \\ -0.54140 89 \\ 0.62285 49 \end{array} \end{pmatrix} \\ &=& \begin{pmatrix} \begin{array}{rrr} 0.56644 42 \\ -0.53781 35 \\ 0.62441 78 \end{array} \end{pmatrix} \end{eqnarray} \]

テストコード matrix-2.f90

program matrix
    implicit none

    double precision :: a(3,3), b(3,1), result(3,1)

    a(1, 1:3) = (/0.9999760, -0.0063567, -0.0027633/)
    a(2, 1:3) = (/0.0063567,  0.9999798, -0.0000088/)
    a(3, 1:3) = (/0.0027633, -0.0000088,  0.9999962/)

    b(1, 1) = 0.5647373
    b(2, 1) = -0.5414089
    b(3, 1) = 0.6228549

    ! 行列の積を求める
    result = matmul(a, b)

    write(*, *)
    write(*, '("L = ", f10.7)') result(1, 1)
    write(*, '("M = ", f10.7)') result(2, 1)
    write(*, '("N = ", f10.7)') result(3, 1)
    write(*, *)

    stop
end program matrix

実行結果

正しく計算されている。

L =  0.5664442
M = -0.5378136
N =  0.6244178

data文を使った実行結果

配列への代入をdata文に置き換える。結果は正しくならない。

data a /0.9999760, -0.0063567, -0.0027633, &
        0.0063567,  0.9999798, -0.0000088, &
        0.0027633, -0.0000088,  0.9999962/
L =  0.5630033
M = -0.5449933
N =  0.6212967

99