Fortranから始めるJulia
はじめに
このウェブサイトは、Fortranで数値計算を行なったことがある人が、Juliaについて勉強するためのサイトです。そのため、Fortranについてはある程度の知識がありコードを書いたことがある人を対象として、Juliaについて解説したいと思います。
目次
- Fortranから始めるJulia
- はじめに
- Juliaはどんなプログラミング言語か
- Fortran使いがJuliaを学ぶメリット
- 1. 実行前にコンパイルが必要ない
- 2. 他の人が作ったライブラリが豊富にあり、かつ導入が簡単
- 3. 最新のモダンなコーディング技術を使えるため、書きやすくメンテナンスのしやすいコードが書ける
- 4. 古（いにしえ）のFORTRANコードから決別できる（かもしれない）
- 5. Fortranと文法が似ている
- 動作確認に使用したJuliaのバージョン
- このウェブサイトの作成者
- 有用な本
- FortranとJuliaの文法比較
- 線形代数計算
- FortranのコードをJuliaへ移植してみる
- マンデルブロ集合
- 万有引力に従う天体の運動
- 基本方程式
- Leap-frog法による二体問題
- Velocity Verlet法によるN体問題
- Velocity Verlet法によるN体問題: Juliaっぽい書き方版
- Velocity Verlet法によるN体問題: さらに早く
- 2次元Ising模型
- FortranコードをJuliaに移植するためのチェックシート
- FortranでできることをJuliaでやってみる
Juliaはどんなプログラミング言語か
Juliaという言語がどのような考えで開発されたか、についてはなぜ僕らはJuliaを作ったかを見てみてください。非常に新しい言語で、数値計算がしやすい言語として開発されています。
Fortran使いがJuliaを学ぶメリット
Fortran使いがJuliaを学ぶメリットについて述べます。 まず、大前提として、本当に本気の大規模な並列計算を伴う数値計算を行う場合には、現状ではFortranを使った方が高速だと思います。多くのスーパーコンピュータではそのスパコンに最適化されたコンパイラが提供されており、それらのコンパイラを用いて適切にチューニングされたFortranのコードは最高速度を達成すると思います。そして、計算期間が数ヶ月以上に渡るような計算の場合には、数パーセントの高速化が非常に重要になってきますが、そのような場合にはJuliaよりFortranの方が向いていると思います。 また、ちゃんと動いている既存の大きめなFortranコードをJuliaで書き換える必要はないと思います。
それでは、Juliaにはどのようなメリットがあるでしょうか。以下の点がFortran使いにとって重要だと思います。
- 実行前にコンパイルが必要ない。コンパイルに四苦八苦することがなくなる
- 他の人が作ったライブラリが豊富にあり、かつ導入が簡単
- 最新のモダンなコーディング技術を使えるため、書きやすくメンテナンスのしやすいコードが書ける
- 古（いにしえ）のFORTRANコードから決別できる（かもしれない）
- Fortranと文法が似ている：配列の添字の始まりが1から(Pythonは0から)、配列のメモリ格納順が列ベース、等、Fortranと似ている
それぞれの項目について説明していきます。
1. 実行前にコンパイルが必要ない
FortranではコードをIntel Fortranコンパイラやgfortranでコンパイルし、実行ファイルを作って実行します。もしLAPACKやBLASなどのライブラリを使いたい場合には、コンパイル時にリンクを行います。複数のファイルからなるソースコードであれば、makefileを用いてコンパイルすることもあるでしょう。 Juliaではこのようなコンパイル作業は必要ありません。Juliaではコードを実行する直前に型推論と最適化を行ない、その後にコードを実行します。ですので、
julia hoge.jlのような形でソースコードをそのまま実行することができます。 自分で書いたコードであればコンパイルはそれほど問題になりませんが、例えば他の人が書いたコードを改造して使いたい場合、コンパイル方法をどのようにすればよいかは計算機環境によって異なるために、コンパイル自体が難しい場合があります。Juliaではコンパイルという作業がありませんので、他の人が書いたコードをそのまますぐに実行することができます。
2. 他の人が作ったライブラリが豊富にあり、かつ導入が簡単
Fortranでは他の人が作ったライブラリを自分のコードに使う場合は、「ソースコードをそのまま使う」か「リンクが可能なライブラリとなっているものをコンパイル時にリンクする」という方法があります。どちらの方法であっても、複数のライブラリを自分のコードにリンクしてコンパイルして実行ファイルを作る、というのはなかなか難しいです。さらに、他の人が作ったライブラリがそもそも自分の計算機環境で動くのか、コンパイルが可能なのか、ということについても気を使う必要があります。
そして、Fortranで書かれたコードというのはウェブ上ではあまり多くありません。よくあるのは、研究室で代々使っているコードや、論文の著者のウェブサイトからダウンロードしたコードや、直接やりとりをしてもらったコードだと思います。あるいは、有償あるいはオープンソースで書かれた20年近く使われ続けている巨大プロジェクトのコードだと思います。Fortranは研究室単位やプロジェクト単位で使われていることが多いため、「自由に使えるパッケージ」という概念がほとんどありません（最近はFortranでもパッケージ管理システムを作ろうという試みがこちらにあるようです）。そのため、必要な機能は自分で実装するか詳しい人にコードのありかを聞くことになります。例えば、非線形な関数$f(x)$の最小値となる$x$を求める場合、どうしますか？ 詳しい方であれば、Intel Math Kernel Library (MKL)に非線形最小二乗問題を解くためのルーチンが含まれていることを知っていると思います。しかし、その存在をどうやって知ればよいでしょうか？ あるいは、2変数関数$f(x,y)$のあるデータ点での値を複数持っているときにスプライン補間でデータを補間したい場合、どうすればいいでしょうか？ Fortranにはこの目的のためのDIERCKXというコードがあるという情報に辿り着けるでしょうか？ そして、どのように使うかわかるでしょうか？
Juliaであれば、非線形関数の最小化のパッケージとしてOptimというパッケージがあり、スプライン補間用パッケージにDierckxがあり（実はDIERCKXのラッパーパッケージ）、どちらも簡単に使うことができます。
3. 最新のモダンなコーディング技術を使えるため、書きやすくメンテナンスのしやすいコードが書ける
Fortranは歴史が長く、古いコードであってもコンパイルが可能で実行が可能であるということが大きな利点にひとつになっていると思います。科学技術分野における数値計算であれば物理現象自体が変わることはそうそうありませんから、非常に良いパフォーマンスを発揮する古いコードも使う価値があります。
一方で、それらのコードを改造して新しい問題に対応させようとすると、よほどのFORTRAN/Fortran経験者でなければ困難になります。そして、FORTRANをいじれる人間はFORTRANをいじったことがある人間だけですから、現代でわざわざFORTRANをいじる新規参入者は少なく、FORTRAN人口は減っていってしまうでしょう。それでも巨大なプロジェクトのコードであれば開発者を確保することができるかも知れませんが、研究室単位、あるいは数人のグループによって開発されたコードの場合、読むのが困難でしょう。
例えば、FORTRAN77時代に変数を複数の箇所で使う場合にはcommon文というものが使われていました。Fortran 90以降においてはcommon文を使うのではなくmoduleを使うことが推奨されるようになりました。しかし、module内に定義された変数の値がどのタイミングでどの関数によって変更されるかは注意深くコードを読んでいかないとわかりません。Fortranしか書いたことがない方にとってはこれは当然かもしれませんが、現代ではこのような「グローバル変数」を使うコーディングは読みにくいコードとして知られています。FortranでもFortran 2008以降のオブジェクト指向のコーディングスタイルを用いればモダンな書き方が可能ですが、現状多くのFortran使用者はオブジェクト指向Fortranを使っていません（英語及び日本語でweb検索した際の情報の少なさがそれを示唆しています。個人的にはオブジェクト指向Fortranは結構好きでよく使っており記事も書いていますが、使用者は実感として少ないです）。
ソースコードは書かれた時間と空間から離れれば離れるほど読みにくくなります。書いた本人が書いた直後であればコードの意図を取ることは簡単で、どんなに読みにくくてもコメントがなくても問題はありません。複数人で開発している場合は書いた人に聞くことでまだコードを理解できます。一方、書いてから10年経ってしまったコードは書いた本人にも意図がわからなくなってしまう場合があります。これを避けるためには、なるべく、書きやすく、読みやすく、メンテナンスしやすい、コードを書く必要があります。Fortranは少数での開発が多かったためかあまりこの視点に立ったコーディングスタイルが広まってきませんでした（気象庁のFortran 標準コーディングルールは有名ですが、今となっては古いと言えるでしょう。現代であればFortranをコーディングする際に気をつけていることが参考になります）。一方、プログラミング業界自体は日進月歩で進歩しているために、この視点に立った多くのノウハウが蓄積されています。一番有名なのはオブジェクト指向における「デザインパターン」でしょうか。 Juliaは2018年にバージョン1になった極めて新しい言語でして、現代のコーディングの良いところをなるべく取り込むようにデザインされています。ですので、Juliaでコードを書くとモダンな読みやすいメンテナンスしやすいコードを書くことができます。
4. 古（いにしえ）のFORTRANコードから決別できる（かもしれない）
2010年代には確実に使われていた古いFORTRANコードを紹介します。そのコードはあるレビュー論文のコード(動的平均場理論の数値計算用のコード)の一部分でして、非線形関数の最小値を求めるコードです。そのコードのコメントに
C========+=========+=========+=========+=========+=========+=========+=$
C PROGRAM: minimize
C TYPE : subroutine
C PURPOSE: conjugent gradient search
C I/O :
C VERSION: 30-Sep-95
C COMMENT: This is a most reliable conjugent gradient routine! It has
C served us well for many years, and is capable to cope with
C a very large number of variables. Unfortunately, we don't
C know who wrote this routine (original name: 'va10a'), and
C we find it very obscure.
C Don't worry, it works just fine.
Cnoprint================================================================とあります。1995年の時点で誰が書いたかわからないのですが、よく動く、という理由で2010年代にも使われていました（今でも使われているかもしれません）。このコードはdoループをdo,enddoで閉じるのではく、行番号とcontinue文で制御しており、ループから抜け出したりgotoで移動したりが多発するコードです。確かによく動き、ちゃんとした解を返します。私は一度このコードをFortran90に書き換えようと試みたことがありますが、私には無理でした。
Juliaの場合、このコードを使わなくてもOptimというパッケージを使えば非線形関数の最小値を求めることができます。
追記
このコードについて調べてみたところ、ついに元のコードの情報が書かれている文献（R.Fletcherの"FORTRAN SUBROUTINES FOR MINIMIZATION BY QUASI-NEWTON METHODS", Research group report(United Kingdom Atomic Energy Authority)）を見つけました。これによると、このコードva10aは1972年の4月に書かれたそうです！
5. Fortranと文法が似ている
配列の添字が1始まりであるため、Fortranのコードを素朴に移植することが容易です。PythonやCは0始まりですから、それらよりも違和感は少ないと思います。 FortranとJuliaがどのくらい文法が似ているかは後述します。
動作確認に使用したJuliaのバージョン
Julia 1.6.1
このウェブサイトの作成者
Yuki Nagai, Ph.D
有用な本
以下の本はJulia言語を理解する上で役に立つと思います。
- 進藤 裕之, 佐藤 建太, "1から始めるJuliaプログラミング", コロナ社(2020)