Fortran (90以上) を書く時に気をつけると良いこと

ここ1年以上メイン言語は Fortran90 です.ここでは,Fortran そのものに関する賛否はさておき,(個人的に) Fortran でこのように書かれていると良いだろうという話をしていきます *1

implicit none を必ず書く

これが何かわからなくてもいいのでとにかく書いてください.
Fortran は恐らく歴史的経緯で暗黙の型宣言がデフォルトで有効になっています.そのため,例えば変数のスペルミスで勝手に何らかの型付きで別の変数として自動宣言されてしまうなど,多くの場合バグの原因になります.

全ての変数を module 経由で渡さない

数値計算コードでよくある話だと思うのですが,入力パラメータが非常に多くなってしまうため,関数の引数として受け取るのではなく,グローバル変数 (Fortran では module を使う) 経由で引数を関数に渡す,という書き方をしている人が少なからずいるようです.これは多くの場合問題になります.例えば,1つ前の計算結果が入力となるような計算の場合,module で受け取ると,ユーザの方で入力変数と出力変数の値を入れ替える必要があります.これが配列と考えると,コピーコストは無視できない可能性があります *2

全てを引数で渡した方が良いという意見もあるかもしれませんが,やはり宣言が長くなってしまうので,記述の面倒さと最適化のし易さを天秤にかけると以下のようになってるとバランスが取れるのではないか,と思います.

  1. スカラの変数・定数は module 経由で渡しても良い
  2. 配列は subroutine/function の引数として渡す
  3. 【例外】subroutine/function 中,配列が read-only (値の書き換えが発生しない) 場合は,下記の通りにする
    • 配列のサイズが定数で決定されている (入力によらず一定) 場合,module 経由で渡しても良い
    • 配列のサイズが変数で決定されている (入力によって変化) 場合,必ず subroutine/function 経由で渡す
  4. 【例外】もし subroutine/function が重要な計算カーネルで,最適化が必要な場合は全て subroutine/function の引数として渡す

最後の例外は,Fortran ではなく C, C++ を用いた最適化 (つまり SIMD 命令を使って最適化する可能性があるとき) のためです.Fortran と C 言語の連携は Fortran2003 から規格されており,それ以前ではなし崩し的というか,各コンパイラ独自でやっていたようです.現在は変数に BIND(C) というアトリビュートをつけると C 言語で見た時の変数名が設定できたりします.

配列は関数の引数として渡す

前述の,入力と出力が入れ替わるような計算のために考えることです.引数経由で渡すことで,入力配列と出力配列を手動で入れ替える必要がなく,subroutine/function の呼び出し時に入力配列と出力配列を入れ替えるだけでよくなります.例えば,配列変数 a, b がそれぞれ入力と出力になる場合,呼び出し元では if 文を使ってこのように書けばいいでしょう.

do i=1,100
  if (i % 2 == 0) then
    call f(a, b)
  else
    call f(b, a)
  end if
end do

subroutine f(a, b)
  implicit none
  integer,intent(in)  :: a(:)
  integer,intent(out) :: b(:)
  ! do something
end subroutine

また,この書き方にはもう1つ利点があります.Fortran の配列には "形状指定配列" という仕様があり,「引数として渡された配列は,配列の元サイズを超えない範囲で,次元およびサイズを自由に変更できる」という仕様です.つまり,元々1次元配列だったものを,サイズを超えない範囲で多次元配列として扱っても良く,その逆も然りということです.
例えば,シミュレーションしている実空間を意識した計算を行う場合には,対象とする実空間 (ex. 2, 3, 4次元) と次元を揃えて計算したいが,それ以外の場合は1次元的に扱って計算する,といったことが可能です.起きそうな事例を考えると,元々は1次元配列として定義した module の変数を,あとから多次元配列に変更した場合,その変数を使っているすべてのコードを修正する必要があります.引数として受け取ると,関数内では多次元配列として扱うことができ,修正範囲が小さくなります.また,Fortran は配列の添字が1から始まりますが,計算によっては0から始まる方が良い場合があります (ex. 周期境界条件を用いたステンシル計算).この場合も,関数として受け取ればその関数内でだけ0始まりで扱うことができます.

ソースファイル外から呼び出されない場合は contains を使う

基本的には,1つの処理をする subroutine/function で1つのソースファイルを構成した方がよいと考えられます.そのため,関数は別のソースファイルから呼び出されます (外部関数とする).しかしながら,入力によっては関数の一部の計算が変わることもあると思います.その場合,外部関数の一部処理を切り出した関数 (内部関数) を作ることになります.この内部関数は基本的に同じソースファイルの,特定の外部関数から呼び出されるため,外部関数と同じように定義すると,どこかのタイミングで同じ名前の内部関数が定義されたとき,コンパイル時に "multiple definition" とエラーが出力されます.また,誤って別の外部関数の内部関数を呼び出してしまう可能性も十分考えられます.これらの問題は,Fortran コンパイラが内部関数と外部関数を理解していれば十分防げるバグです.
contains はこのために,特定の module/subroutine/function の内部関数を定義するためのキーワードです.

例えば,下記のプログラムはコンパイル・リンクすると multiple definition エラーが発生します.

! main.f90
program main
  implicit none

  call sub
  call test
end program

subroutine test
  print *, 'main :: test'
end subroutine
! sub.f90
subroutine sub
  call test
end subroutine

subroutine test
  print *, 'sub :: test'
end subroutine

これは,subroutine test が2つのソースファイルで外部関数として定義されているためです.そのため,これを contains を使ってそれぞれ,program main と subroutine sub の内部関数として定義します.

program main
  implicit none

  call sub
  call test

contains
  subroutine test
    print *, 'main :: test'
  end subroutine
end program
subroutine sub
  call test

contains
  subroutine test
    print *, 'sub :: test'
  end subroutine
end subroutine

これで出力は下記のようになるはずです.

$ ./a.out
 sub :: test
 main :: test

まとめ

個人的に Fortran を書くときに気をつけると良いことを書きました.中には仕様で正確に決まっていない項目や,もっと別の機能があるかもしれません.その際はご指摘いただけると幸いです.

*1:FORTRAN77 以前については,もはや現代においてその書き方などについて考える必要はないと思われます.もちろん過去の遺産で使われているものもあり,知識として持っておく必要はあるかもしれません.

*2:Fortran では動的配列はポインタではないので,C, C++ のようにポインタをスワップして解決というのができません.ただし Fortran でもポインタは使えます.使えますが,これはまた別の型として扱われ,同じ計算をしていても配列型とポインタ型で性能が大きく異なる可能性があります.(大体の場合ポインタ型で計算したほうが性能が悪い)

autotools の出力について

最近,autotools に苦戦しているのでメモ.

autotools を使った時に make のログが以下のように表示される事があるが,これは configure の --enable-silent-rules が有効になっているため.

  CXX hoge.o
  CXX fuga.o
  AR libx.a

それ以外の方法では,configure.in に AM_SILENT_RULES([yes]) をつけるか,make V=0 で可能.
configure によって出力される Makefile を読んでも分かるが,このログは単純に @echo を使った簡単な方法で,自分で作っている Makefile にも使える.もうちょっとちゃんとしてるが,書くとこんな感じ.

DEFAULT_VERBOSE=0

V_CXX=@(V_CXX_at_$(V))
V_CXX_at_=@(V_CXX_at_$(DEFAULT_VERBOSE))
V_CXX_at_0=@echo " CXX " $@;

.cpp.o :
	$(V_CXX) $(CXX) -c -o $@ $<

autotools 上で,カスタムのビルドルール (例えば CUDA カーネルのビルドなど) では,通常この silent rule が動いてくれない.
そのため,automake 側で AM_V_GEN というマクロが用意されている.以下のように使う.

.cxx.o:
	$(AM_V_GEN) $(CXX) -c -o $& $<

とやると,ログには以下のように出力される.

  GEN hoge.o
  ...

カスタムのビルドルールで, silent rule が無効の時は出して欲しいが,有効のときは出力自体して欲しくない場合, AM_V_at で良さそう.AM_V_at は単独の at symbol が silent rule 有効時に設定される.

NVCC=nvcc

.cxx.o :
	$(AM_V_GEN) $(NVCC) -cuda -c -o _$*.cpp $<
	$(AM_V_at) $(CXX) -c -o $@ _$*.cpp

とすると,有効無効でログの出力は以下のようになる.

# silent rule 有効時
  GEN hoge.o

# silent rule 無効時
  nvcc -cuda -c -o _hoge.cpp hoge.cu
  g++ -c -o hoge.o _hoge.cpp

これでスマートな感じの出力が得られる.

カーネルをなんとかして関数テンプレートにしたい

はてなブログに移行しました.そして最初の記事です.

最近 Xeon Phi なるものを触ったりしていますが,かなりとても書きづらいというか pragma ベースで書くのは (OpenMP で楽に並列化できるレベルなら良いですが) 決して楽ではないと僕は感じています.なぜならば pragma ベースで Xeon Phi のコードを実行 (Offload 実行モデル) だとメモリの取り扱いが若干謎で,ここで確保されてるから大丈夫だろうと思ったら SEGV になったりするし,メモリの確保とか解放,コピー操作といったメモリに関する管理をやろうとすると pragma ベースではかなりつらいしめんどくさい,というかそこら中に複雑な pragma が満ちます.それだったらいっその事 OpenCL で書いた方がメモリ管理はこっちでできるし良いと思ったので勉強しています.

幸いにして IntelXeon Phi を動かせる OpenCL SDK を用意していますし,Xeon Phi を使った OpenCL アプリケーションに関するガイドラインもあるので,問題は OpenCL を使った事がない,ということだけです.

ということで若干気に入らないスタイルもあるものの,自分で書くのは面倒なので C++ Wrapper を使いながら,コードを書いているんですが,OpenCL 1.2 以前ではオペレータオーバーロードも,テンプレートもサポートされていないとかいうので辛い*1 *2

私たちは C++ プログラマであるため,同じコードを何度も書くのは全くのナンセンスであり,コピペして型と関数名を変更するなど許されない.それに管理がめんどくさいし,型の数だけ増えればいずれは破綻する.なんとかしてテンプレートみたいな事はしたい.ならばソースファイルで型だけ置換すれば良いではないかという話になるが,それもそれでつらくて困る.

という訳で一時期話題になってその後は追っていない Boost.Compute *3 を調べてみたところ,やはりソースコードを置換しているようで,ただし OpenCL のカーネルビルドオプションで "-DTYPE=float" のように,コンパイル時に型を置換しているっぽい.ただし,これも問題で,TYPE を置換した場合に本来は置換されては困る部分もあるとかそういうこともあるが,カーネルのコンパイル時に分かる事ではあるのでまあ良い事にしましょう.コンパイル万歳.

なので,それを実践するとこうなる.

__kernel
void func(int n, __global TYPE * x) {
  const int tid = get_global_id(0);
  if(tid < n) {
    x[tid] = tid + 1;
  }
}
const std::string kernelSrc = /* std::string でカーネルコードを読み出す */;
cl::Program::Sources source(1, std::make_pair(kernelSrc.c_str(), kernelSrc.size()));
cl::Program program(context, source);
program.build(devices, "-DTYPE=float");
cl::Kernel kernel(program, "func", &err);

これで func は float 型の関数となる.オーバーロードをしたい場合は必要な型の数だけコンパイルして,何らかの方法で (Boost.Compute ではハッシュテーブルか何かで解決していたが) オーバーロード解決してカーネルを呼び出せるようなロジックを作れば良い.

コンパイル時に CPU か CUDA のどちらかをバックエンドにして計算しているコードはあるので,OpenCL のバックエンドを組み立てて上げればパフォーマンスはともかく,そのまま動くはずだ.

*1:OpenCL 2.0 は調べてないので不明.

*2:CUDA ではテクスチャメモリ周りで色々問題が起こったりもするけど使える.CUDA の最新バージョンでは試していない.

*3:OpenCL を用いた Thrust や Intel TBB などと同じ STL-Like なライブラリである.

ログインシェルが bash/sh 以外で MPI のプログラムを動かすときに環境変数が渡らないときは

最近,研究室の環境も zsh で動かしています.

しかし,MPI のプログラムを走らせるときにどうしても環境変数に頼ってしまう部分があって,そのままプログラムを実行すると環境変数が他のノードに渡らなくてプログラムがクラッシュします.(もしかしたら実装の問題かも? MVAPICH2 を動かしていて起きましたが OpenMPI でも同じような問題で ML が流れていた気が)

環境変数に頼らないといけない部分というのは,例えばライブラリの実行パスとかですね.コンパイラ,Boost, MPI とか複数のバージョンの組み合わせを行う必要があるときにどうしても必要でして.僕の場合は普段は変更しないパラメータを環境変数にしてたりして.

mpirun_rsh に環境変数を渡す,という方法はありますがちょっと多すぎて設定が辛い.

色々調べたら MVAPICH2 の場合,mpirun_rsh でラップしたシェルスクリプトを実行すれば良いらしい.

$ cat run.zsh
#! /bin/zsh

# 実行時に必要な環境変数等
LD_LIBRARY_PATH=...:$LD_LIBRARY_PATH
OMP_NUM_THREADS=16

# 第1引数に実行ファイルが指定されているとして
./$@

こんなシェルスクリプト書いて,実行は

$ mpirun_rsh -hostfile ${HOSTFILE} -np ${NHOSTS} run.zsh a.out a b c

とすると run.zsh環境変数が設定されて a, b, c が a.out に引数として渡される.

よく考えると MPI は bash/sh ぐらいでしか実行してない雰囲気は感じる.検索してもあまりそういう情報がなかったのもそういうことな気がする.

MPI ことはじめ資料 (研究室向け)

この間の高速化,並列化みたいなスライドの MPI 版です.OpenMPI を使う必要があるとのことで,OpenMPI を対象として (といってもビルドと実行方法が OpenMPI 用になっているだけですが) 作りました.

C++ と書いておきながらほぼ C です.

# MPI は素で書くのが辛すぎる...

MPI ことはじめ