【Intel VSL】乱数発生時の注意点

　seed（初期条件）を適切に設定すること。2つ以上の乱数を発生させたいときは特に注意。例えば、quasi-random numberの場合、これは確定的な動きをするので、2つの独立するVSLの乱数生成器で(seed = 1として)「乱数発生」させても、両者が全く同じ値をとる。この場合、seed=2としなければならない。

　なお、seedはpseudo-random numberのときは初期値、quasi-random numberのときは次元を表すことに注意。詳しくはVSL Noteを参照。

To obtain a random number sequence from a given basic generator, you should assign initial, or seed values. The assigning procedure is called the generator initialization (the C language function analogous with the initialization function is srand(seed)) in stdlib.h). Different types of basic generators require a different number of initial values. For example, the seed for MCG31m1 is an integral number within the range from 1 to 231–2, the initial values for MRG32k3a are a set of two triples of 32-bit digits, and the seed for MCG59 is an integer within the range from 1 to 259–1. In contrast to the pseudorandom number generators, quasi-random generators require the dimension parameter on input. Thus, each BRNG, including those registered by the user, requires an individual initialization function. However, requiring individual initialization functions within the library interface would limit the versatility of the routines. (VSL Note, p. 31)

【LaTeX】 AucTeX

EmacsでのLaTeX環境を提供してくれるAucTeXをUbuntuにインストール。

sudo apt-get install auctex

(注)個人的に英語でしか使わないので、日本語環境pLaTeXの設定まではわからない。

インストール後、Emacsでfilename.texを編集するには

emacs filename.tex

とすればよい。

Emacs上でfilename.dviへのコンパイルは"Ctl + c”を2回繰り返し、Entキーを押す。

さらに"Ctl + c”をし、Entキーを押すことでfilename.dviを見ることができる。

filename.dviからpdfにコンパイルすには、端末やEmacs（Alt + X）から、

dvipdf filename.dvi

乱数発生

　Intel Math Kernel Library (MKL)に入ってある乱数発生ライブラリーVector Statistical Library (VSL)を使うことにする。
　Fortranの組込みであるrandom_numberはpseudo-random number generatorを基礎としているらしいが、具体的な乱数発生方法はわからない（ちゃんと調べてない）。

　一方、VSLには流行のMersenne Twister (これはpseudo-random number）やquasi-random number generatorを簡単に選択することができる。

　pseudo vs. quasiにはたくさんの議論があるよう。なんとなく、後者が積分を計算する上でパフォーマンスが良い（収束が速い）というイメージ。

link
http://demonstrations.wolfram.com/MersenneTwisterAndFriends/

RCI様様

　Intel MKLはreverse communication interface (RCI)を採用しているので、ユーザーが目的関数やヤコビアンの値を与えるサブルーチンを用意することになる。前者があるのは当たり前だが、後者があるのが面倒、なんて思ったりしたが、実はこのRCIはとてもありがたい。非線形最小二乗法で利用するヤコビアンのルーチン（central difference法）の精度がよろしくなくて、収束しない（プログラムがうまく走るか走らないかはヤコビアンに関する収束条件に大きく依存）なんてことが生じてしまう。もし解析的にヤコビアンを求めることができるなら、それを計算するサブルーチンを用意すればいいのだ。ありがたやRCI。

　自分の場合、ヤコビアンはmaximaを使ってシコシコと。なんとも原始的なアプローチ。コンピュータに詳しいなら、mathematicaとFortranを連携させてこんな無駄な作業をしなくて済むんだろうな。憧れるわー。

（追記）
　Trust Region method（信頼領域法）はglobal convergence（大域的収束）が保証されているという。でも、大域的収束の意味って、単に「少なくとも局所最適解に収束する」ことしか意味してない。もちろん、これは望ましい特徴であるのは間違いない。
　解析的にexactlyな値を与えるJacobian用のサブルーチンを用意しても、Intel MKLの1回のループでは非線形方程式を解いてくれない場合がある。
　ということで、少なくとも２通りのアプローチが考えられる。
　(1) 並列化で遺伝的アルゴリズムを採用
　(2) 収束条件を満たさなければ初期条件を変更するというsequentialな形にする
がある。将来的には別の次元（外側のループ）で並列化することを考えてるから、とりあえず後者の方法を採用してみよう。
　というか、Intel MKLの内容って、RCIと言いながら、この種の反復についてはユーザーが勝手にやってくださいねって感じなんだよな。おかしいよ。

makefileを作るべきかな

　コンパイル&リンクするファイルの数が増えると、いちいち入力するのが面倒になる。emacsでも、起動時は必ず入力しないといけない。

　近いうちに、これをみて勉強してみよう。

【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)を書く必要はない。