fortran コマンドライン引数処理

今日はfortranでコマンドライン引数を処理するサブルーチンを作ってみたので貼っておく。

fortranは言語仕様ではコマンドライン引数へのアクセスができないけども、通常のcompilerではiargc関数、getargサブルーチンが拡張仕様として備わっているので、最近では問題ないのかなと思う。

* Fortran 2003に対応しているcompilerならば、下記の関数を使う。
count = command_argument_count(), CALL get_command_argument(i, arg)

インターフェースはperlの標準モジュール"Getopt::Std"を真似した。使い方も大体同じ。
http://perldoc.perl.org/Getopt/Std.html

例えば、"-f"とか、"-i"とかよくやるやり方でfortranに引数を渡せるというまぁまぁ便利なものになった。

$ ./a.out -f filename -i
 f:filename
 i:1
 s:0

ソースは下記のようなもの。例えば下記のように"f:is"というオプションリストを渡すと、fについては文字列が格納され、i、sにはフラグ（0:off,1:on）がoptsという文字列配列に返される。順番は"f:is"で指定した順に入る。

program getarg
    implicit none

    character(len=256) :: opts(3)

    call getopt("f:is", opts)

    print *, "f:", trim(opts(1))
    print *, "i:", trim(opts(2))
    print *, "s:", trim(opts(3))

    stop
end program getarg

! --------------------------------------------------
subroutine getopt(optlist, opts)
! --------------------------------------------------
! parse command line arguments specified with
! single-character swiches.
!
! Usage: 
! call getopt("f:is", opts)
!
! where, switches followed by ":" takes arguments, 
! and others work as flag.
!
! for example, "./a.out -f file -i" returns
!  opts(1): file
!  opts(2): 1
!  opts(3): 0
! --------------------------------------------------

    type opt
        character(len=1)   :: name
        logical            :: arg  ! flag opt takes arg or not.
        character(len=256) :: val
    end type opt

    character(len=*), intent(in)  :: optlist
    character(len=*), intent(out) :: opts(:)

    integer                :: i, j, nopts
    character(len=256)     :: argv
    character(len=1)       :: optname
    type(opt), allocatable :: opts_tmp(:)

    allocate(opts_tmp( size(opts) ))

    ! parse optlist
    nopts = 0
    do i =1,len_trim(optlist)
        optname = optlist(i:i)
        if (optname /= ":") then
            nopts = nopts+1
            if (nopts > size(opts)) then
                print *, "error: number of options exceeds the size of opts"
            endif

            opts_tmp(nopts)%name = optname
            opts_tmp(nopts)%arg  = .false.
            opts_tmp(nopts)%val  = "0"
        endif

        if (optname == ":") opts_tmp(nopts)%arg = .true.
    enddo

    ! check command line arguments
    do i = 1, iargc()
        call getarg(i, argv)

        if (argv(1:1) == "-") then
            ! serach option list
            do j = 1,nopts
                if (opts_tmp(j)%name == argv(2:2)) then
                    ! found option, set value
                    if (opts_tmp(j)%arg .eqv. .true.) then
                        call getarg(i+1, opts_tmp(j)%val)
                    else
                        opts_tmp(j)%val = "1"
                    endif
                endif
            enddo
        endif
    enddo

    do i = 1, nopts
        opts(i) = opts_tmp(i)%val
    enddo

    deallocate(opts_tmp)

    return
end subroutine getopt

fortran用行列関係の関数いくつかまとめ

fortranで行列関係の関数をいくつか作ったので、まとめておく。といっても、MATLABの良く使いそうなやつをいくつか真似して作ってみただけ。

fortranは数値計算は何も考えなくても早いけども、ちょっとしたことが不便すぎる。いろいろツールを作っておかないとそのちょっとしたことで時間を取られてしまうし、コードも汚くなってしまう。まだfortran77のソースも良く見るけど、fortran90, fortran95の機能を知りながら使う気には到底なれないなあ。

それにしてもMATLABで便利なのって、結局、転置行列X'とか行列積X*Yが表現できることだったり、変数がそのまま残ることだったりするような気がする（速度を出すにはテクニックが必要で、まだ全然使いこなせない）。

あとは、MATLABのsave関数的なものがあると便利そうではある。

! ---------------------------------------------------------
module matrix
! ---------------------------------------------------------
! inv_lapack
!   Ai = inv_lapack(A) : return inverse of A(m,n) using lapack.
!
! whos
!   call whos(X) : shows size, bytes, and type of X.
!                  only real(8) and integer are supported.
!   
! linspace
!   X = linspace(x0,x1,n) : return linearly spaced vector just same as MATLAB.
!
! eye, ones, zeros (same as MATLAB, Octave)
!   X = eye(n,[m])   : get identity matrix
!   X = ones(n,[m])  : get matrix of all ones
!   X = zeros(n,[m]) : get matrix of all zeros
!
! show_dmat
!   call show_dmat(X,[fmt]) : show real(8) X(n,m). 
!
! ---------------------------------------------------------
! Note:
!   build: gfortran matrix.f90 main.f90 -llapack
! ---------------------------------------------------------
! Created by sako  2013/05/14 
! ---------------------------------------------------------
    implicit none

    interface whos
        module procedure whos_d1, whos_d2, whos_i1, whos_i2
    end interface

contains
    ! --------------------------------
    ! estimate inverse of array A(m,n).
    ! Ai = inv_lapack(A)
    ! --------------------------------
    function inv_lapack(A) result (AI)
        real(8), intent(in) :: A(:,:)
        real(8) :: AI(size(A,1),size(A,2))
        integer, parameter :: NB = 64
        integer :: info
        integer :: m, n

        ! for lapack
        integer :: LWORK
        integer, allocatable :: IPIV(:) ! dimension N
        real(8), allocatable :: WORK(:) ! dimension LWORK

        ! get size of A(m,n)
        m = size(A,1)
        n = size(A,2)

        LWORK = n*NB
        allocate(IPIV(n))
        allocate(WORK(LWORK))

        AI = A
        call DGETRF(m, n, AI, m, IPIV(1:min(m,n)), info)  ! factorize
        call DGETRI(n, AI, m, IPIV, WORK, LWORK, info)  ! inverse

        deallocate(IPIV)
        deallocate(WORK)

        if (INFO==0) then
            !write(*,*) "succeeded."
        endif
        if (INFO/=0) then
            print *, 'error: failed to estimate inverse at "inv_lapack".'
        endif

        return
    end function

    ! --------------------------------
    ! show variable size, bytes, type
    ! like whos of MATLAB. (call whos(X))
    ! Size                     Bytes    Kind
    ! ====                     =====    ====
    ! 1x100                      800  double
    ! --------------------------------
    subroutine whos_d1(x)
        real(8), intent(in) :: x(:)
        call show_whos(1, size(x), kind(x), 'double')
        return
    end subroutine whos_d1

    subroutine whos_d2(x)
        real(8), intent(in) :: x(:,:)
        call show_whos(size(x,1), size(x,2), kind(x), 'double')
        return
    end subroutine whos_d2

    subroutine whos_i1(x)
        integer, intent(in) :: x(:)
        call show_whos(1, size(x), kind(x), 'integer')
        return
    end subroutine whos_i1

    subroutine whos_i2(x)
        integer, intent(in) :: x(:,:)
        call show_whos(size(x,1), size(x,2), kind(x), 'integer')
        return
    end subroutine whos_i2

    subroutine show_whos(n, m, b, valtype)
        integer, intent(in) :: n, m, b
        character(len=*), intent(in) :: valtype

        character(len=11) :: strn, strm
        character(len=16) :: strb

        write(strn,'(i11)') n
        write(strm,'(i11)') m
        write(strb,'(i16)') b*n*m  ! bytes

        write(*,'(a48)') "          Size                     Bytes    Kind"
        write(*,'(a48)') "          ====                     =====    ===="
        write(*,'(a11,"x",a11,1x,a16,1x,a7)') adjustr(strn), adjustl(strm), adjustr(strb), adjustr(valtype)
        return
    end subroutine

    ! --------------------------------
    ! return linearly spaced matrix. 
    ! just same as linspace of MATLAB.
    ! --------------------------------
    function linspace(x0,x1,nd) result(X)
        real(8), intent(in) :: x0, x1
        integer, intent(in), optional :: nd
        real(8), allocatable :: X(:)

        real(8) :: dx
        integer :: n, i, ierr
        
        if (present(nd)) then
            n = nd
        else
            n = 100
        endif

        ! error check
        if (n < 1) n = 1  ! force to be 1

        allocate(X(n), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "linspace".'
            stop
        endif

        ! special case for n==1
        if (n == 1) then
            X(1) = x1
            return
        endif

        dx = (x1-x0)/dble(n-1)
        do i = 1,n
            X(i) = x0 + dx*dble(i-1)
        enddo

        return
    end function linspace

    ! --------------------------------
    ! return idendity matrix E(n[,m])
    ! --------------------------------
    function eye(n,m) result(E)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: E(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(E(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_eye".'
            stop
        endif


        ! initialize
        E(1:n,1:n2) = 0.d0
        do i = 1,n
            E(i,i) = 1.d0
        enddo
        return
    end function eye

    ! --------------------------------
    ! return matrix X(n[,m]) of all ones
    ! --------------------------------
    function ones(n,m) result(X)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: X(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(X(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_ones".'
            stop
        endif

        ! initialize
        X(1:n,1:n2) = 1.d0
        return
    end function ones

    ! --------------------------------
    ! return matrix X(n[,m]) of all zeros
    ! --------------------------------
    function zeros(n,m) result(X)
        !implicit none
        integer, intent(in) :: n
        integer, intent(in), optional :: m
        real(8), allocatable :: X(:,:)
        integer :: i, n2, ierr

        if (present(m)) then
            n2 = m
        else
            n2 = n
        endif

        allocate(X(n,n2), stat=ierr)
        if (ierr /= 0) then
            print *, 'error: failed to allocate memory at "matrix_zeros".'
            stop
        endif

        ! initialize
        X(1:n,1:n2) = 0.d0
        return
    end function zeros

    ! --------------------------------
    ! show matrix A(m,n)
    ! --------------------------------
    subroutine show_dmat(A,fmtf)
        !implicit none
        real(8), intent(in) :: A(:,:)
        character(len=*), intent(in), optional :: fmtf
        character(len=25) :: fmtc
        integer :: m, n, i, j
        m = size(A,1); n = size(A,2)

        if (present(fmtf)) then
            fmtc = fmtf
        else
            fmtc = '(f13.4," ")'
        endif
        do i=1,m
            do j=1,n
                write(*,fmtc, advance='no') A(i,j)
            enddo
            write(*,*) ""
        enddo
        return
    end subroutine show_dmat

end module matrix

fortranの副プログラムからの返り値は参照渡しかどうか

ふと、fortranのfunctionとsubroutineではどのように返り値が返ってきているのか気になった。fortranの引数は参照渡しなので、恐らくsubroutineに関しては配列をallocateして返したりしても同じアドレスに違いない。

だけど、functionに関してはどうなのだろう？というのも、function内でallocateして配列を返すようなものを結構書いていたから。fortran95からはスコープから外れるとメモリを解放してくれることになっているので、deallocateなどしていない。

ということで、試してみた。

program test_allocate
    implicit none

    integer, parameter :: n = 1000000
    integer :: array1(n)
    integer, allocatable :: array2(:)
    array1(1:n) = 1

    ! function
    array1 = allocate_array(n)
    print *, loc(array1)

    ! subroutine
    call allocate_array_sub(n, array2)
    print *, loc(array2)
    stop
contains

subroutine allocate_array_sub(n, array)
    integer, intent(in) :: n
    integer, allocatable, intent(out) :: array(:)

    allocate(array(n))
    array(1:n) = 1

    print *, loc(array)

    return
end subroutine allocate_array_sub

function allocate_array(n) result(array)
    integer, intent(in) :: n
    integer, allocatable :: array(:)

    allocate(array(n))
    array(1:n) = 1

    print *, loc(array)

    return
end function allocate_array
end program test_allocate

すると、こんな結果に。

      140120599949328  ! function内でallocate
              6299808  ! functionの外
             39948288  ! subroutine内でallocate
             39948288  ! subroutineの外

ということで、試したのはgfortranだけれども、function内でallocateした配列は値のコピーが返り、subroutineから返る配列は参照渡しで返るみたいだ。本当に言語仕様なのかはよく分からない。

fortranのファイルオープン関数

fortranのファイルオープンをしようとすると、いつも引数の指定方法を忘れてしまう。cとかperlみたいにopenするのに慣れすぎた。

いちいち外部入力装置番号を覚えておかなくてはいけないとか、オプションも自由気ままなんて、あまりにカオスすぎる。

ということで、perlなどのようにファイルを開くことができるようなwrapperを作ってみた。これ以外で使う時は、scratchでファイルを作ってOut of Coreのプログラムを書く時くらいだろうか。バイナリなんかめったにいじらないし。

こんな感じで使える。

    funit = open_file('>', "test.dat")   !書き込み専用
    funit = open_file('<', "test.dat")   !読み込み専用
    funit = open_file('>>', "test.dat")  !追加書き込み専用
    funit = open_file('+<', "test.dat")  !読み書き可能
    funit = open_file('+>', "test.dat")  !読み書き可能（ファイル新規作成）
    funit = open_file('+>>', "test.dat") !読み書き可能追加書き込み

本体のモジュールが次のような感じ。未使用の外部入力装置番号を検索してそれを自動で割り当てて、あとはただのwrapper。未使用の外部入力装置番号のインデックスは、モジュールのprivateメンバ変数にしてしまっても良いかもしれない。

module file
    implicit none

    ! define error codes
    integer, parameter :: ERR_FOPEN = -1
    integer, parameter :: ERR_FMODE = -2
    integer, parameter :: ERR_FUNIT = -8

contains
    ! ----------------------------------------
    ! open file like 'c' and 'perl'.
    ! if success, returns unit number.
    ! if error, returns minus number.
    ! 
    ! write only:
    !   funit = open_file('>', file_name)
    !   funit = open_file('w', file_name)
    ! ----------------------------------------
    integer function open_file(mode, file_name)
        character(len=*), intent(in) :: mode
        character(len=*), intent(in) :: file_name

        integer :: ios
        character(len=7) :: fstatus
        character(len=9) :: faction
        character(len=6) :: fpos
        integer :: funit
        funit = get_file_unit()

        if ((mode == '<') .or. (mode == 'r')) then
            ! read only
            fstatus='old'
            faction='read'
            fpos='rewind'
        elseif ((mode == '>') .or. (mode == 'w')) then
            ! write only
            fstatus='replace'
            faction='write'
            fpos='rewind'
        elseif ((mode == '>>') .or. (mode == 'a')) then
            ! append only
            fstatus='unknown'
            faction='write'
            fpos='append'
        elseif ((mode == '+<') .or. (mode == 'r+')) then
            fstatus='old'
            faction='readwrite'
            fpos='rewind'
        elseif ((mode == '+>') .or. (mode == 'w+')) then
            fstatus='replace'
            faction='readwrite'
            fpos='rewind'
        elseif ((mode == '+>>') .or. (mode == 'a+')) then
            fstatus='unknown'
            faction='readwrite'
            fpos='append'
        else
            ! invalid mode
            open_file = ERR_FMODE
            return
        endif

        ! file open
        open(unit=funit, file=file_name, iostat=ios, status=fstatus, &
             action=faction, position=fpos)

        ! check file open error
        if (ios==0) then
            open_file = funit
        else
            open_file = ERR_FOPEN
        endif
        return
    end function open_file

    ! ----------------------------------------
    ! get file unit which has not opened.
    ! return minus if error found.
    ! ----------------------------------------
    integer function get_file_unit()
        integer :: i
        integer, parameter :: MAX_UNIT = 999
        logical :: fopened
        
        do i = 100, MAX_UNIT
            inquire(unit=i, opened=fopened)
            if (fopened .neqv. .true.) then
                get_file_unit = i
                return
            endif
        enddo

        ! not found unopened file unit.
        get_file_unit = ERR_FUNIT
        return    
    end function get_file_unit
end module file

Okada(1992) DC3Dの使い方

Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040]の式のfortran77プログラムが下記で公開されている。
http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html

ある場所（緯度経度）での変位を出力しようと思ったら、結構頭がこんがらがった。というわけで、wrapperを作ったので、プログラムをメモ代わりにアップしておく。断層パラメータの記載は一般的な定義のままで、位置は断層上端左端にした。ああ、ややこしかった。

注意としては、深さは正の値を入れること、長さLの場合はAL1=0, AL2=Lとすること、幅Wの場合は、AW1=-W, AW2=0とすること。
一般的な断層パラメータの定義が数式と少し違い（走向は北から反時計回りとか）、注意が必要ですね。すぐ忘れそうなので、やはりwrapperを作っておいた方が良いな。もう少しモジュールがたまってきたら、arでライブラリにしてしまいたい。

ちなみに、こんな感じで呼び出す。

test.f90:
program test
    use okada92_util
    implicit none

    type(disp) :: du
    type(flt)  :: f
    real(8) :: lat, lon, h

    integer :: i, j

    f%lambda = 1.d0   ! ずるしている
    f%mu     = 1.d0
    f%lat    = 32.7725d0
    f%lon    = 136.0282d0
    f%depth  = 5000.d0
    f%str    = 261.766d0
    f%dip    = 8.903d0
    f%rake   = 90.d0
    !f%rake   = 0.d0
    f%length = 100000.d0
    f%width  = 25263.533d0
    f%slip   = 5.d0

    do i = -50, 50
        do j = -50, 50
            lat = f%lat + dble(i)*0.1d0
            lon = f%lon + dble(j)*0.1d0
            h   = 0.d0

            du = flt_disp(f, lat, lon, h)
            write(*,'(F7.3,",",F7.3,3(",",F15.7))') lon, lat, du%e, du%n, du%u
        enddo
    enddo

    stop
end program test

モジュールは下記のようなもの。

okada92_util.f90:
module okada92_util
    type disp
        real(8) :: e, n, u
    end type disp

    type flt
        !real(8) :: alpha          ! (lambda+mu)/(lambda+2*mu) !αをインプットにした方が良かったかも
        real(8) :: lambda, mu    ! (lambda+mu)/(lambda+2*mu)
        real(8) :: lat, lon      ! latitude, longitude (deg) at the corner of upper edge
        real(8) :: depth         ! fault depth (positive is ...)
        real(8) :: str           ! strike angle (deg)
        real(8) :: dip           ! dip angle    (deg)
        real(8) :: rake          ! rake angle   (deg)
        real(8) :: length, width ! fault length, width (m)
        real(8) :: slip          ! slip (m)
    end type flt

    ! ellipsoid
    type ellip
        real(8) :: a, b, e2, f, fi
    end type ellip

    ! parameters
    real(8), parameter :: PI = 3.1415926535897932d0
    real(8), parameter :: D2R = PI/180.d0
    real(8), parameter :: R2D = 180.d0/PI
 
contains
    ! ----------------------------------------------------------
    type(disp) function flt_disp(f, lat, lon, h)
    ! ----------------------------------------------------------
    ! Calculates displacement at (lat, lon, h) due to a 
    ! rectangular fault. The equations are given by 
    ! Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040].
    !
    ! This function is a wrapper of 'DC3D' on the website:
    ! http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html
    !
    ! Input:
    !   f : type(flt), includes
    !    f%lambda, mu    ! (lambda+mu)/(lambda+2*mu)
    !    f%lat, lon      ! latitude, longitude (deg) at the corner of upper edge.
    !    f%depth         ! fault depth (positive is downward) 
    !    f%str           ! strike angle (deg)
    !    f%dip           ! dip angle    (deg)
    !    f%rake          ! rake angle   (deg)
    !    f%length, width ! fault length, width (m)
    !    f%slip          ! slip (m)
    !   lat, lon, h : latitude(deg), longitude(deg), h(m)
    !                 positive h is upward.
    ! Output:
    !   disp        : type(disp),
    !     disp%e         ! eastward disp. (m)
    !     disp%n         ! northward disp. (m)
    !     disp%u         ! vertical disp. (m)
    !
    ! Created:
    !   2013/04/29  Sako
    ! ----------------------------------------------------------
        type(flt), intent(in) :: f
        real(8),   intent(in) :: lat, lon, h

        ! inputs of DC3D
        real(8) :: a, x, y, z, d, dip
        real(8) :: al1, al2, aw1, aw2
        real(8) :: disl1, disl2, disl3

        ! output of DC3D
        real(8) :: ux, uy, uz
        real(8) :: uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz
        integer :: iret

        ! tmp
        real(8) :: dx, dy, t

        !real(8) :: alpha          ! (lambda+mu)/(lambda+2*mu)
        a     = (f%lambda+f%mu)/(f%lambda+2.d0*f%mu)
        d     = f%depth   ! 断層上端の深さ。深さは正の値にする。
        dip   = f%dip
        al1   = 0.d0
        al2   = f%length  ! 断層の長さ。
        aw1   = -f%width  ! 断層の幅。図を良く見るとAW1はマイナスということで結構紛らわしい。
        aw2   = 0.d0
        disl1 = f%slip * dcos(f%rake*D2R)  ! すべり量を各軸に分解。
        disl2 = f%slip * dsin(f%rake*D2R)
        disl3 = 0.d0

        ! 目的の緯度経度をローカル座標系のx, yに変換する。走向の定義が北からなので、ここもややこしい。
        call bl2xy(f%lat, f%lon, lat, lon, dx, dy)
        t = (-f%str+90.d0)*D2R;     ! rotate to local coordinate
        x = dx*dcos(t) + dy*dsin(t) ! t = -str+90
        y =-dx*dsin(t) + dy*dcos(t)
        z =-h
        
        !  SUBROUTINE  DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
        !*                AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
        !*                UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET) 
        call DC3D(a, x, y, z, d, dip, al1, al2, aw1, aw2, disl1, disl2, disl3, &
                  ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz, iret)

        if (iret == 0) then
            flt_disp%e = ux*dsin(f%str*D2R) - uy*dcos(f%str*D2R)  ! 忘れずに座標変換。走向の向きは北からだからこうなる。
            flt_disp%n = ux*dcos(f%str*D2R) + uy*dsin(f%str*D2R)
            flt_disp%u = uz
        else
            ! singular location
            flt_disp%e = 0.d0
            flt_disp%n = 0.d0
            flt_disp%u = 0.d0
        endif

        return
    end function flt_disp

    ! ----------------------------------------------------------
    subroutine bl2xy(lat0, lon0, lat, lon, dx, dy)
    ! ----------------------------------------------------------
    ! Calculates distance between reference station(lat0, lon0) and 
    ! station2(lat, lon). 
    ! 
    ! Reference: Ohno S.(1987), Geodetic Methods (in Japanese).
    !   see equations [1.59] and [1.62]. 
    !
    ! Input:
    !   lat0, lon0 : position of station1 (reference)
    !   lat,  lon  : position of station2
    !
    ! Output:
    !   dx, dy     : distance between stations (sta1 -> sta2).
    ! 
    ! Created:
    !   2013/04/29 Sako
    ! ----------------------------------------------------------
        implicit none

        real(8), intent(in)  :: lat0, lon0, lat, lon
        real(8), intent(out) :: dx, dy

        type(ellip) :: ell
        real(8) :: dlat, dlon, Bm
        real(8) :: a, fi, f, b, e2


        ! WGS84
        ell = ellipsoid('WGS84')
        a = ell%a
        fi= ell%fi
        f = ell%f
        b = ell%b
        e2= ell%e2

        dlat = (lat - lat0)*D2R
        dlon = (lon - lon0)*D2R
        Bm = dsin(dlat)
        dx = a*dcos(lat0*D2R) / dsqrt(1.d0-e2*dsin(lat0*D2R)**2.d0)*dlon
        dy = a*(1.d0-e2)*dlat*(1.d0+3.d0/4.d0*e2-3.d0/4.d0*e2*dcos(2.d0*Bm)+1.d0/8.d0*e2*(dlat)**2.d0*dcos(2.d0*Bm))

        return
    end subroutine bl2xy

    ! ----------------------------------------------------------
    type(ellip) function ellipsoid(ellip_name)
    ! ----------------------------------------------------------
        implicit none
        character(len=*), intent(in), optional :: ellip_name
        real(8) :: a, fi, f
        
        ! default is GRS80
        a = 6378137.d0
        fi= 298.257222101d0

        if (present(ellip_name)) then
            if (ellip_name == 'GRS80') then
                a = 6378137.d0
                fi= 298.257222101d0
            elseif (ellip_name == 'WGS84') then
                a = 6378137.d0;
                fi= 298.257223563d0;
            else 
                ! ...
            endif
        endif

        f = 1.d0/fi

        ellipsoid%a = a
        ellipsoid%fi= fi
        ellipsoid%f = f
        ellipsoid%b = a*(1.d0-f)
        ellipsoid%e2= 1.d0 - (1.d0-f)**2.d0

        return
    end function ellipsoid
end module okada92_util

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

逆行列の計算(fortran)

今日は逆行列の計算をしてみたのでメモ。

ソルバはLAPACKのDGETRF,DGETRIで、fortran90。

! --------------------------------
! calculate inverse of matrix A(m,n)
! --------------------------------
function inv_lapack(A) result (AI)
    implicit none

    real(8), intent(in) :: A(:,:)
    real(8) :: AI(size(A,1),size(A,2))
    integer, parameter :: NB = 64
    integer :: info
    integer :: m, n

    ! for lapack
    integer :: LWORK
    integer, allocatable :: IPIV(:) ! dimension N
    real(8), allocatable :: WORK(:) ! dimension LWORK

    ! get size of A(m,n)
    m = size(A,1)
    n = size(A,2)

    LWORK = n*NB
    allocate(IPIV(n))
    allocate(WORK(LWORK))

    AI = A
    call DGETRF(m, n, AI, m, IPIV(1:min(m,n)), info)  ! factorize
    call DGETRI(n, AI, m, IPIV, WORK, LWORK, info)  ! inverse

    deallocate(IPIV)
    deallocate(WORK)

    if (INFO==0) then
        !write(*,*) "succeeded."
    endif

    return
end function

! --------------------------------
! show matrix A(m,n)
! --------------------------------
subroutine show_dmat(A,fmtf)
    implicit none
    real(8), intent(in) :: A(:,:)
    character(len=*), intent(in), optional :: fmtf
    character(len=25) :: fmtc
    integer :: m, n, i, j
    m = size(A,1); n = size(A,2)

    if (present(fmtf)) then
        fmtc = fmtf
    else
        fmtc = '(f13.4," ")'
    endif
    do i=1,m
        do j=1,n
            write(*,fmtc, advance='no') A(i,j)
        enddo
        write(*,*) ""
    enddo
    return
end subroutine show_dmat

呼び出しはこんな感じ。ついでにinterface付き。

program test
    implicit none

    interface
        function inv_lapack(A)
            real(8), intent(in) :: A(:,:)
            real(8) :: inv_lapack(size(A,2),size(A,1))
        end function inv_lapack

        subroutine show_dmat(A,fmtf)
            real(8), intent(in) :: A(:,:)
            character(len=*), intent(in), optional :: fmtf
        end subroutine show_dmat
    end interface

    real(8) :: A(2,2)
    real(8) :: AI(2,2)
    integer :: i, j

    A(1,1) = 1.d0
    A(1,2) = 2.d0
    A(2,1) = 4.5d0
    A(2,2) = 6.0d0

    write(*,'("A=")')
    call show_dmat(A, "(f9.3)")
    AI = inv_lapack(A)
    write(*,'("AI=")')
    call show_dmat(AI, "(f9.3)")

    stop
end program test

実行すると、まあうまくいったぽい。

$ gfortran -o test test.f90 -llapack
$ ./test
A=
    1.000    2.000 
    4.500    6.000 
AI=
   -2.000    0.667 
    1.500   -0.333