[Fortran][ifort] 引数で出力ファイルの名前を指定する

ソース･コードのコンパイル時に出力ファイルの名前を引数として与える方法についてのメモ。

1. (出力ファイル名を1番目の引数として与えると仮定して、)ソース･コードを用意：
program main
implicit none
character(100) filename

call getarg(1,filename)
open (12, file=trim(filename), status='replace')
write (12,*) "Hello!"
close (12)

end program main
• "getarg(1,filename)"はfilenameに引数の1番目を書き込む組み込み型サブルーチン。
2. 端末で次のようにコンパイル・実行する：
   $ ifort -o main.out main.f90
   $ ./main.out ./filename_output
   

このとき、実行ファイルと同じディレクトリにfilename_outputの名前で出力ファイルが作成される。

[Intel MKL] Windowsコマンドプロンプト上でIntel MKLの乱数発生を含むコードをコンパイルする方法

コンパイルしようとしてるコードがIntel MKLの乱数発生を利用している場合のWindowsコマンドプロンプトでのコンパイルの方法についてのメモ（長い！）ただし、Intel Visual Fortranをインストールしていることが前提である。コンパイル・オプションの指定方法はアーキテクチャに依存する。以下ではインテル64アーキテクチャを前提として書くが、その他のアーキテクチャについてはNoteを参照。

手順は以下の通り：

1. ショートカット・アイコンを利用した、自動の環境設定を行うと同時にコマンド・プロンプトを立ち上げる。
2. コマンド・ラインに以下を入力し、ENTERキーを押す：
   "C:\...\mkl\tools\environment\mklvarsem64t.bat"
   
3. ソース・コードfilename.f90を次のようなオプションを付けてコンパイルする：
   ifort filename.f90 mkl_intel_lp64.lib/mkl_intel_ilp64.lib mkl_sequential.lib mkl_core.lib
   
   ※static link + sequentialの例。スレッド化レイヤーに逐次(mkl_sequential.lib)を選択する場合、ランタイム・ライブラリーはしていしなくてよい。

Note: サポートされているアーキテクチャー
IA-32/IA-32 互換 mklのディレクトリ\ia32
インテル64/インテル64 互換 mklのディレクトリ\em64t
IA-64 mklのディレクトリ\ia64

【Intel MKL】 Intel MKL Nonlinear Least Squaresを使う方法

Intel Math Kernel Library (MKL)の中の非線形最小二乗法ルーチンの使い方についてのメモ。ちなみに、自分のマシン環境を前提に書く。具体的には、言語はFortran90、アーキテクチャはIntel(R) 64である。リファレンスはIntel(R) Math Kernel Library for Linux* OS User's Guide January 2010, Doc. N.314774-010US(UG2010)。

• 言語やアーキテクチャによってパスが異なるので注意。

手順は次の通り：

1. MKLのルーチンを呼び出しているソース・コードの頭に、
   program name
   implicit none
   include 'mkl_rci.fi'
   のように、include 'mkl_rci.fi'を付け加える。


2. 次のようにコンパイルする。
   
   (dynamic & threadingの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread
   
   (dynamic & sequentialの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread
   
   (static & threadingの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -Wl,--start-group -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -Wl,--end-group -liomp5 -lpthread
   
   (static & sequentialの場合)
   ifort filename_1.f90 filename_2.f90 ... filename_x.f90 -Wl,--start-group -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread
   • ポイントは、staticの場合、cluster components(この場合使用していない)、interface、threading、そしてcomputational librariesの前後を-Wl,--start-groupと-Wl,--end-groupで囲うこと。
   • sequentialの場合、Compiler Support Run-time librariesは関係ないので、オプション-liomp5(compatibility OpenMP run-time libraryで、US2010が推奨)を外している。

Note.1: Intel MKLの構造
Intel MKLは次の4つの階層を持つレイヤー・モデル：

1. Interface layer
2. Threading layer --- threadingがデフォルト
3. Computational layer
4. Compiler Support Run-time libraries

であり、それぞれのレイヤー毎にユーザーが必要なライブラリの選択を行う。なお、linking modelがstaticかdynamicかで選択するファイル名が異なる(アバウトに言えば、staticなら拡張子が.aで、dynamicなら拡張子が.so)。各レイヤーの内容（選択肢）についてはUG2010 Table 3-6、3-7を参照。

注意点は以下の通り（[x]はレイヤーxに関するもの）：

• [1] 2^31-1(約21億)以上の要素を持つ配列を使用する場合、デフォルト・インターフェイスのLP64ではなく、ILP64を指定する必要がある。これは配列のインデックスを与える上での問題。(UG2010, p.3-6)


• [2] MPIを利用する等、特別な理由がない限りthreadingを選択するべき。(UG2010,p.3-5)


• [2] sequentialを選択する場合、Compiler Support Run-time librariesは関係ない。


• [2] sequentialを選択する場合、*sequential.*libraryを選択する。また、個の場合、link lineにPOSIX threads library(pthread)を追加する。


• [2]Optimization Solver routinesのfunction domainはMKLがデフォで使用するスレッディング機能は使われない(UG2010,p.6-1,Using the Intel(R) MKL Parallelism)。そのため、OpenMPのスレッディングの数を指定に関するオプション設定を気にしなくてもOK。


• [3] 非線形最小二乗法はMKL function domainではOptimization (Trust-Region) Solver routinesに対応。この場合、includeファイルmkl_rci.fiを使う。

また、UG2010は、次のように、dynamic linkを推奨している。

You are strongly encouraged to dynamically link in the compatibility OpenMP* run-time library libiomp or legacy OpenMP* run-time library libguide. Link with libiomp and libguide dynamically even if other libraries are linked statically.

Linking to static OpenMP* run-time library is not recommended because it is very easy with complex software to link in more than one copy of the library. This causes performance problems (too many threads) and may cause correctness problems if more than one copy is initialized. (UG2010,p.5-6)

Note.2: 環境設定
.profileに以下の行を追加。

# setting up MKL environment for bash
. absolute_path_to_installed_MKL/tools/environment/mklvarsem64t.sh

.の後は半角スペースが必要。

Note.3: コンパイル&リンク

ifort filename.f90 -L(MKL path) -I(MKL include)
[-I(MKL include)/{32|em64t|{ilp64|lp64}|64/{ilp64|lp64}}]
[-lmkl_blas{95|95_ilp64|95_lp64}]
[-lmkl_lapack{95|95_ilp64|95_lp64}]
[cluster components]
-lmkl_{intel|intel_ilp64|intel_lp64|intel_sp2dp|gf|gf_ilp64|gf_lp64}
-lmkl_{intel_thread|gnu_thread|pgi_thread|sequential}
[-lmkl_lapack] -lmkl_core
{-liomp5|-lguide} [-lpthread] [-lm]

なお、()はパスを表す。ただし、Note.2で環境設定している場合、-L(MKL path)と-I(MKL include)を書く必要はない。

【Fortran】 Intel MKL Nonlinear Least Squares Problem with Bounded Linear Constraint

　非線形連立方程式を解くためにアルゴリズムを探している。

　いくつかあるだろうけど、ここではIntel Compilerに同梱されているMath Kernel Library (MKL)を使う方法を考える。

　MKLの中では、線形制約下の最小二乗法のルーチンを使うことになる。このルーチンはTrust Region Methodを使っている。Trust Region Methodについては、こちらのレジュメがわかりやすい。

（追記）
　この非線形最小二乗のアルゴリズムは使い物にならないかもしれない。というのも、内部でJacobianを求めるときに引数として用意すべき外部サブルーチンの形式への規制が強すぎるから。具体的には、外部サブルーチンが sub(M,N,X,F) （Mは連立方程式の数、Nは変数の数、Xは変数、Fは関数の値）の形でなければならず、連立方程式に登場するパラメータを引数として定義することができなくなっている。もちろん、サブルーチンの中でテキストファイルから読み込むというアプローチも考えられるけど、I/O時間が余計に掛かってしまって計算が遅くなると思われる。

（追記２）
　上記の（追記）は間違い。グローバル変数を使えば大丈夫。具体的には、moduleファイルを使う。

memo

DQED - Bounded, Constrained Least Squares/Nonlinear Equations

MINPACK Least Squares Minimization

多分これが原因

スプライン補間の近似すげーと思った。
元の離散的なデータから滑らかな曲線で近似できる。

でも、この「滑らか」という部分で勘違いをしてしまっていたようだ。
具体的には、移動平均みたいな、滑らかな「トレンド」を抽出してくれると勘違いしてしまっていた。なんともお恥ずかしい。でもまあ、スプライン近似の定義からありえないのだが（笑）

以下の図は元のデータがギザギザな部分を含む場合。

グリッドの数は400とまあまあ細かい。
これに対してspline.f90の中の3次スプライン補間法を適用。
するとどうだろう。まあ当たり前だけど、ギザギザの部分も見事に再現してくれた。
でも、自分としてはこのギザギザを削ぎ落とした滑らかな曲線がほしいわけだ。まあ考えられる方法はあるのでいいのだ。とりあえず問題発見はできた。これがプログラムの結果が上手くいかない１つの原因と思われる。

あと、もう１つは、外挿の精度の低さ。これも図の右端でわかる。これもプログラムが上手く結果を出さない１つの原因だと思われる。

さて、これからシャワー。今日はテニスだ！

【Fortran90】 3次スプライン補間法

関数を何らかの方法で補間する必要があるので、その方法として3次スプライン補間法を考える。よく使われているらしいから。

でも、これがよく使われる理由とは？データがたくさんあるときでも、ラグランジュ補間のように当てはまりの精度が問題にならないこと、なのかな。なぜ3次かはわからない。

B-スプラインというのも耳にしたことあるけど、これはデータの点を通過しないという点で、個人的には、問題がある。

スプラインのサブルーチンはFortran Source CodesのSPLINEがある。

3次スプライン補間法（Piecewise Cubic Spline）を使うには、次のサブルーチンを使う。

• spline_cubic_val(n,t,y,ypp,tval,yval,ypval,yppval)


• spline_cubic_set(n,t,y,ibcbeg,ybcbeg,ibcend,ybcend,ypp)
• 2階微分の情報を与えるサブルーチン。
• スプライン補間の境界条件はここで設定する。自然境界(natural boundary。2つの境界点における近似関数の2階微分が0となる)を使う場合、ibcbeg=2, ybcbeg=0.0D0, ibcend=2, ybcend=0.0D0とする必要がある。

Codingの手順（自然境界条件の場合）

1. n, t, yにそれぞれデータ・ポイントの数（整数）、argumentの配列（倍精度、n個の要素）、関数値の配列（倍精度、n個の要素）を設定する。

2. ``spline_cubic_set''サブルーチンをコールして2階微分の配列yppを得る。

3. ``spline_cubic_val''サブルーチンをコールして、argumentがtval（倍精度、input）の時の関数の近似値yval（倍精度、output）、1階微分ypval（倍精度、output）、2階微分yppval（倍精度、output）を得る。


Note 1: 微分可能でない関数の近似精度は低い。特に、微分可能でない点の近傍においてそれは顕著。

【ifort】バグ

コーディングに問題があるにも関わらず、コンパイル・エラーが出ないバグについて、発見したものを記録しておく。[ ]の中はifortのヴァージョン。

• [Ver.11.1_072] メインのプログラム・ファイルとモジュール・ファイルで同じ名前の変数を宣言してもエラーが出ずにコンパイルされる。ただし、値は"NaN"となる。

【Fortran】チューニング

並列計算をする前に考えるべきことは、チューニング。
簡単に思いついたものを箇条書きでメモしておく。

• if文はできる限り使わない。
• （例）配列を上手く使う。


• 割り算や累乗をできるだけ避ける。
• （例）ループ内で無駄な処理が重複しないよう、できる限りループの外へ出す。


• メモリへの連続アクセスをできるように配列の添字の使用を行う。
• ただし、キャッシュ・ヒット率とのバランスを考慮する。


• キャッシュ・ヒット率を高める。
• （例）ループ・アンローリング。
• ただし、メモリ連続アクセスとのバランスを考慮する。


• できる限り組み込み関数を使わない。


• I/O処理を少なくするため、インライン展開を行う。
• インライン展開しない場合、コンパイラの最適化が上手く行わず、処理速度が低下することがある。
• もちろん、コードの見やすさとのバランスを考慮すべき。

【Fortran】メモリ連続アクセスとキャッシュヒット

Fortranのコーディングをしていて気づいたのだけど、何が何でもメモリに連続アクセスできるように配列の添字を最優先すればいいというわけでもないっぽい。

同じ変数を使う場合、多少連続アクセスを犠牲にしてもキャッシュ・ヒットを高めることで計算処理を高速化できる場合がある。

要は両者のバランスをうまくとる必要があるということ。この辺はプログラムを走らせながら確認しないとわからない。

【Ubuntu】Ubuntu 10.04にIntel Fortran Compilerをインストールする

Ubuntu 10.04にIntel Fortran Compilerをインストールする。

Last Update: August 13, 2010.

Step.1 まずはNon-Commercial Software Downloadから「Intel® Fortran Compiler Professional Edition for Linux」を選んでダウンロード。

Step.2 次に圧縮ファイルを解凍。

$ tar xvzf l_cprof_p_11.1.072_intel64.tgz

Step.3 rootになってインストールする。

$ sudo su

を入力し、Enterキーを押す。パスワードの入力を求められるので、入力する。これでOK。

# ./install.sh

後はインストーラーにしたがっていけばOK。

Step.4 rootからログアウト。homeディレクトリの".bashrc"の最後に次の一行を追加。

. /opt/intel/Compiler/11.1/072/bin/ifortvars.sh intel64

※初めの"."の後は半角スペースであることに注意。

Note.1 ただ、インストール(Step.3)の途中に

ステップ: 4 / 7 | インストール設定 - 重要な必要条件の不足
--------------------------------------------------------------------------------
インストールに必要な次のコマンドがありません。
g++;libstdc++.so.5 (library)
--------------------------------------------------------------------------------
1. 次の問題を表示する [デフォルト]
2. [必要条件サマリー] ダイアログへ戻る

h. ヘルプ
b. 前のメニューに戻る
q. 中止
--------------------------------------------------------------------------------

と表示され、"g++"と"libstdc++.so.5"をインストールすることが要求される。"g++"については

# aptitude install g++

で解決。しかし、"libstdc++.so.5"は駄目。これについてはNote.2を参照。

Note.2 こちらを参考にして、"libstdc++.so.5"を次の手順でインストールする。

# wget http://ftp.riken.go.jp/pub/Linux/ubuntu/pool/universe/g/gcc-3.3/libstdc++5_3.3.6-17ubuntu1_amd64.deb
# dpkg -i libstdc++5_3.3.6-17ubuntu1_amd64.deb

Windowsが数値計算に向かない１つの理由

PGI® Workstation/ServerWindows®版の詳細からの引用。

# 64ビット Windows版の PGI コンパイラでは、2GB以上の単一オブジェクトを扱うことができる Linux上の -mcmodel=medium オプションと等価な機能を提供しておりません。これは、Microsoft(R) Win64 プログラミング・モデルが、まだ 2GB を超える静的な単一データオブジェクトのハンドリングをサポートしていないためです。従って、2GB を超える静的配列を扱うようなプログラムでは、実行モジュールは生成できません。たとえ Windows 64ビット版であっても、2GBを超える単一配列オブジェクトがプログラム上に存在する場合、静的な配列宣言はできず、動的な配列宣言（Allocatable 配列宣言）を行う必要があります。これは、PGI コンパイラの制約ではなく、Win64 プログラムモデルの問題です。従って、2GB以上のデータオブジェクトを扱うプログラムをコンパイルする場合は、現在のところ、このような制約がない Linux 版の方が適しています。
※　静的な単一配列オブジェクトとは、具体的に言えば、Fortranでは、COMMON 文宣言、ルーチン内でローカルに宣言された配列であっても、その総和が 2GB を超えるもの（メモリマップ上、シングルセクション上に割り付けられる）、C言語では、STATIC 宣言された配列等が相当します。なお、Fortran の allocatable 宣言された配列は 2GB 以上であっても問題なく実行モジュールを生成できます。

Windows x64 (Win64) は、64ビットのアドレス空間は使用できるものの、そのプログラミングモデルにおいては、上記のとおり 単一の静的データあるいはデータセクションのサイズが、2GB 以内に制限されております。従来の静的宣言された配列の多いプログラムでは大きな制約となり、2GB 以上使用するプログラムではリンケージでエラーとなるケースが出てきます。この制約は、コンパイラの問題ではなく、Win64 上の共通の制限ですので、フルに 2GB 以上の個別データ領域を使用するプログラムをコンパイル・実行したい場合は、「Linux」でのコンパイル・実行環境を強く推奨します。

Fortran9xとC++の違い

Fortran9x(以下単純にFortran)とC++の違いについて、知っているものから順に箇条書きにしていこうと思う。

慣れれば気にならなくなる違い

• パラメータの宣言：例えば、「a=0.5」というパラメータaを定義したい場合、Fortranでは宣言部で
  real, parameter :: a = 0.5
  と書けばよい。一方、C++ではプログラムの始めに
  #define a 0.5
  と書く方法がある。ただ、この場合aのデータ型が（普通に宣言したときのデフォルトの）倍精度実数なのか、そもそもデータ型がないのかがよくわからん。


• 倍精度実数の宣言：
  Fortranでは"9.9D0"のように実数9.9を倍精度実数として認識させるのに"D0"をつける必要があるのに対し、C++では"9.9"とすれば倍精度実数として認識される。デフォルト設定について、Fortranが単精度実数になっているのに対して、C++は倍精度実数になっている。


• 多次元配列のメモリ上の記録位置：
  Fortranでは
  integer, dimension(2,2) :: M
  で宣言された配列Mの各要素のメモリ上の配置は、M(1,1)→M(2,1)→M(1,2)→M(2,2)と、左側のindexを優先して並べられる。この場合、ループを書く場合は左側のindexほどループの中に書く必要がある。C++はこれとはまったく逆。右側のindexを優先してメモリ上に記録される。したがって、速いループを書くには、左側のindexほどループの外側に書く必要がある。


• 変数の値の入れ換え：
  メインプログラムにある変数xとy（ここでは整数とする）の値を入れ換えるのに外部手続き（Fortranではサブルーチン、C++では関数）を使いたいと仮定する。このとき、Fortranでは
  subroutine exhange(x,y)
  integer, intent(inout) :: x, y
  integer z
  z = y
  y = x
  x = z
  
  end subroutine
  と書いて"intent"を使えばOK。一方、C++はポインタかエイリアスを使う必要がある。
  • （注）このポインタやエイリアスを使うことによって最適化が不十分になってしまうという問題点があるらしい。この点で計算速度がFortranに比べて劣るか。


• 配列の添え字：Fortranでは配列の添え字が1から始まるが、C++では0から始まる。


• 1次元配列（ベクトル）の書き方：１次元配列をベクトルと捉えて、例えば、v=[1,2,3]を配列vを使って定義する場合、Fortranでは
  v=(/1, 2, 3/)
  であるが、C++では
  v={1, 2, 3}
  と書く。


• 行列(2次元配列)の初期化：
  M = [ 1, 2 ]
  　　　[ 3, 4 ]
  という行列を2次元配列Mで定義するとき、C++は
  int M[2][2] = {1, 2, 3, 4};
  とすればよいが、Fortranでは
  integer, dimension(2,2) :: M = reshape( (/1, 3, 2, 4/), (/2, 2/) )
  とする。いずれもメモリ上に記録されている順に数値が代入されていく。

C++の面倒くさい部分

• 累乗：
  Fortranでは「xのa乗」は"x**a"と書けばいいが、C++は累乗の演算子がなく、面倒らしい。

C++の、便利そうだが実はそうでもなさそうな機能

• 変数の宣言位置：Fortranでは変数の宣言はすべてプログラム冒頭（なんて言ったっけ？）で宣言する必要があるが、C++は変数を使う部分より前であればどこでもOK。
  • 後で参照しやすくするため、一ヶ所でまとめて宣言する方が良いに決まってる。

Intel Fortran Compiler for Linuxを利用する

●利用するときは、

<コンパイラがインストールされているフォルダ> ./bin/ifortvars.sh ia32

と端末で入力。''ia32''の部分はアーキテクチャによって変更する。ただし、これによると、ログアウトするまで有効になるだけで、ログイン時は毎回この設定を行う必要がある。それが面倒な場合は、ホームディレクトリにある''.bashrc''ファイルの一番最後の行に

. /opt/intel/Compiler/11.1/069/bin/ifortvars.sh ia32

を書き加えればOK。


●ソース・コードのコンパイル方法は

$ ifort filename.f90

となる。いろいろとオプションをつけることができるが、それは勉強しよう。

　上のようにコンパイルすると、''a.out''という実行ファイルが作成される。これを実行するには端末で

$ ./a.out

と入力すればよい。

Intel Visual FortranでLAPACKを使う方法

　Intel Visual FortranでMKLのFortran95インターフェイス仕様のLAPACKを使う方法についてメモしておく。マニュアルを見ただけではなかなかうまくいかず、結構試行錯誤した（汗）

　次の簡単な線形連立方程式をIntel MKL LAPACK95のdriver routineである"gesv"を使って解く：

[2 1][x_1] = [4]
[3 1][x_2] = [2]
⇔ Ax = B

（※行列による表現です。わざわざGoogle Mathematical Formulasを使うのは面倒なので。）

　次のコードは、この連立方程式をIntel MKL LAPACK95のdriver routineである"gesv"を使って解くためのコード。

program sample
use mkl95_lapack
use mkl95_precision
implicit none
real, dimension(2,2) :: A
real, dimension(2,1) :: B

a(1,:)=(/2.0, 1.0/)
a(2,:)=(/3.0, 1.0/)
b(:,1)=(/4.0, 2.0/)

call gesv( A, B )

print *, 'solution:'
print *, B(1,1)
print *, B(2,1)

end program sample

　なお、driver routineとはLU分解等を行うcomputational routineを組み合わせたもので、直接解を出力できる。driver routineを使わない場合、いくつかのcomputational routineを自前で接続しないといけない。


（注意点）
・プロジェクトの「ソースファイル」フォルダを右クリックし、「追加」->「既存の項目」を選択。"...\Intel\Compiler\11.1\048\mkl\include"から"lapack.f90"を選び、追加する。

・ソースコードの"implicit none"の前に、"use mkl95_lapack"と"mkl95_precision"（もしくは、"use lapack95"と"use f95_precision"）を書き加える。

・「プロジェクト」タブの「プロパティ」を選択。「リンカー」->「入力」と進み、「追加の依存ファイル」に、アーキテクチャがia32の場合、"mkl_lapack95.lib"を書き加える(アーキテクチャがia64の場合、"mkl_lapack95_ilp64"と"mkl_lapack95_lp64"。ただ、こちらの場合は試していないのでこれは予想)。"mkl_lapack95.lib"にIntelがLAPACKを修正したものが格納されていて、企業秘密ということになっている。コンパイラー利用者はコンパイル時にそのブラックボックスにリンクすることでIntel MKL版のLAPACKを利用する。なお、"mkl_lapack95.lib"は"...\Intel\Compiler\11.1\048\mkl\ia32\lib"にある。

・同じく「プロジェクト」->「プロパティ」->「Fortran」->「ライブラリー」->「インテル(R)マス・カーネル・ライブラリーを使用する」の部分を、「並列」か「シーケンシャル」を選ぶ（「クラスター」はおそらくPCクラスター向け）。

Fortran: 変数のデフォルト値

　Fortranでは、実数（スカラーで、配列ではないもの）のデフォルト値はわけのわからん数値になってしまうので、注意する必要がある。
　その他の実数の配列とスカラーとしての整数、整数の配列についてはデフォルト値は0になっている。