Fortrantips

From Ecal

Revision as of 17:15, 11 September 2012 by Kino (Talk | contribs)
(diff) ← Older revision | Current revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Contents

行番号はもう不要

error処理は大抵iostat,statが用意されている。

open(10,file=filename,iostat=ierr)
if (ierr.ne.0) then
...
endif
read(*,'(i)',iostat=ierr) a
allocate(ia_v(N),stat=ierr)

formatもcharacterとして定義できる。

character(200):: form
form='(a,1x,I4)'
write(*,form) 'fixedform?freeform?',-1
form='*'
write(*,form) 'fixedform',-1 !--- これはNG。 write(*,*)と書く。

loopからも抜けられる。

i_loop: do i=1,n
    j_loop: do j=1,m
        if (j.eq.1) cycle j_loop
        if (func(i,j)>0) exit i_loop
    enddo j_loop
enddo i_loop

直前のloopから抜けるだけならi_loop, j_loopと書いたlabelは不要。

open文

actionでread/writeの制御が可能。

open(unit,file=filename,action='read')

pointerとallocatable

まとめ

繋ぎかた

allocatableだけでなく、 pointerに対してもallocateが可能である。

integer,pointer:: ip_v(:), ip_v2(:)
integer,allocatable:: ia_v(:), ia_v2
integer,allocatable,target:: iat_v(:)
allocate(ip_v(N))
allocate(ip_v2(N))
allocate(ia_v(N))
allocate(ia_v2(N))
allocate(iat_v(N))

pointer=>pointer、pointer=>allocatable,targetは繋げる。それ以外はNG。

pointerをつないだ先にメモリーがあるかは分からない

real(8), pointer:: rp_v(:),rp_v2(:)
allocate(rp_v(100))
rp_v2 => rp_v
write(*,*) associated(rp_v) ! .true.
write(*,*) associated(rp_v2) ! .true.
nullify(rp_v)
write(*,*) associated(rp_v) ! .false.
write(*,*) associated(rp_v2) ! .true.
deallocate(rp_v2) ! --- memoryをfreeして、pointerをnullにする。
write(*,*) associated(rp_v) ! .false.
write(*,*) associated(rp_v2) ! .false.
allocate(rp_v(100))
rp_v2 => rp_v
deallocate(rp_v)
write(*,*) associated(rp_v) ! .false. 
write(*,*) associated(rp_v2) ! .true. pointerの方は相手がallocateされたままなのかどうか分からない

下に示すが、move_alloc()でallocatable変数に対して同じことが可能。

pointerにallocateができてpointer=>pointerと出来るならallocatable,targetは何のためにあるのか。

pointerを繋ぎ変えて

ip_v2 => ip_v
...
ip_v => ip_v2
...

などとやっていくとそのうちわけが分からなくなるので、

allocate(iat_v(N))
ip_v =>iat_v ! pointer=>allocatable,target

とtargetを定義してallocateする変数は必ずallocatableに、 pointerの方はallocatable,targetを指し示すだけにしたほうが誤りが少ない。allocateする変数がallocatableならば変数のlifetimeが終わると自動的にdeallocateしてくれるのでmemory leakしないはず。


pointer,allocatable変数の初期値

  • pointer変数の初期値は他の変数と同様undef。p=>null()もしくはnullify(p)でnull化される。
  • allocatable変数は自動的に初期値を持っており最初は.not.allocated。

p=>null()はnullify(p)と同じ。 歴史的には、最初は初期化するところに書けないnullify(p)だけだったが、後にnull()ができている。

(device driverを書くCなどの言語と違い、fortranの場合は定義した変数、配列は全て使うはずなのでundefという状態がないように初期化されていない変数はcompilerが全て初期化するようf90の段階で規定しておいた方がよかった。C言語、man mallocのovercommit_memory参照)

追記: pointer:: p(:)とlocal変数を定義してもsubroutine中でp=>...と必ずpを何かに繋ぐことがコンパイラにとって明らかな場合はp=>null()と初期化できないことがある。コンパイラのバグだと思う。--Kino 08:48, 17 May 2011 (JST)

pointerとallocatableとの違いは?

  • allocatableは変数を宣言したsubroutineが終了すると自動的にdeallocateされる。(変数がなくなるとmemoryも開放される。)
  • pointerは変数を宣言したsubroutineから抜けても自動的にdeallocateされない。(変数がなくなってもmemoryは保持したまま。)

これらは

  • allocatableはallocatedで
  • pointerはassociatedで

書き込むmomeryがあるか判断する。(同じ名前で良かっただろうに。)

正しく言うと、場合分けがあって、 saveされたallocatable変数は自動的にdeallocateされない。しかし、それはpointerの振舞で、例えば、

integer,allocatable,save:: ia_v(:)
if (.not.allocated(ia_v)) allocate(ia_v(N))
integer,pointer,save:: ip_v(:)=>null()
if (.not.assciated(ip_v)) allocate(ip_v(N))

は同じである。(初期の決めかたがまずかったせいでこのあたりが混乱している。多くの人が分からなくて当然。)


pointerはもともとは多くのexampleがそうなっているようにlist pointerとして用意されたようだ。しかし、実際にはほとんどの人は配列全体を指し示すために使う。

pointerとallocatable、その2

以下はallocatableがtype内に書ける場合の話になる。

現在では、allocatableは「initialize=null化とfinalize=deallocateを自動的に行うpointer」、 つまり、使用範囲に制限を加えたことで有用な機能を得たpointerである。pointer⊃allocatable。

初期の頃はpointerはのろいと言われたそうだが、pointerがarray全体を示している限りそれはもう忘れていい。

少なくともintel fortranはallocatable属性があるとsubroutineが終了するときにallocatedになっていれば自動的にdeallocateされる、と規定している。 deallocateされる時はfortran2003で変数を定義しているsubroutineが終了するとき=変数がなくなるときと規定された。 しかし、fortran2003の準拠が完全かはコンパイラによる。 実際のimplementationから判断するしかない。


type内のallocatableとpointer

昔からtype内にpointerを書けたので配列を定義できていた。allocatableを置けて初めてtype内に配列を定義できるようになったわけではない。 fortran2003規格に沿ったcompilerはtype内でallocatableと書ける。

module m_foo
implicit none
type t_foo
  real(8), allocatable :: ra_v(:)
  real(8), pointer :: rp_v(:)
end type
end module m_foo

subroutine setfoo(s_foo,n)
  use m_foo
  implicit none
  type(t_foo):: s_foo
  integer ::n
!.begin
    allocate(s_foo%ra_v(n))
    allocate(s_foo%rp_v(n))
    s_foo%ra_v=100.0d0
    s_foo%rp_v=200.0d0
!.end
end subroutine setfoo

subroutine main
  use m_foo
  implicit none
  type(t_foo):: s_foo   
!.begin
    call setfoo(s_foo,10) ! type内変数をallocateする場合subroutineをuseしなくていい。
    write(*,*) 'ra_v',s_foo%ra_v(:) ! 100.0d0を10個表示
    write(*,*) 'rp_v',s_foo%rp_v(:) ! 200.0d0を10個表示
!.end
end subroutine main

(allocatableの意味からすると、allocateしたsubroutine内でのみ配列を定義してend subroutineでdeallocateしてもいいが、) 見た目pointerと同じ振る舞いをする。時間に従い変数がどうなるかを表すと

subroutine mainに入る。s_fooの宣言。ra_vのnull化。
subroutine setfooに入る。 s_foo%ra_p,%rp_vのallocation
subroutine mainに戻る。s_foo%ra_p,%rp_vはallocateされたまま。
subroutine mainの最後。s_foo%ra_pはdeallocateされる。%s_fooは消える。%rp_vはdeallocateされないで残る。(%rp_vにアクセスはできなくなる。)
subroutine mainを抜ける。

pointerだと定義したsubroutineの終了などで変数がなくなってもallocateしたmemoryはdeallocateないが、allocatableだと変数がなくなるときにallocateしたmemoryも自動的にdeallocateされる。

なお、subroutineから抜けて変数がundefになることとその変数が指していたメモリーがdeallocateされることとは違うので混同しないよう注意。


最後に二重にallocateした場合の違い

  • allocatableで宣言しallocateした変数をdeallocateなしにもう一度allocateするとどうなるかというと、エラーが起き異常終了する。なお、このエラーはstat=で取れ終了しないようにできる。
  • pointerで宣言しallocateした変数をdeallocateなしにもう一度allocateしてもエラーは起きない。この場合前回allocateした領域のpointerが残っていないと参照不能になりprogramが実行されている限り開放もできなくなる。

type内のallocatableとpointer(2)

不思議なことにtype内変数にtarget属性を書けない。

type foo 
  real(8),allocatable,target:: v(:)
end type

は不可であり

type foo2
  real(8),allocatable:: v2(:)
end type

type(foo2) :: fff
real(8),pointer:: pv(:)

として

pv=> fff%v2

はできない。

subroutineでpointerをallocateするにはuseする

module m_foo
contains
  subroutine foo(n,a)
    integer::n
    real(8),pointer :: a(:)
!...begin
      allocate(a(n))
      a=100.0d0
!...end
  end subroutine foo
end module m_foo
...
  use m_foo, only: foo !--- useしておく
  real(8),pointer :: a(:)
    call foo(10,a) !--- 関数内でallocateする
    write(*,*) a !-- 100.0d0が10個表示される。

最新の規格ではpointerではなくallocatableでも受け付ける。useしないとSEGVが起きる。typeの中にpointer,allocatableを入れてtype定義した変数をprocedureに渡す場合はuseしなくていい。

理由を考えると意味が分かる。 fortranがreference(=C言語のpointer)渡しであることはご存じだと思う。subroutine内でメモリーを取得して親routineへ持っていくためにはreferenceのreferenceを渡さないといけない(単なるreferenceとは違う)ためである。

present、optionalを使用するためにはuseする

--Kino 16:40, 2 September 2010 (JST) 同じファイルに書かれていても以下はNG

     program aaa
     call abc(b)
     end program aaa
     subroutine abc(b,ax)
     integer,intent(inout):: b
     real(8),intent(in),optional::ax
     print *,'pxxx=',present(ax)
     end subroutine abc

objective-cやpythonなどsmalltalkの考え方を取り入れた言語と違い fortranではsubroutineが実行時に自分自身の定義情報をもっておらず、 program aaaとsubroutine abc間の情報のやりとりの仕方はあらかじめユーザーが記述してやらないといけない。 上の例はできそうな気がするが、

     subroutine abc(b,ax,n)
     integer,intent(inout):: b
     real(8),intent(in),optional::ax
     integer,intent(in):: n
     end subroutine abc

こんな風になっているかもしれない。useしないでcall abcとするとだめそうなことが理解できると思う。

pointerを使うにはuseしなければいけなかったのと同じで、 program aaaの方でoptionalとした変数を渡しているという情報をユーザーが付加して(ここがuseやinterfaceの役割)callしないと機能しない。(たまたま動く場合もある。)


moduleでなくてもintereface文を書いてもよい。ただし、書き換えたときにユーザーがinterface文も書き換える必要があり忘れてしまいがちであり,userが作成し引数を書き換える可能性があるroutineにuseでなくinterface文を用いるのはバグの原因となるので勧めない。interfaceによる引数の定義はblas,lapackなど引数の書き換えがないlibraryに留めた無難である。


optional

open文のようにfoo(file=,status=)などと書きたい。

module m_foo
contains
subroutine foo(name, age, school, degree)
  character(*)::name
  integer,optional::age
  character(*),optional:: school
  integer,optional:: degree
    if (present(age)) then
    ...
    endif
end subroutine foo
end module m_foo

use m_foo
  call foo('Taro', age=3)
  call foo('Ichiro', school='Tokyo Univ.', degree=3, age=22)
  call foo('Hanako',20,'Keio Univ.',1) ! 順序通りならば変数名を書かなくてもよい。
  call foo('Tora',90) ! 順序通りならば変数名を書かなくてもよい。

全部optionalも可。使用する際はuseすること。


type変数の初期値

fortranなのでinitialization routineがなく、(C++になれていると奇妙に思えるが)type definition内で初期値を書くことができる。

module m_foo
  type s_foo
    integer,pointer:: ipv(:)=>null() !--- 初期値 =NULL
    integer,pointer:: ipv2(:) !--- undef
    real(8):: r=1000.d0
#ifdef ALLOW_ALLOCATABLE_IN_MODULE
    integer,allocatable:: iav(:) !--- 書けるとしたら自動的に.not.allocatedが初期値。
#endif
  end type
end module m_foo
program test
  use  m_foo
  type(s_foo):: sv
  integer,pointer :: v_ipv(:) !--- undef
#ifdef ALLOW_ALLOCATABLE_IN_MODULE
    write(*,*) allocated(sv%iav) !--- 書けるとしたら .F. 
#endif
    write(*,*) associated(sv%ipv),associated(sv%ipv2),associated(v_ipv),sv%r
    !--- associated(sv%ipv2),associated(v_ipv)の結果は不定。
end program test

初期値があるtypeを使用するためには

type(s_foo),save:: v_foo

とsave属性が必要かもしれない。

pointerのpointer

typeを絡めるとpointerのpointerを作ることも可能

type t_rv1
  real(8),pointer:: v(:)
end type
type t_foo
  type(t_rv1),pointer:: rp_v(:)
end type

  type(t_foo) :: s_foo
  integer,parameter::n=10
  integer:: i,isizeofv(n)
  allocate(s_foo%rp_v(n))
  isizeofv=(/5,5,5,5,5,10,10,10,10,100/)
  do i=1,n
    allocate(s_foo%rp_v(i)%v(isizeofv(i)))
  enddo

これによりindexによりsizeが異なる多次元配列を作ることが可能で、それぞれのsizeがかなり違うことが分かっている場合にmemoryの無駄がなくなる。 今の場合s_foo%rp_v(i)%vは以下のようになる。

i  size(v)
----------
1  xxxxx
2  xxxxx
3  xxxxx
4  xxxxx
5  xxxxx
6  xxxxxxxxxx
7  xxxxxxxxxx
8  xxxxxxxxxx
9  xxxxxxxxxx
10 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx   (x100)

fortranではこれができないと思っている人が多いが、上の通り可能である。

allocatableのallocatableでも可。

move_alloc

fortran2003ではallocatable変数をpointerの様に繋ぎ替えることも可能。 下はintel fortranのpointerを示す拡張機能%locも使っている。

        program test
        integer,allocatable:: a(:),b(:)
        allocate(a(5))
        write(*,*) %loc(a),a
        a=3
        allocate(b(10))
        b=0
        b(1:5)=a
        write(*,*) %loc(b),b
        call move_alloc(from=b,to=a)  ! a=>b相当
        write(*,*) %loc(a),a,allocated(b)
        end

結果(を整形し直したもの)は

 436527152  0  0  0  0  0
 436535680  3  3  3  3  3  0  0  0  0  0
 436535680  3  3  3  3  3  0  0  0  0  0 F

実体はa=>bをしているのがわかる。 (だとするとmove_alloc後はbはallocateされていない状態へ変化するのは当然である。)

pointerの場合、同じmemoryをいくつものpointer変数が指し示すことが可能だが、 allocatableの場合はmemoryとallocatable変数が一対一に対応している。


...どうしてこれだけuserが作った関数のようにcallで呼ぶのだろう?

prototype宣言

interface
 subroutine foox(r,n)
  integer,intent(in):: n
  real(8),intent(inout):: r(n)
 end subroutine foox
end interface

userが作成し引数を書き換える可能性があるroutineにmoduleとuseでなくinterface文を用いるのはバグの原因となるので勧めない。interfaceによる引数の定義はblas,lapackなど引数の書き換えがないlibraryに留めた無難である。

関数の総称、interface & module procedure

同じくinterfaceで書く。 (こんな文法でいいのだろうか。)

こういうものを作る場合、moduleの中に入れるのが普通でしょう。その場合module procedureと書く。

module m_set_func
  interface set_n
    module procedure r_set_n, i_set_n
  end interface
contains
  subroutine r_set_n(n,d)
    integer::n
    real(8):: d(n)
    ...
  end subroutine r_set_n
  subroutine i_set_n(n,d)
    integer::n
    integer:: d(n)
    ...
  end subroutine i_set_n
end module m_set_func

program test
  use m_set_func
  real(8):: dv(10)
  integer:: iv(10)
    call set_n(10,dv) !---
    call set_n(10,iv) !--- real, integerどちらでも受けられる。
end program test

free formatの継続行

free formatの継続行は行末の&だが、なぜかこういう書き方も可能。

      write(*,*) i,j &
     & ,k   ! simply ignore & in this line.

free formatでは2行目の最初の&は認識されない。 PHASEはこういう継続行で書かれている。

implicit none in module

module test
  implicit none
contains
  subroutine foo()
    i=10   ! causes a compiler error 
  end subroutine foo
end module test

上に書くとimplicit noneがmodule内全てに適用される。従ってsubroutine内で型宣言がないとエラーになる。

同様に変数も全体に適用されてしまう。 うっかり

program test
  implicit none
  integer::i
  main_loop: do i=1,...
    call foo(...)
  enddo main_loop
contains
  subroutine foo(...)
     do i=1,...   !--- このiはsubroutine fooのlocal変数iでなくmain_loopのi。
     ...
     enddo 
  end subroutine foo
end program test

と書くと正しくloopが回らないので注意。(処理系によってはloop変数は自動的にlocal変数と解釈してうまく動いてしまうが、うまく動かないのが正しい。)

program test
  implicit none
  call main   !--- ここには変数を書かない方が安全
contains
  subroutine main
    implicit none  !--- 重複していても構わない
    integer::i
    do i=1,...
      call foo(...)
    enddo
  end subroutine main
  subroutine foo(...)
     integer:: i
     do i=1,...   
     ...
     enddo 
  end subroutine foo
end program test

contains前には変数を書かず、変数がある場合はこんな形の別subroutineにした方が安全。

intent

ユーザーは普通procedure側の見方で理解しており、こう思っているはず。procedure内で

  • in: 必ずreferする。redefしない。
  • out: 必ずredefする。始めはundef
  • inout: 値を持ち必ずreferする、必ずredefする。

しかし、これは誤りである。ifort 11.1のmanualのよるとintentはcaller側からの定義になっている。

  • in: procedureに値を渡している。これは変数がundefでないということ。procedure内でredefしてはいけない
  • out: procedure内でredefしないcallerに戻さないといけない。procedure内で始めはundefだからprocedureでreferするにはまずredefする。
  • inout: procedureに値を渡し、redefしてcallerに戻すという両方のことができる。(「できる」というのは何だろうか。)

(undefとは何かというと

program test
integer:: i  !-- initializeしていないのでundef
i=0    !--- ここで値が確定,undefでなくなる。

ということ。)

ユーザーの直観と違うところがいくつかある。 例えば、intent(in)はundefではないということを規定しているだけで、procedure側が値を使うかどうかを規定していないため、次の使い方は正しい。

subourinte foo(i,j)
integer,intent(inout)::i
integer,intent(in)::j
i=1   !--- inoutなのにiをreferしていない。
! jにアクセスしない。
end subroutine foo

caller側から定義を書くことで実装を簡単化しているのだが、まだ複雑でさらに簡単化している。

subourinte foo(j)
integer,intent(out)::j
j=j+1 !---undefのときにjをreferとしている。
end subroutine foo

これはundefのときにreferしていて誤りであるが、これはwarningがでることもなくコンパイルできる。undefの場合にreferする前に必ずredefするかどうかは、例えばif文があるともうコンパイル時にはわからないことで、それなら一切checkはしないというスタンスなのかもしれないが、残念ながらmanualの記述とは違う。

ifortは「-enable-diag sc」がintentを厳密に解釈してwarnigを出すoptionのはずだがversionにより解釈しない。「-check uninit」でもsubroutine内のintent(out)による変数のundef指定などを実行時にcheckしない。

intentに関してcompilerが実際にやっていることは

  • in: procedure内でredefしない
  • out: procedure内のどこかで必ずredefする。
  • inout: procedure内のどこかで必ずredefする。

のみで、ユーザーの期待に反してintentはちゃんとcheckされていない



double precisionとreal(8)

同じ意味であるが、どちらがいいのだろうか。

多くのコンパイラはrealをreal(8)(real(kind=8))と読みかえるoptionを持っている。これを使用するためにはrealと定義しないといけない。

一方、gfortranはreal(10)というintel cpuが持つ80bitの演算が可能であり、real(dp)と精度をinteger::dpで指定したコードを書いておけば対応が容易である。dpはmoduleに書いてuseしておけばいい。 program中では全てreal(8)で書くのが(第一原理計算では)普通であるが、速度や精度を考えるとroutine毎変更した方がいい場合もある。私はreal(dp)という書き方を勧める。--Kino 16:56, 2 September 2010 (JST)

総称関数

歴史的な理由でcomplex<->real,imagの間は対称的な操作になっていない。どういう意味なのか以下に記す。

  • complex(kind=?):: zに対してreal partを出すreal(z),imaginaryを出すimag(z)やaimag(z),conjugateを返すconjg(z),の方は総称(Generic)関数である。つまり型に応じて自動変換してくれる。

しかし、実数から複素数を作るcmplxは総称関数ではない。

real(8):: a,b,z
z=cmplx(a,b)

cmplxはreal(4)(system defaultがreal(4)として)の型であるから、期待した動作にならない。 これは、a,bをreal(4)に直して、complex(4)をつくり、complex(8)に型変換する。 real(8)でcomplex(8)を作りたければ

z=cmplx(a,b,kind=8)
z=dcmplx(a,b)

となる。プログラムするときに注意が必要。

これを見ると、real(kind=dp)という書き方にしてcmplx(a,b,kind=dp)という書き方の方がわかりやすいし、 そもそも実数から複素数変換をする総称関数を作った方がいいだろう。


  • 型変換でなくsqrt,sin,expなどの実際に計算をする関数はだいたい総称関数が存在する。

例えば、complex(8)のsqrtはzsqrtやcdsqrtとなるが、総称関数sqrtがあるのでそれを使えばいい。

しかし実数から複素数を作る場合はまたややこしい。例えば、complex = sqrt( real )は必ずしも保証されていない。

complex(8) :: z
z=sqrt(-1.0d0)

はintel fortranではOKでもgfortranではNGとなる。

gfortranでコンパイルできるようにするには

z= sqrt( dcmplx(-1.0d0,0.0d0) )

と書く。

write(実数,*)

write(6.0,*)'abc'

動きます。

linuxのstacksize

最近のlinuxはstacksizeをunlimitedにできる。process毎でなくsystem全体でstackを共有するそうだ。

できるのだが、http://d.hatena.ne.jp/fd0/20080110/p1 にある通りx86ではstackがheapに衝突してもエラーを出さないそうだ。 またuserがstacksizeをunlimitedにできるかどうかはadministratorが管理することで一般的に許されることではないので、stacksize=unlimitedと思ってprogramを組んではいけない。

また自動的に処理系がつくるarrayをstackに置くか、heapに置くかは処理系に依存する。

  • intel fortranには"-heap-arrays <size>"というstackでなくheapにarrayをとるoptionがある。すべてをheapにとるのではなく<size>以上の配列のみをheapに取ることに注意。
  • gfortranはdefaultでheapに取る(と書いてあるwebpageがある)。


関連して、linuxの最大の欠点はkernelがuserが実際に使用するメモリー量を制限できないことである。 仮想記憶があるのでmallocでもしもの時に備えて巨大なメモリーをheapに取っておいて、実際は少しだけ使う、ということが 許されるが、この実際に使っているメモリー量の制限がkernelからできない。 計算機センターではjob schedulerが大量に実メモリーを使おうとするprocessを停止させているが、あれはkernelの機能をつかっているのではなく 別processが使用メモリー量を監視して停止させているだけである。だからときどき停止が間に合わない。 一度diskとの間で巨大なmemory swapが起きるとノード全体の速度が落ち、速度低下自体の原因を調べることすら困難になる。その場合の対処法はほぼリセットしかない。

ほとんどのデータセンターがlinuxを使っているのでこの実装は急務のはずだがlinux communityとしてはその重要性は低いということらしい。なお、AIXでは実際に使用するメモリー量は制限可能となっている。

stacksizeとsize(, dim)関数

subroutine foo(vr,N,M,...)
real(8):: vr(N,M)
integer:: N,M
...

として

complex(8):: vc(N,M)
call foo(dreal(vc),N,M,...)

という書き方が可能ですが、配列vcに対してdreal(vc)と自動配列をとるのは止めた方がいい。 普通dreal(vc)はstackのとられるかheapに取られるかは処理系による。vrのサイズN*Mが大きくなる可能性がある場合、stackにとられると、heapの方のメモリーが残っていてもstackのメモリーがとれなくなりSEGVになるかもしれない。

なお、allocateすればheapに取られるので

complex(8):: vc(N,M)
real(8),allocatable:: vr(:,:)
allocate(vr(size(vc,dim=1), size(vc,dim=2))
vr=dreal(vc)
call foo(vr,N,M,...)
deallocate(vr)

とでもするしかない。ただ、コンパイラが馬鹿でdreal(vc)の実行でstackに配列を一度とり直してしまうと上と同じ理由でSEGVの可能性がある。

また、arrayのsizeくらいsystem側で知っていて欲しいので

complex(8):: vc(N,M)
real(8),allocatable:: vr(:,:)
allocate(vr(size(vc,dim=1), size(vc,dim=2))
vr=dreal(vc)
call foo(vr,size(vc,dim=1),size(vc,dim=2),...)
deallocate(vr)

としたいところだが、size(vc,dim=1)がinteger(4)なのかinteger(8)なのかは環境(OSとコンパイラとオプション)によりからsubroutine foo側でN,Mが正しく読めない可能性があります。一度同じ型(今の場合はinteger)に落とさないといけません。

complex(8):: vc(N,M)
integer:: isize1_vr, isize2_vr
real(8),allocatable:: vr(:,:)
allocate(vr(size(vc,dim=1), size(vc,dim=2))
vr=dreal(vc)
isize1_vr=size(vc,dim=1)
isize2_vr=size(vc,dim=2)
call foo(vr,isize1_vr, isize2_vr,...)
deallocate(vr)

もちろん、allocateするところでisize?_vrを使えます。

上を全部あわせて、どんな大きなN,M,どのコンパイラ、どの計算機でも動くコードを書くならば

complex(8):: vc(N,M)
integer:: isize1_vr, isize2_vr
real(8),allocatable:: vr(:,:)
isize1_vr=size(vc,dim=1)
isize2_vr=size(vc,dim=2)
allocate(vr(isize1_vr, isize2_vr))
do j=1,isize2_vr; do i=1,isize1_vr;vr(i,j)=dreal(vc(i,j)); enddo; enddo
call foo(vr,isize1_vr, isize2_vr,...)
deallocate(vr)

という面倒なコードにしないといけません。 (integerをinteger(8)にするかどうかはコンパイラオプションで設定できる。)

'なお、size(vc,dim=?,kind=?)という書き方も可能。'

private,public,protected

module foo
integer,private:: i_pri  ! module内でのみアクセス可能。
integer,public:: i_pub   ! module内外でアクセス可能。(default)
integer,protected:: i_pro ! module外ではreadのみ可能。
end module foo

command line argumentと環境変数

fortran2003で規定された。

       program test
       implicit none
       character(200):: name,value
       integer:: n,i
       n= command_argument_count()
       do i=1,n
           call get_command_argument(i,value)
           write(*,*) i,trim(value)
       enddo

       name="SHELL"
       call get_environment_variable(name,value)
       write(*,*) trim(name),"=",trim(value)
       end 

結果

$ ./a.out a  c bc xxx
           1 a
           2 c
           3 bc
           4 xxx
 SHELL=/bin/bash

時間

ようやく実時間を計れるようになった。

        program test
        real(4):: t1,t2,diff
        integer:: i1,i2,i3,irate,imax
        call cpu_time(t1)
        call cpu_time(t2)
        write(*,*) t1,t2  ! cpu time

        call system_clock(i1)
        call system_clock(i2,irate,imax)
        diff = i2-i1
        if (i1>i2) diff=imax-i1+i2
        write(*,*) dble(i2-i1)/dble(irate)  ! real time
        call system_clock(i3,irate,imax)
        if (i2>i3) diff=imax-i2+i3
        write(*,*) dble(i3-i2)/dble(irate)  ! real time

        end

実時間を何に使うのと思うかもしれないが、 例えば、thread並列をするとcpu_timeだとthread全部のcpu時間を返すので、実際にどの程度時間短縮できたのか分からない。 MPIの動的なload balancingをするときに実時間の測定結果を用いるなど。

しかし、system_clockは使いにくい。どうしてrealを返してはいけないのだろうか。 cpu_timeのように差の形で実時間を求めたい場合はCを呼ぶしかない。

double second_()
{

        struct timeval tp;
        struct timezone tzp;
        int i;

        i = gettimeofday(&tp,&tzp);
        return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
}

使っていないfile descriptorが欲しい

file unit番号をuserが指定できるfortranの実装がそもそもおかしいが、それを認めるとして、 人のコードは何番が使われていないのか特に分かりにくい。

logical:: L
integer:: fid
inquire(fid,opened=L)

としてLの戻り値を見れば、fid番のfileがopenされているかわかる。

他にも

character(100):: filename
inquire(file='somefile',number=fid)  ! somefileがopenされていればそのunit番号
inquire(unit=10,name=filename)    ! unit=10でopenされているfile名

などの情報が取れる。

equivalence

fortranは読み手にとって複雑で間違いやすい書き方ができる仕様が多い。これもまたその例。

      integer:: n1,n2,n3,ngabc(3)
      equivalence (n1,ngabc(1)),(n2,ngabc(2)),(n3,ngabc(3))

他人がわからなくなるので無秩序な使用は勧めないが、 abbreviationとしてノートの変数にしたがった利用に止めた方がいい。

しかしこのようなabbreviationならcppで#defineしておいてもいい。

#define n1 ngabc(1)

など

また、これはallocatable or pointerには使えない。

コメント行とopenMP

fixed formのfortranにおいて

C

で始まる行はコメントであることは当然知っていると思う。 openMPを有効にするとopenMPのdirective行である

C$OMP

が有効になるだけでなく

C$<space>

で始まる行も有効なソースになる。fortranソースのlabelを書く場合は<space>無しで数字でも有効になる。

C$       write(*,*)'ok'
C$100    write(*,*)'ok100'

は有効なソースとしてコンパイルされるということだ。 ただ有効になるのはC$を除いてlabel部分に妙な記載がない場合だけで、 例えば、

C$ xxx   write(*,*)'okxxx'

はlabel部分に妙なものが入っていて有効にならない。全部有効にすればいいのにどうしてこうなったのだろうか。 それはともかく、'C$ 'の存在理由はomp_get_num_thread()などopenmpの時だけ使う行はここに書けということだ。 しかし、'C'で始まる行がソースでないのに、'C$ 'だとソースになる可能性があるとはあまりに紛らわしい。 (言わなくてもいいと思うが

C$$
C$$$

はopenMPの時もコメントである。)

free formだと(fixed formでも)

!$OMP

だけでなく、

!$<space>

もopenMPだと有効になるのだが、やはり'!$ 'は紛らわしい。cppの使用を前提として仕様を決めた方がsimpleで良かっただろうに。

昔は

D

の行は「debug」の時はソースとしてコンパイルしたそうだ。ifortでは この「debug」の指定は'-DD' とする。

omp_get_num_threads

よくある間違い。

parallel sectionでないところでは

c$    write(*,*) 'num_threads=',omp_get_num_threads()

の結果は1です。

parallelになってはじめてthreadの数が1でなくなります。

c$    integer:: omp_get_num_threads
c$    integer:: omp_get_thread_num
c$omp parallel
c$omp master
c$    write(*,*) 'num_threads=',omp_get_num_threads()
c$omp end master
c$omp barrier
c$    write(*,*) 'thread_id=',omp_get_thread_num()
c$omp end parallel

(自分で書いていてもc$がコメントなのかわからなくなる。他にCで始まるコメントがあれば尚更。)

ompの型宣言

c$      use omp_lib, only get_thread_num()

cpu_time

cpu_timeはelapsed timeであるのでthread 並列をしているとthread全部足してしまう。

system()

gfortran/fortran 2008 stdでは

implicit none
integer :: ret,system
ret=system(string)

もしくは

call system(string,status=ret)

同じprogram中では関数形かsubroutine形かどちらかしか呼べない。


intel fortranでは

implicit none
integer(4):: ret,system
ret=system(string)

なのでsystem()の戻り値をretの型へ自動型変換してくれるようにifportを呼ぶこと

use ifport
implicit none
 ...
integer:: ret
ret=system(string)

今はifportがなくても、integer(4)でなくても偶然動いているかもしれないが、型の違いはversion upや新しいマシンで突然エラーになることがあるから注意が必要。

型チェック

gfortranですべてを一つのファイルにするとsubroutine/callの型チェックをしてくれる。

unformatted file

      open(fid,file="x.dat",form="unformatted")

      str="data1"
      write(fid) str
      write(fid) n1,n2
      v(1:n1,1:n2)=-999.0
      write(fid) v(1:n1,1:n2)

      write(*,*) v(1,1)

      str="temporary"
      write(fid) str
      n1=2;n2=2
      write(fid) n1,n2
      v(1:n1,1:n2)=10.0
      write(fid) v(1:n1,1:n2)

      write(*,*) v(1,1)

      close(fid)

で書き込んで

      open(fid,file="x.dat",form="unformatted")
      n1=1;n2=1;
      do i=1,2
      read(fid,iostat=ios) str
      if (ios.ne.0) write(*,*) 'error in str',ios
      read(fid,iostat=ios) n1
      if (ios.ne.0) write(*,*) 'error in n1',ios
      read(fid,iostat=ios) v(1:n1,1:n2)
      if (ios.ne.0) write(*,*) 'error in v',ios
      write(*,*) str
      write(*,*) n1,n2
      write(*,*) v(1:n1,1:n2)
      enddo
      close(fid)

で読み込んでもちゃんと読める。

 data1               
          20           1
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000       -999.000000000000     
  -999.000000000000       -999.000000000000     
 temporary           
           2           1
   10.0000000000000        10.0000000000000    

formatted fileと同様に、次のread文は次の行頭から値を読んでくれる。

formatted fileでは「改行」の印を認識していたが、unformatted fileのx.datをhex dumpをするとこうなっていて

0000000 0014 0000 6164 6174 2031 2020 2020 2020
0000020 2020 2020 2020 2020 0014 0000 0008 0000
0000040 0014 0000 000c 0000 0008 0000 0780 0000
0000060 0000 0000 3800 c08f 0000 0000 3800 c08f
*
0003660 0780 0000 0014 0000 6574 706d 726f 7261
0003700 2079 2020 2020 2020 2020 2020 0014 0000
0003720 0008 0000 0002 0000 0002 0000 0008 0000
0003740 0020 0000 0000 0000 0000 4024 0000 0000
0003760 0000 4024 0000 0000 0000 4024 0000 0000
0004000 0000 4024 0020 0000

0014 0000が行頭を示している。(x86)

かなり紛らわしい点、初期を書くとsave属性になる。

      integer:: is=100

と書くとこの変数isは自動的にsave属性がつく。並列化の妨げになる可能性があるため注意。 関数に入ったときに毎回is=100としてくれるわけではない。つまり、上と

      integer:: is
      is=100

とは違う。c言語で期待される動作と違う。

Personal tools