2011-07-30

FortranプログラムをSConsでビルドする

SCons Fortran

前置き

とある方からFortran ソースコードを頂いて、ご丁寧にもMakefileまで付けてくれていたのですが、
まあMakefileなんてもう古いですよね。

そこでmakeからSConsにビルドツールを変更しました。
SCons: A software construction tool - SCons

ビルドツールといえば、OMakeとかCMakeとかもあるのですが
OMakeだとFortranの対応状況が不明。CMakeは使ったこと無いしよく知りませんけどちょっと大げさな印象。configureとか使ってたわけではないし、単純にmakeの代わりが欲しいだけだったのでCMakeは見送り。

SConsはFortran用のモジュールがあるそうですし、Pythonで設定ファイルを書けます。
大規模開発になるとビルドに時間がかかるようになるそうですが、そんな大規模なもの作らないので問題なし。所詮数値計算で使うファイル数なんて大したことないでしょう(少なくとも私の使う範囲では)。

使い方

SConsではSConstructという設定ファイルにビルドルールを書いて行きます。makeでいうMakefileです。
私が使ってる状況は非常に簡単で、ディレクトリ内にある*.f90ソースをすべて使ってコンパイルします。
あと、blasとlapackライブラリをリンクするだけです。

SConstructファイルはこんな感じです。

#SConstruct
import os
env = Environment(ENV=os.environ)
env.Program(target='calculation',source=Glob('*.f90'),LIBS=['blas','lapack'])

calculationが実行ファイル名です(安直)。sourceのところにソースコードのファイル名を指定します。私の場合ディレクトリ内のすべてのf90ファイルだったので、非常に簡単。LIBSでリンクするライブラリを指定。
これで、それぞれの依存関係をすべてチェックしてくれるようになります。
簡単ですね。

コンパイルは

$ scons -Q

とかやれば、勝手にやってくれます。私はvimから:makeと打つのが癖になってるので
sconsを呼び出すMakefileを書きました(無駄..)
いわゆるmake cleanとかは

$ scons -c

でやってくれるそうです。SConsすばらしい。

2011-07-29

Fortranのmaxloc関数でハマった。

Fortran

Fortranのmaxloc関数でハマった。maxlocとは

配列中のすべての要素，配列中の要素の組，または配列の指定された次元中の要素の最大値の位置を返します。
http://www.xlsoft.com/jp/products/intel/cvf/docs/vf-html/az/az08_09.htm

という関数らしいので、便利そうだと思っていたのだが。。。
位置は返すけど、indexは返さないんですね

program mx
  integer,dimension(8   ) :: var1 = (/1,2,3,4,2,8,3,4 /) !  indexは1〜8
  integer,dimension(0:7) :: var2                                       !  indexは0〜7
  var2 = var1
  print *,var1
  print *,var2
  print *,repeat("-",256)
  print *,maxloc(var1,1),var1(maxloc(var1,1))
  print *,maxloc(var2,1),var2(maxloc(var2,1))
end program mx

結果：

           1           2           3           4           2           8           3           4
           1           2           3           4           2           8           3           4
 ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
           6           8
           6           3

添字が1から始まるとは限らないんだから、これはまずいでしょ。
位置という表現は不正確。せめて、注意書きが欲しい。

過去のコードとの互換性もわかりますけど、deprecatedにして新しい関数作ってくださいよ。

2011-07-11

勾配・発散・回転のメモ

Mathematics

ただし直交座標系。写し間違いがあっても知りませんので注意。

grad
- $\Large \nabla f = \sum_i \frac{1}{h_i} \frac{\partial f}{\partial u_i} \mathbf{e_i}$
div
- $\Large \nabla \cdot \mathbf{A} = \frac{1}{h_1h_2h_3}\left[ \frac{\partial}{\partial u_1}(h_2h_3 A_1)+\frac{\partial}{\partial u_2}(h_3h_1 A_2)+\frac{\partial}{\partial u_3}(h_1h_2 A_3) \right]$
rot
- $\Large(\nabla\times\mathbf{A})_{u_1}=\frac{1}{h_2h_3}\left[\frac{\partial}{\partial u_2}(h_3A_3)-\frac{\partial}{\partial u_3} (h_2A_2)\right]$
- $\Large(\nabla\times\mathbf{A})_{u_2}=\frac{1}{h_3h_1}\left[\frac{\partial}{\partial u_3}(h_1A_1)-\frac{\partial}{\partial u_1} (h_3A_3)\right]$
- $\Large(\nabla\times\mathbf{A})_{u_3}=\frac{1}{h_1h_2}\left[\frac{\partial}{\partial u_1}(h_2A_2)-\frac{\partial}{\partial u_2} (h_1A_1)\right]$
Jacobian
- $\Large J=h_1h_2h_3$
Laplacian
- $\Large \nabla^2 f = \frac{1}{h_1h_2h_3}\left[\frac{\partial}{\partial u_1}\left(\frac{h_2h_3}{h_1}\frac{\partial f}{\partial u_1}\right) + \frac{\partial}{\partial u_2}\left(\frac{h_3h_1}{h_2}\frac{\partial f}{\partial u_2}\right) + \frac{\partial}{\partial u_3}\left(\frac{h_1h_2}{h_3}\frac{\partial f}{\partial u_3}\right)\right]$
(h1,h2,h3)
- デカルト: (1,1,1)
- 円筒: (1,r,1)
- 極: (1,r,rsinθ)

2011-07-11

流れ関数による速度の表式

Physics

流れ関数(stream function)を $\psi$ とすると、速度 $\mathbf{V}=(V_x, \hspace{5mm} V_y)$ は

$\Large V_x=\frac{\partial \psi}{\partial y}, \hspace{5mm} V_y=-\frac{\partial \psi}{\partial x}$

と表されるんだけれども、うちらの分野では

$\Large \mathbf{V}=\nabla \times \left( \begin{array}{ccc}0 \\ 0 \\ \psi\end{array} \right)$

とか

$\Large \mathbf{V}=\nabla\psi\times\mathbf{e}_z$

の方が覚えておくと便利っぽい。
ちなみに流れ関数でやると、非圧縮性流体での連続の式

$\Large \nabla \cdot \mathbf{V}=0$

が自動的に満たされるのです。

2011-07-04

札幌C++勉強会 #2

Sapporo.cpp

7/3 札幌市産業振興センターで２回目となる札幌C++勉強会が開催されました。

詳しくはこちらからhttps://sites.google.com/site/sapporocpp/benkyou-kai

私はLTでBoost.Phoenixなんぞを話しましたけど、
まあ５分で話せるライブラリじゃないですね。皆さん結局よくわからなかったんじゃないでしょか。

他の発表者の皆さんが軽めの話題でC++に興味を持ってもらおうと頑張っていらっしゃったのに、
一人だけ空気読まず、Boostとか出しちゃいましたよ。
LTだから！、LTだからいいですよね！？(ごめんなさい...)

資料については上のリンクの勉強会サイトに置いておきます。
LTで話せるような簡単なものなので、ほとんど参考になりませんけどねー。
30分もあれば、proto的なことも含めて話せるかな？でも多分盛り上がらないでしょうねー。

次はAndroidとかやりたいですね。もちろんC++で。

2011-06-15

OSC2011 Hokkaido参加レポート

Sapporo.cpp

OpenSourceConference2011 Hokkaidoに札幌C++勉強会として参加してきました。
ブースへの出展と、セミナー発表(講師は@hotwatermorningさんにお願い)です。

ブース出展では、無料とは思えぬ見事な小冊子を配布しておりました。
12:30頃には全部配布し終わってしまいました。
お客さんの「これ無料なんですか？」の声が嬉しかったですね。
ブースに来ていただいた皆様、どうもありがとうございます。

以下、裏話。

あの小冊子、最初の計画ではもっと簡単なもののはずでした。
それぞれメンバーがコラム作って来て、ホッチキスでとめる、みたいな。
それが計画途中で、「これはマジメに行くしかない」(タイトルがアレなくせに)みたいなことになり、
製本することに。結果的に良かったんじゃないかなーって感じですが。
出費はかかりましたけど。

セミナー発表は私、準備段階で死んでましたね。
実は最初、Google docs上でメンバーが共同で作ってたものがあったのですが、
割とネタ成分は少なめでした。
ネタ的文句をいれて、C++0xの説明に入る、みたいなものが続いていたわけですけど、
量的、完成度的に不安があるものだったので、保険として私が別スライドを急ぎで作ってました。
（最終的にはそっちの保険の方をベースにしたわけですけど）

OSC前日、会場の下見に行った後、プレゼンの作成に移ったわけですが、
作成始める前に、保険の方でプレゼン実演して、
そしたら意外と、メンバーからの反応がいい。
保険で作った方をベースにして作っていこうと。このときすでに20時すぎ？だったかな。
この時の出来具合は、全体の５割程度しか作ってなくて、ああこれは今日中には終わらないなと。
（実は実演したとき、50ページ越えのスライドだったんですが18分で終わっちゃってて、40分の発表だから、これ最終的に100ページ超えんのかなーとか、時間的に結構焦ったんですけどねww)
しかも作っててわかったのですが、後半になるほどネタが入らない。
「regex」と「こんなの絶対おかしいよ！」をどう混ぜろとwww
おかげで、後半まで作り込む時間はありませんでした。

なんとか、プレゼンを形にして、hotさんと打ち合わせしてたんですが、
最初のベースが、私くらいプレゼン慣れしている人を想定した作りになっていて、かなり難しい構成になっていたんですよね。
最終的に、90枚越えでそれを40分でってことなので、長めのプレゼンってこともあり時間配分が難しいのです。
(ちなみに私の最長プレゼンは120分ですねー、ってむしろそれ以上は休憩挟まないと駄目だろ...)
打ち合わせの段階で、マジメに説明してたりすると、時間足りなくなるし、時間かけすぎるとリズム悪くなってダメなことがわかってたので、
「とにかく、時間余ってもいいから次々進めてください！」的な打ち合わせしてました。

まあ、発表自体はなんとか成功したらしく、「hotさんすばらしい」ってところです。

以上、らくがき的レポートでした。

2011-05-01

Fortranの開発環境をWindowsで構築したときのメモ

Fortran

コンパイラ(gfortran)のインストール
参考サイト: MinGW http://www.mingw.org/

http://sourceforge.net/projects/mingw/files/Automated%20MinGW%20Installer/mingw-get-inst/から
mingw-get-inst-20110316をダウンロードし実行
Select Componebtsまでデフォルト設定でOK

Select ComponebtsでFortran Compilerにチェックを入れる
その後インストール作業を続けると、自動でパッケージをダウンロードして
インストールされる。

これでgfortranはOK

Code::Blocks for Fortranのダウンロード
参考：Code::Blocks for Fortran http://darmar.vgtu.lt

http://darmar.vgtu.lt/index.php/cbdownloadから
For Windows 32bit with all plugins (for C/C++ and Fortran)
CodeBlocks_Fortran_v0.6.1-1_Win32_all.zipをダウンロード

codeblocks.exeを実行すると、コンパイラを自動で見つけてくれる
おそらくGNU GCC CompierがDetectedになっているはず。
そのままOKで続行

[確認作業(不要かも)]
SettingsメニューのCompiler and debuggerを選び、
Selected CompilerでGNU Fortran Compilerを選ぶ。
そして、その下のToolcain executablesで
Compiler's installation directoryにC:\MinGWが選ばれていることを確認
その他Program Filesについてもわかる人は確認

使ってみる
テストプロジェクトの作成
File -> New -> Project...と開き、
Fortran applicationを選択してGoで次のダイアログへ

なんかTipsっぽいものが表示されるがそれは無視して、
Project/Directoryの名前を決める。各自適当に
例えば
Project title : test
Folder to create project in: 右のファイル選択ツールをクリックしてマイドキュメントを選択

次のダイアログで、Compilerを
GNU Fortran Compilerに指定しFinish

するとプロジェクトが出来上がる。
左(デフォではね)のManagementパネルの
ProjectタブにFortran Ssourcesというのがあるから、
その中のmain.f90をダブルクリックすると、
Hello Worldを表示させるソースが現れる。

これを編集すればFortranプログラミングができる。

コンパイル＆実行
上部Buildメニューの"Build&run"をクリックすると
コンパイルされ、プログラムが実行される。
コマンドプロンプトが開き、Hello Worldが確認できるはず。

追記：
Intel Fortranを使うときは、iccvars.batだっけ？そんな感じのの環境変数を設定するファイルが付いてくるはずなので、iccvars.bat->codeblocksの順に起動するようにbatファイルを書く。詳細は公式サイトで。