fortran90の構造体へのアクセス速度 メモ

前のブログ(http://sak12.blogspot.com/2013/04/matlab.html)で、MATLABの構造体の扱いには注意しなければいけないということをメモしたけれど、これまでfortranで使っていたものは大丈夫なのか不安になった。

 で、同じように２乗を計算してみた（1000万回）。
どうやら、fortranでは構造体へのアクセス速度は、特に何も考えなくても大丈夫そう。
構造体をよく使っていたのでドキッとしたのだけれど、一安心。

あと、MATLABが速いことを改めて実感。MATLABでドット演算子を使った場合、fortranのforallと速度はほぼ変わらない。
MATLABは遅いと良く言われたりしてるけど、うまく使えばかなり高速になりそうだし、きちんとMATLABの勉強しようかな。

program main
    implicit none

    type struct1
        real(8), allocatable :: x(:)
    end type struct1

    type struct2
        real(8) :: x
    end type struct2

    integer, parameter :: n = 10000000
    type(struct1) :: s1
    type(struct2) :: s2(n)
    real(8)       :: x(n), y(n)
    integer       :: stime(8), etime(8)
    integer :: i

    allocate(s1%x(1:n))

    forall (i=1:n)
        x(i)    = 2.d0
        s1%x(i) = 2.d0
        s2(i)%x = 2.d0
    end forall

    ! x ただの配列同士
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = x(i) * x(i)
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (x):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    ! s1　構造体が配列を持つ
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = s1%x(i) * s1%x(i)
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (s1):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    ! s2　構造体配列
    forall (i=1:n) y(i) = 0.d0
    call date_and_time(values=stime)
    forall (i=1:n)
        y(i) = s2(i)%x * s2(i)%x
    end forall
    call date_and_time(values=etime)
    write(*,'("elapsed time (s2):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0

    stop
end program main

結果はこのとおり。

elapsed time (x):    0.0250(s)
elapsed time (s1):    0.0260(s)
elapsed time (s2):    0.0260(s)

----- おまけ -----
perlでやってみると、こんなに時間がかかる・・・。まあ数値計算はやっちゃだめということですかね。そりゃそうか。

use strict;
use Time::HiRes qw( gettimeofday tv_interval );

my $N = 10_000_000;
my ($i, @x, $y);
foreach $i (1 .. $N) {
    push(@x, 2);
}

# estimate square of @x.
my $t0 = [gettimeofday];
for ($i=0;$i<$N;$i++) {
    $y = $x[$i]*$x[$i];
}
my ($seconds, $microseconds) = gettimeofday;
my $elapsed = tv_interval($t0, [$seconds, $microseconds]);
print "elapsed time: $elapsed\n";

$ perl test.pl 
elapsed time: 1.291174